# Topological Optimization Design for a Multiscale Femoral Prosthesis Model Based on Homogenization Method

Authors: Cheng Cheng; Ning Dai; Jie Huang; Yahong Zhuang; Tao Tang; Mengna Yang
Journal: Scientific Programming
 (2022)

---

## Abstract

Femur is a thigh bone in humans, which is particularly prone to injury and difficult to mend if damaged. The stress shielding effect is caused by a change in bone mechanics transmission route as a result of femur bone damage. The formation of a porous structure provides a better solution to the problem of stress shielding effect reduction. The porous material is basically a solid which has empty space(s) which is not covered by the primary structure of particles that build up the solid’s structure. Traditional porous materials have a uniform density distribution in their microstructure. The density distribution of a material should be established according to different loads in each location of the material to completely represent its characteristics. In this paper, to acquire the ideal density distribution of porous materials and construct high-performance variable-density porous structures, topological optimization is applied into the design of the lattice structure. According to a structure constructed utilising individualised parameters paired with computed tomography (CT) data of the human femur, a topological optimization design of a multiscale femoral prosthesis model based on the homogenization approach is presented. Experiments on different femoral prosthesis architectures revealed that variable-density porous prostheses have superior material and stress distributions than equal-density porous prostheses under the same stress and successfully decreased femoral stress shielding to improve prosthetic stability.


---

## Body

## 1. Introduction

The femur is a bone in the human body that is found in the thigh area. It is the body’s longest and the most powerful bone. It is an important component of the human capacity to stand and move. Many vital muscles, tendons, ligaments, and components of the human circulatory system are also supported by the human femur. The femur bears the majority of the human body weight as one of the primary load-bearing elements of the lower limbs. This bone is very easily injured and difficult to mend after it has been damaged. Furthermore, arthritis, osteoporosis, and other linked disorders have received more attention in aging people [1].Total hip replacement is the most effective way to solve injury and inflammatory pain of the femoral joint. However, many postoperative complications have not been solved completely. The proportion of young patients who considerably exercise and will be rebuilt after the operation has also been increasing [2]. The cause of postoperative adverse reactions is the mismatch between elastic modulus of implant and host femur. Wolff’s law states that the implant with the highest elastic modulus carries a portion of the load that is initially carried by the femur [3]. This syndrome causes bone resorption by altering the direction of mechanical transmission on the femur and lowering the load on the femur [4]. Reducing the stress shielding effect caused by change in the transmission path of bone mechanics is the problem to be solved [5].A porous substance is the one that has pores (voids) in it. The “matrix” or “frame” refers to the skeleton component of the material. A fluid (liquid or gas) is usually injected into the pores. Although the skeleton material is normally solid, structures such as foams might benefit from the notion of porous materials [6]. Along with rapid development of 3D printing technology, the emergence of a porous structure better solves the problem of reducing the stress shielding effect. The introduction of porous structure can greatly reduce the elastic modulus of a metal, which causes the elastic modulus of prosthesis to be similar to that of a natural bone. Stress shielding is reduced as a result, and the implant-bone tissue combination can be improved [7]. Node connectivity, porosity, pore width, and the overall material of unit cells influence the mechanical and biological characteristics of unit cells in completely porous biomaterials [8]. According to the type of unit body structure, porous structures for additive manufacturing can be divided into three types: imitated dot array of metal lattice, curved-surface porous structure generated by controlling the implicit function, also known as a three-period minimum curved-surface porous structure, and the unit body obtained after topological optimization [9].Numerous researches have been undertaken on developing the crystal lattice structure and its performance due to the simplicity of the crystal lattice design and its outstanding qualities, such as being light and having high strength [10, 11]. Arabnejad et al. [12] designed a hip stem filled with conceptualized two-dimensional lattice to optimize the structure with bone resorption and interface failure as conditions to control the density. The design has been proved to result in a more uniform transfer of force from the prosthesis to the femur, reduce bone resorption, and lower the fatigue strength on the stem of the prosthesis. Oladapo et al. [13] designed five kinds of hip bone implants with composite porous cell microstructure. The results showed that controllable homogenization, porosity, and particle size distribution are beneficial in increasing the cell infiltration and biological integration of hip implant composites. Seharing et al. [14] compared and analyzed the mechanical properties of cubic and eight-axis truss gradient lattice structures. The displacement is identified by using finite element analysis based on the simulation of uniaxial compression. The results showed that the cubic gradient lattice structure has the best mechanical properties under appropriate relative density and pore diameter, which is suitable for bone implant. Eldesouky et al. [15] proposed a novel design of tibial implant with a porous rhombic dodecahedron structure. Compared with the solid titanium implant, the clinical performance of such tibial-knee joint implants is improved, which lowers the stress shielding. Zhang et al. [16] proposed a step topology design of functional gradient porous biomaterial with a diamond unit structure to simulate the structure of the femoral shaft. Selective laser melting is selected by taking Ti-6Al-4V powder as the raw material. The defect coupling model predicted the Young’s modulus and yield stress of functional gradient porous biomaterials that revealed a considerable yield at the bearing location.The three-period minimum curved-surface structure is a minimum surface with complex 3D topological spatial structure [17, 18]. This structure has high specific surface area, high porosity, and high long-range ordered structure [19]. Its internal structure is interconnected, and surface is smooth. Such a porous structure is useful for constructing bone implants [20, 21]. Ma et al. [22] designed a gyroid double implicit curved-surface porous structure (double-gyroid) by using the Parametric Technology Corporation (PTC) software and gyroid. They prepared a double-gyroid structure by using laser selective melting technology. From the compression test and hydrodynamics simulation, the mechanical performance and fluid permeability are obtained. The structure is predicted to be suitable for cell culture and medical implant. Corona–Castuera et al. [23] designed and fabricated personalized stainless steel partial hip implants by using tomography data and self-supporting three-period minimum curved-surface structure. The mechanical properties of the implants during compression are adjusted by using an internal rotating element structure. The design and manufacture of implants are developed by considering the clinical conditions of specific patients. The quality of the different types of bone tissues can be adjusted to meet specific clinical requirements. Vijayavenkataraman et al. [24] studied the design of porous bone implant based on the three-period minimum curved surface, and this design is made by using photoetching ceramic manufacturing technology. A total of 12 different initial surface structure elements are considered. The results show that the selection of materials and the design of porous-based three-period minimum curved surface had led to markedly lower compression modulus of the structure than that of the natural bone. Therefore, the design method can be used in the design of bone implant to alleviate the stress shielding effect. Song et al. [25] proposed a design and optimization method for the porous structure of customized root simulated implant. The Procter and Gamble (P and G) structures with four kinds of porosity are designed and prepared by using the three-period minimum curved surface as the cubic samples. The Young’s modulus, Poisson’s ratio, and yield strength of each sample are measured by a compression test. The stress distribution at the interface between the tailored implant and the surrounding bone tissue is studied using finite element analysis under various pore structures and porosities. The results show that the porous implant constructed by using the three-period minimum curved surface can lower the stress shielding effect.The basic principle of designed porous element based on the topological optimization method is the specific load and boundary conditions as well as design and nondesign areas are set in the cubic space by using topological optimization algorithm. In order to create a unit structure with a certain porosity, the relative density is used as the weight reduction goal when creating topological optimization conditions. Simoneau et al. [26] developed a disordered porous-structure filling at the top section of the prosthesis stem, which is proven to better distribute stress and enhance mechanical distribution of the prosthesis. Deering et al. [27] designed a porous scaffold by using selective Voronoi tessellation and priority sowing technology to simulate the natural structure of trabecular bone. During polyhedron expansion, a preference texture is produced in the seed void in the original volume to change the implant’s support direction. Anisotropy is digitally characterized by the mean intercept length and star-shaped volume distribution to determine the similarity with the trabecular direction. Nicali et al. [28] built the bottom structure of the prosthesis to be more efficient than a prosthesis created using standard approaches by taking into account the patients’ masticatory biomechanics. The ultimate structural volume is 2% of the starting model, which differs significantly from the standard model. The material is distributed according to the load point, the direction and modulus of the applied force.A new topological optimization method for variable-density lattice structure is proposed in this paper. A personalized porous prosthesis that can fit the human femur is designed based on the CT data of human femur. Then, based on the homogenization method, the topological optimization design of a multiscale femoral prosthesis is achieved. The numerical analytical results of the femoral stress shielding by three prostheses are compared by combining the mechanical performance of total solid prosthesis, uniform porous prosthesis, and variable-density porous prosthesis.The rest of the paper is structured according to the following section: the topological optimization model of the multiscale lattice structure based on the homogenization method is present in Section2, Parametric modelling of femoral prosthesis is discussed in Section 3, while results and analysis are delivered in Section 4. Finally, the paper is concluded in Section 5.


