Chrysoula Chatzigeorgiou

Boris Piotrowski

Yves Chemisky

Pascal Laheurte

Fodil Meraghni

Numerical investigation of the effective mechanical properties and local stress distributions of TPMS-based and strut-based lattices for biomedical applications

published or not. The documents may come

Numerical investigation of the effective mechanical properties and local stress distributions of TPMS-based and strut-based lattices for biomedical applications

Chrysoula Chatzigeorgiou, Boris Piotrowski, Yves Chemisky, Pascal Laheurte,

Introduction

Over the past decades, the number of orthopedic implant operations has increased [START_REF] Yang | Investigation on the orientation dependence of elastic response in Gyroid cellular structures[END_REF]. This increase is accompanied by the improvement of the patients' postoperative life conditions by solving the biomechanical problems that occur between the implant and the bone and by enhancing the mechanical biocompatibility and durability. Indeed, since the patients are young, they are in demand of orthopedic implants that allow them to continue to practice physical activities [START_REF] Learmonth | The operation of the century: total hip replacement[END_REF]. A commonly observed undesirable phenomenon after an orthopedic implant operation is bone resorption and bone mass reduction (Bono et al., 1999;[START_REF] Fischer | Élaboration in situ d'alliages de titane de structures architecturées par fabrication additive : application aux dispositifs médicaux implantables[END_REF][START_REF] Huiskes | The relationship between stress shielding and bone resorption around total hip stems and the effects of flexible materials[END_REF] which is caused by stress shielding, described by Wolff's law [START_REF] Wolff | The Law of Bone Remodelling[END_REF]. It is mainly induced by the significant stiffness mismatch between the implant and the neighboring bone [START_REF] Wolff | The Law of Bone Remodelling[END_REF]. In more detail, the bone, which has low stiffness that varies in the range from 0.02 (trabecular bone) to 30 GPa (cortical bone), is subjected to a complex dynamic load [START_REF] Al-Tamimi | Topology optimization to reduce the stress shielding effect for orthopedic applications[END_REF]. With the presence of a more rigid implant, the load distribution is perturbated: the implant bears a significant part of the load, contrary to the peri-implant bone which is relieved. Stress shielding occurs in this case, impacting the implant-bone interface durability [START_REF] Al-Tamimi | Topology optimization to reduce the stress shielding effect for orthopedic applications[END_REF]. The bone resorption and mass reduction affects hence the long-term performance of the implants, their stability and induces a pain feeling to the patients [START_REF] Yamako | Improving stress shielding following total hip arthroplasty by using a femoral stem made of β type Ti-33.6Nb-4Sn with a Young's modulus gradation[END_REF].

The orthopedic implants are made of various biocompatible metals and alloys; such as titanium alloys, stainless steel, and cobalt alloys [START_REF] Al-Tamimi | Topology optimization to reduce the stress shielding effect for orthopedic applications[END_REF][START_REF] Schuh | Second generation (low modulus) titanium alloys in total hip arthroplasty[END_REF][START_REF] Wang | Topological design and additive manufacturing of porous metals for bone scaffolds and orthopaedic implants: a review[END_REF]. Among all these materials, titanium and its alloys are mostly utilized because they combine great biocompatibility, high corrosion and fatigue resistance, and very good mechanical properties [START_REF] Didier | Mechanical stability of custommade implants: numerical study of anatomical device and low elastic Young's modulus alloy[END_REF][START_REF] Wang | Topological design and additive manufacturing of porous metals for bone scaffolds and orthopaedic implants: a review[END_REF][START_REF] Yan | Ti-6Al-4V triply periodic minimal surface structures for bone implants fabricated via selective laser melting[END_REF]. However, titanium alloys-based implants are significantly stiffer than the surrounding host bone, and consequently, the stress shielding phenomenon takes place [START_REF] Wang | Topological design and additive manufacturing of porous metals for bone scaffolds and orthopaedic implants: a review[END_REF][START_REF] Yamako | Improving stress shielding following total hip arthroplasty by using a femoral stem made of β type Ti-33.6Nb-4Sn with a Young's modulus gradation[END_REF]. As a result, the reduction of the stress shielding effect is of high importance and it can be achieved in two ways [START_REF] Elmay | Optimisation of mechanical properties of Ti-Nb binary alloys for biomedical applications[END_REF][START_REF] Laheurte | Mechanical properties of low modulus ββ titanium alloys designed from the electronic approach[END_REF]. The first way is the choice of bulk material for the implant, with an apparent elastic modulus as close as possible to the Young's modulus of the natural bone [START_REF] Didier | Mechanical stability of custommade implants: numerical study of anatomical device and low elastic Young's modulus alloy[END_REF][START_REF] Mahmoud | Lattice structures and functionally graded materials applications in additive manufacturing of orthopedic implants: a review[END_REF][START_REF] Schuh | Second generation (low modulus) titanium alloys in total hip arthroplasty[END_REF], making use of solid-solid reversible phase transformations. However, there is a limitation regarding the mechanical properties of these low apparent elastic modulus alloys. Although they are subjected to thermomechanical treatment to reduce their stiffness, the Young's modulus cannot be less than 50 GPa, a value still higher than the stiffness of the natural bone [START_REF] Didier | Mechanical stability of custommade implants: numerical study of anatomical device and low elastic Young's modulus alloy[END_REF]. Furthermore, the second solution consists of the design of a porous structure for the implant instead of a continuum bulk material, and thus the aforementioned material-dependent limitation could be overcome. It has been shown that porous structures, specifically lattice structures, could be tailored to mimic the mechanical properties of the bones and also allow higher tissue regeneration and osseointegration [START_REF] Schuh | Second generation (low modulus) titanium alloys in total hip arthroplasty[END_REF][START_REF] Wang | Topological design and additive manufacturing of porous metals for bone scaffolds and orthopaedic implants: a review[END_REF][START_REF] Yan | Ti-6Al-4V triply periodic minimal surface structures for bone implants fabricated via selective laser melting[END_REF].

