理想的骨折内固定植入物刚度应具有随时间变化的性能，使骨折在不同愈合阶段都能产生合理的力学刺激，可降解材料可以满足这一性能。考虑到材料降解过程的时间依赖性，提出一种具有时变刚度的骨折内固定植入物复合结构拓扑优化设计方法。采用相对密度和降解残留率描述两种具有不同降解速度和弹性模量的材料分布和降解状态，建立降解模拟-力学分析耦合数学模型；基于变密度法设计双材料复合结构，使之具有时变刚度特性。以胫骨骨折治疗用的接骨板为例，采用所提出方法设计具有时变刚度特性的复合结构接骨板，优化结果显示高刚度的材料1形成柱状的支撑结构，低刚度的材料2分布在降解边界和内部。利用骨重塑模拟模型对优化后的接骨板进行评估，经过11个月重塑，使用可降解时变刚度接骨板、钛合金接骨板和不锈钢接骨板的骨痂平均弹性模量分别为8 634、8 521、8 412 MPa，表明使用可降解时变刚度接骨板可使骨痂取得更好的重塑效果。
引用本文： 孙浩, 丁晓红, 徐世鹏, 段朋云, 熊敏, 张横. 具有时变刚度的骨折内固定植入物结构设计及骨重塑效果评估. 生物医学工程学杂志, 2024, 41(3): 595-603. doi: 10.7507/1001-5515.202311037 复制
0 引言
骨折是一种常见的创伤性疾病，内固定术是最为常用且有效的骨折治疗方法。骨折愈合是一个力学调控和生物调控的复杂过程，骨折内固定植入物需要在治疗初期提供足够的支持和稳定性，随着骨愈合的进行，内固定植入物的刚度逐渐降低，实现力流向断骨处转移，促进骨愈合[1-3]。这要求内固定植入物的刚度应随时间发生变化（时变刚度），而使用可降解材料设计的植入物可满足这一要求。参考文献[4]，绘制具有时变刚度的理想可降解植入物的力学性能随骨折愈合过程的变化趋势图（见图1），骨愈合过程可以分为三个阶段：炎症期（阶段1）、骨痂形成期（阶段2）、骨痂重塑期（阶段3）。内固定植入物的时变刚度如虚线所示，在炎症期和骨痂形成期，骨折处血肿逐渐被肉芽组织代替，最终骨化成为编织骨，此时骨痂的弹性模量较低，植入物需要提供良好的固定稳定性；在骨痂重塑期，骨骼逐渐恢复到健康骨骼状态，植入物刚度以较快的速率降低，使应力逐渐向骨痂转移，在骨折完全愈合之后，内固定植入物的刚度接近零。因此如何设计内固定植入物结构，使其刚度特性符合骨愈合规律的要求是一个亟待解决的问题。
图1
理想骨折内固定植入物降解率与力学性能变化趋势
Figure1.
The degradation rate and mechanical trend of ideal fracture internal fixation implants
目前可降解内固定植入物的设计主要通过合金化、冷加工、热处理和表面改性等方法改变可降解材料的降解行为，使其具有与骨愈合相匹配的刚度[5-7]。结构拓扑优化设计方法通过设计结构的材料分布使结构得到所需的力学性能，Wu等[8]通过考虑骨重塑过程的时间依赖性，推导了时域的灵敏度，开发了一种新的下颌骨接骨板拓扑优化方法，但未考虑骨固定板降解；Chandra等[9]考虑到生物降解过程中钢板的尺寸、均匀的降解率和整个过程中钢板的固定稳定性，设计了一种可生物降解的植入钢板，却忽视了材料降解对结构的影响。由于内固定植入物的力学性能随降解时间发生变化，而力学性能的变化会影响骨愈合效果，因此如何模拟降解过程，并量化材料降解对内固定植入物结构力学性能的影响，设计具有时变刚度的内固定植入物使之符合骨愈合规律的要求，具有重要的意义。目前，通过拓扑优化设计方法对可降解内固定植入物进行结构设计以改变可降解材料降解行为的研究仍较少[10-11]。
本文提出一种具有时变刚度的骨折内固定植入物复合结构拓扑优化设计方法，建立考虑结构降解的有限元模型，获得降解过程中中间结构的力学性能，通过结构拓扑优化实现具有不同降解速度和弹性模量的双材料可降解植入物，并利用骨重塑模拟模型对设计结果的时变刚度性能和骨重塑效果进行评估。
1 可降解结构的力学分析
1.1 具有时变刚度特性的复合结构设计构思
目前，单一材料的可降解内固定植入物的刚度很难与骨愈合周期相匹配，会出现降解过快发生固定失效和降解过慢导致应力屏蔽[12-14]，因此考虑采用两种不同的可降解材料。如图2所示，生物降解材料1降解速度慢、刚度大，材料2性能与其完全相反。在给定的边界条件下，在设计域内确定两种材料的最优布局，以实现在早期结构刚度缓慢降低，而后期结构刚度能够以相对较快速率降低。
图2
考虑时变刚度特性的复合结构设计
Figure2.
Design of composite structures considering time-varying stiffness characteristics
1.2 问题描述
进行内固定植入物的拓扑优化设计，模拟结构的降解过程并量化材料降解对结构性能的影响是研究的关键。分别使用相对密度(xe)和降解残留率(ye)描述两种材料的分布和降解状态，如式(1)～(2)所示。通过xe决定单元e为材料1还是材料2；通过ye表示降解状态，ye = 1，代表单元e完全未降解，ye = ymin代表单元e完全降解。
 |
 |
