Orthopedic implant structure optimization method giving consideration to bone mass stability and fatigue life

Abstract

The invention discloses an orthopedic implant structure optimization method giving consideration to bone mass stability and fatigue life. The method mainly comprises porous structure topological optimization for realizing target domain bone mass maintenance and a porous structure generation strategy based on design variable iteration updating. Calculating the sensitivity of a target function and a constraint condition to a design variable for optimal iteration of implant density distribution; by measuring the function relation between the elastic modulus of the porous structure and the pseudo-density of the design domain unit, pseudo-density distribution of the design domain structure is converted into equivalent elastic modulus distribution of the porous structure, and then the porous structure is generated through parameterization. According to the method, a mechanical environment is constructed by optimizing a porous structure, loads are guided to be transmitted to bones, the bone mass of a target area is improved and kept through force stimulation, the sensitivity of a target function and constraint conditions to design variables is calculated, iteration is carried out through a moving asymptote optimization algorithm, and the bone mass is obtained. And finally obtaining the porous structure capable of keeping the bone mass of the target area and ensuring the fatigue life of the structure.

Description

Technical Field

[0001] The present invention belongs to the technical field of computational structure optimization and intelligent design algorithms, and specifically relates to an orthopedic implant structure optimization method that takes into account both bone mass stability and fatigue life. Background Art

[0002] Intervertebral fusion has become a key technology for treating various spinal lesions. Typically, the intervertebral disc is removed and a porous structure is placed between the vertebrae to promote spinal fusion. To avoid fusion settlement and stress shielding, the elastic modulus of the porous structure and bone must be matched. This allows the porous material to be sparsely distributed within the design space, ensuring that the porous structure's performance meets structural design requirements while maintaining low cost and modulus. Therefore, optimizing the material distribution of the porous structure is a challenge that requires particular attention.

[0003] Topology optimization is a numerical method that iteratively optimizes material placement within a design domain to optimally achieve specific functional objectives. Because it can effectively control manufacturing costs and improve structural performance, it has become one of the most important tools in structural design. Topology optimization allows for the design of structures using optimization methods given constraints and an objective function. However, in the field of bone implant design, the objective function and constraints often encompass both the bone and the scaffold, and there is a time dimension to the continuous remodeling of bone under force. Existing research primarily targets minimizing scaffold strain energy and employs volume as a design constraint, lacking consideration for the long-term stability of the bone-implant complex.

[0004] Unlike topology optimization design methods that target the initial strength or stiffness of porous structures after implantation, the topology optimization design method involved in the present invention can maintain the time dependence of the bone density change process, aiming to reshape the bone in a manner that adapts to the changing biomechanical environment. At the same time, the Tsai-Hill fatigue criterion, a representative unit of the porous structure, is used as a constraint to achieve the effect of increasing and maintaining bone mass in the target area and extending the fatigue life of the porous structure. Summary of the Invention

[0005] Purpose of the invention: The purpose of the present invention is to provide an orthopedic implant structure optimization method that takes into account both bone mass stability and fatigue life. The method can find the optimal layout of the porous structure by calculating the sensitivity of the objective function and constraints to the design variables and iterating using the moving asymptote algorithm. Then, based on the mechanical properties of the porous structure, after certain post-processing, the porous structure is finally parameterized and generated.

[0006] Technical solution: The method for optimizing the structure of an orthopedic implant that takes into account both bone mass stability and fatigue life according to the present invention comprises the following steps:

[0007] (1) Based on the topology optimization theory of variable density method, a topology optimization model is established with the maximum bone mass of the bone bridge as the objective function in the target domain and the Tsai-Hill fatigue criterion of the global P-norm porous structure in the design domain as the constraint condition;

[0008] (2) Calculate the displacement response of the structure under the load boundary conditions of the host bone, and calculate the sensitivity of the objective function and constraints to the design variables based on the displacement response;

[0009] (3) Using the MMA method to optimize the design variables until convergence, the pseudo-density distribution of the design variables is obtained;

