盧子興 裴 鶴

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

引言

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

1 松質(zhì)骨的有限元模型

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

式中，骨小梁長(zhǎng)度 L沿 z方向，r(z)是 z處的半徑，Tb.Th是最大厚度(一般是兩終端處厚度)與最小厚度(即通常是小梁中間部位厚度)的平均值，而經(jīng)驗(yàn)常數(shù)取 α=1.373 6，β=40.9 μm，由多次實(shí)驗(yàn)的均值確定［6］。

圖1 模擬松質(zhì)骨微結(jié)構(gòu)的14面體單胞實(shí)體模型。(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(橫向骨小梁長(zhǎng)度)和l(縱向骨小梁長(zhǎng)度)。模型的各向異性比定義為14面體單胞的高度H和寬度D的比值，即λ=H/D，可以證明它僅與l/b有關(guān)。在有限元計(jì)算時(shí)，骨小梁組織材料的性能參數(shù)采用文獻(xiàn)［8］中得出的皮質(zhì)骨的平均壓縮力學(xué)性能，即彈性模量E0為12 GPa，泊松比為0.3，極限壓縮強(qiáng)度為107 MPa。文中模型全部采用四面體實(shí)體單元SOLID45劃分。為使有限元模型在邊界面節(jié)點(diǎn)位置上符合周期性邊界條件［9］的要求，劃分單元時(shí)首先劃分邊界面處單元，相對(duì)的邊界面單元通過(guò)復(fù)制保持相同，使得節(jié)點(diǎn)能順利耦合，劃分完邊界面后再劃分體內(nèi)單元，通過(guò)控制單元尺寸使計(jì)算達(dá)到一定的精度并具有較高的效率。

2 邊界條件及計(jì)算方法

2.1 邊界條件

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

2.2 計(jì)算方法

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

另外，由文獻(xiàn)［11－12］可知，松質(zhì)骨的表觀模量和表觀強(qiáng)度之間有極強(qiáng)的相關(guān)性。當(dāng)預(yù)測(cè)了松質(zhì)骨的表觀模量后，可根據(jù) Hou等提出的線性關(guān)系［12］預(yù)測(cè)其表觀極限壓縮強(qiáng)度，即

圖2 有限元計(jì)算不同載荷步時(shí)模型的應(yīng)力云圖及網(wǎng)格圖。(a)Mises應(yīng)力分布;(b)“殺死”單元;(c)“生”單元;(d)繼續(xù)加載后的Mises應(yīng)力分布;(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等則由實(shí)驗(yàn)確定了松質(zhì)骨極限壓縮強(qiáng)度與屈服強(qiáng)度之間的關(guān)系［13］，即

因此，經(jīng)過(guò)有限元計(jì)算，可以得到松質(zhì)骨的極限壓縮強(qiáng)度σult;通過(guò)式(3)，便可得到松質(zhì)骨的屈服強(qiáng)度。

3 數(shù)值結(jié)果及討論

3.1 骨體積分?jǐn)?shù)與名義骨小梁厚度的關(guān)系

圖3給出了BV/TV(骨體積分?jǐn)?shù))與Tb.Th(名義骨小梁厚度)之間的關(guān)系。從中可以看出:BV/TV隨著Tb.Th的增加而明顯增大，故通過(guò)控制Tb.Th可以控制BV/TV的變化;l/b對(duì)BV/TV有較大影響，在Tb.Th相同的情況下，隨 l/b增大，BV/TV逐漸減小;此外，通過(guò)改變 l/b，還可以控制 λ的大小。

3.2 壓縮模量預(yù)測(cè)及同實(shí)驗(yàn)結(jié)果的對(duì)比

圖3 BV/TV與Tb.Th的關(guān)系(b=1.2 mm)Fig.3 The relations between BV/TV and Tb.Th

圖4給出了相對(duì)壓縮模量Eyy/E0與BV/TV的關(guān)系及同實(shí)驗(yàn)結(jié)果的對(duì)比。從中可以看出:一是隨著BV/TV的增加，壓縮模量隨之單調(diào)增加;與前人所得腰椎松質(zhì)骨壓縮實(shí)驗(yàn)［14］(取相近Tb.Th及b的數(shù)據(jù))相比，實(shí)驗(yàn)值略高于計(jì)算值，說(shuō)明了此模型的有效性;二是壓縮模量與BV/TV間表現(xiàn)較好的平方律關(guān)系，與 Carter和 Hayes給出的規(guī)律［1］相符。因此，該模型能為分析個(gè)體骨質(zhì)疏松患者的骨骼質(zhì)量提供一定的理論依據(jù)。需指出，有限元模型計(jì)算得到的模量值基本上不隨骨小梁長(zhǎng)度變化，而真實(shí)松質(zhì)骨的力學(xué)性能與骨小梁長(zhǎng)度并不是無(wú)關(guān)的。為此，Grant等在大量實(shí)驗(yàn)數(shù)據(jù)的基礎(chǔ)上，提出了針對(duì)椎體松質(zhì)骨的模量修正因子αE及屈服強(qiáng)度修正因子 ασ的概念［14］，即

