Homogenization of trabecular structures

Numerical modeling of implants and specimens made from trabecular structures can be difficult and time-consuming. Trabecular structures are characterized as spatial truss structures composed of beams. A detailed discretization using the finite element method usually leads to a large number of degrees of freedom. It is attributed to the effort of creating a very fine mesh to capture the geometry of beams of the structure as accurately as possible. This contribution presents a numerical homogenization as one of the possible methods of trabecular structures modeling. The proposed approach is based on a multi-scale analysis, where the whole specimen is assumed to be homogeneous at a macro-level with assigned effective properties derived from an independent homogenization problem at a meso-level. Therein, the trabecular structure is seen as a porous or two-component medium with the metal structure and voids filled with the air or bone tissue at the meso-level. This corresponds to a twolevel finite element homogenization scheme. The specimen is discretized by a reasonable coarse mesh at the macro-level, called the macro-scale problem, while the actual microstructure represented by a periodic unit cell is discretized with sufficient accuracy, called the meso-scale problem. Such a procedure was already applied to modeling of composite materials or masonry structures. The application of this multi-scale analysis is illustrated by a numerical simulation of laboratory compression tests of trabecular specimens.


Introduction
The scientific and technical development of new materials and structures is significantly supported by additive manufacturing methods rapidly developing in recent years, which allow for manufacturing products of complex shapes and geometry. For example, 3D printed metallic trabecular structures find application in modern implantology as external layers of implants, improving mainly the osteointegration of bone tissue. The main advantages of trabecular structures are seen in the enlargement of the contact area between implants and bone tissue and in the possibility to modify mechanical properties, which helps to eliminate some problems of conventional implants, e.g., the stress shielding effect. Furthermore, this complex morphology allows the ingrowth of bone cells, thus guarantees the better stability of implants in bones [1].
Numerical modeling of specimens of complex geometry such as implants made from trabecular structures can be difficult and time-consuming. A detailed discretization using the finite element method leads very often to a large number of degrees of freedom. This is attributed to the effort of creating a very fine mesh to capture the geometry of beams of the structure as accurately as possible. One possible solution can provide the application of the homogenization method in the framework of a multi-scale analysis. In this method, the whole specimen is assumed to be homogeneous at a macro-level with assigned effective properties derived from an independent homogenization problem at a meso-level. Therein, the trabecular structure is seen as a porous or two-component medium with the metal structure and voids filled with the air or bone tissue at the meso-level. This corresponds to a two-level finite element homogenization scheme, where the specimen is discretized by a reasonable coarse mesh at the macro-level (called the macro-scale problem), while the actual microstructure in terms of a certain representative volume element (RVE) is discretized with sufficient accuracy (called the meso-scale problem) [2]. Such an RVE is usually presented in the form of a periodic unit cell (PUC) or a statistically equivalent periodic unit cell (SEPUC) constructed, such as to match the real meso-structure as close as possible [3].
The aim of this paper is a numerical homogenization of elastic properties of trabecular structures. Basic principles of numerical homogenization are briefly introduced in Section 2, and a practical application of this approach is described in following Section 3, where the main results of a simulation of laboratory compression tests are presented. Recapitulation and main benefits are summarized in Conclusions.

First order homogenization method
In the present study, the homogenized elastic properties are obtained by the strain-based first-order homogenization method, where the representative volume element (RVE) is considered as a heterogeneous periodic unit cell (PUC) loaded on its outer boundary by a displacement field compatible with a macroscopically uniform strain E in an equivalent homogeneous medium. Then, the local displacements and strains are split as are the corresponding strains, and x is the spatial coordinate [4]. The computational model is defined in terms of PUC taking the trabecular structure into account. Because the present formulation assumes the macroscopic strain E in (1) and (2) be prescribed, the solution is searched in terms of the fluctuation rather than total displacements. The stepping stone in deriving  u is the Hill lemma written as where the left-hand side represents the volume average of virtual work done by local stress and strain fields. In the framework of finite element discretization, we approximate ()  ux in terms of their nodal values r to write Substituting from (4) and (5) where where L is the material stiffness matrix. It is perhaps evident that in order to satisfy the volume average of the fluctuation part of the displacement field must disappear. This is achieved by adopting the periodic boundary conditions, which in this particular example of a rectangular PUC amounts to the same displacements  u on the opposite boundaries of PUC. The system of equations (6) is used to provide the coefficients of the effective stiffness matrix D as volume averages of the local fields derived from the solution of six consecutive elasticity problems for this 3D geometry, when the periodic unit cell is loaded, in turn, by each of the six components of E . The volume stress averages normalized with respect to E then furnish individual columns of effective elastic stiffness matrix D [2].