[0010] (4) Based on the asymptotic homogenization theory, the relationship between the elastic modulus of the Thiessen polygon and the number of axial seed points is fitted, and the elastic modulus distribution is calculated according to the pseudo-density distribution and mapped to the seed point distribution of the Thiessen polygon;

[0011] (5) Expand the seed points to construct a porous prototype structure of Thiessen polygons;

[0012] (6) The obtained porous structure is placed in the non-design domain through translation, rotation and Boolean operation to obtain the final implant structure.

[0013] Furthermore, the topology optimization model in step (1) is:

[0014] find x

[0015]

[0016] stf1(x)=σ T-H -1≤0

[0017] f2(x)=F-KU=0

[0018]

[0019] Where x=[x1,x2,…,x N ] T , x i ∈[0,1]i=1,2,…,N, is the design variable of the variable density method, It represents the bone density of bone unit i in the target area within the mth year after surgery, and its value is the minimum value ρ min to 1, f1(x) is the inequality constraint of the fatigue life of the porous structure, σ T-H represents the global P-norm constraint of the Tsai-Hill fatigue criterion of the representative unit of the porous structure, f2(x) represents the equilibrium equation of the applied load, F represents the load boundary condition of the host bone, K is the global stiffness matrix of the porous structure, and U is the global displacement.

[0020] Furthermore, in step (1), the target domain is the area where the objective function acts, and the design domain is the area where the design variables act; the target domain and the design domain are in a spatially separated state.

[0021] Furthermore, the change in bone mass in the target area is the process of bone adapting to changes in external biomechanical stimulation, and the stimulation is expressed as:

[0022]

[0023] Wherein, the superscript m represents the mth year after surgery, i represents the i-th unit, Q represents the strain energy, and ρ represents the bone density. Within a certain time interval Δt, the bone density change of the i-th unit is It is expressed as follows:

[0024]

[0025] The time step of each iteration of the topology optimization algorithm is set to 1 year to simulate the bone remodeling process in the target area. Q1 and Q2 are the bone density change rate constants, η1 and η2 are the lower and upper limits of the “dead zone” stimulation threshold during bone remodeling, 2l represents the width of the “dead zone”, K is the reference value constant of the mechanical steady state, and η F is the damage limit of overload absorption, higher than η F Mechanical stimulation will cause bone resorption at a constant rate P.

[0026] Furthermore, the mathematical form of the Tsai-Hill fatigue criterion for the porous structure is expressed as:

[0027] During the iterative optimization process of the porous structure layout, the Tsai-Hill fatigue criterion FI of the porous structure representative unit i is i Expressed as:

[0028]

[0029] Among them, σ1, σ2, σ3 represent the principal stresses in three directions respectively, τ 12 ,τ 23 ,τ 13 Respectively represent the corresponding shear stress, X, Y, Z represent the fatigue strength of the unit in the three axial directions, R, S, T represent the shear strength of the unit in the three directions, and V is the auxiliary matrix, which is expressed as:

[0030]

[0031] The global P norm is used to condense the Tsai-Hill fatigue stress of all elements in the porous structure into a global constraint that is approximately equal to the maximum Tsai-Hill fatigue stress of the structure, expressed as:

[0032]

[0033] Among them, P is the norm parameter, Ω e for the entire design domain.

[0034] Furthermore, the process of calculating the sensitivity of the objective function to the design variables in step (2) is as follows:

[0035] Sensitivity of total bone density in the target area It is expressed as follows:

[0036]

[0037] Among them, V i is the volume of unit i;

[0038] Introducing the adjoint problem for the mth year:

[0039]

[0040] Among them, μ m is the global adjoint displacement vector, and the sensitivity of the total bone density in the target domain is simplified to:

[0041]

[0042] Global stiffness matrix K m About the design variable x e The partial derivative consists of two parts: one is the stiffness matrix k of the bone unit in the target domain i,m , the other part is the stiffness matrix K of the intervertebral fusion cage unit in the design domain e , the stiffness of the intervertebral fusion unit is independent of the time of bone remodeling, that is, the global stiffness matrix K m About the design variable x e The partial derivative of is expressed as:

[0043]

[0044] Among them, v i,m and v e,m are the local adjoint displacements of bone unit i in the target domain in the mth year, respectively, represented by the corresponding index vector L i and L e It is expressed as follows:

[0045] v i,m =L i μ m ,v e,m =L e μ m

[0046] Combined with the chain derivation rule, formula (19) is transformed into:

[0047]

[0048] The above formula contains arrive The derivation can be recursively deduced from m, m-1... to 0, and About x e The partial derivative is 0, and the final sensitivity relationship is obtained:

[0049]

[0050] The corresponding adjoint problem is expressed as:

[0051]

[0052] According to the above recursive relation, the solution of the adjoint problem of step m depends on step m+1, thus obtaining an explicit expression of the sensitivity of the objective function.

[0053] Furthermore, the sensitivity of the calculation constraint conditions to the design variables in step (2) is implemented as follows:

[0054] In finite element analysis, the stress of entity element i is expressed as:

[0055] σ i =D i B i u i,m =[σ1,σ2,σ3,τ 12 ,τ 23 ,τ 13 ] T

[0056] B i and D i represents the strain interpolation matrix and elastic tensor of element i, u i,m represents the nodal displacement of element i;

[0057] Tsai-Hill fatigue constraint with global P norm about x e The derivative of is:

[0058]

[0059] The first term on the right side of the equation (31) is expanded to:

[0060]

[0061] The second term on the right side of the equation (31) is the partial derivative of the stress of the Tsai-Hill fatigue constraint with respect to the solid element i, which is expressed as a combination of the six component partial derivatives:

[0062]

[0063] Combined with the auxiliary matrix V, it can be expressed as:

[0064]

[0065] Introducing the adjoint problem, we obtain the global P-norm Tsai-Hill fatigue constraint with respect to x e The partial derivative of is:

[0066]

[0067] where λ e The adjoint displacements calculated for the following adjoint problem:

[0068] K m λ e =L i T B i T D i Vσ i .

[0069] Beneficial effects: Compared with the prior art, the beneficial effects of the present invention are as follows: the present invention takes into account the continuous changes of the target domain bone under biomechanical stimulation, and realizes the growth and long-term maintenance of bone mass in the target area after the implantation of the porous structure, thereby improving the therapeutic effect from the initial to long-term after the implantation of the porous structure; compared with the optimization design method targeting the overall stiffness or strength of the structure, the present invention more flexibly considers the continuous changes in the stiffness or strength of the local area of ​​the porous structure, and realizes the topological optimization design of the non-uniform mechanical properties of the porous structure; the present invention adapts to the porosity process requirements of the porous structure in additive manufacturing, and can achieve the manufacturing requirements of the optimal configuration that meets the topological optimization constraints and different additive manufacturing processes by adjusting the upper and lower limits of the design variable density. BRIEF DESCRIPTION OF THE DRAWINGS

[0070] Figure 1 is a flow chart of the present invention;

[0071] Figure 2 is a schematic diagram of a finite element model of an intervertebral fusion cage calculation example according to an embodiment of the present invention;

[0072] Figure 3 is an iterative history statistical graph of the objective function of an embodiment of the present invention;

[0073] Figure 4 It is a schematic diagram of the optimization results and verification of the intervertebral fusion cage calculation example according to the embodiment of the present invention. DETAILED DESCRIPTION

[0074] The present invention is described in further detail below with reference to the accompanying drawings:

[0075] The present invention proposes an orthopedic implant structure optimization method that takes into account both bone mass stability and fatigue life, which mainly includes a porous structure topology optimization to achieve target domain bone mass retention and a porous structure generation strategy based on iterative update of design variables; the porous structure topology optimization calculates the sensitivity of the objective function and constraint conditions to the design variables, and is used for optimization iteration of the implant density distribution; the porous structure generation strategy based on iterative update of design variables measures the functional relationship between the elastic modulus of the porous structure and the pseudo-density of the design domain unit, and converts the pseudo-density distribution of the design domain structure into the equivalent elastic modulus distribution of the porous structure, thereby parameterizing the generation of the porous structure. Figure 1 As shown, the specific steps include:

[0076] Model initialization, the present invention adopts Figure 2 The finite element model shown in the figure reconstructs the three-dimensional model of the intervertebral fusion cage and adjacent vertebrae based on the anatomical characteristics of the intervertebral space. The model is divided into design domain, non-design domain, and target domain. The resolution of the model is 400 μm, and all areas are divided into C3D8 meshes with a total number of approximately 122W elements. The length, height, and thickness of the intervertebral fusion cage are L = 25 mm, H = 10 mm, and W = 10 mm, respectively. Regions 1 and 2 belong to the non-design domain, representing cortical bone and cancellous bone, respectively, with elastic moduli and Poisson's ratios of 16,000 MPa and 0.3 and 1,200 MPa and 0.3, respectively. Their bone structures are reconstructed from CT data of L3 and L4. Regions 3 and 4 constitute the intervertebral fusion cage, which are the intervertebral fusion cage frame and solid part, respectively. The elastic modulus and Poisson's ratio are both 110,000 MPa and 0.3. Region 3 belongs to the non-design domain, and region 4 represents the design domain. Region 5 belongs to the target domain, representing cancellous bone with an elastic modulus and Poisson's ratio of 1,200 MPa and 0.3. A composite load of a positive pressure of F = 1,200 N in the negative z direction, a flexion moment of 7.5 Nm, and a lateral bending moment of 7.5 Nm is applied to the upper surface of the L3 vertebra. The lower surface of the L4 vertebra constrains the degrees of freedom in six directions.

[0077] Step 2: Based on the topology optimization theory of the variable density method, a topology optimization model is established with the maximum bone bridge volume as the objective function and the global P-norm Tsai-Hill fatigue criterion of the design domain as the constraint condition, which is expressed as:

[0078]

[0079] in, It represents the bone density of bone unit i in the target area within the mth year after surgery, and its value is the minimum value ρ min to 1, f1(x) is the inequality constraint of the fatigue life of porous structure, where σ T-Hrepresents the global P-norm constraint of the Tsai-Hill fatigue criterion of the representative unit of the porous structure, f2(x) represents the equilibrium equation of the applied load, F represents the load boundary condition of the host bone, K is the global stiffness matrix of the porous structure, and U is the global displacement.

[0080] During the iterative optimization process of the porous structure layout, the Tsai-Hill fatigue criterion EI of the porous structure representative unit i is i Expressed as:

[0081]

[0082] Among them, σ1, σ2, σ3 represent the principal stresses in three directions respectively, τ 12 ,τ 23 ,τ 13 Respectively represent the corresponding shear stress, X, Y, Z represent the fatigue strength of the unit in the three axial directions, R, S, T represent the shear strength of the unit in the three directions, and V is the auxiliary matrix, which is expressed as:

[0083]

[0084] The global P norm is used to condense the Tsai-Hill fatigue stress of all elements in the porous structure into a global constraint that is approximately equal to the maximum Tsai-Hill fatigue stress of the structure, expressed as:

[0085]

[0086] Among them, P is the norm parameter, Ω e for the entire design domain.

[0087] Step 3: Finite element analysis, calculate the displacement response of the structure under the boundary condition F, and calculate the sensitivity based on the displacement response.

[0088] (3.1) The sensitivity analysis of the objective function f0(x) adopts the following:

[0089] After m+1 years of remodeling, the bone density of unit i in the target domain is It can be calculated according to the bone remodeling control equation as follows:

[0090]

[0091] The value of bone density in the m+1th year is not only related to the bone density in the mth year, but also depends on the physiological curve during bone remodeling. To obtain an explicit expression for the partial derivative of the bone density in the m+1th year with respect to the design variables, a step-by-step recursive calculation is required. The value of is determined by the mechanical stimulation of the bone unit through the physiological curve of bone remodeling. The mechanical stimulation M of the i-th unit in the target domain of the m-th year isi,m Determined by the following relationship:

