各向異性松質骨壓縮模量和強度的數(shù)值預測

盧子興 裴 鶴

(北京航空航天大學 航空科學與工程學院，北京 100191)

引言

松質骨的力學性質主要體現(xiàn)為表觀強度和剛度，而這些性質又與松質骨的微結構特性密切相關，如骨小梁形狀、生長方向、密度和連通性等;并且對剛度而言，還與組織密度有關［1］。已有的實驗結果表明，骨小梁的細化和連通性的損失均會造成整骨強度和剛度在準靜態(tài)和沖擊載荷下的退化［2－4］。因此，研究微結構和組織密度特性隨骨小梁尺度的變化及與松質骨表觀力學特性之間的關系具有重要的意義。然而，這種關系卻具有高度的非線性，且隨病患個體微結構的變化而改變［1，5］，很難用數(shù)學公式表達出來，故對其行為進行有限元分析成為一種必要的手段。例如，Guo等人應用由梁單元構成的14面體模型對松質骨微結構進行了模擬［5］，計算結果比較符合松質骨的宏觀力學性能，但使用梁單元并不能很好地反映骨小梁的形態(tài)學特征;而Diamant等利用變截面骨小梁模型組成的簡單立方結構對松質骨進行研究，得出受力情況下松質骨微結構細微處的應力狀態(tài)，但簡單立方結構顯然不能很好地模擬松質骨的微結構［6－7］。本研究發(fā)展了一種有限元參數(shù)化建模方法，用14面體模型來模擬棒狀骨小梁構成的松質骨微結構，對各向異性松質骨的彈性性能和強度進行了有限元計算分析，給出了壓縮模量與實驗結果的比較，并對松質骨的表觀壓縮強度進行了數(shù)值預測。

1 松質骨的有限元模型

為評估松質骨力學行為，參考前人的相關工作，建立了對應不同形態(tài)學參數(shù)的松質骨有限元單胞模型，如圖1所示。其中，骨小梁使用Gefen等人基于大量實驗數(shù)據(jù)提出的單根骨小梁模型［7］，該模型比簡單圓棒形骨小梁模型能更好地表達骨小梁的形態(tài)學特性，即Tb.Th(名義骨小梁厚度)和Tb.Sp(用骨小梁長度L表征)。下式給出建立這種骨小梁模型的形狀控制方程，即

式中，骨小梁長度 L沿 z方向，r(z)是 z處的半徑，Tb.Th是最大厚度(一般是兩終端處厚度)與最小厚度(即通常是小梁中間部位厚度)的平均值，而經驗常數(shù)取 α=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

有限元模型的形狀控制參數(shù)為:Tb.Th、b(橫向骨小梁長度)和l(縱向骨小梁長度)。模型的各向異性比定義為14面體單胞的高度H和寬度D的比值，即λ=H/D，可以證明它僅與l/b有關。在有限元計算時，骨小梁組織材料的性能參數(shù)采用文獻［8］中得出的皮質骨的平均壓縮力學性能，即彈性模量E0為12 GPa，泊松比為0.3，極限壓縮強度為107 MPa。文中模型全部采用四面體實體單元SOLID45劃分。為使有限元模型在邊界面節(jié)點位置上符合周期性邊界條件［9］的要求，劃分單元時首先劃分邊界面處單元，相對的邊界面單元通過復制保持相同，使得節(jié)點能順利耦合，劃分完邊界面后再劃分體內單元，通過控制單元尺寸使計算達到一定的精度并具有較高的效率。

2 邊界條件及計算方法

2.1 邊界條件

針對14面體模型能周期性排滿整個空間的特點，在計算時施加周期性邊界條件［9］，這樣能夠保證模型在變形后仍然保持周期性的空間排列，并允許邊界上的面發(fā)生翹曲變形。按照施加周期性邊界條件的要求，模型相對表面上對應節(jié)點在該平面的法線上以同樣的位移膨脹或收縮，在其他方向上則保持一致的平動位移，且在所有方向上均具有相同的轉角。

2.2 計算方法