Numerical study
The usability and efficiency of the proposed homogenization procedure are demonstrated through the numerical simulation of laboratory compressive tests of selected trabecular structures. Trabecular structures are usually designed as spatial truss structures composed of beams. Models of the shape of Dode Thick, Diamond and Rhombic dodecahedron represent the primary three different morphologies, see Fig. 1.

Experiments
For the numerical simulation of laboratory tests, the samples were selected with basic cells in the form of Dode Thick. The samples are trabecular cubes 14x14x14 mm connected to the homogeneous parts from two opposite sides. The shapes of the homogeneous parts were adopted particularly for the compressive and tension tests. All samples are made of the titanium alloy Ti6Al4V, called Rematitan (Concept Laser, Germany), in the form of powder with the maximal granularity of 63 mm. They were produced by the additive SLM method (Selective Laser Melting). The 3D printing was performed by the machine M2 Cusing (Concept Laser, Germany) in an argon atmosphere. Additionally, the samples were subjected to a special heat treatment to eliminate inhomogeneities and internal tensions, which came from the metal powder sintering. Further details regarding the trabecular structures can be found in [1,5].

Finite element model
Because of symmetry, the numerical simulation on the macro-scale considered only oneeight of the sample. For the verification of the homogenization method, two models were created. The first model discretizes the shape of the trabecular structure almost perfectly by a very fine finite element mesh, see Fig. 3 a), where the mesh consists of 1,880,711 nodes and 8,134,568 elements. On the other hand, the second model uses two-level model. The coarse mesh discretizes the sample at the meso-level, see Fig. 3 b). Each element of the trabecular structure takes effective elastic properties derived from an RVE analysis performed at the meso-level displayed in Fig. 4. The finite element mesh of the macro-level consists of 3,461 nodes and 17,285 elements only and the meso-level mesh has 25,061 nodes and 110,427 elements if the voids are also discretized.

Results
The elastic material properties of the titanium alloy were obtained from nanoindentation to get Young's modulus E = 118 GPa assuming the Poisson's ratio n = 0.3.
However, the components of the effective stiffness matrix D show nearly an isotropic response of the representative volume element, which is also expected according to its geometry (Fig. 4). The calculated effective Young's modulus is Eeff = 3.563 GPa. Stressstrain diagrams of compression tests depicted in Fig. 5 show non-linear behavior of trabecular samples in the beginning of the loading process and mainly in post-peak parts, where the jumps and waves of softening and hardening correspond to the collapse of beams rows observed during the tests. Comparison of the elastic part of force-displacement diagrams for Dode Thick samples measurements and the numerical simulation shows Fig.  6, where the slopes of the curves are comparable just in the elastic part. The value of effective Young's modulus calculated from the measurements according to the international standard ISO 13314:2011 is E = 3.74 GPa. The overall shortening determined for the loading force F = 15 kN is summarized in Table 1. It must be noted that the computed displacements are one half of the overall shortening of the whole specimen. The computational time of the analysis is listed in Table 2., where the significant speedup is evident in the homogenization method.

Conclusions
The contribution presented an application of the two-level homogenization method of trabecular structures. In this method, the sample was discretized by a coarse mesh at the macro-level, adopting the effective elastic properties found from the solution of a mesoscale problem. The method seems to be a useful tool when modeling such structures. It showed a perfect match in results when compared to both detailed finite element simulation and experimental measurements. The main benefit of the multi-scale approach is a significant saving of computational time. The next step in the development of the method will be a simulation of new structures and materials designed, e.g., gyroid structures, gradually substituting trabecular structures [5]. Moreover, this method can be used for modelling of bone-implant interaction.
Financial support for this work was provided by GAČR project No. 18-24867S. This financial support is gratefully acknowledged.