[0092]

[0093] Where, the superscript m represents the mth year after surgery, i represents the i-th unit, Q represents the strain energy, and ρ represents the bone density. Within a certain time interval Δt, the bone density change of the i-th unit is It is expressed as follows:

[0094]

[0095] The time step of each iteration of the topology optimization algorithm is set to 1 year to simulate the bone remodeling process in the target area. Q1 and Q2 are the bone density change rate constants, η1 and η2 are the lower and upper limits of the “dead zone” stimulation threshold during bone remodeling, 2l represents the width of the “dead zone”, K is the reference value constant of the mechanical steady state, and η F is the damage limit of overload absorption, higher than η F Mechanical stimulation will cause bone resorption at a constant rate P.

[0096] Q i,m is the strain energy of the i-th unit in the target domain in the m-th year, Q i,m It is obtained by the following formula:

[0097]

[0098] Among them, u i,m is the displacement of the i-th unit node in the target domain in the m-th year, k i,m is the stiffness matrix of the i-th unit in the target domain in the m-th year. According to the partial derivative of the mechanical stimulus with respect to the design variables and the bone remodeling control equation, we can get Partial derivative of:

[0099]

[0100] Among them, V i is the volume of the i-th unit, and the above formula is transposed:

[0101]

[0102] Expanding the partial derivative of the strain energy in the above formula with respect to the design variables yields:

[0103]

[0104] Substitute it into formula (6):

[0105]

[0106] Since the stiffness of unit i has a certain functional relationship with bone density, the third term on the right side of the equal sign in the above formula, the partial derivative of the stiffness with respect to the design variable, can be replaced by the partial derivative of the bone density with respect to the design variable:

[0107]

[0108] Partial derivatives of combined bone density with respect to design variables:

[0109]

[0110] Next, we need to process the second term on the right side of the equal sign to simplify it into an explicit expression, u i,m Represents the displacement of the i-th unit node, which is a local displacement. It has the following relationship with the global displacement:

[0111] u i,m =L i U m #(15)

[0112] U m represents the global displacement in the mth year, L i Represents the index from the global displacement to the corresponding local displacement, and the global displacement is related to the global external force vector F and the global stiffness matrix K in the mth year m The relationship between them is as follows:

[0113] K m U m =F#(16)

[0114] The partial derivatives of both sides with respect to the design variables are taken simultaneously, noting that the external force vector is assumed to be constant during the entire optimization process:

[0115]

[0116] Combine equations (15) and (17) and substitute them into equation (14):

[0117]

[0118] Therefore, the sensitivity of the total bone density of the target domain can be expressed as:

[0119]

[0120] According to the above formula, the adjoint problem of the mth year is introduced:

[0121]

[0122] where μ m is the global adjoint displacement vector, then Equation (19) can be simplified as:

[0123]

[0124] Global stiffness matrix K m About the design variable x e The partial derivative consists of two parts: one is the stiffness matrix k of the bone unit in the target domain i,m , the other part is the stiffness matrix k of the intervertebral fusion cage unit in the design domain e , the stiffness of the intervertebral fusion unit is independent of the time of bone remodeling, that is, the global stiffness matrix K m About the design variable x e The partial derivative can be expressed as:

[0125]

[0126] Among them, v i,m and v e,m are the local adjoint displacement of bone unit i in the target domain at the mth year, and the local adjoint displacement of the intervertebral fusion cage unit in the design domain at the mth year, which are expressed as follows:

[0127] v i,m =L i μ m ,v e,m =L e μ m #(twenty three)

[0128] Combined with the chain rule, formula (21) can be transformed into:

[0129]

[0130] The above formula contains arrive From the derivation of formula (18), we can get About x e The partial derivative of is:

[0131]

[0132] Substitute into the first term on the right side of the equal sign in formula (24):

[0133]