## 2. Topological Optimization Model of Multiscale Lattice Structure Using Homogenization Method

This section is divided into the following sections.


### 2.1. Calculation of Mechanical Properties for a Lattice Structure

The lattice microstructure with uniform periodic array is classified as the continuum. The homogenization process might be used to produce identical mechanical characteristics. The microstructure lattice cell is a rhombic dodecahedron. The volume percent of the lattice cell are stated as follows:(1)xe=VsV,where, xe is the relative density of lattice cell (volume fraction), Vs is the volume of lattice cell’s solid component, and V is the volume of lattice cell.The equivalent elastic modulus of the lattice structureEijklH can be calculated as follows:(2)EijklH=1V∫VsEijmnMmnkldV,where Eijkl is the local elastic modulus, this is the elastic modulus for the solid component of the lattice cell’s material, and Mijkl refers to ε¯ij, the local structure tensor εij related to the macro-strains and micro-strains as follows:(3)εij=Mijklε¯ij.Mijkl can be calculated as follows:(4)Mijkl=12δikδjl+δilδjk−εij∗kl.In (4), δij is the Kronecker delta and εij∗kl is the microstructure strain corresponding to ε¯kl, the component kl of the macro-strain tensor.Micro-strain can be expressed as follows:(5)σij=EijklMklmnε¯mn,=EijklMklmnEpqmnH−1σ¯pq.To avoid the repeated use of the homogenization method and to speed up the iteration rate, this paper is aimed to calculate the equivalent mechanical properties of the lattice cells at different sampling sites by sampling the density and fit the functional relationship between the equivalent mechanical properties and the relative density of the lattice cells. The samples are collected in equal density intervals of the relative density in the range of 0–1. The corresponding equivalent mechanical properties are calculated by the homogenization method. Based on the polynomial interpolation formula, the relationship between the properties of the lattice cell and relative density is fitted into the function formula. Figure1 shows the functional relationship between the equivalent elastic modulus E¯ii, shear modulus G¯ii, Poisson’s ratio ν¯ij, and relative density of cells xe. Among these variables, Es is the elastic modulus of the solid material, while νs refers to the Poisson’s ratio of the solid materials. The cube in Figure 1 is a rhombic dodecahedron structure called a lattice structure cell. The function curve in Figure 1 is obtained by the least squares fitting. Quadratic polynomial fitting provides the elastic and shear modulus, whereas cubic polynomial fitting provides Poisson’s ratio. The parameters related to the fitting function are shown in Table 1, in which R2 represents the fitting accuracy. If R2 is closer to 1, the fitting accuracy is higher. If R2 is closer to 0, the fitting accuracy is lower. The fitting function could be directly used in the iterative process of topological optimization to avoid monotony caused by the repeated use of homogenization in each iteration and speed up the iteration.Figure 1
Equivalent mechanical properties of the lattice cell.Table 1
Functional relationship between equivalent mechanical properties and relative densityxe of lattice cells.
Equivalent mechanical propertiesFitting functionFitting accuracyR2E¯xx/Es=E¯yy/Es=E¯zz/Es0.745xe2+0.101xe+0.0220.998 57G¯xy/Es=G¯yz/Es=G¯xz/Es0.515xe2−0.098xe+0.0310.997 79ν¯xy/νs=ν¯yz/νs=ν¯xz/νs0.449xe3−0.644xe2+0.292xe+0.7060.989 65


### 2.2. Multiscale Topological Optimization