當14面體按體心立方排列時，能周期性地排滿整個空間，故采用單倍胞體計算能提高效率。筆者對文獻［10］提出的計算周期性排列單胞力學性能的邊界條件進行適當變化，得到了計算壓縮模量的方法(限于篇幅，這里不再贅述)。下面簡要介紹壓縮極限強度的計算方法。圖2給出了對應不同載荷步時模型的Von Mises應力云圖及其網(wǎng)格圖。從中可以看出，當模型承受y向壓縮時(椎體松質骨通常主要承受縱向壓縮載荷，其他方向承受載荷較小，故主要討論模型沿y向的承載能力)，縱向骨小梁連接處發(fā)生了應力集中。在進行應力判斷后，這個位置網(wǎng)格的Mises應力超過了材料的極限壓縮應力，從而將此處的單元“殺死”，得到如圖2中所示的“生”單元模型。在此基礎上增加位移載荷，進行下一步計算，又一次可得到模型的Mises應力云圖，再繼續(xù)進行下一次的判斷，如此往復，最終可以得到模型的壓縮應力-應變曲線和對應的極限壓縮強度 σult。

另外，由文獻［11－12］可知，松質骨的表觀模量和表觀強度之間有極強的相關性。當預測了松質骨的表觀模量后，可根據(jù) Hou等提出的線性關系［12］預測其表觀極限壓縮強度，即

圖2 有限元計算不同載荷步時模型的應力云圖及網(wǎng)格圖。(a)Mises應力分布;(b)“殺死”單元;(c)“生”單元;(d)繼續(xù)加載后的Mises應力分布;(e)繼續(xù)“殺死”單元;(f)繼續(xù)“生”單元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 數(shù)值結果及討論

3.1 骨體積分數(shù)與名義骨小梁厚度的關系

圖3給出了BV/TV(骨體積分數(shù))與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的數(shù)據(jù))相比，實驗值略高于計算值，說明了此模型的有效性;二是壓縮模量與BV/TV間表現(xiàn)較好的平方律關系，與 Carter和 Hayes給出的規(guī)律［1］相符。因此，該模型能為分析個體骨質疏松患者的骨骼質量提供一定的理論依據(jù)。需指出，有限元模型計算得到的模量值基本上不隨骨小梁長度變化，而真實松質骨的力學性能與骨小梁長度并不是無關的。為此，Grant等在大量實驗數(shù)據(jù)的基礎上，提出了針對椎體松質骨的模量修正因子αE及屈服強度修正因子 ασ的概念［14］，即

式中Et和Em以及σyt和σym分別代表模量和屈服強度的真值和測量值。研究發(fā)現(xiàn)，α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時，松質骨呈現(xiàn)棒狀骨小梁結構，縱向骨小梁是主要承載者。二是當BV/TV繼續(xù)增大時，對應不同l/b值的曲線逐漸匯聚;當BV/TV增大到一定程度后，將出現(xiàn)板狀骨小梁(這里不予討論)。三是l/b值接近1且BV/TV較低時，預測的σult能較好地與實驗結果相符合;當BV/TV增大時，松質骨內出現(xiàn)板狀骨小梁，使松質骨抗壓能力得到增強，因此實驗測得的 σult比使用棒狀骨小梁模型預測的 σult偏大，這符合真實情況。

圖7 σult與BV/TV的關系及同實驗值的比較Fig.7 The relations between σultand BV/TV and the comparison with the experimental values

3.5 極限壓縮強度預測與經驗公式預測的比較

Hou等在大量實驗數(shù)據(jù)的基礎上，提出了預測σult的經驗公式［12］。利用不同 l/b情況下縱向表觀壓縮模量Eyy的預測結果，代入式(2)，得到極限應力值(記為經驗公式預測)與上述有限元計算的結果進行對比得到圖8。從中可以看出:一是有限元計算和經驗公式均能對松質骨的σult做出較好的預測;二是經驗公式預測不能反映各向異性參數(shù) l/b對強度預測的影響，而有限元計算則可反映各向異性參數(shù)l/b對強度預測的影響，l/b的影響從一定程度上說明了σult實驗數(shù)據(jù)分散的合理性。

圖8 σult有限元計算結果與經驗公式對比Fig.8 The comparison between the FE results of σult and the empirical formula

4 結論

1)14面體單胞模型能有效地用于松質骨壓縮模量和強度的數(shù)值預測，預測結果與實驗數(shù)據(jù)及經驗公式相吻合;

2)相對于利用模量對強度進行預測的經驗公式而言，本研究的有限元方法能夠反映各向異性對松質骨強度的影響，從而在一定程度上解釋了實驗中強度數(shù)據(jù)的分散性;

3)由模型受力云圖可以判斷破壞發(fā)生的位置，從而為今后進行斷裂方面的分析奠定基礎。

［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.