Finite element analysis aimed at characterizing the mechanical properties of a milky gyroid structure. The milky gyroid structure has already been presented in the article “Gyroid and minimal surfaces for CAD“.

The following article is an extract from a project carried out for the course of Advanced Techniques for the Design of Prosthetic Devices. Held for Medical Engineering at the University of Rome Tor Vergata.

Authors: Mastrofini Alessandro & Muscedere Erica

FEA targets

The following analysis aims to identify the overall stiffness of a cubic cell with side 2 elementary cells characterized by a gyroid geometry. The stiffness is calculated as a function of the thickness of the gyroid surface as well as as a function of the relative density of the cell.


Several considerations are made to set up a sequence of simulations such as to guarantee reliable results with a low computational cost.
They are considered thin shell elements. The relationship between the thickness 𝑑 and a characteristic length of the cell lies in the range:

\frac{t}{L} \in[0.01 \div 0.05]

A choice of a fine mesh was made in order to consider a total of 17078 elements.

mesh giroide
Fig. 1

To create a simulation that takes into account the real behavior of a compression test, a displacement is applied and the resulting force is measured. We are interested in characterizing the cell stiffness. This approach also allows you to make sure that all points of the top face (subjected to displacement) move uniformly, best simulating a compression test. An approach to tensions, on the other hand, would instead have to deal with the problem of ensuring uniformity of displacements or, using spider elements, with the singularity caused by the single point of one end (fig. 2a).
The lower face is constrained with a carriage type constraint such as to allow only the translations on the plane perpendicular to the compression axis (fig. 2b)

Fig. 2: Application of forces with spider element (the node on which it is impossible to apply the force is also visible) (a); carriage type constraint on the bottom face (blocks only the normal translation to the base plane) (b)

The cell has maximum dimensions of 𝐿 = 20 mm and five different thicknesses are considered

t=\{0.2,0.4,0.6,0.8,1.0\} \mathrm{mm}

Hence the relative porosity:

\frac{\rho}{\rho_{s}}=\{6.44,12.89,19.33,25.78,32.23\} \%

By applying a displacement and measuring the force it is possible to characterize the stiffness:

K=\frac{F}{\Delta x}

The results are in fig. 5b.

Ultimate stress in gyroid

A different approach was used to characterize the yield stress. An analysis was carried out by applying a load 𝐹 = 1 N on the upper surface by measuring its displacements and stress. This approach is forced to collide with the limitations of the software that does not allow the application of forces (distributing the total load through spider elements, fig. 2a) on edge elements of different types, edges and vertices. The extremal single point is therefore excluded. In estimating stiffness, this approach results in a maximum error of 11% compared to the displacement approach. The displacement approach was used to estimate the stiffness since it provides more consistent results.

Instead, by measuring the stress state and applying the Von Mises criterion, an estimate of the maximum stress within the structure is obtained. By making the ratio between the maximum reachable voltage and that applied (𝐹 / A), it is possible to calculate how much it will be amplified. Therefore, it is possible to identify the yield stress of the structure as the stress that will bring within it a maximum stress state such as to make the material overcome the yield point.

The results are reported in fig. 3b and table 2.

weight stifness ratio
Fig. 3: Stiffness with varying thickness (a) and weight / stiffness ratio (b).


From the analysis of the stiffness, a link with the Gibson-Ashby model was sought.
In the light of the available data, there is neither a completely linear model (cells dominated by compression / traction stresses) nor a quadratic trend (cells dominated by bending moments), so we sought an optimal fit to least squares according to the law generic:

\frac{K}{K_{s}} &=C_{1}\left(\frac{\rho}{\rho_{s}}\right)^{n} \\
\frac{\sigma}{\sigma_{0 S}} &=C_{2}\left(\frac{\rho}{\rho_{s}}\right)^{m}

Where, the stiffness 𝐾𝑠 = 𝐸𝐴 / L e refers to the full cube and the term 𝜌 / 𝜌s represents the porosity. For the different thicknesses:

\frac{\rho}{\rho_{s}}=\{6.44 \%, 12.89 \%, 19.33 \%, 25.78 \%, 32.23 \%\}

Searching for the optimal least squares fitting we obtain:

  • C1 = 0.414
  • n = 1.386

Results are in fig. 5a.

For the yield stress data we obtain:

  • C2 = 0.222
  • m = 1
Fig. 4: False color diagram of deformation (a) and stresses according to the Von Mises criterion (b) in the case of application of a vertical force with spider elements

