基于超彈性模型的牙周膜生物力學響應

黃輝祥 湯文成 吳 斌 嚴 斌

(1東南大學機械工程學院, 南京 211189)

(2南京林業大學機電工程學院, 南京 210037)

(3南京醫科大學口腔醫學院, 南京 210029)

牙周膜是連接牙齒與牙槽骨之間的纖維結蒂組織,其具有支持牙齒,傳遞、吸收和分散咬合力的作用,在口腔正畸及其生物力學研究中扮演著重要的角色．正畸矯治中,正畸力通過牙齒作用在牙周組織上,進而直接或間接地作用在頜骨上,從而引起一系列的生物反應,致使牙槽骨發生改建,最終使牙齒的位置發生變化,產生移動．牙周膜的應力、應變是正畸牙齒移動的始動因素,參與牙槽骨的重建過程[1-2]．由于牙周膜的復雜性及實驗力學的局限性,當前對牙周膜的生物力學研究多采用有限元方法．要獲得準確的有限元分析結果,合適的牙周膜材料本構模型是關鍵,牙周膜的生物力學響應可以通過本構模型來進行預測,從而為牙齒動態仿真及正畸力量化奠定理論基礎．以前的有限元分析中牙周膜多采用線彈性本構模型[3-4],但是牙周膜表現出明顯的非線性,其生物力學特性也更適于用非線性本構模型來描述[5],其本構關系曲線具有指數形式,當前研究局限于采用經典常規的超彈性模型來研究牙周膜的力學性能,低階模型描述能力差,高階模型參數多且參數因牙周膜實驗數據缺少而不易確定,因此,用盡可能少的參數描述牙周膜材料的實際力學性能且無顯著偏差的超彈性模型是必要的．但利用指數形式的超彈性模型開展牙周膜生物力學響應的研究鮮有報道．本文基于指數形式的超彈性模型,借助于有限元模型子程序對正畸力下牙周膜的生物力學響應進行了分析研究．

1 超彈性本構模型

超彈性模型的本構關系可以用應變能的形式表示[6]:

(1)

(2)

(3)

(4)

第二Piola-Kirchhoff應力張量為

(5)

則不可壓縮材料的Cauchy應力為

(6)

式中，p為靜水壓力.

超彈性本構模型較多,目前得到廣泛應用的主要有Mooney-Rivlin模型、Neo-Hookean模型、Yeoh模型和Ogden模型等[2]．對于生物軟組織來說,其應力-應變關系曲線表現為指數形式[7],即初始應力隨應變平緩變化,當應變達到一定程度時,應力隨應變迅速增加,而Mooney-Rivlin模型、Neo-Hookean模型及低階Ogden模型等的指數性較差．為了能夠更好地描述生物軟組織的力學特性,一些學者在原有經典模型的基礎上提出了指數形式的修正模型,主要有Mooney-Rivlin指數模型[8]、Ogden指數模型[7]等．

Mooney-Rivlin指數模型為

W=c1(exp(c2(I1-3))-1)+c3(I2-3)

(7)

式中,c1,c2,c3為材料參數．

由式(6)、(7)及單軸拉伸條件(σ22=σ33=0),可得

p=2λ-1[c1c2exp(c2(λ2+2λ-1-3))+c3(λ2+λ-1)]

(8)

σ11=2c1c2exp(c2(λ2+2λ-1-3))·
(λ2-λ-1)+2c3(λ-λ-2)

(9)

Ogden指數模型為

(10)

式中,α1,α2為材料參數．

(11)

σ11=2c1c2exp(c2(λα1+2λ-α1/2-3))·

(λα1-λ-α1/2)+2c3α2(λα2-λ-α2/2)

(12)

在實際工程應用中,通常可直接測量獲得工程應力τ(即第一Piola-Kirchhoff應力)和工程應變ε.工程應力τ與Cauchy應力σ(真應力)的關系為

σ=τFT

(13)

利用伸長率與工程應變關系λi=1+εi(i=1,2,3),式(9)與式(12)可改寫為拉伸方向上的工程應力τ與工程應變ε的關系．

2 本構模型參數確定