Porous materials, and particularly lattice structures, present considerably complex geometries at two-scales (the lattice and the implant shape). However, thanks to the recent improvements in additive manufacturing (AM) technologies, the fabrication of complex geometries, such as lattices, is nowadays easily achievable [START_REF] Van Bael | The effect of pore geometry on the in vitro biological behavior of human periosteum-derived cells seeded on selective lasermelted Ti6Al4V bone scaffolds[END_REF]. They can be fabricated by combining both optimized topological design and the right additive manufacturing process concerning the application [START_REF] Wang | Topological design and additive manufacturing of porous metals for bone scaffolds and orthopaedic implants: a review[END_REF]. For example, additively manufactured lattices and porous structures have been studied for various medical applications, such as tissue engineering [START_REF] Hollister | Porous scaffold design for tissue engineering[END_REF][START_REF] Valainis | Integrated additive design and manufacturing approach for the bioengineering of bone scaffolds for favorable mechanical and biological properties[END_REF], implants (Ataee et al., 2018b), human phantoms in radiotherapy [START_REF] Okkalidis | Assessment of 11 available materials with custom three-dimensional-printing patterns for the simulation of muscle, fat, and lung hounsfield units in patient-specific phantoms[END_REF].

The first step in designing a lattice structure is the particular emphasis on the topological aspect of its architecture. It is important to notice that the topology of a lattice greatly impacts its mechanical properties. This can be viewed as an opportunity to tune these mechanical properties locally to retrieve the mechanical response of the bone to be replaced [START_REF] Wang | Topological design and additive manufacturing of porous metals for bone scaffolds and orthopaedic implants: a review[END_REF][START_REF] Xu | Design of lattice structures with controlled anisotropy[END_REF]. Additionally, the relative density or the volume fraction of the lattices plays an important role in the control of their effective properties [START_REF] Al-Ketan | Functionally graded and multi-morphology sheet TPMS lattices: design, manufacturing, and mechanical properties[END_REF]. However, stress concentrations may be observed in various areas of the lattice structures, due to the nature of their topology [START_REF] Lohmuller | Stress concentration and mechanical strength of cubic lattice architectures[END_REF]. Hence, it is also essential to minimize the stress concentration on the lattices to improve mechanical performances such as fatigue response. The stress concentration observed in several patterns [START_REF] Al-Ketan | Functionally graded and multi-morphology sheet TPMS lattices: design, manufacturing, and mechanical properties[END_REF] should be avoided to maximize the durability properties of the implant [START_REF] Soro | Quasi-static and fatigue properties of graded Ti-6Al-4V lattices produced by Laser Powder Bed Fusion (LPBF)[END_REF].

More specifically, the lattices or the architectured materials are a category of porous materials that possess a defined and controlled topology. There are two distinct lattice topology families, the strut-based and the surface-based topologies [START_REF] Al-Ketan | Topology-mechanical property relationship of 3D printed strut, skeletal, and sheet based periodic metallic cellular materials[END_REF]. The first family includes periodic topologies based on cylindrical beams/struts, for instance, octet-truss [START_REF] Deshpande | Effective properties of the octet-truss lattice material[END_REF], BCC [START_REF] Maskery | An investigation into reinforced and functionally graded lattice structures[END_REF], the cubic Bravais lattices, and some more complex combinations of the beams [START_REF] Favre | A continuous crystallographic approach to generate cubic lattices and its effect on relative stiffness of architectured materials[END_REF][START_REF] Lohmuller | Stress concentration and mechanical strength of cubic lattice architectures[END_REF]. The second family includes topologies based on the Triply Periodic Minimal Surfaces (TPMS) concept whose lattices exhibit smoother connections. Some TPMS-based lattices are the Schoen's Gyroid, Schoen's IWP, Schwartz's Primitive, and Neovius [START_REF] Schoen | Infinite Periodic Minimal Surfaces without Self-Intersections[END_REF]. Furthermore, the TPMS-based lattices can be divided into two subfamilies, the sheet and the skeletal, based on their geometry [START_REF] Al-Ketan | Topology-mechanical property relationship of 3D printed strut, skeletal, and sheet based periodic metallic cellular materials[END_REF][START_REF] Li | Comparison of mechanical properties and energy absorption of sheet-based and strut-based gyroid cellular structures with graded densities[END_REF]. The difference between the sheet and skeletal types is explained in detail in the materials and numerical methods section.

Currently, the investigation of the mechanical properties of the lattices, and especially the TPMS-based lattices, is a topic of high interest in the whole scientific community and not only in biomechanical applications. For example, [START_REF] Maskery | Effective design and simulation of surface-based lattice structures featuring volume fraction and cell type grading[END_REF] have numerically investigated skeletal TPMS-based lattices with a constant volume fraction as part of their work. Moreover, [START_REF] Li | Comparison of mechanical properties and energy absorption of sheet-based and strut-based gyroid cellular structures with graded densities[END_REF] have implemented a homogenization method to study the mechanical properties of both sheet and skeletal Gyroid lattices and also investigated their anisotropic behavior. Furthermore, [START_REF] Lu | The anisotropic elastic behavior of the widely-used triply-periodic minimal surface based scaffolds[END_REF] have studied the mechanical behavior and the anisotropy of five skeletal TPMS-based lattices using a homogenization method and periodic boundary conditions in a FE software for a biomedical application. [START_REF] Chen | On hybrid cellular materials based on triply periodic minimal surfaces with extreme mechanical properties[END_REF] have investigated the effective elastic properties of four sheet TPMS lattices using a Fast-Fourier Transform (FFT) full-field simulation method adapted for such periodic cells. They have also combined two different types of TPMS to create a hybrid structure and to improve their mechanical properties.

Also, [START_REF] Al-Ketan | Topology-mechanical property relationship of 3D printed strut, skeletal, and sheet based periodic metallic cellular materials[END_REF] have fabricated steel lattices with strut-based, sheet TPMS-based, and skeletal TPMS-based topologies and have compared their mechanical properties experimentally while Ataee et al. (Ataee et al., 2018a, p.) have fabricated skeletal Gyroid scaffolds with commercially pure titanium employing a Selective Laser Melting (SLM) machine to determine their mechanical properties, which demonstrated promising results in terms of biomechanical applications. [START_REF] Yan | Ti-6Al-4V triply periodic minimal surface structures for bone implants fabricated via selective laser melting[END_REF] have fabricated skeletal Gyroid and Diamond lattices, using Ti-6Al-4V, to mimic the porosity and the stiffness of the natural bone. Several numerical studies in this field are summerized in Table 1 for comparison.