[0134] Formula (26) can be recursively extended from m, m-1... to 0, and About x e The partial derivative value is 0, and the final sensitivity relationship can be obtained:

[0135]

[0136] The corresponding adjoint problem can be expressed as:

[0137]

[0138] The Φ update method is as follows:

[0139]

[0140] According to the above recursive relationship, the solution of the adjoint problem of step m depends on step m+1, from which an explicit expression of the sensitivity of the objective function can be obtained.

[0141] (3.2) The sensitivity analysis of the constraint f1(x) adopts the following:

[0142] In finite element analysis, the stress of entity element i can be expressed as:

[0143] σ i =D i B i u i,m =[σ1,σ2,σ3,τ 12 ,τ 23 ,τ 13 ] T #(30)

[0144] B i and D i represents the strain interpolation matrix and elastic tensor of element i, u i,m represents the nodal displacement of element i;

[0145] Tsai-Hill fatigue constraint with global P norm about x e The derivative of is:

[0146]

[0147] The first term on the right side of the equal sign in formula (31) can be expanded as:

[0148]

[0149] The second term on the right side of the equation (31) is the partial derivative of the stress of the Tsai-Hill fatigue constraint with respect to the solid element i, which can be expressed as a combination of six component partial derivatives:

[0150]

[0151] Combined with the auxiliary matrix V, it can also be expressed as:

[0152]

[0153] Global displacement, global external force vector F and global stiffness matrix K m The relationship between them is as follows:

[0154] Km U m =F#(35)

[0155] u i,m Represents the displacement of the i-th unit node, which is a local displacement. It has the following relationship with the global displacement:

[0156] u i,m =L i U m #(36)

[0157] The partial derivatives of the design variables on both sides of Equation (35) are calculated simultaneously, noting that the external force vector is assumed to be constant during the entire optimization process:

[0158]

[0159] Shift the items to get:

[0160]

[0161] The partial derivatives of the design variables on both sides of formula (30) are calculated:

[0162]

[0163] Substituting formula (38) into formula (39) yields:

[0164]

[0165] Substituting Equations (32), (34) and (40) into Equation (31), we obtain the Tsai-Hill fatigue constraint of the global P norm with respect to x e The partial derivative of is:

[0166]

[0167] According to the above formula, the adjoint problem is introduced accordingly:

[0168] K m λ e =L i T B i T D i Vσ i #(42)

[0169] Finally, the global P-norm Tsai-Hill fatigue constraint with respect to x e The partial derivative of is expressed as:

[0170]

[0171] Step 4: Based on the above recursive relationship, the explicit expression of the constraint sensitivity is obtained. Combined with the sensitivity of the objective function, the MMA optimization problem is constructed. The optimization problem is solved and the design variables are updated using the Lagrangian dual algorithm. The iteration is stopped when the bone mass change in the target domain is less than 0.1%. The optimization result is obtained and the pseudo-density distribution of the design variables is obtained.

[0172] Step 5: Porous structure generation strategy based on iterative update of design variables.

[0173] Thiessen polygons with different numbers of seed points were constructed for porous structures. The elastic moduli of Thiessen polygons with different porosities were calculated using finite element analysis. The relationship between the elastic modulus of the Thiessen polygons and the number of axial seed points was fitted based on the asymptotic homogenization theory. The elastic modulus distribution was calculated based on the pseudo-density distribution obtained from topology optimization and mapped to the seed point distribution of the Thiessen polygons. The seed points were then expanded to construct the porous prototype structure of the Thiessen polygons. Finally, the resulting porous structure was combined with the non-design domain through operations such as translation, rotation, and Boolean operations to obtain the final implant structure.

[0174] Combining the measured relationship between the elastic modulus and pseudo-density (porosity) of the Thiessen polygon structure, the pseudo-density distribution of the design domain structure was converted into an elastic modulus distribution. The relationship between the measured elastic modulus of the Thiessen polygon and the number of axial seed points used to construct the Thiessen polygon was mapped to a seed point distribution based on the elastic modulus distribution of the design domain structure. A Thiessen polygon structure with a rod diameter of 400 μm was generated. The rods on the porous Thiessen polygon surface were connected to obtain a porous Thiessen polygon structure for the solid portion.

