ABAQUS二次开发在骨重建仿真中的应用

井 野，吕林蔚*，宋 阳，张春秋

（天津理工大学a.天津市先进机电系统设计与智能控制重点实验室，b.机电工程国家级实验教学示范中心，天津300384）

骨骼是支撑人体结构的重要器官，其主要功能是为人体承担各种外力载荷[1]。WOLFF和JOUOS[2]认为骨具有力学功能适应性，应力分布与骨骼的结构之间存在一定联系。骨组织可以通过重建作用，根据外力载荷的不同，改变其自身结构，以适应新的力学环境。根据WOLFF定律[3-4]，骨重建是通过骨形成和骨吸收调节骨组织的含量并试图用最少的骨量来获得最大的结构刚度，以满足力学性能要求的过程。通过这种方式，对骨组织结构不断进行优化，防止骨骼中微损伤的积累。并在一定程度上也会引起骨宏观尺度力学特性（例如密度和刚度）的变化。

自从WOLFF定律被提出以来，通过建立骨重建的数学模型，定量地模拟力学环境对骨量调节的影响过程，逐渐成为了研究重点[5]。现有的骨重建数值仿真研究，主要是以MULLENDER等[6]根据骨功能适应性理论提出的自适应模型以及WEINANS等[7]改进Fyhrie和FYHRIE[8]所建立的骨自优化方程为基础。其中自适应模型是一个可以依据所受力学激励与稳态之间的差值而进行调整和改变的反馈系统。而骨自优化方程实际上是根据外部载荷情况，通过调整密度，以获得单位骨量应变能密度（strain energy density，SED）的恒定预设值的目标函数。GONG等[9]通过高阶非线性方程来调控骨重建过程，将定量骨重建理论拓展到椎体边缘骨质增生的形成，证明了骨质增生等骨形态异常也符合自适应重建过程，骨形态异常与力学环境有关。此外，在许多研究中[10-11]，表观密度被认为是确定骨组织力学性能及影响骨重建过程的重要变量，可以用其描述及评价骨强度和刚度的变化。并且由于以前的研究已经证明了维持骨量与骨组织中的局部应变值之间有密切的联系，因此也有学者[12]利用最小化局部应变能密度的求解，来描述骨重建过程。

随着计算机技术的发展，有限元法（finite element method，FEM）成为骨生物力学研究的重要手段[13]。KUMAR等[14]基于ANSYS软件，分析了正常行走、站立、跑步和跳跃活动中的股骨应力情况，用来预测骨折风险及选择人造骨材料。YANG等[15]采用不同的弹性模量与密度关系，研究二维股骨近端的骨重建情况，发现两者参数的选择会对骨重建结果产生影响。

ABAQUS软件是功能强大的有限元分析工具，具有极强的通用性和仿真计算能力[16]。本研究通过Python脚本语言对ABAQUS软件进行二次开发，把基于应变能密度理论的骨重建自优化方程进行编程，建立骨组织重建模型，来分析外部载荷作用下骨结构的变化，预测骨组织的密度分布。希望证明基于Python脚本语言的骨重建程序算法是一种可靠的有限元计算手段，为今后的骨重建数值仿真研究提供新的方法。

1 理论与方法

1.1 骨自优化理论

本研究的骨重建数值仿真是以有限元技术为手段，以应变能准则为基础的骨自优化理论，模拟骨组织受力重建过程。总体思路是通过骨组织表观密度的变化率控制其重建，利用有限元法迭代计算出模型各个节点力学激励数值，与参考值比较，判断是否达到平衡态。骨重建数值模拟流程如图1所示。

图1 骨重建数值模拟流程图Fig.1 Flow chart of numerical simulation of bone remodeling

根据FROST[17]提出的骨组织中存在着“力学调控稳态机制”，当施加外部载荷使骨组织内部的力学环境发生改变时，骨细胞原本处于的最佳生理平衡状态即被破坏，而应力、应变及应变能密度等力学特性因素会受到其影响，与平衡时的稳态值存在差值，此时骨组织将启动重建程序，使其再次重新回到平衡状态。FYHRIE和CARTER[18]利用数学方法将Frost的力学稳态理论参数化，提出了骨自优化控制方程，用来描述骨组织表观密度变化率为：

式中，ρ为骨组织表观密度，可以表征骨组织内部的结构特性；B为骨重建变化率的时间常数，与骨组织本身有关；S为力学激励，可以选择包括等效应力、等效应变及应变能密度在内的力学特性因素；K为力学稳态的参考值常数。