To the best of our knowledge, no previous study has investigated and compared both the effective mechanical properties and local stress distribution at the same time of a large range of sheet, skeletal TPMS-based, and strut-based lattice structures. This enables the comparison of isoelasticity lattice structures leading to an optimal one according to a specific application.

The present study proposes an integrated approach aimed at systematically investigating and comparing various TPMS-based and strutbased lattices with respect to elastic properties, stiffness anisotropy, and local stress distribution for biomedical implants. The numerical simulation procedure consists of various steps that include the geometrical design and FE-adapted mesh generation procedures of the lattices, the application of periodic boundary conditions, and implementation of a homogenization method towards the computation of their effective mechanical properties. The second section presents all the aforementioned steps of the methodology. The following section starts with the numerical simulation and comparison of the effective mechanical properties of the lattices based on the topology and the volume fraction.

Next, the investigation of stiffness anisotropy is performed and finally, the equivalent Von Mises stresses distribution is analyzed for several TPMS-based and strut-based lattices. The obtained numerical results are discussed in detail in the relevant section and subsequently, they are specifically detailed in two cases demonstrating the capabilities of the developed approach. The paper ends with some concluding remarks that summarize the main outcomes of the work.

Materials and numerical methods

Lattice structures; design and generation

A lattice structure is defined by a periodic unit cell considering V solid as the volume of the solid part and V total as the total volume of the bounding box such that:

Volume fraction = V solid V total (1)
The geometry of the lattices necessarily controls their volume fraction.

In this work, the surface-based lattices that have been investigated are the Schoen's Gyroid, Schwartz's Primitive, Neovius, and IWP topologies, while the strut-based lattices are the Octet-truss and Diagonal topologies. The Diagonal topology may also be referred to as bodycentered-cubic (BCC) topology. The term Diagonal is used in this paper. The numerical study of a single unit cell is sufficient for the study of the mechanical properties of various lattice topologies in the elastic regime. The implemented method for FE analysis is described in 2.2. and 2.3. sections.

The design and the generation of the strut-based lattices are performed using the computational environment of the software Abaqus/ CAE 2018. A Python script is developed to enable and automate the design procedure of the strut-based lattices. The strut-based lattices consist of symmetric and periodic struts. The imported data to the script are the parametric shape and the number of the basic solids that constitute the unit cell, their orientation, and their geometric characteristics, such as the length and the radius. Then, the appropriate periodic boundary conditions are applied to the lattice structure through a homogenization method, which is described in detail in the Numerical Periodic Homogenization section.

In Fig. 1 the two strut-based unit cells that are among the investigated lattice types in this study, i.e. Diagonal and Octet-truss unit cells, are presented.

The lattice structures based on the TPMS are designed with a different approach. The TPMS are periodic surfaces in all three directions that have a zero mean curvature at every point and each one divides the 3D space into two non-interconnected sub-spaces. The mathematical equations used for their definition offer a pretty accurate approximation of the surfaces [START_REF] Yan | Ti-6Al-4V triply periodic minimal surface structures for bone implants fabricated via selective laser melting[END_REF]. The sinusoid level-set equations have the form of f (x, y, z) = t , where t is a constant value [START_REF] Maskery | Effective design and simulation of surface-based lattice structures featuring volume fraction and cell type grading[END_REF]. The value of the t parameter controls the offset of the minimal surfaces and thus the volume ratio between the two sub-spaces. In the case of t = 0 each volume fraction of the sub-spaces is 50%. When the parameter t has a positive or negative value, the volume ratio between the two sub-spaces changes. Therefore, the volume fraction of the obtained lattice structures depends on the value of the design parameter t.

The investigated TPMS and their corresponding mathematical t-parametrized equations are presented in what follows:

(a) Schoen's Gyroid: Furthermore, there are two different types of TPMS-based lattices, the skeletal, and the sheet. The skeletal type can be obtained when one of the two sub-spaces is the solid part, while the other is the void, i.e., f (x, y, z) ≥ t or f (x, y, z) < t. The sheet type can be obtained when a prescribed thickness is given to the surface generated by an equation with t = 0, or mathematically, it can be obtained when the generated solid part satisfies the double inequalityt ≤ f(x, y, z) ≤ t (Al-Ketan et al., 2020). There is a relation between the design parameter t and the volume fraction of the unit cells.

cos(x) * sin(y) + cos(y) * sin(z) + cos(z) * sin(x) = t (2) (b) Schwartz's Primitive: cos(x) + cos(y) + cos(z) = t (3) (c) Schoen's IWP: cos(x)*cos(y) + cos(y)*cos(z) + cos(z)*cos(x) = t (4) (d) Neovius: 3 * (cos(x) + cos(z) + cos(y)) + 4 * cos(x) * cos(z) * cos(y) = t (5)
The methodology adopted to generate TPMS-based unit cells is the design of the two opposite surfaces with the offset of -t and t, respectively. The simultaneous design of a sheet and two skeletal unit cells within a bounding box is obtained. The summation of the three volume fractions is equal to the volume fraction of the bounding box, i.e. 100%, because their volumes are complementary. The relative density of each skeletal cell is less or equal to the value of 50%.

Moreover, the nature of particular TPMS limits the volume fraction of the generated unit cells. The t parameter must be ranging in an allowable interval to avoid a lack of surface connectivity that generates several disconnected surfaces. Indeed, when the values of the t param-eter exceed an allowable value, the surface presents discontinuities, such that a structure cannot be fabricated. The allowable value is different for each topology and can be founded by a CAD analysis of the sets of closed domains. As a consequence, the feasible range of the effective elastic properties varies for the TMPS-based unit cells. For example, the Neovius sheet unit cell can exhibit a volume fraction in the range of 5% -31%. Nevertheless, this limitation is not important in mechanical applications since the lattice structures are usually utilized at low volume fractions [START_REF] Lohmuller | Architectural effect on 3D elastic properties and anisotropy of cubic lattice structures[END_REF].