The homogenization approach is used to propose the density method. Essentially, the classic density technique generates a distribution of 0 and 1, which is related to macrostructure optimization. The topological optimization design of the lattice structure included the microscale, which involved the qualities of a microstructure, in addition to the macroscale.The lattice cell corresponds to the finite element one by one in the iteration of the density method. The equivalent elastic modulusEe (i.e., unit elastic modulus) of the lattice structure cell is calculated rapidly according to the fitting formula in Table 1, and then, the numerical solution is performed as follows:(6)x=x1,x2,…,xnT,cx=UTKU=∑e=1NEexeueTkeue,s.t.VxV0=f,F=KU,0<xmin≤xe≤1,where, xe refers to the design variable (unit relative density), and its value is in the range 0,1,xmin refers to the minimum relative density to avoid occurrence of singularity, x refers to the vector of design variables, N refers to the number of design variables, c refers to the overall flexibility, Ee refers to the unit elastic modulus (equivalent to the elastic modulus of the cell), U refers to the global displacement matrix, F refers to the global stress matrix, K refers to the global stiffness matrix, ke refers to the unit stiffness matrix, ue refers to the unit displacement matrix, Vx and V0 refer to the solid and total volumes of the design domain, respectively, and f refers to the volume fraction.Different methods are used to solve the optimization results. The optimization criterion algorithm is common. The design variables could be iterated constantly in the solution process. The iterative process depended on the following heuristic algorithm(7)xenew=maxxmin,xe−m,ifxeBeη≤maxxmin,xe−m,xeBeη,ifmaxxmin,xe−m<xeBeη<min1,xe+m,min1,xe+m,ifxeBeη≥min1,xe+m,where, m represents the maximum variation; and η refers to the numerical damping coefficient. Be could be solved as follows:(8)Be=−∂c/∂xeλ∂V/∂xe,where λ refers to the Lagrangian function.Contrary to the traditional density method, the new technique does not have a penalty factor. The sensitivity of the objective function is derived on the basis of the fitting formula in Table1 as follows:(9)∂c∂xe=1.098xe+0.102ueTkeue.When the unit volume is used in each element, the sensitivity of the volume is as follows:(10)∂V∂xe=1.To summarise, the above represented the essential principle of topological optimization of the lattice structure using the homogenization method: The fitting function of comparable mechanical characteristics of the rhombic dodecahedron lattice cell is calculated and fitted using the homogenization approach. By merging the density approach with the structure topological optimization model of the minimal compliance issue under volume constraints, a structural topological optimization model of the minimum compliance problem under volume constraints is created. The unit mechanical characteristics of the new technique are derived using the homogenization method, rather of the previous density method. There is no penalty factor employed. The fitting function equation is used to calculate sensitivity. As a result, the optimization results comprise a series of gray units with the relative densities ranging from 0 to 1. A variable-density lattice structure would eventually emerge from these gray components.


### 2.3. Parametric Modelling of Variable-Density Lattice Structure

To construct the lattice structural morphology based on discrete elements, extracting the information of key position is necessary to assist the modelling. As shown in Figure2, the spatial position of each lattice cell could be positioned by extracting the node coordinates of the element. In terms of equal-density lattice structure, array modelling could be directly performed according to the geometric parameters of the cell, as shown in Figure 2(b). If the node position coordinates are retrieved without addressing the transition between cells with varying densities for a variable-density lattice structure, the problem of unequal connectivity between cells will arise. In the case of the X-shaped cell, the connection’s strength is decided by the weakest cell at the time shown in Figure 2(c). To solve a nonsmooth connection, extracting the adjacency of each unit and the connection between unit nodes, in addition to the node information, is required to transform the unit density into the node density and construct the smooth variable-density lattice model based on the node density shown in Figure 2(d).Figure 2
Example of the lattice model based on the discrete element: (a) unit node model; (b) equal-density lattice model; (c) nonsmooth variable-density lattice model; (d) smooth variable-density lattice model.
(a)(b)(c)(d)The autogeneration of the geometric model for the variable-density lattice mainly included four parts as shown in Figure3: (1) acquisition of unit density distribution after optimization; (2) calculation of the node density; (3) construction of the cell geometric model, and (4) generation of the lattice geometric model.Figure 3
Schematic diagram of the variable-density lattice model.After optimization, the density of the unit density can be retrieved from the topological optimization result. The following weighted average equation is used to determine node density:(11)ρnj=∑ixeiVei∑iVei,where, ρnj refers to density of the jth node, xei and Vei refer to density and volume of the j element adjacent to node i, respectively.Nonuniform rational B-spline is used to describe the model in the existing mainstream computer-aided design (CAD) system. Its essence is an element expressed in the form of the tensor product. As expressing models with complex geometry, this element is usually expressed topologically by boundary representation (B-Rep) and format of the constructed solid geometry (CSG).The geometric model of the cell is constructed based on the B-Rep method. The geometric and topological information of the model are constructed parametrically according to the spot, line, and surface from bottom to top. The expression of B-Rep is depicted in Figure4 using the X-shaped cell in Figure 3 as an example. The parametric construction and operation of various geometric models based on the spot, line, and surface can be achieved from bottom to top according to the topological information of the shape in the B-Rep structure once the B-Rep structure is established.Figure 4
Schematic diagram of the B-Rep expression in the X-shaped cell.For the variable-density lattice structure, adjusting the coordinate position on the top of the cell according to the node density is necessary. Then, the geometric model of the deformed cell is constructed by constructing the line and surface based on the B-Rep relation. After the deformed cells corresponding to various elements are constructed, the parametric model of complex lattice structure could be constructed by CSG through intersection, union, difference, and other operations. In CAD system, geometric models are expressed by boundaries, while 3-dimensional (3D) solid models are described through their surfaces. Therefore, the above methods can be used in the surface parametric construction of 3D models to realize the autogeneration of complex 3D variable-density lattice models.


### 2.4. Algorithm Implementation

The basic flow of topological optimization algorithm for multiscale lattice structure based on the homogenization method is shown in Figure5.Figure 5
Basic flowchart.The proposed method is based on the homogenization and density methods. First, equivalent mechanical properties of the lattice cell are determined by using the homogenization method and fitted as functions with relative density of the lattice cell. Then, using the density approach in MATLAB, iteration is optimized using the relative density of microstructure cells as the design value and the minimal flexibility as the goal function. The fitting formula in Table1 is used to compute the equivalent mechanical characteristics of cells. The equivalent volume element approach is used to equate the lattice cells to dense solid components. The finite element problem is solved via Analysis System (ANSYS). The relative density distribution of the element, element nodes, and other information are produced after interactive MATLAB and ANSYS calculation. Finally, the relative density is used to construct a one-to-one connection between the element and the lattice structure cell. Thus, the variable-density lattice structure is parameterized and generated based on the element node information.