式中Et和Em以及σyt和σym分別代表模量和屈服強(qiáng)度的真值和測(cè)量值。研究發(fā)現(xiàn)，αE對(duì)Tb.Th和BV/TV并不敏感，而對(duì)Tb.Sp(文中為l+b/2)存在一定的線性相關(guān)性，且 αE與 ασ之間也存在線性關(guān)系?；谏鲜鲂拚蜃拥母拍詈皖愃频年P(guān)系式，筆者對(duì)有限元計(jì)算得到的壓縮模量進(jìn)行了修正，給出的預(yù)測(cè)結(jié)果如圖4所示。

圖4 Eyy/E0隨BV/TV的變化Fig.4 The variation of Eyy/E0with BV/TV

3.3 壓縮應(yīng)力-應(yīng)變曲線

圖5 壓縮應(yīng)力-應(yīng)變曲線Fig.5 The curves of compressive stress and strain

圖5給出了 b=0.9 mm、l/b=1.5、BV/TV=10.3%時(shí)模型在單向受壓情況下的應(yīng)力-應(yīng)變曲線，同時(shí)給出了一個(gè)由壓縮實(shí)驗(yàn)得到的椎體松質(zhì)骨樣本的壓縮應(yīng)力-應(yīng)變曲線。可見，兩者在彈性階段基本一致，但實(shí)驗(yàn)測(cè)得的強(qiáng)度要大于有限元計(jì)算的結(jié)果(計(jì)算極限壓縮強(qiáng)度為1.9 MPa)。這是因?yàn)檎鎸?shí)松質(zhì)骨的骨組織中含有豐富的結(jié)締組織，而不僅僅是無(wú)機(jī)物，從而增強(qiáng)了松質(zhì)骨的韌性，使得屈服強(qiáng)度與極限強(qiáng)度較之純骨組織有所提高。顯然，對(duì)于特定的 l/b(如取1.0，1.25，1.50和1.75)，通過(guò)改變 BV/TV，可以得到一系列的壓縮應(yīng)力-應(yīng)變曲線(限于篇幅，這里僅給出 l/b=1.50對(duì)應(yīng)的結(jié)果)，如圖6所示。從中可以看出，隨著BV/TV的增加，極限壓縮強(qiáng)度逐漸增大，但所對(duì)應(yīng)的應(yīng)變卻逐漸減小。

圖6 不同 BV/TV下的壓縮應(yīng)力-應(yīng)變曲線(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 極限壓縮強(qiáng)度預(yù)測(cè)與實(shí)驗(yàn)值的比較

如上所述，針對(duì) l/b=1.0，1.25，1.50和 1.75幾種情況，再通過(guò)改變 BV/TV，可以得到一系列的壓縮應(yīng)力-應(yīng)變曲線。從中提取對(duì)應(yīng)每種情況的極限壓縮強(qiáng)度σult，并用強(qiáng)度修正因子對(duì)其修正，獲得的 σult與 BV/TV 的關(guān)系及同實(shí)驗(yàn)值［12，15］的比較，如圖7所示。從中可以看出:一是當(dāng) BV/TV較低時(shí)(BV/TV<0.15)，l/b值對(duì) σult存在明顯影響，相同BV/TV時(shí)，σult隨著l/b的增加而增大，這表明，較低BV/TV時(shí)，松質(zhì)骨呈現(xiàn)棒狀骨小梁結(jié)構(gòu)，縱向骨小梁是主要承載者。二是當(dāng)BV/TV繼續(xù)增大時(shí)，對(duì)應(yīng)不同l/b值的曲線逐漸匯聚;當(dāng)BV/TV增大到一定程度后，將出現(xiàn)板狀骨小梁(這里不予討論)。三是l/b值接近1且BV/TV較低時(shí)，預(yù)測(cè)的σult能較好地與實(shí)驗(yàn)結(jié)果相符合;當(dāng)BV/TV增大時(shí)，松質(zhì)骨內(nèi)出現(xiàn)板狀骨小梁，使松質(zhì)骨抗壓能力得到增強(qiáng)，因此實(shí)驗(yàn)測(cè)得的 σult比使用棒狀骨小梁模型預(yù)測(cè)的 σult偏大，這符合真實(shí)情況。

圖7 σult與BV/TV的關(guān)系及同實(shí)驗(yàn)值的比較Fig.7 The relations between σultand BV/TV and the comparison with the experimental values

3.5 極限壓縮強(qiáng)度預(yù)測(cè)與經(jīng)驗(yàn)公式預(yù)測(cè)的比較

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

圖8 σult有限元計(jì)算結(jié)果與經(jīng)驗(yàn)公式對(duì)比Fig.8 The comparison between the FE results of σult and the empirical formula

4 結(jié)論

1)14面體單胞模型能有效地用于松質(zhì)骨壓縮模量和強(qiáng)度的數(shù)值預(yù)測(cè)，預(yù)測(cè)結(jié)果與實(shí)驗(yàn)數(shù)據(jù)及經(jīng)驗(yàn)公式相吻合;

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

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

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