WEINANS等[19]将前面学者的相关研究统一为骨自优化理论，把单位质量的应变能密度Un/ρ作为式（1）中的力学激励S为：

式中，Un为应变能密度，E为骨的表观弹性模量，C为模量-密度常数，r为影响指数，ρcb为皮质骨的表观密度，假定为最大骨密度。

上述骨自优化理论的数学表达是本文研究的编程依据，是定量研究骨重建过程的理论基础。

1.2 ABAQUS的二次开发

本研究对ABAQUS进行二次开发，编写骨重建算法程序，为提高有限元仿真的效率。

1.2.1 前处理部分

前处理部分的创建部件、定义装配件、设置分析功能、定义载荷及边界条件和划分网格在ABAQUS主界面中完成。将材料属性参数编译成程序，创建截面属性，并将其赋予给部件。根据骨自优化理论数学表达式（3），在每次骨重建迭代计算之前设定密度、弹性模量和泊松比等材料属性。将初始密度假设为1.74 g/cm3，弹性模量与密度有关，取泊松比为0.3。

1.2.2 计算过程

前处理过程完成后，执行作业模块，创建作业、提交运算及等待计算结束。根据骨自优化控制方程，将力学激励S设置为Un/ρ，即骨单位质量的应变能密度，取参考值K为0.153 7，时间常数B设置为1[20]，则目标函数转化为：

密度的变化率随骨组织内部应力的分布情况而改变，最终目标是令密度的变化率dρ/dt=0，即调整单位大小的应变能密度达到平衡态预设的K值。满足平衡条件，迭代优化结束，计算完成。

当模型的单元数量较多时，有限元计算量大，计算时间较长，为了提高工作效率，可以采用并行计算方式，在ABAQUS软件作业管理器中选择使用多个处理器进行计算。本研究采用双核处理器进行计算。

1.2.3 后处理部分

打开得到结果的数据文件，通过访问对象数据库（objectdatabase，ODB）文件的子对象分析步及增量步，查寻并提取应变能密度数据计算其结果的改变量。将密度添加至结果变量，以场变量的形式输出。

2 有限元模型

2.1 二维悬臂梁模型

探究初始密度值对骨自优化过程的影响与基于Wolff定律的骨重建研究中均认为二维悬臂梁模型在一定程度上可以模拟长骨结构，利用其受力后密度分布情况来预测骨自优化的最终结构形态是一种可行的方法[20-21]。

本研究选取悬臂梁模型，尺寸参数为长80 mm，宽10 mm，外力载荷取5.0 N的集中力，施加在悬臂梁竖直方向右侧中点，并在其左端施加完全约束[20]，二维悬臂梁结构模型如图2所示。通过与文献[20-21]实验结果的对比，以验证本算法的正确性。

图2 二维悬臂梁结构模型Fig.2 Two-dimensional cantilever beam structuremodel

2.2 股骨近端模型

对股骨近端区域进行定量计算机断层扫描（quantitative computerized tomography，QCT），扫描参数为像素值0.937 5 mm，分辨率512×512，层厚2.500 0 mm，电压80 kV，电流260 mAs，获取56张股骨近端QCT影像。利用QCT影像在（materialise′s interactivemedical imagecontrol system，MIMICS）软件中建模，并用MAGICS切取一个股骨冠状面，将其导入ABAQUS，得到股骨近端有限元模型如图3所示。

图3 股骨近端有限元模型Fig.3 Finiteelement model of proximal femur

本研究假定骨组织是各向同性材料，密度与弹性模量之间的关系式为E=3 790ρ3，泊松比为0.3。最大表观密度ρmax设为1.74 g/cm3，最小表观密度ρmin设为0.01 g/cm3，所有元素的初始密度值为1.74 g/cm3。分别在股骨关节处及外展肌的附着区域施加4种外力，股骨近端模型载荷条件如表1所示，模型底部所有节点施加完全约束。

表1 股骨近端模型载荷条件Tab.1 Loading conditionsof theproximal femur model

3 结果

3.1 二维悬臂梁模型

悬臂梁模型采用带缩减积分的四节点平面应力（CPS4R）单元，网格尺寸为1.00 mm，划分800个单元，进行有限元计算。规定初始密度为1.74 g/cm3，提取密度变量输出，迭代200次后悬臂梁密度分布结果如图4所示。

图4 迭代200次后悬臂梁密度分布结果Fig.4 Density distribution results of cantilever after 200 iterations