## 2.1. Calculation of Mechanical Properties for a Lattice Structure

The lattice microstructure with uniform periodic array is classified as the continuum. The homogenization process might be used to produce identical mechanical characteristics. The microstructure lattice cell is a rhombic dodecahedron. The volume percent of the lattice cell are stated as follows:(1)xe=VsV,where, xe is the relative density of lattice cell (volume fraction), Vs is the volume of lattice cell’s solid component, and V is the volume of lattice cell.The equivalent elastic modulus of the lattice structureEijklH can be calculated as follows:(2)EijklH=1V∫VsEijmnMmnkldV,where Eijkl is the local elastic modulus, this is the elastic modulus for the solid component of the lattice cell’s material, and Mijkl refers to ε¯ij, the local structure tensor εij related to the macro-strains and micro-strains as follows:(3)εij=Mijklε¯ij.Mijkl can be calculated as follows:(4)Mijkl=12δikδjl+δilδjk−εij∗kl.In (4), δij is the Kronecker delta and εij∗kl is the microstructure strain corresponding to ε¯kl, the component kl of the macro-strain tensor.Micro-strain can be expressed as follows:(5)σij=EijklMklmnε¯mn,=EijklMklmnEpqmnH−1σ¯pq.To avoid the repeated use of the homogenization method and to speed up the iteration rate, this paper is aimed to calculate the equivalent mechanical properties of the lattice cells at different sampling sites by sampling the density and fit the functional relationship between the equivalent mechanical properties and the relative density of the lattice cells. The samples are collected in equal density intervals of the relative density in the range of 0–1. The corresponding equivalent mechanical properties are calculated by the homogenization method. Based on the polynomial interpolation formula, the relationship between the properties of the lattice cell and relative density is fitted into the function formula. Figure1 shows the functional relationship between the equivalent elastic modulus E¯ii, shear modulus G¯ii, Poisson’s ratio ν¯ij, and relative density of cells xe. Among these variables, Es is the elastic modulus of the solid material, while νs refers to the Poisson’s ratio of the solid materials. The cube in Figure 1 is a rhombic dodecahedron structure called a lattice structure cell. The function curve in Figure 1 is obtained by the least squares fitting. Quadratic polynomial fitting provides the elastic and shear modulus, whereas cubic polynomial fitting provides Poisson’s ratio. The parameters related to the fitting function are shown in Table 1, in which R2 represents the fitting accuracy. If R2 is closer to 1, the fitting accuracy is higher. If R2 is closer to 0, the fitting accuracy is lower. The fitting function could be directly used in the iterative process of topological optimization to avoid monotony caused by the repeated use of homogenization in each iteration and speed up the iteration.Figure 1
Equivalent mechanical properties of the lattice cell.Table 1
Functional relationship between equivalent mechanical properties and relative densityxe of lattice cells.
Equivalent mechanical propertiesFitting functionFitting accuracyR2E¯xx/Es=E¯yy/Es=E¯zz/Es0.745xe2+0.101xe+0.0220.998 57G¯xy/Es=G¯yz/Es=G¯xz/Es0.515xe2−0.098xe+0.0310.997 79ν¯xy/νs=ν¯yz/νs=ν¯xz/νs0.449xe3−0.644xe2+0.292xe+0.7060.989 65


## 2.2. Multiscale Topological Optimization

The homogenization approach is used to propose the density method. Essentially, the classic density technique generates a distribution of 0 and 1, which is related to macrostructure optimization. The topological optimization design of the lattice structure included the microscale, which involved the qualities of a microstructure, in addition to the macroscale.The lattice cell corresponds to the finite element one by one in the iteration of the density method. The equivalent elastic modulusEe (i.e., unit elastic modulus) of the lattice structure cell is calculated rapidly according to the fitting formula in Table 1, and then, the numerical solution is performed as follows:(6)x=x1,x2,…,xnT,cx=UTKU=∑e=1NEexeueTkeue,s.t.VxV0=f,F=KU,0<xmin≤xe≤1,where, xe refers to the design variable (unit relative density), and its value is in the range 0,1,xmin refers to the minimum relative density to avoid occurrence of singularity, x refers to the vector of design variables, N refers to the number of design variables, c refers to the overall flexibility, Ee refers to the unit elastic modulus (equivalent to the elastic modulus of the cell), U refers to the global displacement matrix, F refers to the global stress matrix, K refers to the global stiffness matrix, ke refers to the unit stiffness matrix, ue refers to the unit displacement matrix, Vx and V0 refer to the solid and total volumes of the design domain, respectively, and f refers to the volume fraction.Different methods are used to solve the optimization results. The optimization criterion algorithm is common. The design variables could be iterated constantly in the solution process. The iterative process depended on the following heuristic algorithm(7)xenew=maxxmin,xe−m,ifxeBeη≤maxxmin,xe−m,xeBeη,ifmaxxmin,xe−m<xeBeη<min1,xe+m,min1,xe+m,ifxeBeη≥min1,xe+m,where, m represents the maximum variation; and η refers to the numerical damping coefficient. Be could be solved as follows:(8)Be=−∂c/∂xeλ∂V/∂xe,where λ refers to the Lagrangian function.Contrary to the traditional density method, the new technique does not have a penalty factor. The sensitivity of the objective function is derived on the basis of the fitting formula in Table1 as follows:(9)∂c∂xe=1.098xe+0.102ueTkeue.When the unit volume is used in each element, the sensitivity of the volume is as follows:(10)∂V∂xe=1.To summarise, the above represented the essential principle of topological optimization of the lattice structure using the homogenization method: The fitting function of comparable mechanical characteristics of the rhombic dodecahedron lattice cell is calculated and fitted using the homogenization approach. By merging the density approach with the structure topological optimization model of the minimal compliance issue under volume constraints, a structural topological optimization model of the minimum compliance problem under volume constraints is created. The unit mechanical characteristics of the new technique are derived using the homogenization method, rather of the previous density method. There is no penalty factor employed. The fitting function equation is used to calculate sensitivity. As a result, the optimization results comprise a series of gray units with the relative densities ranging from 0 to 1. A variable-density lattice structure would eventually emerge from these gray components.