In Fig. 2., the sheet and skeletal parts of a unit cell for Schoen's Gyroid and Schwartz's Primitive TPMS are presented. It is obvious that the two Gyroid skeletal parts are identical, i.e., the red and the blue parts, with the same volume fraction. However, it is not the case for the Schwartz's Primitive, where the two skeletal parts have different shapes. This property depends on the type of the TPMS. The IWP and Neovius topologies present two non-identical skeletal parts, as well.

Once the unit cells have been designed, a FE-adapted mesh is generated and application of the periodic boundary conditions is carried out. This step is described in section 2.2.

Numerical Periodic Homogenization

A FE-compatible periodic homogenization method [START_REF] Suquet | Elements of homogenization for inelastic solid mechanics[END_REF] is implemented on the lattice structures for accurately estimating their effective properties. The lattices are periodic structures and they represent the macroscopic scale of the problem. The unit cell is the periodic part and henceforth, it is referred to as Representative Volume Element (RVE), which represents the microscopic scale of the problem.

For the application of the periodic boundary conditions (PBC) the mesh of the RVE should be periodic which means that for every node i on a boundary face of the RVE, another node j exists on the opposite face with the same relative position. The constraint driver method is used and a displacement gradient is applied to every pair of opposite nodes, node i, and node j. This displacement gradient is connected to the macroscopic infinitesimal strain tensor, which is related to the constraint driver. This is expressed in the following kinematic equation (Eq. ( 6)) and it is presented in Fig. 3,

u ′ i = u ′ j ⇔ u i -u j = ε * ( x i -x j ) ∀x ∈ V (6)
where i and j are the opposite nodes, u ′ i = u ′ j is the displacement perturbation, which is periodic for the two nodes, u i -u j is the displacement vector difference between the nodes, x i -x j is the position vector difference between the nodes, and ε is the macroscopic strain.

More details about this method are available in [START_REF] Praud | Fully integrated multi-scale modelling of damage and time-dependency in thermoplastic-based woven composites[END_REF].

According to Hooke's law σ kl = C klmn ε mn , the C klmn is the homogenized effective stiffness matrix and the σ kl and ε mn are the macroscopic effective stress tensor and strain tensor, respectively. Using the Voigt notation for the stress and strain tensors, as well as for the elastic stiffness matrix, one could write: σ K = C KL ε L . Six loading cases are implemented for every RVE. For every loading case, a component of the strain vector had the value 1, while the remaining five are equal to 0, according to the Еq. 7. The values of the stress vector are calculated from the reaction forces (RF) of each face in Abaqus and thus, the connection between the stiffness matrix and the stress tensor is obtained. The outcome of each loading case is one column of the stiffness matrix C KL . The complete stiffness matrix is the assembly of the six computed columns.

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ε 11 ε 22 ε 33 2ε 12 2ε 13 2ε 23 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 1 0 0 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ → ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ σ 11 σ 22 σ 33 σ 12 σ 13 σ 23 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ≡ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ C 11 C 21 C 21 0 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ε 11 ε 22 ε 33 2ε 12 2ε 13 2ε 23 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 1 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ → ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ σ 11 σ 22 σ 33 σ 12 σ 13 σ 23 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ≡ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 0 0 0 C 44 0 0 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ . (7)

Effective mechanical properties of the lattices

All the unit cell topologies, strut-based and TPMS-based, that have been investigated, with various volume fractions, present a stiffness matrix with cubic symmetry. Therefore, the elastic constitutive law can be described by three independent parameters, C 11 , C 12 and C 44 . The cubic stiffness matrix form is the following:

C ij = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ C 11 C 12 C 12 0 0 0 C 12 C 11 C 12 0 0 0 C 12 C 12 C 11 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 44 0 0 0 0 0 0 C 44 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (8)
Components C 11 and C 12 of the stiffness matrix are related to the tensile/compressive behavior of the materials while the C 44 component is related to the shear behavior. The stiffness matrix components are related to the computation of the effective Young's modulus, E ′ , effective Shear modulus, G ′ , and effective Poisson's ratio, v ′ , for the cubic materials as follows [START_REF] Lohmuller | Architectural effect on 3D elastic properties and anisotropy of cubic lattice structures[END_REF]: Considering cubic symmetry, the degree of anisotropy of the unit cells can be evaluated using an anisotropic ratio, the Zener ratio, Z [START_REF] Li | Comparison of mechanical properties and energy absorption of sheet-based and strut-based gyroid cellular structures with graded densities[END_REF][START_REF] Lohmuller | Architectural effect on 3D elastic properties and anisotropy of cubic lattice structures[END_REF], related to C 11 , C 12 and C 44 with:

E ′ = C 2 11 + C 11 C 12 -2C 2 12 C 11 + C 12 (9) G ′ = C 44 (
Z = 2C 44 C 11 -C 12 (12)
A value Z = 1 corresponds to isotropic behavior, while Z > 1 means that the predominant contribution of the C 44 component leads the structure to be suitable for shear. Conversely, Z < 1 means that the structure possesses a predominant tensile/compression behavior [START_REF] Al-Ketan | Functionally graded and multi-morphology sheet TPMS lattices: design, manufacturing, and mechanical properties[END_REF][START_REF] Lohmuller | Architectural effect on 3D elastic properties and anisotropy of cubic lattice structures[END_REF].

All the aforementioned effective mechanical properties, the Zener ratio, and local stress distribution are computed for all the investigated topologies and volume fractions. The results are presented in detail in the Results section.

Meshing the unit cell for finite elements analysis

In both cases of the TPMS-based and the strut-based lattices, the mesh of the unit cells is generated in Abaqus. Three-dimensional tetrahedral quadratic elements C3D10 are systematically adopted for the sake of comparing the final numerical results for the different types of studied unit cells. Fig. 4 shows two meshed unit cells in Abaqus, namely Gyroid sheet and Octet-truss unit cells.

A mesh size investigation is carried out for all the lattice types to ensure the numerical convergence results, their reliability, and mesh- Moreover, a mesh sensitivity analysis is carried out for the local stress distribution of the studied topologies and the convergence is reached. The mesh size does not influence on the elastic constants even though the stress values vary slightly. However, the local stress distributions present the same shape when the mesh size changes, and this is what matters in this study.

Results and comparison

This section presents and compares all the numerical results of all the lattice types. It is organized into three parts: first, the presentation of the computed effective elastic properties of all the lattices at various volume fractions followed by the study of the anisotropic behavior of the lattices. Finally, the quantitative analysis of the local stress distribution on the unit cells is presented. The combination of the aforementioned analyses shall lead to the design of a suitable lattice with tailored mechanical properties for a hip implant application. The results are discussed in the Discussion section.

Mechanical variation by topology and density

The mechanical response of various TPMS-based and strut-based topologies is investigated, namely Gyroid sheet, Gyroid skeletal, Primitive sheet, Primitive skeletal internal and external, Neovius sheet, IWP sheet, IWP skeletal internal, Diagonal, and Octet-truss. Fig. 5 illustrates the relationship between the normalized mechanical properties and the volume fraction of the studied unit cells. The elastic properties of the unit cells have been normalized to focus only on the influence of the geometry, regardless of the material properties. Some numerical results taken from the literature are also added to Fig. 5 a) to compare to and validate our results.

The direct effect of the TPMS design parameter t on all the elastic properties of the unit cells is shown in Fig. 6. This way of presenting the results accelerates the design procedure. The design parameter t is also directly linked to the C 11 and C 12 elastic components of the stiffness matrix and the effective values of the Young's modulus and shear modulus. The C 44 component is equal to G ′ as determined by Eq. ( 10)

and its relation to the volume fraction of the unit cells is presented in Fig. 6 b). By using the results above and by applying the mechanical properties of Ti-6Al-4V to the normalized values, we compute the apparent elastic modulus and the density of the unit cells [START_REF] Favre | A continuous crystallographic approach to generate cubic lattices and its effect on relative stiffness of architectured materials[END_REF]. The computed values are placed on an Ashby with reference to other materials. In Fig. 7 the Ashby chart with the titanium-based lattices is presented.

Comparison of the anisotropy of the lattices

Topology-induced anisotropy is an important feature in the effort of designing mechanically compatible implants. The Zener ratio for all the investigated structures is computed and compared to highlight the impact of one topology selection on the anisotropic response. Fig. 8 shows the Zener ratio evolution for all the investigated types of unit cells with various volume fractions.

Local stress distributions of the lattices

The local stress distribution is another important feature to take into consideration when designing architectured implants, which are required to present good high cycle fatigue. The local stress distribution in the elastic region of all the investigated unit cells is computed and compared statistically to identify the most homogenized and suitable one.

Fig. 9 shows the local stress distributions of the studied unit cells at the same volume fraction of 30%. The stress at every material point is computed for the unit cells topologies and the stress values are distributed over the range 0 -1. Axis y denotes the number of the elements that are subjected to particular stress values over the population of elements for each model (N stress /N total ). The local stress distribution is computed for all the implemented loading cases. The unit cells are triply periodic, consequently, the stress distributions are identical for the tensile loads. The curves for the shear loading cases are identical as well. Fig. 9 shows the results for a tensile loading case.

For the sake of clarity, the results are split and presented in two separate graphs, i.e. Fig. 9a) and b). The Gyroid skeletal structure is set as the reference lattice type and it is presented on both graphs to compare all the unit cells among them. Then, a statistical analysis of the results is carried out to highlight the suitable topologies by using the weighted arithmetic mean and weighted standard deviation of the sample.

Table 2 presents the weighted arithmetic mean and weighted standard deviation of the investigated normalized local stress distributions to identify the optimal structure with the most homogenized distribution. The weighted arithmetic mean is an indicator of the apparent elastic modulus of the unit cells; the higher the arithmetic mean value the higher the elastic modulus of the unit cell. Moreover, the standard deviation is an indicator of the apparent σ max of the unit cells; taking into consideration that the bulk material is the same for all the topologies, some elements of the unit cells with a higher standard deviation are subjected to higher stress values compared to those of the unit cells with a low standard deviation. This means that their σ max is closer to their σ yield . Fig. 10 shows the local stress distribution for both loading cases, i.e. tensile and shear loads, on several unit cell types. The stress distributions are different for the two loading cases per topology. Fig. 10 a),c), e) and g) present the distributions for a tensile load, while Fig. 10 b),d),f) and h) present the distributions for a shear load.

While the topologies presented in Fig. 9 have the same volume fraction but not the same apparent elastic modulus, a different comparison can be adapted among the topologies, which is portrayed in

Case study: application of the methodology to replace cortical bone according to three objectives

Fig. 12 presents the range of the elastic modulus of the cortical bone and its density as they are reported in [START_REF] Grimal | Quantitative ultrasound assessment of cortical bone properties beyond bone mineral density[END_REF]. The target range of the normalized elastic modulus is portrayed by the horizontal area and the target of density is portrayed by the vertical highlighted area to limit the topologies that satisfy these two requirements.

Then, by applying the material properties of the biocompatible Ti-6Al-4V alloy, which are E Bulk = 110 GPa and ν = 0.3, to the selected topologies and by setting the effective elastic modulus target as E ′ =

18500 MPa, we can study their local stress distribution. The local stress distributions of several unit cells with the same effective modulus are presented in Fig. 13. Table 4 summarizes the weighted arithmetic mean and weighted standard deviation of the investigated local stress distributions.

Discussion and analysis

A unit cell specific to medical applications should be generated by taking into consideration several parameters, namely the topology type, the volume fraction, and the material properties. The problem is governed by the aforementioned parameters that can be optimized. These parameters induce different effective elastic properties, anisotropic response, and local stress distributions for every unit cell topology and consequently for the implants. The following chart (Fig. 14) presents the collected constraints and objectives of the present problem. The study is carried out in the elastic regime and the parameters taken into consideration are those in the boxes with the bold outline.

