范若尋,劉杰,甘樹坤,張敏
(1.吉林化工學院汽車工程學院,吉林 132102;2. 吉林大學機械科學與工程學院,長春 130025)
隨著醫療技術的發展,骨性能測試變得簡單可行,其中骨表觀力學參數常被用來評判人體骨質量及預測骨折風險[1-2]。但由于骨結構分層水平特性,其表觀力學性能往往由不同水平骨材料性能所決定[3-5]。因此,骨微觀力學性能對表觀性能,甚至是骨質量與骨折風險均有影響,尤其是骨在微觀水平的失效參數,對骨在不同水平力學性能有直接影響:在表觀水平,骨微觀失效應變在一定程度上決定了骨結構強度極限;在組織水平,骨微觀失效應變則影響著松質骨小梁的損傷與重建過程[2,6-9]。因此,為深入了解骨疾病的發病機理,以及預測人體骨質量與骨折風險,除了掌握表觀力學性能,更需對骨的微觀力學性能參數進行預測與研究。
有限元分析被認為是解決生物醫學工程,特別是骨力學問題的一種有效方法[2-3,9]。然而,有限元模型材料參數的選取卻是一個影響分析精度的重要問題[6,10-11]。一些仿真參數,如彈性模量,隨著納米壓痕等實驗技術的發展,已能準確測得[2,8]。而當對骨結構進行斷裂模擬時,需輸入骨在微觀水平的材料參數,例如骨組織屈服與斷裂應變等,由于這類參數難以通過實驗測得,所以針對骨結構的斷裂仿真模擬,一般是依據文獻或經驗公式賦予微觀水平材料參數,這在一定程度上降低了有限元模型的生物逼真度與分析準確性[3-4,9]。
既然骨在微觀水平的力學性能參數對臨床醫學及有限元分析均有重要意義,如何測得該水平上的力學參數即成為一個重要問題。因此,本研究的目的為開發一種骨微觀力學性能參數預測方法,以大鼠股骨皮質骨為研究對象,以表觀壓縮實驗數據為研究依據,通過結合實驗與有限元仿真技術,預測骨的一種微觀力學性能參數——微觀失效應變。
選取四只7月齡雄性大鼠,獲取其右側股骨,在股骨中部截取5 mm長皮質骨,作為實驗樣本,置于試驗機中進行壓縮實驗,壓頭壓縮速度以位移量控制。將皮質骨上表面與壓頭下表面的接觸面積作為計算表觀強度極限的橫截面積,四個樣本的橫截面積范圍處在26~29 mm2,依據樣本的力-位移曲線計算其表觀應力-應變曲線。
基于大鼠皮質骨微觀影像,利用逆向建模軟件重構皮質骨外表面,通過有限元軟件應用三維四面體單元建立皮質骨有限元模型,最后在皮質骨模型上、下方分別加入剛性壓板,其中下壓板作為固定端,約束所有自由度,同時將0.5 mm的壓縮位移施加在上壓板,以模擬實驗條件,見圖1。由于網格的疏密程度對有限元仿真結果有較大影響,所以針對同一樣本分別建立了網格尺寸為5、10、15、20、25、30、35、40 μm的有限元模型,進行網格敏感性分析,以確定網格尺寸。大鼠股骨皮質骨的彈性模量由課題組前期所做納米壓痕實驗測得,7月齡大鼠股骨皮質骨的彈性模量平均值為34 000 MPa,泊松比為0.3[2]。

圖1 大鼠股骨皮質骨有限元模型示意圖
Fig1Schematicdrawingoftheratfemoralcorticalbonefiniteelementmodel
大鼠股骨皮質骨有限元模型的斷裂由內部等效應變控制[6-7]。這里采用單元刪除技術:當單元達到失效條件時,即等效應變超過皮質骨材料的微觀失效應變時,該單元即視為被刪除,將被刪除單元的彈性模量設定為1 MPa[2,4]。隨著壓縮載荷持續增加,越來越多的單元被刪除,當被刪除單元達到一定數量時,有限元模型即會發生表觀斷裂失效。
為模擬上述過程,本研究在有限元軟件中編寫了子程序USDFLD,將上述分析思路與單元失效判定條件寫入子程序。由于其它仿真輸入參數均為已知,僅有微觀失效應變這一未知輸入參數。因此,只需在子程序中嘗試賦予微觀失效應變數值即可模擬出皮質骨在壓縮載荷作用下的斷裂。將仿真得到的表觀應力-應變曲線與樣本壓縮實驗測得的應力-應變曲線相對比,觀測二者擬合情況,如擬合不成功,則需調節預設微觀失效應變,直至曲線擬合成功(擬合成功的標準為仿真與實驗測出的表觀斷裂應變差值小于5%),即認定此時預設的微觀失效應變為擬獲取的大鼠股骨皮質骨微觀失效應變,具體仿真流程見圖2。
圖3顯示了同一樣本采用不同網格所建立的有限元模型在壓縮載荷作用下得到的表觀斷裂應變??梢钥闯觯畲笈c最小網格尺寸模型所得到的數值差異達到了25.977%,同時表觀斷裂應變隨著網格尺寸增大顯著減小,但當網格尺寸小于20 μm 時,不同有限元模型得到的表觀斷裂應變數值差異卻都小于5%。因此,考慮到仿真準確性與計算成本,大鼠股骨皮質骨有限元模型的網格尺寸均選取為20 μm。

圖2大鼠股骨皮質骨微觀失效應變預測流程
Fig2Simulatedprocessonthemicro-levelfailurestrainofratfemoralcorticalbone

