Multiscale Homogenization Techniques for TPMS Foam Material for Biomedical Structural Applications

Multiscale techniques, namely homogenization, result in significant computational time savings in the analysis of complex structures such as lattice structures, as in many cases it is inefficient to model a periodic structure in full detail in its entire domain. The elastic and plastic properties of two TPMS-based cellular structures, the gyroid, and the primitive surface are studied in this work through numerical homogenization. The study enabled the development of material laws for the homogenized Young’s modulus and homogenized yield stress, which correlated well with experimental data from the literature. It is possible to use the developed material laws to run optimization analyses and develop optimized functionally graded structures for structural applications or reduced stress shielding in bio-applications. Thus, this work presents a study case of a functionally graded optimized femoral stem where it was shown that the porous femoral stem built with Ti-6Al-4V can minimize stress shielding while maintaining the necessary load-bearing capacity. It was shown that the stiffness of cementless femoral stem implant with a graded gyroid foam presents stiffness that is comparable to that of trabecular bone. Moreover, the maximum stress in the implant is lower than the maximum stress in trabecular bone.


Introduction
Foam-like materials based on triply periodic minimal surfaces (TPMS) have a high potential in biomedical applications. Its highly porous configuration is beneficial for cell growth and proliferation, and thus, it becomes a valid alternative to design scaffolds for biomedical applications, such as bone tissue growth. The modification of the cell density allows adjusting the mechanical properties of the implant so as to avoid unwanted phenomena such as stress shielding, which occurs when the stiffness of the implant is much higher than the stiffness of the tissue it is replacing, as a consequence of the Wolff's Law [1]. In comparison to other scaffold shapes, gyroid-based scaffolds present higher permeability than other triply periodic minimal surfaces (TPMS) [2]. Additionally, TPMS-based foams built using additive manufacturing technology have been shown to be able to achieve the necessary mechanical properties to replace human bone [3,4].
TPMS foams are the adaptation of the surface equation into a solid. The surface equation for different TPMS is the following, namely the gyroid (G) and the primitive (P): where L is the length of the unit cell. If t is equal to zero then the gyroid surface divides the space into two equal parts. In this work, the term TPMS foam refers to the sheet foam variation of the TPMS. The surface may be adapted into a solid in two different manners, sheet foam or truss-like foam. The difference between sheet foams and truss-like foams consists in how the surface is converted to a solid. The surface divides the space into two parts, and if the space enclosed by one of the sides of the surface is chosen, a solid network is obtained. On the other hand, if the surface is "thickened", a sheet network is obtained.
According to Al-Ketan et al. [5] sheet solid foams exhibit higher mechanical properties for the same apparent density as truss-like solids. In terms of surface area, which is also very relevant to the study of scaffold viability, sheet-like foams present higher surface area [6], because the surface is projected twice while the truss-like foam consists only of one repetition of the surface. It is also due to this that the gyroid sheet solid is also named double-gyroid [7].Addionally, TPMS based foam materials have been known to improve thermal properties such as conduction and convection [8,9].
An example of the two foams (gyroid and primitive) is shown in Figure 1.