We also underline that the data on yield stresses are very close to each other and very difficult to fit with a trend of this type. The fitting with a Gibson-Ashby model has the lowest error in the linear case but remains very high in order to consider the model representative of the experimental data (Fig. 5b).

Fitting with GIbson-Ashby model

For the relationship between the stiffnesses, the model shows a trend between the linear and quadratic cases. For the coefficient 𝐢1 it is possible to impose the dependence of (𝜌 / 𝜌s) and obtain the fitting:

  • 𝑛 = 1.386, best fitting con RMSE = 2.14 β‹… 10βˆ’4
  • 𝑛 = 1, linear fitting, RMSE = 2.28 β‹… 10βˆ’2
  • 𝑛 = 2, quadratic fitting, RMSE = 3.41 β‹… 10βˆ’2
Fig. 5: Comparison between porosity and stiffness (a) and yield stress (b) ratios. There is also the fitting of the Gisbon-Ashby model which is reflected in the case of stiffness with 𝑛 = 1.386 (a) but not in the case of yield stress (b). In the case of yield stress, the fitting is better using a straight line or a parabola, limited to the reference domain.

For the coefficient 𝐢2 the estimate is more complex due to the small numerical difference between the available data. However, it is also possible to obtain a good fitting (RMSE <10βˆ’3) with a straight line or a parabola, valid only in the porosity interval considered and such as not to respect the Gibson-Ashby model (fig. 5b)
In particular,

  • π‘š = 1, cell dominated by compression/tensile stress, RMSE = 2.2 β‹… 10βˆ’2
  • π‘š = 3βˆ•2, cell dominated by momentum, RMSE = 2.9 β‹… 10βˆ’2


Account to the following data, the cell appears to have a trend more similar to a cell dominated by compression / tension stresses.
For stiffness, the data match what one would expect.
Furthermore, the cells considered having a side of 2 cm, would resist loads of 30 Γ· 170 N, depending on the thickness considered, before reaching yield.

By also observing the information obtained from the buckling analysis, it is possible to see how both the cell with thickness 𝑑 = 0.2 mm and the one with 𝑑 = 1 mm will yield sooner due to the achievement of the yield point and due to instability. This is also due to the particular geometry of the gyroid which tends to distribute the tensions uniformly unless there is a slight increase in the points of contact between the octants of the unit cell. This point has a slight build-up of tension. This accumulation could also be due to a lower accuracy and precision of the starting model in that surrounding.

It should also be noted that the analysis was carried out considering ABS material for 3D printing. This would mean that the final structure could yield before the estimated limit due to the anisotropy resulting from 3D printing techniques. The mechanical properties, at most transversely isotropic in the printing surface, of the printed structures could lead to an overestimated resistance in the points where tensions accumulate.

Buckling phenomenon in the gyroid

The buckling analysis was carried out in two different ways. The results are reported in table 1. The two limit thicknesses were analyzed, 𝑑 = 0.2 mm and 𝑑 = 1 mm. A first analysis was carried out by applying a vertical load by means of spider elements leaving the lower surface constrained with trolleys (as previously described). The results show, contrary to what one would expect, a higher limit for the thinner cell.

However, by analyzing the deformation and buckling modes, this analysis shows how the structure tended to screw. This phenomenon is not very representative of what would happen to a cell that is part of a cellular structure subjected to a compressive load. Therefore, a second analysis was carried out by constraining the action of the load on the axis perpendicular to the base plane, or by preventing the loaded points from rotating with respect to the compression plane. The results show a much higher load factor than free load. Also, the thicker cell shows a higher limit.

Loading conditionsThicken π‘‘ [mm]Load factor
Tab 1: Buckling analysis results

It should be noted that, in the case of an unconstrained load, the cells would yield due to instability. In the case where the constraints on the upper surface are considered, the cells would both fail sooner due to exceeding the elastic limit.

Fig. 6: Results of the buckling analysis. Load constrained with 𝑑 = 0.2 mm (a) and 𝑑 = 1 mm (b); unconstrained load with 𝑑 = 0.2 mm (c) and 𝑑 = 1 mm (d).

Learn more about simulations

The results for the above simulations are reported.

Thicken π‘‘ [mm]Stiffness πΎ [N/mm]πœŽπ‘¦ [Mpa]Volume [mm3]πœŒβˆ•πœŒπ‘  %πΎβˆ•πΎπ‘  (10βˆ’3)𝜎0βˆ•πœŽ0𝑆 (10βˆ’3)
Tab. 2: Risultati FEA