从密度分布结果来看，当迭代200次时，沿长轴方向最外侧单元密度分布均最高，而对称轴线两侧向内单元密度按一定比例逐渐降低，悬臂梁被优化成明显的桥梁工程中常见的桁架结构，与利用ANAYS算法探究初始密度对骨自优化结果的影响和基于Wolff定律的骨重建研究的计算结果一致[20-21]，证明了程序算法的正确性。

3.2 股骨近端模型

二维股骨有限元模型采用CPS4R单元，网格尺寸为2.00 mm，划分3 492个单元，施加载荷，进行有限元计算。求解过程迭代50次，模拟50天骨重建过程。在每次迭代中计算每个单元的骨组织密度，提取密度数据并得到股骨近端密度分布结果，如图5所示。可以看出，在股骨表面施加外力载荷组，随着迭代进行，股骨近端密度分布呈现出3个特征。

图5 股骨近端密度分布结果Fig.5 Density distribution result of proximal femur

1）外层的股骨干表面骨密度逐渐增强，形成皮质骨区域，最大密度出现在该处。

2）股骨头及大转子部位密度减小，形成松质骨区域，密度值在0.442～1.308 g/cm3。

3）股骨干髓腔趋近于0，呈现为“空腔”，Ward三角区也出现了密度最低值，因为该处密度值低，导致出现股骨颈骨折的风险较高。

股骨近端应力分布结果如图6所示。根据WOLFF定律，骨组织在需要的部位就生长，在不需要的区域就吸收[2]，对比图5和图6可以发现，当骨承受较大的载荷时，该区域的应力会较大，骨密度需要增大，反之，应力小的部位，骨密度则较小。

图6 股骨近端应力分布结果Fig.6 Stressdistribution result of proximal femur

4 讨论

本研究利用ABAQUS软件的二次开发环境，结合现有骨自优化理论，编写骨程序算法，建立二维悬臂梁模型和股骨近端模型，验证骨重建程序算法的正确性及准确性。

将初始密度为1.74 g/cm3悬臂梁模型进行200次迭代计算，与宫赫和朱兴华[20]研究初始密度值对骨自优化结果的影响、马宗民等[21]基于WOLFF定律的骨重建研究，所建立悬臂梁模型的计算结果一致，两者均采用ANSYS软件进行仿真。从密度分布结果来看，沿长轴方向最外侧单元密度分布均最高，沿其轴线两侧向内单元的密度成比例逐渐降低。显然，这种情况符合骨自优化理论及Wolff定律，试图用最少的骨量来获得最大的结构刚度以满足力学要求，说明本研究编写的程序算法是正确可靠的。

股骨近端模型的密度分布表明，当关节力和外展肌载荷组施加在股骨上，外力由股骨上端传递到股骨干时，股骨干与股骨头内部的密度大小表现出了明显差异。股骨干表面逐渐变“硬”，呈现为皮质骨，最大密度主要出现于该处，股骨头、转子部位的密度较小，呈现为松质骨，而髓腔内及Ward三角区的骨密度趋近于0，是易引起骨折的脆弱区。所得的密度分布情况，可以表现真实股骨的结构特征，载荷分配较多的区域，需要的骨组织更致密，密度较高，载荷分配较少的区域，用较少的骨量维持结构，说明模拟结果遵循了WOLFF定律。以应变能密度、等效应力作为激励变量，利用ANSYS软件数值模拟骨重建过程也得到了与本研究结果相一致的密度分布[22]。结果表明，股骨近端会受所处力学环境的影响，通过骨重建作用，调整、优化其自身结构，证明了所编写骨重建算法的准确性。

本研究存在一定的局限性。首先是二维模型的局限性，使用二维模型减少了计算成本，但平面模型不能完全地模拟骨组织的真实结构及受力情况。其次，由于实际中骨骼是各向异性材料，且人体真实受力情况十分复杂，有限元分析的工况与真实力学环境存在差异。下一步将把算法拓展，进行三维模型以及不同加载条件下的骨重建研究。但总体上看，股骨近端模型的仿真结果得到的骨组织密度分布情况，与其他学者使用ANSYS软件的计算结果一致，符合骨组织调节自身结构适应力学环境的规律。

5 结论

综上所述，本研究利用Python语言对ABAQUS进行二次开发，依据骨自优化理论，以单位质量的应变能密度作为参考变量，建立以密度值变化率为目标函数的程序算法，对模型进行有限元仿真。通过悬臂梁模型计算得到的桁架结构，证明程序算法可以满足WOLFF定律。股骨近端模型的计算结果能够反映出股骨内部骨组织的主要结构特征，说明程序算法模拟骨重建过程符合股骨实际情况。希望本研究的方法及结论可以为骨生物力学仿真研究提供新思路。