式中，xmin和ymin取0.001；m为单元个数。
1.3 结构降解模拟-力学分析耦合模型
采用均匀腐蚀方法，以平均降解率模拟结构的降解[15-16]。为简化模拟过程，假设：① 降解只发生在材料的固液交界处且均匀进行，降解边界随材料的降解而移动；② 材料的降解具有均一的降解速率，不考虑其他环境参数（如温度、pH值等）；③ 腐蚀环境相同，忽略贴近骨面和贴近组织的两面对降解速率的影响。需要说明的是，均匀腐蚀与实际腐蚀存在一定的差异，实际腐蚀是一个随机过程，不仅和腐蚀环境有关，而且存在一定的不确定性。但降解实验表明，以均匀腐蚀模拟结构的降解与降解实验中随机腐蚀降解的质量损失率差距在0.1%以内[17]，表明均匀降解模拟方法是可行的。
降解模拟过程如图3所示，单元颜色由深到浅表示材料从未开始降解到完全降解，可降解边界使用黑色粗虚线表示。图3a为结构网格，图3b为降解过程，降解过程示意如图3b1～3b5所示。其中从b1到b3可降解边界处的一层单元格颜色逐渐由深变浅，表示该层单元逐渐降解；从b3到b4，可降解边界从外层单元移动到相邻的内层单元；然后从b4到b5内层单元开始降解。在一个降解时间步完成后，对其降解结果进行结构性能分析，获得结构性能随时间的变化关系。单元降解的数学表达为：
图3
降解过程示意图
a. 结构网格；b. 降解过程；c. 单元
a. meshes of the structure; b. simulation of structure degradation; c. the location mark of element
 |
式中，i表示降解时间步，de表示单元e的降解速率，t为单位降解时间步的时间。
单元e和相邻单元位置标记如图3c所示，单元e的位置为(k, l)。通过判断降解单元相邻单元的状态，来决定该单元是否发生降解，为了保证数值计算的稳定，设定yb=0.1，当相邻单元的残留率小于yb时可沿该单元方向进行降解，如当yk − 1,l < yb时，认为单元e能够从左侧开始降解，当yk + 1,l < yb时，认为单元e可以从右侧开始降解，上下方向类似[18]。通过对相邻单元降解状态的判断，单元e的降解速率可以表示为：
 |
式中，Ω表示e相邻单元的集合；α表示Ω中的第α个单元；d为材料的降解速率；Nα由式(5)得到：
 |
为了便于求导，使用Heaviside方程将式(4)近似为连续方程，单元的降解方程可以改写为：
 |
阶跃函数H为：
 |
式中，γ是一个常数。
由式(3)和式(6)，单元e残留率的更新方程为：
 |
使用Kreisselmeier-Steinhauser函数将式(8)最大函数问题近似为连续方程[19-20]：
 |
式中，γ = 6，η = 6。
结构平衡方程如式(10)所示：
 |
式中，K是结构的整体刚度矩阵，由于材料的降解，K与时间相关；U为位移向量；F为力向量。
对于可降解结构，K由单元刚度矩阵ke与ye组装得到，ye是关于降解速率d和降解时间i的函数。将材料降解对结构力学性能的影响耦合到有限元分析过程中，可降解结构的刚度矩阵为：
 |
式中，B表示几何矩阵；D为弹性系数矩阵；Ωe为单元的体积。
降解模拟算法包括以下步骤。步骤1：初始化残留率ye，使用k和l对单元e进行标记；步骤2：判断相邻单元的降解状态；步骤3：使用式(6)计算单元e的降解速率de；步骤4：使用式(9)计算单元e的残留率ye；步骤5：根据ye，使用式(10)进行有限元分析；重复步骤1～5，直到指定的时间步td。
2 双材料可降解结构的拓扑优化
优化设计的目的是通过确定两种具有不同可降解性能的材料的最佳分布，使复合结构在指定降解时间步中间结构的刚度之和最大。使用结构柔顺度表征结构刚度，可降解复合结构柔顺度最小的优化数学模型可以表示为：
 |
式中，X为设计变量的集合；C为结构柔顺度，ʌ为降解时间步的集合；v1为复合结构中材料1的最终体积，f是给定的体积分数，v0为总体积。
为获得轮廓清晰的优化结构，使用基于固体各向同性材料惩罚模型（solid isotropic material with penalization，SIMP）的材料插值机制：
 |
 |
式中，E为弹性模量；p为惩罚系数，取3.0；下标1、2分别表示可降解材料1、2。
2.1 灵敏度分析
在每个降解时间步对结构进行有限元分析，计算结构柔顺度，其灵敏度可以通过求导得到。目标函数对于设计变量X的灵敏度可以表示为：
 |
结构柔顺度C是关于xe、ye和i的函数：
 |