牙周膜本構模型研究在國內較少,缺少相關的實驗數據,而在國外是研究熱點,但由于人類牙周膜結構的復雜性和特殊性,較多的是采用動物實驗.文獻[9]數據來源于對2名女性患者的體內拉伸實驗,雖然獲得的數據較少,但可信度較高[5]．因此本文采用文獻[9]中牙周膜的單軸拉伸實驗數據分別對Ogden指數模型和Mooney-Rivlin指數模型以最小二乘法進行擬合,擬合曲線如圖1所示,擬合計算獲得的模型參數見表1．數據擬合的相對誤差為

(14)

圖1 Ogden與Mooney-Rivlin指數模型的擬合曲線

表1 超彈性指數模型參數

從圖1可見,Ogden指數模型和Mooney-Rivlin指數模型與單軸拉伸實驗數據的擬合效果較好,其數據擬合的相對誤差分別為0.1639和0.1589．雖然兩模型數據擬合的相對誤差較大,但其擬合曲線在整個應變區域是穩定的,擬合誤差較大可能是文獻[9]中的牙周膜單軸拉伸實驗數據的分散性大所致．從計算的擬合相對誤差來說,Mooney-Rivlin指數模型更適合描述牙周膜的瞬時超彈性．

3 有限元模擬與分析

3.1 模型的有限元實現及驗證

本文采用Mooney-Rivlin指數模型對牙周膜生物力學響應進行分析,由于通用有限元軟件ABAQUS材料庫中缺少Mooney-Rivlin指數材料本構模型,需要對其用戶子程序UHYPER或UMAT進行二次開發．

為驗證Mooney-Rivlin指數模型UHYPER子程序的有效性,對牙周膜拉伸實驗情況進行模擬計算,施加的拉伸載荷為1MPa,利用UHYPER子程序進行計算獲得的應力分布云圖如圖2所示．選取其中一個單元,提取該單元一個積分點處的應力和應變值,得到應力-應變關系曲線(即數值仿真曲線),與理論模型曲線進行對比,如圖3所示．從圖3可看出,應用UHYPER子程序后ABAQUS計算結果與實驗數據擬合結果具有良好的一致性,表明了本文開發的Mooney-Rivlin指數模型UHYPER子程序的正確性及有效性．

圖2 牙周膜拉伸模擬應力云圖

圖3 應力-應變關系曲線對比

3.2 有限元建模及分析

采用CT對正畸患者的上中切牙進行掃描,借助于醫學圖像建模軟件Mimics和逆向工程軟件Catia建立牙齒、牙周膜和牙槽骨的三維幾何模型,如圖4所示,牙周膜厚度為0.2mm,位于牙齒和牙槽骨之間．由于牙齒和牙槽骨的彈性模量是牙周膜的20000～30000倍,因此在有限元模型中將牙齒和牙槽骨作為剛性體以節省計算資源,它們和牙周膜之間的接觸約束采用綁定連接,牙周膜采用C3D8H單元,在牙冠上唇舌方向施加大小為1N的正畸力,利用開發超彈性模型的ABAQUS子程序UHYPER及線彈性模型分別對牙周膜在正畸力作用下的生物力學響應進行模擬分析研究．線彈性模型中設置牙周膜彈性模量為0.667MPa,泊松比為0.49[11]．

圖4 中切牙、牙周膜和牙槽骨CAD模型

圖5分別為超彈性模型和線彈性模型下的牙周膜的應力響應云圖．從圖中可看出兩者的應力分布存在較大的差異,線彈性模型牙周膜的應力分布相對于超彈性模型的較均勻,且最大應力大于超彈性模型的最大應力．然而,由于牙周膜幾何結構的復雜性致使應力分布情況并不均勻,這表明采用非線性超彈性模型來模擬計算牙周膜的生物力學響應比使用線彈性模型更為合理．兩者應力主要集中分布在牙周膜牙頸處和根尖部位,其最大應力出現在舌側牙周膜牙頸處,為壓應力,同一側根尖處產生拉應力,而唇側的應力情況則相反,研究結果與文獻[2,12]一致．此外,在超彈性牙周膜根尖和牙頸之間的部分存在局部的應力集中．2個模型的最大主應力值分別為2.7190MPa和0.1451MPa．

圖5 超彈性和線彈性模型牙周膜應力對比

圖6為正畸力作用下的超彈性牙周膜的應變響應,其分布情況與應力相似．

