基于超彈性模型的牙周膜生物力學(xué)響應(yīng)

黃輝祥 湯文成 吳 斌 嚴(yán) 斌

(1東南大學(xué)機(jī)械工程學(xué)院, 南京 211189)

(2南京林業(yè)大學(xué)機(jī)電工程學(xué)院, 南京 210037)

(3南京醫(yī)科大學(xué)口腔醫(yī)學(xué)院, 南京 210029)

牙周膜是連接牙齒與牙槽骨之間的纖維結(jié)蒂組織,其具有支持牙齒,傳遞、吸收和分散咬合力的作用,在口腔正畸及其生物力學(xué)研究中扮演著重要的角色．正畸矯治中,正畸力通過(guò)牙齒作用在牙周組織上,進(jìn)而直接或間接地作用在頜骨上,從而引起一系列的生物反應(yīng),致使牙槽骨發(fā)生改建,最終使牙齒的位置發(fā)生變化,產(chǎn)生移動(dòng)．牙周膜的應(yīng)力、應(yīng)變是正畸牙齒移動(dòng)的始動(dòng)因素,參與牙槽骨的重建過(guò)程[1-2]．由于牙周膜的復(fù)雜性及實(shí)驗(yàn)力學(xué)的局限性,當(dāng)前對(duì)牙周膜的生物力學(xué)研究多采用有限元方法．要獲得準(zhǔn)確的有限元分析結(jié)果,合適的牙周膜材料本構(gòu)模型是關(guān)鍵,牙周膜的生物力學(xué)響應(yīng)可以通過(guò)本構(gòu)模型來(lái)進(jìn)行預(yù)測(cè),從而為牙齒動(dòng)態(tài)仿真及正畸力量化奠定理論基礎(chǔ)．以前的有限元分析中牙周膜多采用線彈性本構(gòu)模型[3-4],但是牙周膜表現(xiàn)出明顯的非線性,其生物力學(xué)特性也更適于用非線性本構(gòu)模型來(lái)描述[5],其本構(gòu)關(guān)系曲線具有指數(shù)形式,當(dāng)前研究局限于采用經(jīng)典常規(guī)的超彈性模型來(lái)研究牙周膜的力學(xué)性能,低階模型描述能力差,高階模型參數(shù)多且參數(shù)因牙周膜實(shí)驗(yàn)數(shù)據(jù)缺少而不易確定,因此,用盡可能少的參數(shù)描述牙周膜材料的實(shí)際力學(xué)性能且無(wú)顯著偏差的超彈性模型是必要的．但利用指數(shù)形式的超彈性模型開(kāi)展牙周膜生物力學(xué)響應(yīng)的研究鮮有報(bào)道．本文基于指數(shù)形式的超彈性模型,借助于有限元模型子程序?qū)φο卵乐苣さ纳锪W(xué)響應(yīng)進(jìn)行了分析研究．

1 超彈性本構(gòu)模型

超彈性模型的本構(gòu)關(guān)系可以用應(yīng)變能的形式表示[6]:

(1)

(2)

(3)

(4)

第二Piola-Kirchhoff應(yīng)力張量為

(5)

則不可壓縮材料的Cauchy應(yīng)力為

(6)

式中，p為靜水壓力.

超彈性本構(gòu)模型較多,目前得到廣泛應(yīng)用的主要有Mooney-Rivlin模型、Neo-Hookean模型、Yeoh模型和Ogden模型等[2]．對(duì)于生物軟組織來(lái)說(shuō),其應(yīng)力-應(yīng)變關(guān)系曲線表現(xiàn)為指數(shù)形式[7],即初始應(yīng)力隨應(yīng)變平緩變化,當(dāng)應(yīng)變達(dá)到一定程度時(shí),應(yīng)力隨應(yīng)變迅速增加,而Mooney-Rivlin模型、Neo-Hookean模型及低階Ogden模型等的指數(shù)性較差．為了能夠更好地描述生物軟組織的力學(xué)特性,一些學(xué)者在原有經(jīng)典模型的基礎(chǔ)上提出了指數(shù)形式的修正模型,主要有Mooney-Rivlin指數(shù)模型[8]、Ogden指數(shù)模型[7]等．

Mooney-Rivlin指數(shù)模型為

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

(7)

式中,c1,c2,c3為材料參數(shù)．

由式(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指數(shù)模型為

(10)

式中,α1,α2為材料參數(shù)．

(11)

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

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

(12)

在實(shí)際工程應(yīng)用中,通?？芍苯訙y(cè)量獲得工程應(yīng)力τ(即第一Piola-Kirchhoff應(yīng)力)和工程應(yīng)變?chǔ)?工程應(yīng)力τ與Cauchy應(yīng)力σ(真應(yīng)力)的關(guān)系為

σ=τFT

(13)

利用伸長(zhǎng)率與工程應(yīng)變關(guān)系λi=1+εi(i=1,2,3),式(9)與式(12)可改寫(xiě)為拉伸方向上的工程應(yīng)力τ與工程應(yīng)變?chǔ)诺年P(guān)系．

2 本構(gòu)模型參數(shù)確定