## 2.3. Parametric Modelling of Variable-Density Lattice Structure

To construct the lattice structural morphology based on discrete elements, extracting the information of key position is necessary to assist the modelling. As shown in Figure2, the spatial position of each lattice cell could be positioned by extracting the node coordinates of the element. In terms of equal-density lattice structure, array modelling could be directly performed according to the geometric parameters of the cell, as shown in Figure 2(b). If the node position coordinates are retrieved without addressing the transition between cells with varying densities for a variable-density lattice structure, the problem of unequal connectivity between cells will arise. In the case of the X-shaped cell, the connection’s strength is decided by the weakest cell at the time shown in Figure 2(c). To solve a nonsmooth connection, extracting the adjacency of each unit and the connection between unit nodes, in addition to the node information, is required to transform the unit density into the node density and construct the smooth variable-density lattice model based on the node density shown in Figure 2(d).Figure 2
Example of the lattice model based on the discrete element: (a) unit node model; (b) equal-density lattice model; (c) nonsmooth variable-density lattice model; (d) smooth variable-density lattice model.
(a)(b)(c)(d)The autogeneration of the geometric model for the variable-density lattice mainly included four parts as shown in Figure3: (1) acquisition of unit density distribution after optimization; (2) calculation of the node density; (3) construction of the cell geometric model, and (4) generation of the lattice geometric model.Figure 3
Schematic diagram of the variable-density lattice model.After optimization, the density of the unit density can be retrieved from the topological optimization result. The following weighted average equation is used to determine node density:(11)ρnj=∑ixeiVei∑iVei,where, ρnj refers to density of the jth node, xei and Vei refer to density and volume of the j element adjacent to node i, respectively.Nonuniform rational B-spline is used to describe the model in the existing mainstream computer-aided design (CAD) system. Its essence is an element expressed in the form of the tensor product. As expressing models with complex geometry, this element is usually expressed topologically by boundary representation (B-Rep) and format of the constructed solid geometry (CSG).The geometric model of the cell is constructed based on the B-Rep method. The geometric and topological information of the model are constructed parametrically according to the spot, line, and surface from bottom to top. The expression of B-Rep is depicted in Figure4 using the X-shaped cell in Figure 3 as an example. The parametric construction and operation of various geometric models based on the spot, line, and surface can be achieved from bottom to top according to the topological information of the shape in the B-Rep structure once the B-Rep structure is established.Figure 4
Schematic diagram of the B-Rep expression in the X-shaped cell.For the variable-density lattice structure, adjusting the coordinate position on the top of the cell according to the node density is necessary. Then, the geometric model of the deformed cell is constructed by constructing the line and surface based on the B-Rep relation. After the deformed cells corresponding to various elements are constructed, the parametric model of complex lattice structure could be constructed by CSG through intersection, union, difference, and other operations. In CAD system, geometric models are expressed by boundaries, while 3-dimensional (3D) solid models are described through their surfaces. Therefore, the above methods can be used in the surface parametric construction of 3D models to realize the autogeneration of complex 3D variable-density lattice models.


## 2.4. Algorithm Implementation

The basic flow of topological optimization algorithm for multiscale lattice structure based on the homogenization method is shown in Figure5.Figure 5
Basic flowchart.The proposed method is based on the homogenization and density methods. First, equivalent mechanical properties of the lattice cell are determined by using the homogenization method and fitted as functions with relative density of the lattice cell. Then, using the density approach in MATLAB, iteration is optimized using the relative density of microstructure cells as the design value and the minimal flexibility as the goal function. The fitting formula in Table1 is used to compute the equivalent mechanical characteristics of cells. The equivalent volume element approach is used to equate the lattice cells to dense solid components. The finite element problem is solved via Analysis System (ANSYS). The relative density distribution of the element, element nodes, and other information are produced after interactive MATLAB and ANSYS calculation. Finally, the relative density is used to construct a one-to-one connection between the element and the lattice structure cell. Thus, the variable-density lattice structure is parameterized and generated based on the element node information.


## 3. Parametric Modelling of Femoral Prosthesis

The ilium and femur of participants were scanned using CT for the model used in this experiment. The elastic modulus distribution in the appropriate region can be determined using the corresponding relationship equation, which can be reflected in the CT data. In the human body, the density of blood is approximately equal to that of water, which isρwater=1.0g/cm3, the corresponding apparent density is ρAPP_0=0, and gray value is 0. The bone density in the area with the largest bone density is ρcort, and the gray value is 1613. Given that the hardest compact bone contains almost no liquid, its apparent density can be considered equal to the actual density. Therefore, the linear relationship between its apparent density ρAPP and gray value Ht is as follows:(12)ρAPP=Htρcort1613.The cortical bone is greatly different from the cancellous bone in density. The mapping relation of the elastic moduli could be expressed as piece-wise functions as follows:Cortical bone:(13)ECortical=15010ρAPP2.18ν=0.3ρAPP>0.28.Cancellous bone:(14)ECancellous=6850ρAPP1.49ν=0.3ρAPP≤0.28,where, ECortical and ECancellous refer to the elastic moduli of the cortical and cancellous bones, respectively, and ν represents Poisson’s ratio. (11) and (12) are utilised to get the distribution of elastic moduli in the femur, which is then used as the foundation for further microstructure design, modelling, and simulation experiments.The neck plane of the femur (T20), lesser trochanter plane of femur (T0), and isthmus plane of femur (TN) at the top of the femur’s CT image is shown in Figure 6. These planes could reflect the data of the patient’s femur and are suitable for designing prosthesis. The long and short diameters of the pulp cavity for different sections are measured and obtained on the abovementioned three planes and are recorded as LL20, LB20, LL0, LB0, LLN, and LBN. The distance from the neck plane of femur T20 to the isthmus plane of femur TN is measured and recorded as L. The frontal projection data of CT data is used to pick CT data that might convey the morphology of the femur for the basis of individualised modelling. A section perpendicular to the center line of the neck of the femur is established at the intersection between the facial femur T20 and pulp cavity. The distance from the vertical point to center of the femur’s head is measured and recorded as LD. The angle between the central line of the neck and the central line of pulp cavity is measured and recorded as α. The radius of the head of femur is measured and recorded as Rc. For the prosthetic design, this data would serve as the primary reference sizes.Figure 6
Extraction of the feature parameters of the femoral microstructure: (a)TN face; (b) T0 face; (c) T20 face; (d) axial size of the femur.
(a)(b)(c)(d)The parametric modelling system of the femoral prosthesis is built based on the open-cascade open-source modelling platform in this work to better adapt the prosthesis shape of the femur. The parametric modelling of the prosthesis is accomplished utilising the software based on the structural characteristics of the femur measured. Simultaneously, CT scans are utilised to simulate the ilium and femur in contact with the femoral prosthesis in order to verify and investigate the stress on the femoral prosthesis. Figure7 depicts the completed model.Figure 7
An example of the parametric modelling: (a) cortical bone of the femur; (b) cancellous bone of the femur; (c) femoral stem prosthesis; (d) femoral head prosthesis; (e) acetabular cartilage; (f) cancellous bone of ilium; (g) cortical bone of ilium; (h) assembled model.