圖6 超彈性模型牙周膜應變

4 結語

由于牙周膜拉應力-應變關系呈現非線性的指數形式,本文將Mooney-Rivlin指數模型引入到牙周膜生物力學響應的模擬分析中,利用人體牙周膜拉伸實驗數據,擬合出模型的材料參數.因實驗數據較少且發散性大,導致數據擬合的相對誤差較大,會影響模型對牙周膜力學特性的描述能力,因此,較為精確可靠的力學實驗數據是牙周膜力學模型研究亟待解決的問題．在超彈性指數模型本構關系的基礎上,開發了模型的ABAQUS用戶子程序UHYPER,通過拉伸實驗測試數據的模擬計算驗證了子程序開發的正確性,同時該子程序具有比較理想的計算精度,從而為自定義的本構模型開發提供了借鑒．Mooney-Rivlin指數模型能夠較好地描述牙周膜的瞬時彈性響應．

牙周膜不僅表現出瞬時的超彈性,還表現為與時間歷程相關的黏彈性,后續工作重點是如何將牙周膜的超彈性與黏彈性結合起來,建立一個能夠全面描述牙周膜生物力學特性的黏-超彈性本構模型,同時開展人體體內實驗測試研究來進一步驗證和評價所構建的本構模型．

)

[1]Kawarizadeh A, Bourauel C, Jager A. Experimental and numerical determination of initial tooth mobility and material properties of the periodontal ligament in rat molar specimens [J].EurJOrthod,2003,25(6): 569-578.

[2]Wu Bin, Tang Wencheng, Yan Bin. Study on stress distribution in periodontal ligament of impacted tooth based on hyperelastic model [C]//IEEEInformationEngineeringandComputerScience.Wuhan, China,2009: 1550-1553.

[3]Lu Hongfei, Mai Zhihui, Chen Qi, et al. Initial stress distribution of the maxillary anterior teeth, periodontal ligament and alveolar bone by different intruding loadings: a three-dimensional finite element analysis [J].JournalofClinicalRehabilitativeTissueEngineeringResearch,2011,15(48): 8964-8967.

[4]Hohmann A, Kober C, Young P, et al. Influence of different modelling strategies for the periodontal ligament on finite element simulation results [J].AmJOrthodDentofacialOrthop, 2011,139(6):775-783.

[5]魏志剛,湯文成,嚴斌,等.基于黏彈性模型的牙周膜生物力學研究[J].東南大學學報:自然科學版,2009,39(3):484-489.

Wei Zhigang,Tang Wencheng,Yan Bin,et al. Biomechanical analysis of periodontal ligament based on viscoelastic model [J].JournalofSoutheastUniversity:NaturalScienceEdition,2009,39(3): 484-489.(in Chinese)

[6]Natali A N, Carniel E L, Pacan P G, et al. Experimental-numerical analysis of minipig’s multi-rooted teeth [J].JournalofBiomechanics,2007,40(8):1701-1708.

[7]Zhan Gao, Lister K, Desai J P. Constitutive modelling of liver tissue: experiment and theory [J].AnnalsofBiomedicalEngineering,2010,38(2):505-516.

[8]Chui C, Kobayashi E, Chen X,et al. Compression and elongation experiments and non-linear modelling of liver tissue for surgical simulation [J].MedBiolEngComput,2004,42(6):787-798.

[9]Yoshida Noriaki, Koga Yoshiyuki, Peng Chien-Lun, et al. In vivo measurement of the elastic modulus of the human periodontal ligament [J].MedEngPhy,2001,23(8):567-572.

[10]Hibbitt D, Karlsson B, Sorensen P.ABAQUS/standarduser’smanual[M]. Pawtucket, RI,USA: Dassault Systèmes Simulia Corp, 2010.

[11]Tanne K, Yoshida S, Kawata T, et al. An evaluation of the biomechanical response of the tooth and periodontium to orthodontic forces in adolescent and adult subjects [J].BritishJournalofOrthodontics, 1998,25(2):109-115.

[12]Cattaneo P M, Dalstra M, Melsen B. Strain in periodontal ligament and alveolar bone associated with orthodontic tooth movement analyzed by finite element [J].OrthodCraniofacRes,2009,12(2):120-128.