牙周膜本構(gòu)模型研究在國(guó)內(nèi)較少,缺少相關(guān)的實(shí)驗(yàn)數(shù)據(jù),而在國(guó)外是研究熱點(diǎn),但由于人類(lèi)牙周膜結(jié)構(gòu)的復(fù)雜性和特殊性,較多的是采用動(dòng)物實(shí)驗(yàn).文獻(xiàn)[9]數(shù)據(jù)來(lái)源于對(duì)2名女性患者的體內(nèi)拉伸實(shí)驗(yàn),雖然獲得的數(shù)據(jù)較少,但可信度較高[5]．因此本文采用文獻(xiàn)[9]中牙周膜的單軸拉伸實(shí)驗(yàn)數(shù)據(jù)分別對(duì)Ogden指數(shù)模型和Mooney-Rivlin指數(shù)模型以最小二乘法進(jìn)行擬合,擬合曲線如圖1所示,擬合計(jì)算獲得的模型參數(shù)見(jiàn)表1．?dāng)?shù)據(jù)擬合的相對(duì)誤差為

(14)

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

表1 超彈性指數(shù)模型參數(shù)

從圖1可見(jiàn),Ogden指數(shù)模型和Mooney-Rivlin指數(shù)模型與單軸拉伸實(shí)驗(yàn)數(shù)據(jù)的擬合效果較好,其數(shù)據(jù)擬合的相對(duì)誤差分別為0.1639和0.1589．雖然兩模型數(shù)據(jù)擬合的相對(duì)誤差較大,但其擬合曲線在整個(gè)應(yīng)變區(qū)域是穩(wěn)定的,擬合誤差較大可能是文獻(xiàn)[9]中的牙周膜單軸拉伸實(shí)驗(yàn)數(shù)據(jù)的分散性大所致．從計(jì)算的擬合相對(duì)誤差來(lái)說(shuō),Mooney-Rivlin指數(shù)模型更適合描述牙周膜的瞬時(shí)超彈性．

3 有限元模擬與分析

3.1 模型的有限元實(shí)現(xiàn)及驗(yàn)證

本文采用Mooney-Rivlin指數(shù)模型對(duì)牙周膜生物力學(xué)響應(yīng)進(jìn)行分析,由于通用有限元軟件ABAQUS材料庫(kù)中缺少M(fèi)ooney-Rivlin指數(shù)材料本構(gòu)模型,需要對(duì)其用戶(hù)子程序UHYPER或UMAT進(jìn)行二次開(kāi)發(fā)．

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

圖2 牙周膜拉伸模擬應(yīng)力云圖

圖3 應(yīng)力-應(yīng)變關(guān)系曲線對(duì)比

3.2 有限元建模及分析

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

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

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

圖5 超彈性和線彈性模型牙周膜應(yīng)力對(duì)比

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

圖6 超彈性模型牙周膜應(yīng)變

4 結(jié)語(yǔ)

由于牙周膜拉應(yīng)力-應(yīng)變關(guān)系呈現(xiàn)非線性的指數(shù)形式,本文將Mooney-Rivlin指數(shù)模型引入到牙周膜生物力學(xué)響應(yīng)的模擬分析中,利用人體牙周膜拉伸實(shí)驗(yàn)數(shù)據(jù),擬合出模型的材料參數(shù).因?qū)嶒?yàn)數(shù)據(jù)較少且發(fā)散性大,導(dǎo)致數(shù)據(jù)擬合的相對(duì)誤差較大,會(huì)影響模型對(duì)牙周膜力學(xué)特性的描述能力,因此,較為精確可靠的力學(xué)實(shí)驗(yàn)數(shù)據(jù)是牙周膜力學(xué)模型研究亟待解決的問(wèn)題．在超彈性指數(shù)模型本構(gòu)關(guān)系的基礎(chǔ)上,開(kāi)發(fā)了模型的ABAQUS用戶(hù)子程序UHYPER,通過(guò)拉伸實(shí)驗(yàn)測(cè)試數(shù)據(jù)的模擬計(jì)算驗(yàn)證了子程序開(kāi)發(fā)的正確性,同時(shí)該子程序具有比較理想的計(jì)算精度,從而為自定義的本構(gòu)模型開(kāi)發(fā)提供了借鑒．Mooney-Rivlin指數(shù)模型能夠較好地描述牙周膜的瞬時(shí)彈性響應(yīng)．

牙周膜不僅表現(xiàn)出瞬時(shí)的超彈性,還表現(xiàn)為與時(shí)間歷程相關(guān)的黏彈性,后續(xù)工作重點(diǎn)是如何將牙周膜的超彈性與黏彈性結(jié)合起來(lái),建立一個(gè)能夠全面描述牙周膜生物力學(xué)特性的黏-超彈性本構(gòu)模型,同時(shí)開(kāi)展人體體內(nèi)實(shí)驗(yàn)測(cè)試研究來(lái)進(jìn)一步驗(yàn)證和評(píng)價(jià)所構(gòu)建的本構(gòu)模型．

)

[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]魏志剛,湯文成,嚴(yán)斌,等.基于黏彈性模型的牙周膜生物力學(xué)研究[J].東南大學(xué)學(xué)報(bào):自然科學(xué)版,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.