Effective mechanical properties related to topology, volume fraction, and design parameter t

In section 3.1, the normalized effective elastic properties of various unit cell topologies are presented either versus the volume fraction of the unit cells as shown in Fig. 5 or versus the design parameter t as shown in Fig. 6. In Fig. 5, one can observe the relation between the topologies, both strut-based and TPMS-based, and the volume fraction of the unit cells. The evolution of the effective Young's modulus in Fig. 5 a) and shear modulus in Fig. 5 b) is different for each type of topology. The slopes of the curves for E * (Fig. 5 a)) and G * (Fig. 5 b)) do not however exhibit large differences. The ones for E * present slightly higher difference than those for G * ; Gyroid skeletal curve in Fig. 5 a) exhibits a steeper slope of the E * compared to Octet-truss. Moreover, the G * -values gap is narrower compared to the E * -values gap. The maximum G * -values gap is 0.12 approximately, while the maximum E * -value gap is 0.17, which means that a larger range of elastic modulus can be satisfied by a particular volume fraction of various topologies compared to the range of the shear modulus. Fig. 5 c) shows the computed Poisson's ratio of the unit cells and how is it affected by the topology. Some topologies, such as Gyroid skeletal and Primitive sheet, are highly sensitive to the

Table 3

The table summarizes the weighted arithmetic mean and weighted standard deviation of the investigated normalized local stress distributions. In Fig. 6, the normalized effective values for the effective Young's modulus, shear modulus, C 11 , and C 12 components are directly related to the design parameter t, without the computation of the volume fraction. These graphs can be helpful to design the suitable lattice unit cells according to effective mechanical requirements when their volume fraction is not required to be known. The design parameter t determines the thickness of the sheet TPMS-based topologies and the diameter of the skeletal TPMS-based topologies and thus it is limited by the fabrication process.

The anisotropic ratio of the unit cells

In section 3.2, the anisotropic behavior is computed and presented. Lattice structures present anisotropic behavior unlike the structures with a continuous distribution of material (Bono et al., 1999). From the evaluation of the stiffness matrix, the cubic symmetry of all the mesostructures numerically simulated is verified and the three independent elastic parameters extracted. The cubic anisotropy is further investigated using the Zener ratio and the results are illustrated in Fig. 8.

Firstly, the Zener ratio for the majority of the topologies tends to unity with the increase in the volume fraction of the unit cells. Indeed, the lattices tend to be bulk cubes when their volume fraction increases and so their behavior is getting closer to isotropy. Secondly, the Gyroid sheet and IWP sheet structures hardly present any anisotropic behavior. Their Zener ratio values are almost constant and close to one. Therefore, these two structures can be used in an application that requires isotropic behavior.

A third observation has to do with the most anisotropic unit cell type, which is the Diagonal structure. Fig. 8 demonstrates that the anisotropy of the Diagonal unit cell is highly dependent on its volume fraction. At low volume fractions, the anisotropy of the unit cell is considerably larger than at high volume fractions. A fourth observation is a high difference between the levels of anisotropy for the Diagonal and IWP skeletal internal structures, even if the two structures possess the same BCC configuration. Their topologies are similar but the IWP skeletal internal structure has smoother connections than the Diagonal one.

Local stress distribution: comparison of topologies under two loading cases (tensile and shear) and with various volume fractions

Furthermore, since fatigue resistance and the related durability are necessary for medical implants to avoid a postoperative failure, especially for young patients (Bono et al., 1999), the impact of stress concentrations is investigated in the lattices. This is investigated with the aid of the local stress distribution of the unit cells. The stress values could be compared with the yield stress of the unit cells and thus, they are related to the fatigue life. All the results have been normalized to be adapted to all material types allowing the comparison of several topologies.

The local stress distributions of various unit cells, volume fractions, and loading cases are investigated and compared. As the distributions are topology-dependent, they are different for each topology. First, investigations to check any volume-fraction-and loading-casedependency have been carried out. Then, a comparison of several topologies has been carried out to investigate which topologies present the narrowest/most homogenized distributions. The narrowness of the distributions indicates that there are neither unstressed parts of the unit cells nor overstressed ones. The unstressed parts point out that a unit cell with a lower volume fractions could be used and the overstressed parts point out stress concentrations and low fatigue resistance. Both are undesirable results.

First, Fig. 10 shows the local stress distribution for two loading cases, tensile and shear loads, for several studied unit cells. Fig. 10 a),c), e), g) exhibit the distribution for tensile load per topology with various volume fractions, while Fig. 10 b), d), f) and h) exhibit the distributions for a shear load. The most suitable structure would be the one that exhibits a homogenized local stress distribution on all the elements. The results would be all the elements to be subjected to a single stress value in both tensile and shear loading cases with SD = 0. We observe that each topology presents a different stress distribution for each loading case. Moreover, the majority of the distributions of each topology are not sensitive to the changes in the volume fraction. Primitive sheet topology is slightly sensitive to volume fraction changes for the tensile loading case and Neovius sheet topology for the shear load as it is shown in Fig. 10 c) and f) respectively. Moreover, the majority of the local stress distribution exhibit approximately a normal distribution, apart from Neovius unit cell for the tensile load, Fig. 10 e), and IWP skeletal internal unit cell for the tensile load as well in Fig. 10 g). The latter means that in the last two cases many elements are subjected to two particular ranges of stresses. Since there is no significant difference in the local stress distributions at various volume fractions, only one volume fraction is selected to study the distribution of different topologies at the same time. In Fig. 9 the investigated volume fraction is set as 30% for a tensile loading case. The Octet-truss curve exhibits a peak at very low-stress values, which means that a great number of elements are subjected to very low stresses, with values close to zero. Therefore, a large amount of material in this unit cell could be removed by performing topology optimization without modifying significantly its apparent mechanical properties by achieving a lower specific apparent elastic modulus. Neovius sheet and IWP skeletal internal unit cells exhibit a similar distribution but their elements are subjected to slightly higher stresses.

