各向异性松质骨压缩模量和强度的数值预测

卢子兴 裴 鹤

(北京航空航天大学 航空科学与工程学院，北京 100191)

引言

松质骨的力学性质主要体现为表观强度和刚度，而这些性质又与松质骨的微结构特性密切相关，如骨小梁形状、生长方向、密度和连通性等;并且对刚度而言，还与组织密度有关［1］。已有的实验结果表明，骨小梁的细化和连通性的损失均会造成整骨强度和刚度在准静态和冲击载荷下的退化［2－4］。因此，研究微结构和组织密度特性随骨小梁尺度的变化及与松质骨表观力学特性之间的关系具有重要的意义。然而，这种关系却具有高度的非线性，且随病患个体微结构的变化而改变［1，5］，很难用数学公式表达出来，故对其行为进行有限元分析成为一种必要的手段。例如，Guo等人应用由梁单元构成的14面体模型对松质骨微结构进行了模拟［5］，计算结果比较符合松质骨的宏观力学性能，但使用梁单元并不能很好地反映骨小梁的形态学特征;而Diamant等利用变截面骨小梁模型组成的简单立方结构对松质骨进行研究，得出受力情况下松质骨微结构细微处的应力状态，但简单立方结构显然不能很好地模拟松质骨的微结构［6－7］。本研究发展了一种有限元参数化建模方法，用14面体模型来模拟棒状骨小梁构成的松质骨微结构，对各向异性松质骨的弹性性能和强度进行了有限元计算分析，给出了压缩模量与实验结果的比较，并对松质骨的表观压缩强度进行了数值预测。

1 松质骨的有限元模型

为评估松质骨力学行为，参考前人的相关工作，建立了对应不同形态学参数的松质骨有限元单胞模型，如图1所示。其中，骨小梁使用Gefen等人基于大量实验数据提出的单根骨小梁模型［7］，该模型比简单圆棒形骨小梁模型能更好地表达骨小梁的形态学特性，即Tb.Th(名义骨小梁厚度)和Tb.Sp(用骨小梁长度L表征)。下式给出建立这种骨小梁模型的形状控制方程，即

式中，骨小梁长度 L沿 z方向，r(z)是 z处的半径，Tb.Th是最大厚度(一般是两终端处厚度)与最小厚度(即通常是小梁中间部位厚度)的平均值，而经验常数取 α=1.373 6，β=40.9 μm，由多次实验的均值确定［6］。

图1 模拟松质骨微结构的14面体单胞实体模型。(a)立体图;(b)平面图Fig.1 The tetrakaidecahedral unit cell model for simulating the microstructure of cancellous bone. (a)threedimensional sketch; (b)twodimensional sketch

有限元模型的形状控制参数为:Tb.Th、b(横向骨小梁长度)和l(纵向骨小梁长度)。模型的各向异性比定义为14面体单胞的高度H和宽度D的比值，即λ=H/D，可以证明它仅与l/b有关。在有限元计算时，骨小梁组织材料的性能参数采用文献［8］中得出的皮质骨的平均压缩力学性能，即弹性模量E0为12 GPa，泊松比为0.3，极限压缩强度为107 MPa。文中模型全部采用四面体实体单元SOLID45划分。为使有限元模型在边界面节点位置上符合周期性边界条件［9］的要求，划分单元时首先划分边界面处单元，相对的边界面单元通过复制保持相同，使得节点能顺利耦合，划分完边界面后再划分体内单元，通过控制单元尺寸使计算达到一定的精度并具有较高的效率。

2 边界条件及计算方法

2.1 边界条件

针对14面体模型能周期性排满整个空间的特点，在计算时施加周期性边界条件［9］，这样能够保证模型在变形后仍然保持周期性的空间排列，并允许边界上的面发生翘曲变形。按照施加周期性边界条件的要求，模型相对表面上对应节点在该平面的法线上以同样的位移膨胀或收缩，在其他方向上则保持一致的平动位移，且在所有方向上均具有相同的转角。

2.2 计算方法

当14面体按体心立方排列时，能周期性地排满整个空间，故采用单倍胞体计算能提高效率。笔者对文献［10］提出的计算周期性排列单胞力学性能的边界条件进行适当变化，得到了计算压缩模量的方法(限于篇幅，这里不再赘述)。下面简要介绍压缩极限强度的计算方法。图2给出了对应不同载荷步时模型的Von Mises应力云图及其网格图。从中可以看出，当模型承受y向压缩时(椎体松质骨通常主要承受纵向压缩载荷，其他方向承受载荷较小，故主要讨论模型沿y向的承载能力)，纵向骨小梁连接处发生了应力集中。在进行应力判断后，这个位置网格的Mises应力超过了材料的极限压缩应力，从而将此处的单元“杀死”，得到如图2中所示的“生”单元模型。在此基础上增加位移载荷，进行下一步计算，又一次可得到模型的Mises应力云图，再继续进行下一次的判断，如此往复，最终可以得到模型的压缩应力-应变曲线和对应的极限压缩强度 σult。