## 4. Results and Analysis

The results and analysis of the proposed work consist of the following parts:


### 4.1. Mechanical Model Establishment

In this experiment, the weight delivered to the femur is the same as when a human walks normally. The peak approach is used to model the walking process of an 80 kg human, which is broken down into three distinct moments: early (toe off), middle (single-leg support), and late (heel striking). In each typical instant, a load is applied, as shown in Tables2, 3, and 4, correspondingly. A person’s average walking pace is 4 km/h, with a maximum peak value of 3.8 times their body weight. Additionally, complicated miscellaneous loads with many functions are delivered to the upper femur during this phase, which are summarised in this experiment as joint force, abductor force, lateral femur muscle strength, and iliopsoas muscle strength. Figure 8 illustrates the loading technique.Table 2
Joint force and muscle strength on the upper femur during walking at the early stage.
DirectionJoint force (N)Abductor force (N)Lateral femur muscle strength (N)Iliopsoas muscle strength (N)X1160−5040−17.6Y33600−131.2Z2800−13521400−123.2Table 3
Joint force and muscle strength on the upper femur during walking at the middle stage.
DirectionJoint force (N)Abductor force (N)Lateral femur muscle strength (N)Iliopsoas muscle strength (N)X290−5040−17.6Y8000−131.2Z700−13521120−123.2Table 4
Joint force and muscle strength on the upper femur during walking at the late stage.
DirectionJoint force (N)Abductor force (N)Lateral femur muscle strength (N)Iliopsoas muscle strength (N)X290−1000−88Y8000−565Z700−254700−616Figure 8
loading of the femur:F1r, joint force; F2, abductor force; F3, lateral femur muscle strength; F4, iliopsoas muscle strength.


### 4.2. Optimization Analysis