圖3 單元尺寸對表觀斷裂應變的敏感性分析
Fig3Meshsensitivityanalysisfortheeffectofelementsizesontheapparentfailurestrain
通過與表觀壓縮實驗數據進行反演分析,反復調節預設微觀失效應變,四個皮質骨模型均成功擬合出實驗測得的表觀應力-應變曲線,其中預測出的微觀失效應變以及模型發生表觀斷裂時失效單元數量所占百分比見表1??梢钥闯?,同一月齡的大鼠股骨皮質骨的微觀失效應變數值差異較小,而且當失效單元數量占所有單元比例達到6%~7%時,皮質骨有限元模型即可能發生表觀斷裂失效。
表1大鼠股骨皮質骨的微觀失效應變數值與失效單元百分比
Table1Thepredictedmicro-levelfailurestrainsforratfemoralcorticalbonesandthefailedelementpercentage

樣本微觀失效應變(%)失效單元百分比(%)樣本14.626.35樣本24.536.21樣本34.656.43樣本44.757.22
圖4顯示了四個大鼠股骨皮質骨樣本通過有限元分析與表觀壓縮實驗獲得的表觀應力-應變曲線對比。可以看出仿真所得的曲線形狀與實驗較為吻合,同時仿真預測出的表觀強度極限和斷裂應變數值也與實驗結果相差無幾,差異均保持在5%以內。
本研究通過采用一種將表觀壓縮實驗與有限元仿真相結合的方法,預測出大鼠股骨皮質骨的一種微觀力學性能參數。長期以來,骨裂與骨折等骨科疾病一直困擾著人們,特別是老年人,一旦發生骨折將導致一系列嚴重后果[5,11]。然而研究表明,骨折的發生除了外力作用,并非僅與骨密度的降低或骨結構的退變等表觀因素相關,在很大程度上也取決于微觀力學性能的退變[3,12]。因此,研究骨微觀力學性能參數的預測方法,對骨疾病的診斷與預防都將起到重要作用。

圖4四個皮質骨樣本有限元分析與壓縮實驗獲得的表觀應力-應變曲線對比
Fig4Comparisonoftheapparentstress-straincurvesbetweenfiniteelementanalysisandcompressiveexperimentforfourcorticalbonesamples
骨微觀力學性能參數的測試一般需借助不同水平的實驗手段完成[1]。相比表觀實驗,微觀實驗不但操作復雜,而且對樣本精度也有較高要求,所以,針對微觀力學性能參數的測試長期以來都是一個難點[6,13]?;诖?,本研究通過結合表觀壓縮實驗與有限元仿真技術,開發了一種骨微觀力學性能參數預測方法。該方法不但避免了難度較大的微觀實驗,而且通過與實驗應力-應變曲線對比發現仿真得出的皮質骨表觀斷裂應變數值與實驗數值差異均小于5%,同時通過觀察發現皮質骨樣本與皮質骨有限元模型在壓縮載荷作用下的斷裂模式也較為相似,這些均說明了本研究所開發預測方法的可行性以及獲取皮質骨微觀失效應變參數的準確性。因此,本預測過程為骨微觀力學性能參數的獲取提供一定幫助。同時,該方法不僅可以預測微觀失效應變,也可結合相關表觀實驗預測其它微觀力學性能參數,例如骨材料的彈性模量或松質骨的微觀屈服應變等微觀力學性能參數。
斷裂模擬過程中,為保證仿真分析精度,首先要設置足夠小的載荷增量步,以保證當預設微觀失效應變有微小調整時,仿真得出的應力-應變曲線會有相應變化;同時在與實驗曲線進行反演分析時發現,當對預設微觀失效應變進行調節時,相比表觀斷裂應變,仿真得出的表觀強度極限的變化幅度更大,即當對預設值進行微小調整時,表觀斷裂應變的預測值不會有太大波動,但表觀強度極限的變化卻較大。結合該特點,在模擬過程中,首先確保表觀強度極限的預測值與實驗結果相符,以大致確定預設微觀失效應變的范圍,然后在此范圍內小幅度調整,以找到與實驗所測數值最為吻合的預設失效應變。
本研究在斷裂模擬過程中應用了單元刪除技術,該技術能夠模擬出單元從開始變形至最后失效的力學性能變化過程,但無法表達出單元的屈服過程,因此,該技術目前更適用于模擬材料的脆性斷裂。依據實驗結果顯示,大鼠皮質骨呈現出脆性斷裂特征,所以應用單元刪除技術能夠準確模擬出皮質骨模型的斷裂過程。而人體松質骨會表現出韌性斷裂特征,因此為了預測松質骨的微觀力學參數,日后準備基于ABAQUS平臺,通過開發UMAT子程序,應用連續損傷技術,以模擬松質骨單元的屈服與斷裂過程,從而對松質骨與皮質骨結構進行全面的斷裂模擬。
本研究提出了一種骨微觀力學性能參數預測方法,以表觀壓縮實驗為研究依據,通過模仿實驗條件,應用USDFLD子程序模擬出大鼠股骨皮質骨有限元模型在壓縮載荷作用下的斷裂過程,并通過與實驗所測數據進行對比、反演,預測出大鼠股骨皮質骨的微觀失效應變。結果表明四只7月齡大鼠,其股骨皮質骨的微觀失效應變數值處于4.53%~4.75%之間。經實驗對比驗證顯示,該方法能夠準確預測出骨結構在微觀水平的力學性能參數。