另外，由文献［11－12］可知，松质骨的表观模量和表观强度之间有极强的相关性。当预测了松质骨的表观模量后，可根据 Hou等提出的线性关系［12］预测其表观极限压缩强度，即

图2 有限元计算不同载荷步时模型的应力云图及网格图。(a)Mises应力分布;(b)“杀死”单元;(c)“生”单元;(d)继续加载后的Mises应力分布;(e)继续“杀死”单元;(f)继续“生”单元Fig.2 The stress nephogram and gridding map of mode by the FEM under different loading steps.(a)Mises stress under one loading step;(b)killing elements;(c)reactivating elements;(d)Mises stress under another loading step;(e)killing elements again;(f)reactivating elements again

而Crawford等则由实验确定了松质骨极限压缩强度与屈服强度之间的关系［13］，即

因此，经过有限元计算，可以得到松质骨的极限压缩强度σult;通过式(3)，便可得到松质骨的屈服强度。

3 数值结果及讨论

3.1 骨体积分数与名义骨小梁厚度的关系

图3给出了BV/TV(骨体积分数)与Tb.Th(名义骨小梁厚度)之间的关系。从中可以看出:BV/TV随着Tb.Th的增加而明显增大，故通过控制Tb.Th可以控制BV/TV的变化;l/b对BV/TV有较大影响，在Tb.Th相同的情况下，随 l/b增大，BV/TV逐渐减小;此外，通过改变 l/b，还可以控制 λ的大小。

3.2 压缩模量预测及同实验结果的对比

图3 BV/TV与Tb.Th的关系(b=1.2 mm)Fig.3 The relations between BV/TV and Tb.Th

图4给出了相对压缩模量Eyy/E0与BV/TV的关系及同实验结果的对比。从中可以看出:一是随着BV/TV的增加，压缩模量随之单调增加;与前人所得腰椎松质骨压缩实验［14］(取相近Tb.Th及b的数据)相比，实验值略高于计算值，说明了此模型的有效性;二是压缩模量与BV/TV间表现较好的平方律关系，与 Carter和 Hayes给出的规律［1］相符。因此，该模型能为分析个体骨质疏松患者的骨骼质量提供一定的理论依据。需指出，有限元模型计算得到的模量值基本上不随骨小梁长度变化，而真实松质骨的力学性能与骨小梁长度并不是无关的。为此，Grant等在大量实验数据的基础上，提出了针对椎体松质骨的模量修正因子αE及屈服强度修正因子 ασ的概念［14］，即

式中Et和Em以及σyt和σym分别代表模量和屈服强度的真值和测量值。研究发现，αE对Tb.Th和BV/TV并不敏感，而对Tb.Sp(文中为l+b/2)存在一定的线性相关性，且 αE与 ασ之间也存在线性关系。基于上述修正因子的概念和类似的关系式，笔者对有限元计算得到的压缩模量进行了修正，给出的预测结果如图4所示。

图4 Eyy/E0随BV/TV的变化Fig.4 The variation of Eyy/E0with BV/TV

3.3 压缩应力-应变曲线

图5 压缩应力-应变曲线Fig.5 The curves of compressive stress and strain

图5给出了 b=0.9 mm、l/b=1.5、BV/TV=10.3%时模型在单向受压情况下的应力-应变曲线，同时给出了一个由压缩实验得到的椎体松质骨样本的压缩应力-应变曲线。可见，两者在弹性阶段基本一致，但实验测得的强度要大于有限元计算的结果(计算极限压缩强度为1.9 MPa)。这是因为真实松质骨的骨组织中含有丰富的结缔组织，而不仅仅是无机物，从而增强了松质骨的韧性，使得屈服强度与极限强度较之纯骨组织有所提高。显然，对于特定的 l/b(如取1.0，1.25，1.50和1.75)，通过改变 BV/TV，可以得到一系列的压缩应力-应变曲线(限于篇幅，这里仅给出 l/b=1.50对应的结果)，如图6所示。从中可以看出，随着BV/TV的增加，极限压缩强度逐渐增大，但所对应的应变却逐渐减小。

图6 不同 BV/TV下的压缩应力-应变曲线(l/b=1.50)Fig.6 The curves of compressive stress and strain under the different values of BV/TV(l/b=1.50)

3.4 极限压缩强度预测与实验值的比较