式中，u表示单元位移矩阵；k0为单元刚度矩阵。
根据式(16)定义的柔顺度，在降解时间步i对设计变量xe的导数为：
 |
其中：
 |
 |
 |
材料1的体积为：
 |
则材料1的体积对于设计变量xe的导数为：
 |
为克服棋盘格效应，采用基于Heaviside的灵敏度过滤方法使密度变量的值趋近0或1。使用如式(23)所示的Heaviside方程对中间密度场xei进行过滤，xei由式(24)计算得到。
 |
 |
式中，β为常数，控制函数的平滑度，如果β = 0，则是原始的密度过滤器，如果β接近于无穷大，则Heaviside为阶跃函数。本文取β的初始值为1，每100次迭代对β的值加倍，最大值为512。
2.2 优化设计流程
双材料可降解结构的拓扑优化设计流程如图4所示，主要步骤如下。
图4
双材料可降解结构的拓扑优化设计流程
Figure4.
Topological optimization design process for dual material degradable structures
第一步：构造初始几何模型并进行离散化。定义数学模型的参数，给定体积分数和过滤半径，设置xe和ye的初始值，并定义可降解边界、载荷和约束等边界条件；第二步：根据式(9)进行降解模拟，计算得到残留率；第三步：根据式(10)和(11)对考虑材料降解的复合结构进行有限元分析；第四步：根据式(15)进行灵敏度分析。重复第二步到第四步，直到指定的时间步。第五步：利用MMA优化算法更新设计变量；第六步：重复第二步到第五步，直到满足约束条件或目标函数收敛。当满足以下两个条件之一时，迭代终止：① 连续两次迭代的结构柔顺度变化小于收敛判据；② 迭代步数达到指定值，本文中取100。
3 胫骨接骨板设计
3.1 设计模型及材料性能
使用如图5a所示的胫骨-接骨板系统简化模型，接骨板的尺寸如图5b所示，接骨板中间区域（橘黄色）为可降解结构，两侧（灰色）为不可降解结构。接骨板植入人体后，可降解结构垂直于x轴和z轴的四个面都能够发生降解，其中垂直于x方向两个面的面积较大，这两个面的降解对结构刚度的影响也较大。为了将设计问题简化，我们假定可降解结构仅能从垂直于x方向的两个面发生降解，忽略垂直于z方向两个面的降解。如图5c所示，将垂直于z方向的面作为设计平面，利用考虑时变刚度特性的可降解复合结构拓扑优化设计方法对设计平面进行优化设计，然后沿z方向拉伸得到三维可降解复合结构，实现可降解接骨板的拓扑优化设计。
图5
胫骨-接骨板固定系统
a. 设计模型；b. 接骨板模型；c. 设计区域；d. 接骨板力学边界条件；e. 接骨板的简化模型；f. 降解边界条件
Figure5. Bone and bone plate fixation systema. design model; b. bone plate model; c. design area; d. mechanical boundary conditions of bone plates; e. simplified model of bone plate; f. degradation boundary conditions
如图5d所示，胫骨平台的力学边界条件可以简化为底端全约束，顶端皮质骨受到垂直向下的载荷。骨愈合过程中，骨痂的刚度逐渐增大，但都远小于接骨板的刚度，由于载荷分流，接骨板承担大部分的载荷，设计平面的力学边界条件简化为图5e所示，中心点全约束，中心线其余节点约束竖直方向自由度，上下边受到垂直方向的压力，F1 = 1 600 N。可降解边界条件如图5f所示，上下侧为不可降解边界，左右侧为可降解边界。根据文献[12]，材料属性如表1所示。
3.2 拓扑优化结果及其时变刚度特性分析
采用2.2小节的设计流程，使用尺寸为0.1 mm的四节点四边形平面单元对设计域进行网格划分，拓扑优化过滤半径设置为单元尺寸的两倍，材料1的体积分数为0.5，优化目标为第1和第90天中间结构柔顺度之和最小。优化迭代历程如图6a所示，随着迭代的进行逐渐获得具有清晰边界的材料分布轮廓，优化目标达到了收敛状态。拓扑优化设计结果如图6b所示，高刚度材料1（黑色）形成了柱状的支撑结构，在可降解边界处外层为材料2，结构左右呈现对称性，高刚度材料1（黑色）主要分布在左右侧可降解边界处，在左右侧中部形成V型开口，并在靠近加载侧形成了类似X型的交叉结构，连接左右两侧的高刚度材料，提高了结构整体的稳定性。时变刚度接骨板设计结果如图6c所示，两侧为惰性金属，中间区域为可降解复合结构。需要说明的是，结构的最优拓扑形态主要与所受的载荷相关，本文采用了承受偏心压力的圆柱模拟胫骨，考虑了实际胫骨承受的主要载荷为压力和弯矩，因此图6的设计结果具有合理性。
图6
拓扑优化设计结果
a. 拓扑优化迭代历程；b. 二维设计结果；c. 可降解时变刚度接骨板
Figure6. Topology optimization design resulta. iteration history with two degradable interfaces; b. 2D design results; c. degradable time-varying stiffness plate
考虑时变刚度特性拓扑优化设计结果的降解过程如图7所示。由于材料1刚度较大，其降解对结构刚度的影响程度较大。降解开始时，可降解边界附近的材料2先逐渐降解，结构刚度缓慢下降。当外层的材料2降解完，结构内部的材料1开始降解，结构刚度下降的速度开始增大。
图7
考虑时变特性拓扑优化设计结果的降解过程
Figure7.
Degradation process of topology optimization design results considering time-varying characteristics
3.3 骨重塑效果评估
骨重塑是骨骼受到外部生物力学刺激时产生的一种自适应过程。骨重塑模拟模型能够描述在生物力学刺激作用下骨密度随时间的变化过程。利用应变能密度预测骨重塑过程[21-22]，力学刺激Ξ可以表示为：
 |