[0175] like Figure 3 As shown in the figure, with the increase of the number of iterations, the bone volume of the bone bridge in the intervertebral fusion cage gradually increased and stabilized in the third year, indicating that the model has found an optimal solution under the given conditions.

[0176] The intervertebral fusion cage structure of the 179th topology optimization iteration is as follows Figure 4 As shown in (a), the depth of the unit color indicates the level of pseudo-density. The porous result of the solid part is combined with the frame part through translation, rotation and Boolean operation to obtain the final intervertebral fusion cage. Figure 4 (b) shows the stress distribution cloud diagram of the final porous structure of the intervertebral fusion cage when the lower surface nodes are fixed and the upper surface nodes are loaded with a load of 1200N in the negative Z direction, which illustrates the effectiveness of this method for intervertebral fusion cage structure optimization.

[0177] The foregoing is merely a preferred embodiment of the present invention. It should be noted that those skilled in the art may make various improvements and modifications without departing from the principles of the present invention, and such improvements and modifications are also within the scope of protection of the present invention. Structures, devices, and operating methods not specifically described or explained herein shall, unless otherwise specified or limited, be implemented in accordance with conventional means in the art.

Claims

1. A method for optimizing the structure of an orthopedic implant that takes into account both bone mass stability and fatigue life, characterized in that: The following steps are involved: (1) Based on the topology optimization theory of variable density method, a topology optimization model is established with the maximum bone mass of the bone bridge as the objective function in the target domain and the Tsai-Hill fatigue criterion of the global P-norm porous structure in the design domain as the constraint condition; (2) Calculate the displacement response of the structure under the load boundary conditions of the host bone, and calculate the sensitivity of the objective function and constraints to the design variables based on the displacement response; (3) Using the MMA method to optimize the design variables until convergence, the pseudo-density distribution of the design variables is obtained; (4) Based on the asymptotic homogenization theory, the relationship between the elastic modulus of the Thiessen polygon and the number of axial seed points is fitted, and the elastic modulus distribution is calculated according to the pseudo-density distribution and mapped to the seed point distribution of the Thiessen polygon; (5) Expand the seed points to construct a porous prototype structure of Thiessen polygons; (6) The obtained porous structure is placed in the non-design domain through translation, rotation and Boolean operation to obtain the final implant structure.

2. The orthopedic implant structure optimization method according to claim 1, wherein: The topology optimization model in step (1) is: find x s.t.f1(x)=σ T-H -1≤0 f2(x)=F-KU=0 Where x=[x1,x2,…,x N ] T , x i ∈[0,1]i=1,2,…,N, is the design variable of the variable density method, It represents the bone density of bone unit i in the target area within the mth year after surgery, and its value is the minimum value ρ min to 1, f1(x) is the inequality constraint of the fatigue life of the porous structure, σ T-H represents the global P-norm constraint of the Tsai-Hill fatigue criterion of the representative unit of the porous structure, f2(x) represents the equilibrium equation of the applied load, F represents the load boundary condition of the host bone, K is the global stiffness matrix of the porous structure, and U is the global displacement.

3. The orthopedic implant structure optimization method according to claim 1, wherein: In step (1), the target domain is the area where the objective function acts, and the design domain is the area where the design variables act; the target domain and the design domain are in a spatially separated state.

4. The orthopedic implant structure optimization method according to claim 2, wherein: The change in bone mass in the target area is the process of bone adapting to changes in external biomechanical stimuli. The stimulus is expressed as: Wherein, the superscript m represents the mth year after surgery, i represents the i-th unit, and Q represents the strain energy. ρ represents bone density. Within a certain time interval Δt, the bone density change of the i-th unit is It is expressed as follows: η1=K(1-l) η2=K(1+l) The time step of each iteration of the topology optimization algorithm is set to 1 year to simulate the bone remodeling process in the target area. Q1 and Q2 are the bone density change rate constants, and η1 and η2 are the "dead zone" in the bone remodeling process. The lower and upper limits of the stimulation threshold, 2l represents the width of the "dead zone", K is the reference value constant of the mechanical steady state, η F is the damage limit of overload absorption, higher than η F Mechanical stimulation will cause bone resorption at a constant rate P.