The femoral prosthesis is analyzed using finite element analysis. Then, the prosthesis is optimized topologically by using the topological optimization model of the multiscale variable-density lattice structure based on the homogenization method. The prosthetic material selected is Ti-6Al-4V alloy, with elastic modulus of 110 gigapascals (GPa), Poisson’s rate of 0.3, and density of 4.5 × 10−6 kg/mm3.The load is delivered to the node in this experiment, and the force on the femoral prosthesis is uniformly distributed to the place where the femur is positioned. Changing the node location subject to the load adjusts the position of the load in the prosthesis. There are a lot of burdens to deal with. The computation must converge once the load is applied, which necessitates the creation of three analytical phases. The peak load on the top of the femur in the early stages of walking is applied in the load module in the first analytical phase, as indicated in Table2. The peak load at the middle stage of walking is applied in second analysis as shown in Table 3. The load on the upper part of the femoral prosthesis in the late stage of walking is applied in the third analytical step as shown in Table 4. At the same time, fixed constraints are applied to the lower part of the parametric model of the upper femur.The prosthesis of the femoral head is the mainly stressed part, the whole solid structure is adopted. The optimization analysis of the femoral stem prosthesis is studied. After optimization, the finite element analysis is used to examine the complete solid femoral prosthesis model, uniform porous femoral prosthesis model, and the variable-density porous femoral prosthesis model. Figure9 shows the analytical results of the global deformation and equivalent stress for all-metal femoral stem prosthesis. Figure 10 shows the uniform porous femoral prosthesis model (porosity, 50%) as well as analytical results of the global deformation and equivalent stress. Figure 11 shows the optimized variable-density porous femoral prosthesis model and the analytical results of the global deformation and equivalent stress.Figure 9
Numerical analysis of the all-metal femoral stem prosthesis: cloud charts of the (a) equivalent stress and (b) global deformation.
(a)(b)Figure 10
Numerical analysis of the uniform porous femoral prosthesis: (a) uniform prosthesis model and the cloud charts of the (b) equivalent stress, and (c) global deformation.
(a)(b)(c)Figure 11
Numerical analysis of the variable-density porous femoral prosthesis: (a) variable-density prosthesis model and the cloud charts of the (b) equivalent stress and (c) global deformation.
(a)(b)(c)As shown in Figure9(a), the solid femoral stem prosthesis weighed approximately 151.9 g. The cloud chart of the equivalent stress showed a certain stress concentration at the neck of the femur, and the maximum stress is approximately 117.9 megapascal (MPa). The axial direction of the femoral prosthesis transmits the stress, while the stress on the lower surface is significant. This phenomenon is due to the fixed support constraint in this study. The cloud chart of the solid structure deformation is shown in Figure 9(b). The maximum deformation is located near the neck of the femur, and the maximum value is approximately 4.98 mm. The deformation gradually decreased towards the tail of the femoral stem prosthesis. The whole deformation tended to rotate around the tail of the prosthesis.The uniform porous structure treatment is performed below the neck of the femur by using the rhombic dodecahedral element with porosity of 50% as shown in Figure10(a) and weight of approximately 75.74 g, which is ∼50.14% lower than the whole solid structure. The cloud chart of the equivalent stress as shown in Figure 10(b), indicated a certain stress concentration at the neck of the femur and the transition between the solid and porous parts. The normal area of the load in the prosthesis decreased, and the stress increased accordingly. The maximum stress is approximately 459.75 MPa, considerably more than that of the whole solid femoral stem prosthesis. However, the stress concentration slightly improved. The stress concentration at the neck of the femur had slowed down, while the stress concentration at the lower end of the femoral prosthesis has significantly decreased. The cloud chart of the uniform density structure deformation is shown in Figure 10(c). With a maximum value of 20.87 mm, the largest distortion is seen around the femur’s neck. The greatest deformation achieved by employing the porous structure is 4.19 times more than the total deformation of the solid structure. The distortion is mostly located above the femur’s neck. From the transition between solid and porous portions to the tail, the deformation is minimal.The femoral prosthesis is optimized topologically by using the topological optimization of multiscale variable-density lattice structure based on the homogenization method. The optimized femoral prosthesis structure is shown in Figure11(a). The structure weighed ∼87.23 g, which is ∼42.57% lower than the whole solid structure material and ∼15.17% more than the uniform porous structure material. The stress condition is shown in Figure 11(b). The cloud chart of the equivalent stress showed that the variable-density porous structure is similar to the whole solid structure. The maximum stress is ∼150.44 MPa, which is ∼27.6% more than that of the whole solid structure and ∼67.28% lower than that of uniform porous structure. The stress concentration is significantly improved, and the stress concentration in the neck of the femur is slightly slowed down. The cloud chart of the variable-density porous structure deformation is shown in Figure 11(c). Similarly, the maximum deformation is located near the neck of the femur, and the maximum value is ∼6.61 mm. The variable-density porous structure is adopted. Its deformation is similar to that of the whole solid structure, which is ∼32.73% more than the whole solid structure deformation and ∼68.33% lower than that of uniform porous structure. The deformation gradually decreased towards the tail of the femoral stem prosthesis. The whole deformation tended to rotate around the tail of the prosthesis.The region most prone to bone resorption in hip prosthesis replacement is shown in Figure12. This study is performed to analyze the bone resorption in the feature region. To better evaluate the bone resorption, the stress shielding rate is used as a measurement of bone resorption. Stress shielding rate is calculated as follows:(15)ψ=1−σpσo×%,where ψ refers to the stress shielding rate; σp refers to the stress after prosthesis implantation; and σo refers to the stress before prosthesis implantation.Figure 12
Featured regions of the femur: (a) regions prone to bone resorption and (b) stress distribution.
(a)(b)The mean stresses in the highlighted region are 6.82, 19.49, and 16.25 MPa for the complete solid prosthesis, uniform porous prosthesis, and variable-density porous prosthesis, respectively. Figure12 depicts the stress distributions in the feature areas of the three prostheses. When compared to a solid construction, the porous structure significantly improves stress shielding on the femur. The gravitational shielding caused by all-metal and porous stems differs significantly. This significant difference represented the fact that low-stretch implants are superior than high-stretch implants.Topology optimization of the multiscale lattice structure proposed in this paper has achieved variable-density topological optimization. This approach produced not only an excellent lightweight effect, but also a superior stress distribution, demonstrating the lattice structure’s capability and obviating the need for stress shielding.


## 4.1. Mechanical Model Establishment

In this experiment, the weight delivered to the femur is the same as when a human walks normally. The peak approach is used to model the walking process of an 80 kg human, which is broken down into three distinct moments: early (toe off), middle (single-leg support), and late (heel striking). In each typical instant, a load is applied, as shown in Tables2, 3, and 4, correspondingly. A person’s average walking pace is 4 km/h, with a maximum peak value of 3.8 times their body weight. Additionally, complicated miscellaneous loads with many functions are delivered to the upper femur during this phase, which are summarised in this experiment as joint force, abductor force, lateral femur muscle strength, and iliopsoas muscle strength. Figure 8 illustrates the loading technique.Table 2
Joint force and muscle strength on the upper femur during walking at the early stage.
DirectionJoint force (N)Abductor force (N)Lateral femur muscle strength (N)Iliopsoas muscle strength (N)X1160−5040−17.6Y33600−131.2Z2800−13521400−123.2Table 3
Joint force and muscle strength on the upper femur during walking at the middle stage.
DirectionJoint force (N)Abductor force (N)Lateral femur muscle strength (N)Iliopsoas muscle strength (N)X290−5040−17.6Y8000−131.2Z700−13521120−123.2Table 4
Joint force and muscle strength on the upper femur during walking at the late stage.
DirectionJoint force (N)Abductor force (N)Lateral femur muscle strength (N)Iliopsoas muscle strength (N)X290−1000−88Y8000−565Z700−254700−616Figure 8
loading of the femur:F1r, joint force; F2, abductor force; F3, lateral femur muscle strength; F4, iliopsoas muscle strength.


## 4.2. Optimization Analysis