The computed weighted arithmetic mean and weighted standard deviation of the stresses are presented and compared in Table 2. The weighted arithmetic mean indicates the apparent elastic modulus of the unit cells while the weighted SD indicates the apparent σ yield of the topologies. Indeed, the topologies with high SD present a broad stress distribution with a great number of elements subjected to higher stress values than the mean stress. For the same applied mechanical load on two different topologies, some elements of the topology with the highest SD are subjected to high stresses that may be close or higher than the σ yield (critical areas). It is important the σ max always to be lower than the σ yield and thus, the most homogenized distribution is the most suitable.

The highest weighted arithmetic mean belongs to Neovius Sheet topology, which means that this unit cell at a volume fraction of 30% presents the highest apparent elastic modulus. This can also be observed in Fig. 5 a). The weighted mean stress of the Diagonal structure is the lowest and so does its apparent elastic modulus. The sheet TPMS-based topologies present higher mean stress than the skeletal TPMS-based and strut-based topologies which means that their apparent elastic modulus is also higher than the elastic modulus of the skeletal TPMS-based and strut-based topologies.

Neovius Sheet and Diagonal topologies present the highest and lowest weighted SD values, respectively. Neovius Sheet topology exhibits the broadest local stress distribution, while Diagonal structure exhibits the narrowest and the most homogenized one. This means that the density of Neovius Sheet could be optimized in terms of local stress distribution by suppressing the elements subjected to low-stress values. The majority of the elements of Diagonal topology are subjected to stresses that are closer to the weighted arithmetic mean compared to the other unit cells. Primitive sheet topology presents a slightly higher standard deviation which means that the majority of its elements are subjected to stresses that are very close to the weighted arithmetic mean as well.

Taking into consideration this first analysis with the statistical results, Diagonal topology can be considered as the most suitable one when the objective of the application is both the minimization of the elastic modulus and mass, with the drawback of presenting the highest anisotropic ratio. A deeper investigation based on further criteria, such as damage criteria, is necessary to select a topology that exhibits the best behavior in high cycle fatigue though. unit cells exhibits two peaks. The majority of their elements are subjected to two ranges of stresses. Their volume fractions differ significantly which are 32, 8% and 40% for the IWP skeletal internal and Diagonal, respectively. This is validated by the different weighted mean stress values. Diagonal structure has the narrowest local stress distribution, according to the lowest SD and it could be considered as the most suitable one when the objective is a homogenized distribution. On the contrary to the former advantages of the Diagonal topology, IWP skeletal internal presents a larger number of elements that are subjected to low-stress values and this fact allows volume and mass optimization.

Local stress distribution: comparison at an iso-effective elastic modulus

4.6. Case study: application of the methodology to replace cortical bone according to three objectives Cortical bone exhibits an elastic modulus in the range of 3 -30 GPa (Al-Tamimi et al., 2017) and a density that varies between 1.6 and 2 g* cm -3 [START_REF] Grimal | Quantitative ultrasound assessment of cortical bone properties beyond bone mineral density[END_REF]. All the objectives have to be reached by a topology to replace it successfully. Using Fig. 5 a) and the material properties of Ti-6Al-4V alloy, i.e. E = 110 GPa and ν = 0.3, we define the range of the targets as shown in Fig. 12. The first objective is the architectured structure to exhibit the same apparent elastic properties as the cortical bone, being portrayed by the horizontal highlighted area. The second objective is the density of the cortical bone to be reached, which is portrayed by the vertical highlighted area. Six topologies satisfy these objectives; the IWP sheet (red), Primitive skeletal external (black), Primitive skeletal internal (dark green), Gyroid sheet (light green), Primitive sheet (light blue), and Gyroid skeletal (dark blue).

The third objective is the selection of a lattice structure that exhibits the most homogenized normalized local stress distribution. The local stress distribution of the unit cells that satisfy the first and second objectives are presented in Fig. 13. Indeed, in Fig. 13 a particular apparent elastic modulus value is chosen to compare the stress distribution of the five topologies and it is set as E ′ = 18500 MPa. A first observation is the identical volume fractions and local stress distributions for the two types of Primitive skeletal unit cells, even though their shapes are not the same. Their volume fraction is the lowest among the volume fractions of all the studied structures. In Table 4 the weighted arithmetic mean and standard deviation are presented to analyze statistically the distributions presented in Fig. 13. Gyroid sheet topology exhibits the lowest standard deviation which means that its normalized local stress distribution is the most homogenized among all the others. There are no elements that are subjected to higher normalized stress values than 0.72 as shown in the graph. According to this first analysis, this topology could be considered as the optimal one for the replacement of the cortical bone.

Conclusions

In the present study, a comprehensive qualitative and quantitative approach is performed to numerically design, analyze, and finally select the suitable lattice type as the internal architecture of an implant. A broad database of TPMS-based and two strut-based unit cells with various volume fractions is designed for this investigation, which is carried out with the aid of finite element analysis and the application of a periodic homogenization method in the elastic regime.

The following conclusions have been reached:

• The direct relation of various normalized effective mechanical properties to the design parameter t of the TPMS equations makes straightforward the design procedure. • The most isotropic behavior is exhibited by Gyroid sheet and IWP sheet topologies. Their Zener ratios are quasi-constant and insensitive to volume fraction variation compared to other studied topologies.

• For all the studied topologies, the distributions of the normalized local stresses computed for the unit cells are not sensitive to volume fraction changes, in both tensile and shear loading cases. • According to the statistical analysis carried out, it is shown that the Diagonal structure exhibits a narrower stress distribution than the IWP skeletal internal topology at the same apparent elastic modulus. • Gyroid sheet topology satisfies several objectives of the problem presented (cortical bone replacement), such as the effective elastic properties and the density of the cortical bone and it exhibits a homogenized local stress distribution. This topology is selected as the most suitable one among the five topologies, based on the statistical analysis carried out.

The same methodology can be employed by replacing the third criterion of local stress distribution with a criterion based on anisotropy and Zener ratio evolution, which is presented in Fig. 8. The reader should keep in mind that this analysis has been performed in the elastic regime. It is then necessary but not sufficient to reach the final selection of the most suitable topology. Further criteria, such as nonlinear behavior or ductile damage criteria should be taken into consideration in the design procedure.