Lattice Materials for Tissue Engineering
Treatments for large bone defects in the future show significant potential when using scaffold-based bone tissue engineering.
Bone implants are one of the most commonly used devices in healthcare and also one of the riskiest, as revision surgeries are often necessary. The most common cause for revision surgery is the loosening of the implant due to bone resorption, which occurs as a consequence of the stress shielding phenomenon [10].
The 3D matrix known as a bone scaffold is what enables and promotes osteoinducible cells to connect to and proliferate on its surfaces [11]. Some of the main concerns in designing scaffold are its [11] proper architecture in terms of porosity and pore sizes; 5.
controlled deliverability of bioactive molecules or drugs; Regarding mechanical properties, one can summarize the design of the bone scaffold into two main aspects: the stiffness of the scaffold, which must be low enough to allow for adequate load transfer to the bone; and strength, to withstand the naturally occurring loads in the human body. While bone presents anisotropy, there are no extremely weak directions and therefore, it is sufficient that it matches the mechanical properties of the surrounding bones as an average [3,4,12,13].
Another important characteristic that must be taken into account for implant design is porosity control. Porosity, pore size and pore interconnectivity are vital for the viability of the implant. Higher porosity facilitates the recruitment and penetration of cells from the surrounding tissue and vascularization and thus, benefits bone ingrowth [12,14]. The influence of pore size is not as consensual among researchers as several references determine different pore sizes as optimal for bone ingrowth [12,[15][16][17]. The specific surface area of scaffold is also highly relevant, as high surface areas imply a larger area for bone tissue ingrowth, and scaffolds with smaller pores will present higher surface area [12,18]. However, a higher surface area implies larger frictional forces, which is an obstacle to permeability [3]. The transportation of cells, nutrients, and growth factors require blood flood in the scaffold and thus, permeability is relevant in the design of scaffold [3,12].
In 2007, Harrysson et al. [19] mentioned in their work the possibility of using optimization algorithms to develop non-stochastic porous implants to improve the properties of the femoral stem. Since then, several works followed, with some examples of functionally graded porous cementless femoral stems. In the work of Brian et al. [20], the functional gradient was prescribed axially and radially instead of being calculated from optimization analysis.
Based on the body-centered cubic (BCC) structure, Alkhatib and Tarlochan [21] also developed femoral stem implants with graded density. In their work, the gradient consisted of layers with different densities.
Functionally graded scaffolds for stress shielding minimization may be obtained using optimization algorithms based on bone remodeling. Munteanu et al. [22], Orellana et al. [23] used structural optimization to reduce the stiffness of the femoral stem. However, in their work, the implant is not porous.
In the works of Pais et al. [24,25], which consisted on the development of porous implants, it was shown that the bio-inspired algorithm combined with an experimental material law for the gyroid infill obtained through 3D printing (FFF) [26,27] achieved structural configuration with adequate stiffness to replace human bone.

Homogenization Techniques
Multiscale modeling refers to an approach in which analysis of the material is conducted at one length scale but the outcomes of the analysis are referent to several properties of the material at another length scale [28].
The use of numerical homogenization techniques allows for significant savings in computational time. Often in composites, it is not necessary (and inefficient) to model the entire structure of the composite. Instead, only a representative region is chosen to model all the constituents of the composite [29]. This approach can be extended to lattice materials, by simplifying the assumption, where the composite presents two or more phases (fiber and matrix), and the lattice will only present one phase, being the rest a void phase. In summary, the porous material is transformed into an equivalent solid material with homogenized properties [30].
The mechanical properties of the obtained scaffold are given as a function of their relative density. The relative density ρ * of the scaffold is given by where ρ is the density of the lattice and ρ m denotes the density of the constituent material. Therefore, from now on, whichever property is being studied, it being stiffness or strength, it is presented as a dependency on this parameter. When the lattice unit cell presents cubic symmetry, which occurs for most cases indicated in the previous section, the stiffness tensor C will only present three independent constants, C 11 , C 12 and C 44 . Therefore, the stress-strain relation can be given as: where σ ij are the components of the macroscopic stress tensor and ε ij are the macroscopic strain tensor components.

Boundary Conditions
The most correct way to evaluate the mechanical properties of the lattice material by analyzing one unit cell is through the application of periodic boundary conditions, instead of uniform tensile or linear displacement boundary conditions. The formulations for some possible boundary conditions are presented next.

Periodic Boundary Conditions
Considering a periodic RVE (Figure 2) Ω, the boundary Γ of the RVE can be decomposed into two parts Γ + and Γ − . Each point x + on Γ + is connected to just one unique point x − on Γ − and the normal vectors to each point are symmetrical so that n + = −n − [29]. The local displacement field u can be decomposed as where the macroscopic displacement gradient H is equal to the macroscopic strain ε up to a rotation,û = u 0 + H · x is the macrodisplacement that corresponds to the applied strain, andũ(x) is a micro-displacement. The unit cell simulations assume that in the bulk of the materialũ (x + a · λ i · e i ) =ũ(x) (6) for any integer λ i and all positions of x, and where a denotes the unit-cell length and e i denotes the principal directions. Thus, it is assumed that the micro-displacement field shares the periodicity of the lattice [31].
The periodicity of micro-displacements is enforced by kinematically constraining the difference in the displacements of paired nodes, and setting this difference as equal to the displacement deduced from the macroscopic strain [31] so that

Linear Displacement and Uniform Traction Boundary Condition
Linear displacement boundary conditions consist of applying to the boundary of the RVE the displacement field that would occur if the strain were uniform inside the RVE where ε(u) is the micro strain and x is the position vector on the boundary ∂Ω. This method presents the advantage of having no restriction to its application except that no rigid part must intersect the boundary [32]. Uniform traction boundary conditions consist of applying on the boundary of the RVE the stress vector field that would occur if the stress were uniform inside the RVE where σ is the Cauchy stress tensor. Using this method, there is also no restrictions except that no holes must intersect the boundary [32].

FE Homogenization
The finite element (FE) homogenization technique can be used in order to predict the effective elastic, as well as the elastic-plastic properties of the material. This technique excels at obtaining the properties of materials with complex microstructures, such as composites [29], even though in this case the focus is on lattice materials.
Considering the RVE Ω, any micro-field f, such as stress or strain within the RVE can have the following average functions defined where V is the volume of the lattice. The effective mechanical properties do not depend on the body forces or boundary conditions. Thus, to predict the properties of the material, the following weak-form quasistatic equilibrium differential equation is considered The boundary conditions do not affect the material properties but they must satisfy Hill's energy law, which states that the energy on the micro-level has to be the same as the effective energy for the homogenized material For any point of the RVE, the constitutive model is given as Based on this last relationship, as well as the quasi-static equilibrium differential equation and the boundary condition which satisfies the Hill energy principle, the stress σ(x) and strain ε(x) can be obtained through FE analysis and then the average values for stress and strain can be computed through where n e is the number of elements in the model, n q is the number of integration points in the element e and W(y I ) is the weight of the integration point. σ ij (y I ) and ε ij (y I ) are the stress and strain respectively, evaluated at the integration point y I . The average values, which consist of a volume average of the properties along the material, can then be used to calculate the effective elastic stiffness tensor C and the effective stress tensor for elastic-plastic analysis.

Mechanical Properties Elastic Properties
The elastic constants are obtained from the following relations from the elastic constants C 11 and C 12 where each component of the constitutive matrix is obtained from the linear elastic analysis. By enforcing a unit strain, the stress tensor will correspond to one line in the constitutive matrix C.
The isotropy of the surface may be more easily perceived using Young's modulus surface, which is a representation of Young's modulus in every direction. For a perfectly isotropic material, where Young's modulus will be the same in all directions, the surface should have the appearance of a sphere. The isotropy of the structure is proved by the Zener ratio A H , which allows the near equality in Equation (18). If A H is close to 1, the structure is isotropic [33].

Plastic Properties
Unlike elastic deformation, plastic deformation is irreversible when the loading is removed. When the material reaches its yield point, the material starts displaying plastic behavior and lower material stiffness. After reaching this point, plastic behavior may follow several different models.
If the material is assumed to display linear strain-hardening characterized by the tangential modulus E T . The total deformation caused by a stress increase dσ is given by and so a strain-hardening parameter H is defined as which can also be given by which corresponds to the slope of the strain-hardening portion of the stress-strain curve [34].
In the elastic domain, the stiffness exhibited by the material will be given by where F is the applied force, δ is the shown displacement, E is the Young's modulus, A is the cross-sectional area and L is the length. Thus, the stiffness matrix for an element is If the material reaches its yield stress due to an increase in F, that increase dF will generate a dδ of dδ = dε e + dε p L Considering that the increase in force is given by the tangential stiffness of the material will be which leads to an element stiffness of Finally, the equivalent yield strength of the lattice is given by the apparent stress-strain curve, by retrieving the stress level that leads to 0.2% plastic strain. The apparent stress (macro-stress) σ 11 is obtained as while the apparent strain (macro-strain) is given by:

Scaling Laws
Cellular materials are usually described by scaling laws which correlate the homogenized mechanical properties to the volume fraction of the cell [35]. The volume fraction of the cell is given by The homogenized Young's modulus is given by and the homogenized yield stress is given by The n coefficient depends on whether the cell exhibits stretching-dominated or bendingdominated behavior. The stretching/bending dominated designation is referent to the deformation mode of cell walls. Usually, higher density foams will tend to show near stretching-dominated behavior while at lower densities these will tend to show bendingdominated behavior. Table 1, adapted from the review paper by [36] summarizes the n coefficient in Equations (32) and (33), while the D 1 coefficient will range between [0. 1,4] and the D 2 coefficient will range between [0.1, 1] [36]. Table 1. Exponent coefficient n in the Gibson-Ashby model for stretching-dominated (SD) and bending-dominated (BD) structures.

SD BD
Homogenized Young's modulus 1 2 Homogenized yield stress 1 1.5 In the literature, however, it is common to find modifications to the Gibson-Ashby model, where none of the values shown in Table 1 are necessary. It's also relevant to mention that this model will change depending on the manufacturing process, material, and geometry of the cells. For some complex materials characterized as a function of their relative density, the fitting to the model allows determining if the structure exhibits bending or stretching-dominated behavior.
In this work, the correlations will describe the homogenized mechanical properties using polynomial models as shown in (34) and (35)

Materials
For the homogenization analysis, each unit cell model presents a Young's modulus equal to 3000 MPa and a Poisson ratio equal to 0.3. For the plasticity analysis, the material is considered to be elastic-perfectly-plastic, being that three different yield stress values were tested, 10 MPa, 30 MPa and 50 MPa.

Models
The STL models of each surface were created in MATLAB by computing the equation with an isolevel value t equal to zero. Each facet in the triangular mesh was projected in its normal direction by half of the thickness value, in the positive and negative directions. The solid model was meshed into quadratic 10-node tetrahedral elements from the STL of the outer surface. For both the gyroid and the primitive three models were obtained, each corresponding to a low density, medium density, and high density. Each density level corresponds to a different wall thickness to RVE edge length ratio, as can be visualized in Table 2.

Boundary Conditions
For the elastic analysis, two different types of boundary conditions were used, namely periodic boundary conditions and homogenization boundary conditions.

Periodic Boundary Conditions
The enforcement of periodic boundary conditions required the definition of paired nodes in the faces. Each node on one face is paired with the closest node on the opposite face. Due to the fact that both faces do not present the same number of nodes, some nodes do not have any equation constraint associated with them. The primitive structure unit cell does not have edges or vertices, while the gyroid structure has edges and vertices. The nodes on the edges and vertices in the gyroid model were ignored for the enforcement of periodic boundary conditions.

Homogenization Boundary Conditions
These boundary conditions, adapted from [37] were used both in the elastic and plastic analysis. If a normal strain is applied, for example ε 11 , the following displacements are imposed: whereas if a shear strain is applied, for example ε 12 , the following displacements are imposed The application of these boundary conditions ensures that the prescribed displacement in the top of the RVE is a unit displacement and all other components of the strain tensor remain zero.

Bio-Inspired Remodelling Algorithm
In biomedical implants, stress shielding can be minimized as the stiffness of the implant is optimized according to the bone stiffness [13].
The bio-inspired remodelling, taking into consideration that bone remodelling phenomena acts as structural optimization [38], updates the density of the structure iteratively. Therefore, the bio-inspired algorithm, based on [39] can be summarized according to Figure 3.
In summary, first, the mechanical properties of a point are determined from E[MPa] = f E ρ g/cm 3 so the constitutive matrix is obtained and then, linear-elastic analysis is performed. Using scaling laws, the correlations are adjusted to the material of the implant. Then, the critical variable is used in order to determine which points must be remodelled, meaning, having its density increased or decreased. Points with higher stress or strain energy density can have its density increased and points with lower stress or strain energy density may have its density decreased. As verified in bone remodelling phenomenon (in which unsolicited bone will decay while a stress stimulus promotes bone growth), in tihs bio-inspired algorithm points with the highest values of the critical variable will increase its density and the points with the lowest values of the critical variable will have its density decreased. This density update is enforced using σ[MPa] = f σ ρ g/cm 3 . The remodelling process occurs until convergence, or, until the control density defined by the user is achieved. Figure 3 shows a schematic flowchart of the remodelling algorithm.

Stress Shielding Evaluation
Stress shielding is usually evaluated as a function of the mismatch between the stresses in the bone and stresses in the implant [40].
In order to avoid stress shielding, the stiffness of the implant must be lower than the stiffness of bone. Therefore, the loads transferred to the bone are higher, resulting in a stress stimulus leading to bone growth. In order to establish a term of comparison, a model with the geometry of the area being remodelled in the implant is modeled as bone, through a bone remodelling algorithm [39]. The stiffness of that bone, is then compared to the stiffness of the implant.
The stiffness of the implant is evaluated as follows: where and d = ∑ d t z n t (39) where t refers to the nodes where the load is applied and n t is the number of nodes where the load is applied. Because at each iteration the variables are scaled to the elastic limit of the porous material, the applied force is calculated from the stress field as: (40) where B is the deformability matrix and σ is the stress tensor in Voigt notation. Figure 4 shows the plot of the relative density of the cell as a function of its thickness and unit cell length. For the same wall thickness, the gyroid presents a higher relative density than the primitive structure.

Elastic Properties
The elastic properties of both structures can be consulted in Table 3. The gyroid structure presents higher Young's modulus and higher isotropy than the primitive structure. Additionally, the use of periodic boundary conditions in the gyroid structure led to higher Young's modulus values than the primitive structure.

Plastic Properties
For the gyroid, the apparent stress-strain plots are shown in Figure 5 for RVEs consisting of 1, 2, and 3 unit cells on each side, while for the primitive structure, only 2 by 2 cells were tested, where the results are shown in Figure 6. It can be concluded that the gyroid structure presents higher strength than the primitive structure. For both geometries, it can be seen that higher-density foams exhibit hardening (near positive slope in the plastic domain of the stress-strain curve) while medium-density foams present near perfectly plastic behavior (near zero slopes in the plastic domain of the stress-strain curve) and lower density foams present a yield point.
Because the aim of the plastic analysis is to obtain the elastic limit of the material, these analyses only consisted of a maximum strain of 10%. In order to develop material laws that can be used in the bio-inspired algorithm, the apparent yield stress was obtained from the apparent stress-strain curves, resulting in the data shown in Tables 4-6 for the gyroid structure and Table 7 for the primitive structure. Table 4. Summary of plastic properties of the gyroid structure-low density.  Table 5. Summary of plastic properties of the gyroid structure-medium density.  Table 6. Summary of plastic properties of the gyroid structure-high density.

Material Law and Comparison to Literature Data
By adjusting the obtained points to a third degree polynomial function, the coefficients in Equations (34) and (35) are obtained, which can be consulted in Table 8. Both structures exhibit bending-dominated behavior instead of stretching-dominated behavior, meaning that the cell walls tend to bend instead of being compressed. Stretchingdominated structures tend to be stiffer which can be an issue regarding the minimization of stress shielding. This is in agreement with other studies found in the literature, such as the work of Abueidda et al. [41] where gyroid sheet foams were modeled to fit the Gibson-Ashby model and it was concluded from the coefficients of that model that the structure exhibits bending-dominated behavior. Al-Ketan et al. [5] also concluded that these sheet-like structures exhibit behavior between stretching and bending dominated.
The developed material law assumes that the material is isotropic and perfectly plastic. From the analysis of elastic properties in Table 3 it was concluded that the primitive structure is less isotropic than the gyroid structure. Figure 7 shows the developed laws in comparison to the several works from experimental testing in the literature. Most experimental works use models with a higher number of unit cells than this work. It can be seen from the figures that the error in homogenized modulus is, for most cases, low, and thus, the material law will accurately describe the mechanical properties of gyroid foams at several densities. Moreover, this highlights the importance of experimental testing in the validation and calibration of the material law. The geometry presented by these complex materials is almost always unattainable through conventional manufacturing, leaving additive manufacturing as the only solution to obtain physical parts or prototypes. However, additively manufactured parts tend to be highly anisotropic, showing different mechanical properties in the building direction. The developed correlations do not account for the anisotropy in the manufacturing process. The anisotropy can be minimized by choosing a process with higher isotropy, for example, SLM or SLS if dealing with metals or polymers respectively, instead of extrusion-based processes, where the interlayer bonding effect highly influences the mechanical properties of printed parts.
Additionally, the plot shows, through the shaded area the upper and lower bounds of the Gibson-Ashby law, according to the work of Maconachie et al. [36]. The Gibson-Ashby model, which predicts the coefficients for metallic foams, was developed based on analytical models, extensive testing on polymeric foams, and empirical fits to experimental data [36,49]. Finally, in order to evaluate the accuracy of the developed models in comparison to the literature, the mean absolute percentage error was calculated for each of the models (Young's modulus and yield stress for both unit cells). These values are presented in Table 9.

Study Case-Femoral Stem
In order to establish conclusions regarding the suitability of the gyroid material implant, remodeling analysis was performed on a cementless femoral stem model. These implants are prone to stress shielding, and therefore, the use of a porous structure can reduce stiffness and enhances bone ingrowth [40]. Previously, it was established that the gyroid foam is more isotropic than the primitive foam. Thus, as the remodeling model assumes an isotropic material, the gyroid material model is used.

Implant Model
The femoral stem was modeled according to geometries found in the literature using CAD software and then meshed into linear hexahedral elements. The final mesh has 2808 nodes and 2020 elements.
Only the proximal area of the stem implant is remodeled, the neck and taper area as well as the distal part of the stem remains as solid material (Figure 8). The chosen base material is Ti-6Al-4V alloy whose properties are shown in Table 10, as well as a comparison to femoral trabecular bone, which is used later as reference. The base of the implant was fixed while a vertical force was applied to the top of the taper. This force was set as unitary as at each iteration of the remodeling algorithm the variables are brought to the elastic limit.  Figure 8 shows the evolution in stiffness of the implant and average relative density along the iterations of the remodeling algorithm. The final density distribution is also presented at the central cross-section.

Optimized Implant
Finally, in order to evaluate the suitability of the implant, its strength must also be enough to support the loads in the bone. Figure 9 shows that the load that leads the trabecular to its elastic limit is lower than the load that leads the implant to its elastic limit. The force that takes the implant to the elastic limit is 4248.78N, while the force that can be applied to the bone while remaining below the elastic limit is 1280.85N. Thus, the implant is also suitable in terms of strength.

Conclusions
With this work it was possible to correlate the mechanical properties of two different TPMS, the gyroid and the primitive with the apparent density. The simulations were run considering isotropic and perfectly plastic properties of the solid material and it was verified that this assumption works as a reasonable approximation in comparison to literature data from experimental tests. Additionally, it was shown that the remodelling algorithm provided an optimal configuration for a cementless femoral stem implant based on the gyroid structure. Future directions for the study should include the experimental testing of implant prototypes for final validation of the material law and remodelling algorithm.