如上所述，针对 l/b=1.0，1.25，1.50和 1.75几种情况，再通过改变 BV/TV，可以得到一系列的压缩应力-应变曲线。从中提取对应每种情况的极限压缩强度σult，并用强度修正因子对其修正，获得的 σult与 BV/TV 的关系及同实验值［12，15］的比较，如图7所示。从中可以看出:一是当 BV/TV较低时(BV/TV<0.15)，l/b值对 σult存在明显影响，相同BV/TV时，σult随着l/b的增加而增大，这表明，较低BV/TV时，松质骨呈现棒状骨小梁结构，纵向骨小梁是主要承载者。二是当BV/TV继续增大时，对应不同l/b值的曲线逐渐汇聚;当BV/TV增大到一定程度后，将出现板状骨小梁(这里不予讨论)。三是l/b值接近1且BV/TV较低时，预测的σult能较好地与实验结果相符合;当BV/TV增大时，松质骨内出现板状骨小梁，使松质骨抗压能力得到增强，因此实验测得的 σult比使用棒状骨小梁模型预测的 σult偏大，这符合真实情况。

图7 σult与BV/TV的关系及同实验值的比较Fig.7 The relations between σultand BV/TV and the comparison with the experimental values

3.5 极限压缩强度预测与经验公式预测的比较

Hou等在大量实验数据的基础上，提出了预测σult的经验公式［12］。利用不同 l/b情况下纵向表观压缩模量Eyy的预测结果，代入式(2)，得到极限应力值(记为经验公式预测)与上述有限元计算的结果进行对比得到图8。从中可以看出:一是有限元计算和经验公式均能对松质骨的σult做出较好的预测;二是经验公式预测不能反映各向异性参数 l/b对强度预测的影响，而有限元计算则可反映各向异性参数l/b对强度预测的影响，l/b的影响从一定程度上说明了σult实验数据分散的合理性。

图8 σult有限元计算结果与经验公式对比Fig.8 The comparison between the FE results of σult and the empirical formula

4 结论

1)14面体单胞模型能有效地用于松质骨压缩模量和强度的数值预测，预测结果与实验数据及经验公式相吻合;

2)相对于利用模量对强度进行预测的经验公式而言，本研究的有限元方法能够反映各向异性对松质骨强度的影响，从而在一定程度上解释了实验中强度数据的分散性;

3)由模型受力云图可以判断破坏发生的位置，从而为今后进行断裂方面的分析奠定基础。

［1］Carter DR，Hayes WC.The compressive behavior of bone as a two-phase porous structure ［J］.Bone Joint Surg，1977，59(2):954－962.

［2］Passi N，Gefen A.Trabecular bone contributes to strength of the proximal femur under mediolateral impact in the avian ［J］.Biomech Eng，2005，127(1):198－203.

［3］Gibson LJ.The mechanical behavior of cancellous bone［J］.Biomech，1985，18(5):317－328.

［4］Teo JC，Si-Hoe KM，Keh JE，et al.Relationship between CT intensity，micro-architecture and mechanical properties of porcine vertebral cancellous bone［J］.Clin Biomech，2006，21(3):235－244.

［5］Guo XE，Kim C H.Mechanical consequence of trabecular bone loss and its treatment:a three-dimensional model simulation［J］.Bone，2002，30(2):404－411.

［6］Diamant I，Shahar R，Masharawi Y，et al.A method for patient-specific evaluation of vertebral cancellous bone strength:In vitro validation［J］.Clinical Biomechanics.2007，22(3):282－291.

［7］Gefen A，Portnoy S，Diamant I.Inhomogeneity of tissue-level strain distributions in individual trabeculae:Mathematical model studiesofnormaland osteoporosis cases ［J］. Medical Engineering＆ Physics，2008，30(5):624－630.

［8］Duchemin L，Bousson V，Raossanaly C，et al.Prediction of mechanical properties of cortical bone by quantitative computed tomography［J］.Medical Engineering ＆ Physics，2008，30(3):321－328.

［9］Li Shuguang.Boundary conditions for unit cells from periodical microstructures and their implications［J］.Compos Sci Tech，2008，68(9):1962－1974.

［10］Zihui Xia，Chuwei Zhou，Qiaoling Yong，et al.On selection of repeated unit cell model and application of unified periodic boundary conditions in micro-mechanical analysis of composites［J］.International Journal of Solids and Structures.2006，43(2):266－278

［11］Goulet RW，Goldstein SA，Ciarelli MJ，et al.The relationship between the structural andorthogonal compressive properties of trabecular bone［J］.Biomech，1994，27(1):375 －389.

［12］Hou FJ，Lang SM，Hoshaw SJ，et al.Human vertebral body apparent and hard tissue stiffness［J］.Biomech，1998，31(1):1009－1015.

［13］Crawford RP，Cann CE，Keaveny TM.Finite element models predict in vitro vertebral body compressive strength better than quantitative computed tomography［J］.Bone，2003，33(1):744－750.

［14］Grant B，Sarah K，Easleya，Tony M，et al.Side-artifact errors in yield strength and elastic modulus for human trabecular bone and their dependence on bone volume fraction and anatomic site［J］.Journal of Biomechanics，2007，40(1):3381－3388.

［15］Chamith S，Rajapaksea J S，Thomsen J S，et al.An expression relating breaking stress and density of trabecular bone［J］.Journal of Biomechanics，2004，37(2):1241－1249.