The femoral prosthesis is analyzed using finite element analysis. Then, the prosthesis is optimized topologically by using the topological optimization model of the multiscale variable-density lattice structure based on the homogenization method. The prosthetic material selected is Ti-6Al-4V alloy, with elastic modulus of 110 gigapascals (GPa), Poisson’s rate of 0.3, and density of 4.5 × 10−6 kg/mm3.The load is delivered to the node in this experiment, and the force on the femoral prosthesis is uniformly distributed to the place where the femur is positioned. Changing the node location subject to the load adjusts the position of the load in the prosthesis. There are a lot of burdens to deal with. The computation must converge once the load is applied, which necessitates the creation of three analytical phases. The peak load on the top of the femur in the early stages of walking is applied in the load module in the first analytical phase, as indicated in Table2. The peak load at the middle stage of walking is applied in second analysis as shown in Table 3. The load on the upper part of the femoral prosthesis in the late stage of walking is applied in the third analytical step as shown in Table 4. At the same time, fixed constraints are applied to the lower part of the parametric model of the upper femur.The prosthesis of the femoral head is the mainly stressed part, the whole solid structure is adopted. The optimization analysis of the femoral stem prosthesis is studied. After optimization, the finite element analysis is used to examine the complete solid femoral prosthesis model, uniform porous femoral prosthesis model, and the variable-density porous femoral prosthesis model. Figure9 shows the analytical results of the global deformation and equivalent stress for all-metal femoral stem prosthesis. Figure 10 shows the uniform porous femoral prosthesis model (porosity, 50%) as well as analytical results of the global deformation and equivalent stress. Figure 11 shows the optimized variable-density porous femoral prosthesis model and the analytical results of the global deformation and equivalent stress.Figure 9
Numerical analysis of the all-metal femoral stem prosthesis: cloud charts of the (a) equivalent stress and (b) global deformation.
(a)(b)Figure 10
Numerical analysis of the uniform porous femoral prosthesis: (a) uniform prosthesis model and the cloud charts of the (b) equivalent stress, and (c) global deformation.
(a)(b)(c)Figure 11
Numerical analysis of the variable-density porous femoral prosthesis: (a) variable-density prosthesis model and the cloud charts of the (b) equivalent stress and (c) global deformation.
(a)(b)(c)As shown in Figure9(a), the solid femoral stem prosthesis weighed approximately 151.9 g. The cloud chart of the equivalent stress showed a certain stress concentration at the neck of the femur, and the maximum stress is approximately 117.9 megapascal (MPa). The axial direction of the femoral prosthesis transmits the stress, while the stress on the lower surface is significant. This phenomenon is due to the fixed support constraint in this study. The cloud chart of the solid structure deformation is shown in Figure 9(b). The maximum deformation is located near the neck of the femur, and the maximum value is approximately 4.98 mm. The deformation gradually decreased towards the tail of the femoral stem prosthesis. The whole deformation tended to rotate around the tail of the prosthesis.The uniform porous structure treatment is performed below the neck of the femur by using the rhombic dodecahedral element with porosity of 50% as shown in Figure10(a) and weight of approximately 75.74 g, which is ∼50.14% lower than the whole solid structure. The cloud chart of the equivalent stress as shown in Figure 10(b), indicated a certain stress concentration at the neck of the femur and the transition between the solid and porous parts. The normal area of the load in the prosthesis decreased, and the stress increased accordingly. The maximum stress is approximately 459.75 MPa, considerably more than that of the whole solid femoral stem prosthesis. However, the stress concentration slightly improved. The stress concentration at the neck of the femur had slowed down, while the stress concentration at the lower end of the femoral prosthesis has significantly decreased. The cloud chart of the uniform density structure deformation is shown in Figure 10(c). With a maximum value of 20.87 mm, the largest distortion is seen around the femur’s neck. The greatest deformation achieved by employing the porous structure is 4.19 times more than the total deformation of the solid structure. The distortion is mostly located above the femur’s neck. From the transition between solid and porous portions to the tail, the deformation is minimal.The femoral prosthesis is optimized topologically by using the topological optimization of multiscale variable-density lattice structure based on the homogenization method. The optimized femoral prosthesis structure is shown in Figure11(a). The structure weighed ∼87.23 g, which is ∼42.57% lower than the whole solid structure material and ∼15.17% more than the uniform porous structure material. The stress condition is shown in Figure 11(b). The cloud chart of the equivalent stress showed that the variable-density porous structure is similar to the whole solid structure. The maximum stress is ∼150.44 MPa, which is ∼27.6% more than that of the whole solid structure and ∼67.28% lower than that of uniform porous structure. The stress concentration is significantly improved, and the stress concentration in the neck of the femur is slightly slowed down. The cloud chart of the variable-density porous structure deformation is shown in Figure 11(c). Similarly, the maximum deformation is located near the neck of the femur, and the maximum value is ∼6.61 mm. The variable-density porous structure is adopted. Its deformation is similar to that of the whole solid structure, which is ∼32.73% more than the whole solid structure deformation and ∼68.33% lower than that of uniform porous structure. The deformation gradually decreased towards the tail of the femoral stem prosthesis. The whole deformation tended to rotate around the tail of the prosthesis.The region most prone to bone resorption in hip prosthesis replacement is shown in Figure12. This study is performed to analyze the bone resorption in the feature region. To better evaluate the bone resorption, the stress shielding rate is used as a measurement of bone resorption. Stress shielding rate is calculated as follows:(15)ψ=1−σpσo×%,where ψ refers to the stress shielding rate; σp refers to the stress after prosthesis implantation; and σo refers to the stress before prosthesis implantation.Figure 12
Featured regions of the femur: (a) regions prone to bone resorption and (b) stress distribution.
(a)(b)The mean stresses in the highlighted region are 6.82, 19.49, and 16.25 MPa for the complete solid prosthesis, uniform porous prosthesis, and variable-density porous prosthesis, respectively. Figure12 depicts the stress distributions in the feature areas of the three prostheses. When compared to a solid construction, the porous structure significantly improves stress shielding on the femur. The gravitational shielding caused by all-metal and porous stems differs significantly. This significant difference represented the fact that low-stretch implants are superior than high-stretch implants.Topology optimization of the multiscale lattice structure proposed in this paper has achieved variable-density topological optimization. This approach produced not only an excellent lightweight effect, but also a superior stress distribution, demonstrating the lattice structure’s capability and obviating the need for stress shielding.


## 5. Conclusion

This paper puts forward the topological optimization of a multiscale lattice structure based on the homogenization method. The following conclusions could be drawn: First, the homogenization approach is used to calculate and fit the function for calculating the equivalent mechanical performance of a cell. The results can speed up the computation and eliminate the need for the homogenization approach to be used many times. Second, a variable-density lattice structure design based on topological optimization is effectively implemented, leading to a novel lattice structure design idea. As a result, the characteristics of lattice materials are improved, resulting in even better results. Third, the proposed technology could well be employed in conjunction with 3D printing production. The technology can be used in a variety of sectors, including aerospace and tailored medical treatment, to produce small-batch high-performance lattice structures. The approach is also used to produce high-performance automotive lattice components, reducing the research and development cycle and saving money. Finally, the lattice structure had some special properties, such as high heat dissipation and high energy absorption. The proposed method can be further applied to heat dissipation, energy absorption, and other specific applications to realize the optimal design of high-performance lattice structures.



---

*Source: 4243476-2022-07-13.xml*