Fig. 1 .

 1 Fig. 1. A Diagonal unit cell (left) and an Octet-truss unit cell (right).

Fig. 2 .

 2 Fig. 2. Cubic unit cell partitioned by the a) Gyroid and b) Schwartz's Primitive TPMS.

Fig. 3 .

 3 Fig. 3. Constraint driver connection with the cubic unit cell (Praud et al., 2020).

10 )Fig. 4 .

 104 Fig. 4. Meshed Gyroid sheet and Octet-truss structures with tetrahedral quadratic elements in Abaqus software.

Fig. 5 .

 5 Fig. 5. The effective normalized a) Young's modulus, b) Shear modulus, and c) Poisson's ratio versus the volume fraction of the unit cells are presented for all the studied TPMS-based and strut-based lattices indicated in the legend.

Fig. 6 .

 6 Fig. 6. The design parameter t directly linked to the normalized effective values of the a) Young's modulus, b) Shear modulus, c) C * 11 and d) C * 12 elastic components.

Fig. 7 .Fig. 8 .

 78 Fig. 7. An Ashby chart with the Young's modulus and density of the studied topologies by applying the mechanical properties of Ti-6Al-4V (pink color) and the target are; the cortical bone (yellow color). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 9 .

 9 Fig. 9. Normalized local stress distribution of various unit cells at the same volume fraction of 30%. The results are split into two separate graphs for easier comparison.

Fig. 11 .

 11 Fig.10shows the local stress distribution for both loading cases, i.e. tensile and shear loads, on several unit cell types. The stress distributions are different for the two loading cases per topology. Fig.10 a), c), e) and g) present the distributions for a tensile load, while Fig.10 b), d), f) and h) present the distributions for a shear load.While the topologies presented in Fig.9have the same volume fraction but not the same apparent elastic modulus, a different comparison can be adapted among the topologies, which is portrayed in Fig.11. Fig. 11 presents the local stress distributions for three unit cells with the same effective elastic modulus, E ′ = 7230 Mpa, but at different volume fractions, with bulk material properties the following: E Bulk = 110 GPa and ν = 0.3. The topology of the unit cells induces different local stress distributions which play an important role in the high cycle fatigue of the lattices and the implants.Table 3 summarizes the weighted arithmetic mean and weighted standard deviation of the investigated local stress distributions.

Fig. 10 .

 10 Fig. 10. Distribution of the normalized local stresses over the range 0 -1 of a)-b) Gyroid skeletal, c)-d) Primitive sheet, e)-f) Neovius sheet and g)-h) IWP skeletal internal unit cells for both loading cases, tensile and shear loads.

Fig. 11 .

 11 Fig. 11. Normalized local stress distribution of Diagonal, IWP skeletal internal, and Octet-truss unit cells over the range 0 -1 at the same effective Young's modulus, 7230 MPa. (Iso-effective elastic modulus comparison).

  of the volume fraction. On the contrary, Neovius sheet and IWP sheet topologies exhibit a quasi-constant Poisson's ratio.

Fig. 12 .

 12 Fig. 12. Mechanical requirements to be satisfied by the topologies to replace cortical bone successfully.

Fig. 13 .

 13 Fig. 13. Normalized local stress distribution of Gyroid sheet, Gyroid skeletal, Primitive skeletal internal, Primitive skeletal external, and IWP sheet unit cells over the range 0 -1 at the same effective Young's modulus, E ′ = 18500 MPa. (Iso-effective elastic modulus comparison).

Fig. 14 .

 14 Fig. 14. Constraints, in orange, and objectives, in blue color, of the problem to be addressed. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.)

Fig. 11

 11 Fig. 11 presents a case study; the comparison of the local stress distributions of various unit cells with the same effective elastic properties. The elastic modulus objective is set as E ′ = 7230 MPa, with a maximum error of 0.3%, for the three studied topologies, which present different distributions and have different volume fractions. The IWP skeletal internal and Diagonal structures have the same BCC topology, but the Diagonal is a strut-based topology whereas IWP skeletal internal is a TPMS-based one with smoother surfaces. The stress distribution of both

Table 1

 1 Summary of several numerical studies of TPMS-based lattices.

	Authors, Year of	Investigated lattices	Aim of the study
	publication
	Lu et al. (2019)	5 skeletal TPMS-	Numerical investigation of the elastic
		based topologies	anisotropy by using a homogenization
			method
	Li et al. (2019)	Skeletal and sheet	Anisotropy, effective elastic modulus
		Gyroid	by using a homogenization method
	Al-Ketan et al.	IWP sheet, skeletal	Effect of topology on the effective
	(2018a)	internal and Diagonal	elastic modulus
	Abueidda et al.	6 sheet TPMS-based	Elastic properties and anisotropy
	(2016)	topologies
	Maskery et al.	6 skeletal TPMS-	Effect of orientation, volume fraction,
	(2018)	based topologies	and topology on the elastic modulus

Table 2

 2 The table summarizes the weighted arithmetic mean and weighted standard deviation of the investigated normalized local stress distributions.

		Weighted arithmetic mean	Weighted SD
	Gyroid sheet	0.358	0.105
	Gyroid skeletal	0.289	0.127
	Primitive sheet	0.329	0.095
	IWP sheet	0.335	0.157
	Diagonal	0.245	0.093
	Neovius sheet	0.370	0.209
	IWP skeletal Int	0.246	0.152
	Octet truss	0.260	0.154

Table 4

 4 The table summarizes the weighted arithmetic mean and weighted standard deviation of the investigated normalized local stress distributions.

		Weighted arithmetic mean	Weighted SD
	Gyroid sheet	0.368	0.111
	Gyroid skeletal	0.322	0.150
	Primitive skeletal Int	0.305	0.190
	Primitive skeletal Ext	0.306	0.190
	IWP sheet	0.347	0.160

Acknowledgment

The authors acknowledge Arts et Metiers Foundation for the financial support in the framework of ENSAM Health Network.

Declaration of competing interest

The declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. . (continued).

Appendix