5. The orthopedic implant structure optimization method according to claim 1, wherein: The mathematical form of the Tsai-Hill fatigue criterion for porous structures is expressed as: During the iterative optimization process of the porous structure layout, the Tsai-Hill fatigue criterion FI of the porous structure representative unit i is i Expressed as: FI i =F(σ2-σ3) 2 +G(σ3-σ1) 2 +H(σ1-σ2) 2 +2Lτ 12 2 +2Mt 23 2 +2Nτ 13 2 =s i T Vs i Among them, σ1, σ2, σ3 represent the principal stresses in three directions respectively, τ 12 ,τ 23 ,τ 13 Respectively represent the corresponding shear stress, X, Y, Z represent the fatigue strength of the unit in the three axial directions, R, S, T represent the shear strength of the unit in the three directions, and V is the auxiliary matrix, which is expressed as: The global P norm is used to condense the Tsai-Hill fatigue stress of all elements in the porous structure into a global constraint that is approximately equal to the maximum Tsai-Hill fatigue stress of the structure, expressed as: Among them, P is the norm parameter, Ω e for the entire design domain.

6. The orthopedic implant structure optimization method according to claim 1, wherein: The process of calculating the sensitivity of the objective function to the design variables in step (2) is as follows: Sensitivity of total bone density in the target area It is expressed as follows: Among them, V i is the volume of unit i; Introducing the adjoint problem for the mth year: Among them, μ m is the global adjoint displacement vector, and the sensitivity of the total bone density in the target domain is simplified to: Global stiffness matrix K m About the design variable x e The partial derivative consists of two parts: one is the stiffness matrix k of the bone unit in the target domain i,m , the other part is the stiffness matrix K of the intervertebral fusion cage unit in the design domain e , the stiffness of the intervertebral fusion unit is independent of the time of bone remodeling, that is, the global stiffness matrix K m About the design variable x e The partial derivative of is expressed as: Among them, v i,m and v e,n are the local adjoint displacements of bone unit i in the target domain in the mth year, respectively, represented by the corresponding index vector L i and L e It is expressed as follows: v i,m =L i μ m ,v e,m =L e μ m Combined with the chain derivation rule, formula (19) is transformed into: The above formula contains arrive The derivation can be recursively deduced from m, m-1... to 0, and About x e The partial derivative is 0, and the final sensitivity relationship is obtained: The corresponding adjoint problem is expressed as: According to the above recursive relation, the solution of the adjoint problem of step m depends on step m+1, thus obtaining an explicit expression of the sensitivity of the objective function.

7. The orthopedic implant structure optimization method according to claim 1, wherein: The process of calculating the sensitivity of the constraints to the design variables in step (2) is as follows: In finite element analysis, the stress of entity element i is expressed as: s i =D i B i you i,m =[σ1,σ2,σ3,τ 12 ,t 23 ,t 13 ] T B i and D i represents the strain interpolation matrix and elastic tensor of element i, u i,m represents the nodal displacement of element i; Tsai-Hill fatigue constraint with global P norm about x e The derivative of is: The first term on the right side of the equation (31) is expanded to: The second term on the right side of the equation (31) is the partial derivative of the stress of the Tsai-Hill fatigue constraint with respect to the solid element i, which is expressed as a combination of the six component partial derivatives: Combined with the auxiliary matrix V, it can be expressed as: Introducing the adjoint problem, we obtain the global P-norm Tsai-Hill fatigue constraint with respect to x e The partial derivative of is: where λ e The adjoint displacements calculated for the following adjoint problem: K m λ e =L i T B i T D i Vσ i 。

Images