式中，U为应变能密度，ρ为骨密度。骨密度在一定时间间隔
t的变化量
ρ可以表示为：
 |
 |
 |
式中，t为时间；R1和R2为常数，分别取18 (g/cm3)2/(J/月)和600 (g/cm3)3/(J2/月)；θ1和θ2是力学刺激的阈值，分别由式(26)和(27)计算得到；s和δ为常数，分别取0.000 36 J/(g*cm3)和0.1。骨密度变化量与力学刺激之间的曲线描述如图8所示，可以分为低载退化、惰性、骨重塑和过载退化四个区间。在低载退化和过载退化区间骨密度减小，在惰性区间骨密度不变，在骨重塑区间骨密度增大。
图8
骨密度变化率与力学刺激关系示意图
Figure8.
The schematic of relationship between the rate of change in bone density and mechanical stimulus
假设发生骨重塑的骨材料属性为各向同性，其弹性模量可以通过密度进行量化，经验公式为：
 |
式中ζ为弹性模量；τ为常数，取3 790 MPa/(g/cm3)3。
使用图5a所示的简化模型模拟骨痂的重塑过程。边界条件如图5d所示，顶端受到垂直向下的载荷，底端全约束。载荷为普通人体重（70 kg）的三倍，即F = 2 100 N。骨痂的初始密度为1.0 g/cm3，最大为1.74 g/cm3，最小为0.1 g/cm3。接骨板除设计域外结构和螺钉为钛合金，弹性模量为110 000 MPa。皮质骨和松质骨的弹性模量分别为15 750 MPa和1 100 MPa。不锈钢的弹性模量为210 000 MPa。
采用弹性模量表征骨重塑过程的骨痂单元变化，平均弹性模量一定程度上能够反映骨痂整体结构刚度的变化。骨痂的尺寸和初始状态如图9a所示，图9b、c、d分别为使用可降解时变刚度接骨板、相同结构尺寸的钛合金接骨板和不锈钢接骨板固定时，经过6个月和11个月时的骨痂弹性模量分布云图。接骨板位于骨痂的左侧，由于接骨板造成的应力遮挡效应，靠近接骨板侧的骨痂受到的力学刺激较小，因此弹性模量较低。与使用钛合金接骨板和不锈钢接骨板相比，使用可降解时变刚度接骨板时，骨痂重塑过程中左侧弹性模量小于3 790 MPa的区域明显较小。经过11个月重塑之后，使用可降解时变刚度接骨板时，骨痂弹性模量基本都大于3 790 MPa；而使用钛合金接骨板和不锈钢接骨板时，骨痂左侧仍有部分区域弹性模量小于3 790 MPa。结果表明，使用可降解时变刚度接骨板可以使靠近接骨板侧的骨痂在重塑过程中受到更多的力学刺激，有效地减缓应力遮挡效应，取得更好的重塑效果。对比相同重塑时间下使用三种接骨板的骨痂平均弹性模量的差值，结果如表2所示，随着重塑的进行，三种骨板对应的骨痂的平均弹性模量均逐渐增大，使用可降解时变刚度接骨板时的骨痂平均弹性模量均高于钛合金接骨板和不锈钢接骨板，对比不同骨板的骨痂平均弹性模量的差值，随着时间的增加差值逐渐变大，并且差值增大的速率也逐渐增大，说明可降解时变刚度接骨板的结构刚度减小的速率逐渐增大。
图9
骨痂弹性模量分布云图
a. 初始状态；b. 可降解时变刚度接骨板；c. 钛合金接骨；d. 不锈钢接骨板
Figure9. Distribution cloud map of elastic modulus of callusa. initial state; b. degradable time-varying stiffness bone plate; c. titanium alloy bone plate; d. stainless steel bone plate
表2
使用不同接骨板时重塑过程中骨痂平均弹性模量
Table2.
The average elastic modulus of callus during the reshaping process using different bone plates
4 结论
本文根据理想骨折内固定植入物的设计要求，提出了一种时变刚度骨折内固定植入物复合结构拓扑优化设计方法。通过设计可降解复合结构的构型，实现了对其时变刚度特性的调控设计，并利用骨重塑模拟模型进行了评估。优化的结果显示，高刚度的材料1形成柱状的支撑结构，低刚度的材料2分布在降解边界和内部。降解初期分布在两侧降解边界的材料2优先降解，结构刚度缓慢下降。当两侧的材料2降解完后，材料1开始降解，由于材料1作为主要的支撑，其降解导致结构刚度下降速度加快。因此设计出的可降解时变刚度接骨板可以在骨愈合前期提供良好的固定稳定性，而在愈合后期可以快速降解给予骨痂更多的生物力学刺激。与相同结构钛合金接骨板和不锈钢接骨板相比，可降解时变刚度接骨板能够有效减缓应力遮挡效应，使骨痂取得更好的重塑效果，为基于生物力学的骨折内固定植入物的设计提供了参考。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：孙浩、丁晓红负责算法程序设计、结果分析以及论文的撰写和修改；徐世鹏、段朋云负责数据分析指导及论文审阅修订；熊敏、张横负责算法程序设计的指导。
0 引言
骨折是一种常见的创伤性疾病，内固定术是最为常用且有效的骨折治疗方法。骨折愈合是一个力学调控和生物调控的复杂过程，骨折内固定植入物需要在治疗初期提供足够的支持和稳定性，随着骨愈合的进行，内固定植入物的刚度逐渐降低，实现力流向断骨处转移，促进骨愈合[1-3]。这要求内固定植入物的刚度应随时间发生变化（时变刚度），而使用可降解材料设计的植入物可满足这一要求。参考文献[4]，绘制具有时变刚度的理想可降解植入物的力学性能随骨折愈合过程的变化趋势图（见图1），骨愈合过程可以分为三个阶段：炎症期（阶段1）、骨痂形成期（阶段2）、骨痂重塑期（阶段3）。内固定植入物的时变刚度如虚线所示，在炎症期和骨痂形成期，骨折处血肿逐渐被肉芽组织代替，最终骨化成为编织骨，此时骨痂的弹性模量较低，植入物需要提供良好的固定稳定性；在骨痂重塑期，骨骼逐渐恢复到健康骨骼状态，植入物刚度以较快的速率降低，使应力逐渐向骨痂转移，在骨折完全愈合之后，内固定植入物的刚度接近零。因此如何设计内固定植入物结构，使其刚度特性符合骨愈合规律的要求是一个亟待解决的问题。
图1
理想骨折内固定植入物降解率与力学性能变化趋势
Figure1.
The degradation rate and mechanical trend of ideal fracture internal fixation implants
目前可降解内固定植入物的设计主要通过合金化、冷加工、热处理和表面改性等方法改变可降解材料的降解行为，使其具有与骨愈合相匹配的刚度[5-7]。结构拓扑优化设计方法通过设计结构的材料分布使结构得到所需的力学性能，Wu等[8]通过考虑骨重塑过程的时间依赖性，推导了时域的灵敏度，开发了一种新的下颌骨接骨板拓扑优化方法，但未考虑骨固定板降解；Chandra等[9]考虑到生物降解过程中钢板的尺寸、均匀的降解率和整个过程中钢板的固定稳定性，设计了一种可生物降解的植入钢板，却忽视了材料降解对结构的影响。由于内固定植入物的力学性能随降解时间发生变化，而力学性能的变化会影响骨愈合效果，因此如何模拟降解过程，并量化材料降解对内固定植入物结构力学性能的影响，设计具有时变刚度的内固定植入物使之符合骨愈合规律的要求，具有重要的意义。目前，通过拓扑优化设计方法对可降解内固定植入物进行结构设计以改变可降解材料降解行为的研究仍较少[10-11]。
本文提出一种具有时变刚度的骨折内固定植入物复合结构拓扑优化设计方法，建立考虑结构降解的有限元模型，获得降解过程中中间结构的力学性能，通过结构拓扑优化实现具有不同降解速度和弹性模量的双材料可降解植入物，并利用骨重塑模拟模型对设计结果的时变刚度性能和骨重塑效果进行评估。
1 可降解结构的力学分析
1.1 具有时变刚度特性的复合结构设计构思
目前，单一材料的可降解内固定植入物的刚度很难与骨愈合周期相匹配，会出现降解过快发生固定失效和降解过慢导致应力屏蔽[12-14]，因此考虑采用两种不同的可降解材料。如图2所示，生物降解材料1降解速度慢、刚度大，材料2性能与其完全相反。在给定的边界条件下，在设计域内确定两种材料的最优布局，以实现在早期结构刚度缓慢降低，而后期结构刚度能够以相对较快速率降低。
图2
考虑时变刚度特性的复合结构设计
Figure2.
Design of composite structures considering time-varying stiffness characteristics
1.2 问题描述
进行内固定植入物的拓扑优化设计，模拟结构的降解过程并量化材料降解对结构性能的影响是研究的关键。分别使用相对密度(xe)和降解残留率(ye)描述两种材料的分布和降解状态，如式(1)～(2)所示。通过xe决定单元e为材料1还是材料2；通过ye表示降解状态，ye = 1，代表单元e完全未降解，ye = ymin代表单元e完全降解。
|
|
式中，xmin和ymin取0.001；m为单元个数。
1.3 结构降解模拟-力学分析耦合模型
采用均匀腐蚀方法，以平均降解率模拟结构的降解[15-16]。为简化模拟过程，假设：① 降解只发生在材料的固液交界处且均匀进行，降解边界随材料的降解而移动；② 材料的降解具有均一的降解速率，不考虑其他环境参数（如温度、pH值等）；③ 腐蚀环境相同，忽略贴近骨面和贴近组织的两面对降解速率的影响。需要说明的是，均匀腐蚀与实际腐蚀存在一定的差异，实际腐蚀是一个随机过程，不仅和腐蚀环境有关，而且存在一定的不确定性。但降解实验表明，以均匀腐蚀模拟结构的降解与降解实验中随机腐蚀降解的质量损失率差距在0.1%以内[17]，表明均匀降解模拟方法是可行的。
降解模拟过程如图3所示，单元颜色由深到浅表示材料从未开始降解到完全降解，可降解边界使用黑色粗虚线表示。图3a为结构网格，图3b为降解过程，降解过程示意如图3b1～3b5所示。其中从b1到b3可降解边界处的一层单元格颜色逐渐由深变浅，表示该层单元逐渐降解；从b3到b4，可降解边界从外层单元移动到相邻的内层单元；然后从b4到b5内层单元开始降解。在一个降解时间步完成后，对其降解结果进行结构性能分析，获得结构性能随时间的变化关系。单元降解的数学表达为：
图3
降解过程示意图
a. 结构网格；b. 降解过程；c. 单元
a. meshes of the structure; b. simulation of structure degradation; c. the location mark of element
|
式中，i表示降解时间步，de表示单元e的降解速率，t为单位降解时间步的时间。
单元e和相邻单元位置标记如图3c所示，单元e的位置为(k, l)。通过判断降解单元相邻单元的状态，来决定该单元是否发生降解，为了保证数值计算的稳定，设定yb=0.1，当相邻单元的残留率小于yb时可沿该单元方向进行降解，如当yk − 1,l < yb时，认为单元e能够从左侧开始降解，当yk + 1,l < yb时，认为单元e可以从右侧开始降解，上下方向类似[18]。通过对相邻单元降解状态的判断，单元e的降解速率可以表示为：
|
式中，Ω表示e相邻单元的集合；α表示Ω中的第α个单元；d为材料的降解速率；Nα由式(5)得到：
|
为了便于求导，使用Heaviside方程将式(4)近似为连续方程，单元的降解方程可以改写为：
|
阶跃函数H为：
|
式中，γ是一个常数。
由式(3)和式(6)，单元e残留率的更新方程为：
|
使用Kreisselmeier-Steinhauser函数将式(8)最大函数问题近似为连续方程[19-20]：
|
式中，γ = 6，η = 6。
结构平衡方程如式(10)所示：
|
式中，K是结构的整体刚度矩阵，由于材料的降解，K与时间相关；U为位移向量；F为力向量。
对于可降解结构，K由单元刚度矩阵ke与ye组装得到，ye是关于降解速率d和降解时间i的函数。将材料降解对结构力学性能的影响耦合到有限元分析过程中，可降解结构的刚度矩阵为：
|
式中，B表示几何矩阵；D为弹性系数矩阵；Ωe为单元的体积。
降解模拟算法包括以下步骤。步骤1：初始化残留率ye，使用k和l对单元e进行标记；步骤2：判断相邻单元的降解状态；步骤3：使用式(6)计算单元e的降解速率de；步骤4：使用式(9)计算单元e的残留率ye；步骤5：根据ye，使用式(10)进行有限元分析；重复步骤1～5，直到指定的时间步td。
2 双材料可降解结构的拓扑优化
优化设计的目的是通过确定两种具有不同可降解性能的材料的最佳分布，使复合结构在指定降解时间步中间结构的刚度之和最大。使用结构柔顺度表征结构刚度，可降解复合结构柔顺度最小的优化数学模型可以表示为：
|
式中，X为设计变量的集合；C为结构柔顺度，ʌ为降解时间步的集合；v1为复合结构中材料1的最终体积，f是给定的体积分数，v0为总体积。
为获得轮廓清晰的优化结构，使用基于固体各向同性材料惩罚模型（solid isotropic material with penalization，SIMP）的材料插值机制：
|
|
式中，E为弹性模量；p为惩罚系数，取3.0；下标1、2分别表示可降解材料1、2。
2.1 灵敏度分析
在每个降解时间步对结构进行有限元分析，计算结构柔顺度，其灵敏度可以通过求导得到。目标函数对于设计变量X的灵敏度可以表示为：
|
结构柔顺度C是关于xe、ye和i的函数：
|
式中，u表示单元位移矩阵；k0为单元刚度矩阵。
根据式(16)定义的柔顺度，在降解时间步i对设计变量xe的导数为：
|
其中：
|
|
|
材料1的体积为：
|
则材料1的体积对于设计变量xe的导数为：
|
为克服棋盘格效应，采用基于Heaviside的灵敏度过滤方法使密度变量的值趋近0或1。使用如式(23)所示的Heaviside方程对中间密度场xei进行过滤，xei由式(24)计算得到。
|
|
式中，β为常数，控制函数的平滑度，如果β = 0，则是原始的密度过滤器，如果β接近于无穷大，则Heaviside为阶跃函数。本文取β的初始值为1，每100次迭代对β的值加倍，最大值为512。
2.2 优化设计流程
双材料可降解结构的拓扑优化设计流程如图4所示，主要步骤如下。
图4
双材料可降解结构的拓扑优化设计流程
Figure4.
Topological optimization design process for dual material degradable structures
第一步：构造初始几何模型并进行离散化。定义数学模型的参数，给定体积分数和过滤半径，设置xe和ye的初始值，并定义可降解边界、载荷和约束等边界条件；第二步：根据式(9)进行降解模拟，计算得到残留率；第三步：根据式(10)和(11)对考虑材料降解的复合结构进行有限元分析；第四步：根据式(15)进行灵敏度分析。重复第二步到第四步，直到指定的时间步。第五步：利用MMA优化算法更新设计变量；第六步：重复第二步到第五步，直到满足约束条件或目标函数收敛。当满足以下两个条件之一时，迭代终止：① 连续两次迭代的结构柔顺度变化小于收敛判据；② 迭代步数达到指定值，本文中取100。
3 胫骨接骨板设计
3.1 设计模型及材料性能
使用如图5a所示的胫骨-接骨板系统简化模型，接骨板的尺寸如图5b所示，接骨板中间区域（橘黄色）为可降解结构，两侧（灰色）为不可降解结构。接骨板植入人体后，可降解结构垂直于x轴和z轴的四个面都能够发生降解，其中垂直于x方向两个面的面积较大，这两个面的降解对结构刚度的影响也较大。为了将设计问题简化，我们假定可降解结构仅能从垂直于x方向的两个面发生降解，忽略垂直于z方向两个面的降解。如图5c所示，将垂直于z方向的面作为设计平面，利用考虑时变刚度特性的可降解复合结构拓扑优化设计方法对设计平面进行优化设计，然后沿z方向拉伸得到三维可降解复合结构，实现可降解接骨板的拓扑优化设计。
图5
胫骨-接骨板固定系统
a. 设计模型；b. 接骨板模型；c. 设计区域；d. 接骨板力学边界条件；e. 接骨板的简化模型；f. 降解边界条件
Figure5. Bone and bone plate fixation systema. design model; b. bone plate model; c. design area; d. mechanical boundary conditions of bone plates; e. simplified model of bone plate; f. degradation boundary conditions
如图5d所示，胫骨平台的力学边界条件可以简化为底端全约束，顶端皮质骨受到垂直向下的载荷。骨愈合过程中，骨痂的刚度逐渐增大，但都远小于接骨板的刚度，由于载荷分流，接骨板承担大部分的载荷，设计平面的力学边界条件简化为图5e所示，中心点全约束，中心线其余节点约束竖直方向自由度，上下边受到垂直方向的压力，F1 = 1 600 N。可降解边界条件如图5f所示，上下侧为不可降解边界，左右侧为可降解边界。根据文献[12]，材料属性如表1所示。
3.2 拓扑优化结果及其时变刚度特性分析
采用2.2小节的设计流程，使用尺寸为0.1 mm的四节点四边形平面单元对设计域进行网格划分，拓扑优化过滤半径设置为单元尺寸的两倍，材料1的体积分数为0.5，优化目标为第1和第90天中间结构柔顺度之和最小。优化迭代历程如图6a所示，随着迭代的进行逐渐获得具有清晰边界的材料分布轮廓，优化目标达到了收敛状态。拓扑优化设计结果如图6b所示，高刚度材料1（黑色）形成了柱状的支撑结构，在可降解边界处外层为材料2，结构左右呈现对称性，高刚度材料1（黑色）主要分布在左右侧可降解边界处，在左右侧中部形成V型开口，并在靠近加载侧形成了类似X型的交叉结构，连接左右两侧的高刚度材料，提高了结构整体的稳定性。时变刚度接骨板设计结果如图6c所示，两侧为惰性金属，中间区域为可降解复合结构。需要说明的是，结构的最优拓扑形态主要与所受的载荷相关，本文采用了承受偏心压力的圆柱模拟胫骨，考虑了实际胫骨承受的主要载荷为压力和弯矩，因此图6的设计结果具有合理性。
图6
拓扑优化设计结果
a. 拓扑优化迭代历程；b. 二维设计结果；c. 可降解时变刚度接骨板
Figure6. Topology optimization design resulta. iteration history with two degradable interfaces; b. 2D design results; c. degradable time-varying stiffness plate
考虑时变刚度特性拓扑优化设计结果的降解过程如图7所示。由于材料1刚度较大，其降解对结构刚度的影响程度较大。降解开始时，可降解边界附近的材料2先逐渐降解，结构刚度缓慢下降。当外层的材料2降解完，结构内部的材料1开始降解，结构刚度下降的速度开始增大。
图7
考虑时变特性拓扑优化设计结果的降解过程
Figure7.
Degradation process of topology optimization design results considering time-varying characteristics
3.3 骨重塑效果评估
骨重塑是骨骼受到外部生物力学刺激时产生的一种自适应过程。骨重塑模拟模型能够描述在生物力学刺激作用下骨密度随时间的变化过程。利用应变能密度预测骨重塑过程[21-22]，力学刺激Ξ可以表示为：
|
式中，U为应变能密度，ρ为骨密度。骨密度在一定时间间隔t的变化量ρ可以表示为：
|
|
|
式中，t为时间；R1和R2为常数，分别取18 (g/cm3)2/(J/月)和600 (g/cm3)3/(J2/月)；θ1和θ2是力学刺激的阈值，分别由式(26)和(27)计算得到；s和δ为常数，分别取0.000 36 J/(g*cm3)和0.1。骨密度变化量与力学刺激之间的曲线描述如图8所示，可以分为低载退化、惰性、骨重塑和过载退化四个区间。在低载退化和过载退化区间骨密度减小，在惰性区间骨密度不变，在骨重塑区间骨密度增大。
图8
骨密度变化率与力学刺激关系示意图
Figure8.
The schematic of relationship between the rate of change in bone density and mechanical stimulus
假设发生骨重塑的骨材料属性为各向同性，其弹性模量可以通过密度进行量化，经验公式为：
|
式中ζ为弹性模量；τ为常数，取3 790 MPa/(g/cm3)3。
使用图5a所示的简化模型模拟骨痂的重塑过程。边界条件如图5d所示，顶端受到垂直向下的载荷，底端全约束。载荷为普通人体重（70 kg）的三倍，即F = 2 100 N。骨痂的初始密度为1.0 g/cm3，最大为1.74 g/cm3，最小为0.1 g/cm3。接骨板除设计域外结构和螺钉为钛合金，弹性模量为110 000 MPa。皮质骨和松质骨的弹性模量分别为15 750 MPa和1 100 MPa。不锈钢的弹性模量为210 000 MPa。
采用弹性模量表征骨重塑过程的骨痂单元变化，平均弹性模量一定程度上能够反映骨痂整体结构刚度的变化。骨痂的尺寸和初始状态如图9a所示，图9b、c、d分别为使用可降解时变刚度接骨板、相同结构尺寸的钛合金接骨板和不锈钢接骨板固定时，经过6个月和11个月时的骨痂弹性模量分布云图。接骨板位于骨痂的左侧，由于接骨板造成的应力遮挡效应，靠近接骨板侧的骨痂受到的力学刺激较小，因此弹性模量较低。与使用钛合金接骨板和不锈钢接骨板相比，使用可降解时变刚度接骨板时，骨痂重塑过程中左侧弹性模量小于3 790 MPa的区域明显较小。经过11个月重塑之后，使用可降解时变刚度接骨板时，骨痂弹性模量基本都大于3 790 MPa；而使用钛合金接骨板和不锈钢接骨板时，骨痂左侧仍有部分区域弹性模量小于3 790 MPa。结果表明，使用可降解时变刚度接骨板可以使靠近接骨板侧的骨痂在重塑过程中受到更多的力学刺激，有效地减缓应力遮挡效应，取得更好的重塑效果。对比相同重塑时间下使用三种接骨板的骨痂平均弹性模量的差值，结果如表2所示，随着重塑的进行，三种骨板对应的骨痂的平均弹性模量均逐渐增大，使用可降解时变刚度接骨板时的骨痂平均弹性模量均高于钛合金接骨板和不锈钢接骨板，对比不同骨板的骨痂平均弹性模量的差值，随着时间的增加差值逐渐变大，并且差值增大的速率也逐渐增大，说明可降解时变刚度接骨板的结构刚度减小的速率逐渐增大。
图9
骨痂弹性模量分布云图
a. 初始状态；b. 可降解时变刚度接骨板；c. 钛合金接骨；d. 不锈钢接骨板
Figure9. Distribution cloud map of elastic modulus of callusa. initial state; b. degradable time-varying stiffness bone plate; c. titanium alloy bone plate; d. stainless steel bone plate
表2
使用不同接骨板时重塑过程中骨痂平均弹性模量
Table2.
The average elastic modulus of callus during the reshaping process using different bone plates
4 结论
本文根据理想骨折内固定植入物的设计要求，提出了一种时变刚度骨折内固定植入物复合结构拓扑优化设计方法。通过设计可降解复合结构的构型，实现了对其时变刚度特性的调控设计，并利用骨重塑模拟模型进行了评估。优化的结果显示，高刚度的材料1形成柱状的支撑结构，低刚度的材料2分布在降解边界和内部。降解初期分布在两侧降解边界的材料2优先降解，结构刚度缓慢下降。当两侧的材料2降解完后，材料1开始降解，由于材料1作为主要的支撑，其降解导致结构刚度下降速度加快。因此设计出的可降解时变刚度接骨板可以在骨愈合前期提供良好的固定稳定性，而在愈合后期可以快速降解给予骨痂更多的生物力学刺激。与相同结构钛合金接骨板和不锈钢接骨板相比，可降解时变刚度接骨板能够有效减缓应力遮挡效应，使骨痂取得更好的重塑效果，为基于生物力学的骨折内固定植入物的设计提供了参考。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：孙浩、丁晓红负责算法程序设计、结果分析以及论文的撰写和修改；徐世鹏、段朋云负责数据分析指导及论文审阅修订；熊敏、张横负责算法程序设计的指导。
