牙齒正畸過程中，矯牙托槽的介入和滑動容易導致唇頰軟組織出現較大反應，矯治初期常見軟組織損傷和潰瘍。口腔正畸醫療領域采用臨床案例統計方法進行定性分析，缺乏生物力學機制的定量解釋。為此，開展唇頰—托槽—牙齒的三維模型有限元分析，計算托槽引起的唇頰軟組織的力學反應，其中涉及復雜耦合的接觸非線性、材料非線性和幾何非線性。首先，根據唇頰生物組成特點，優化選取二階奧格登（Ogden）超彈性本構模型，對類脂肪材料的唇頰軟組織進行表征。其次，根據口腔活動特點，建立托槽介入和正交滑動的兩階段仿真模型，并對關鍵接觸參數進行優化設置。最終，采用整體模型—子模型的兩層次分析方法，基于整體模型計算得到的位移邊界，實現子模型高精度應變的高效求解。針對正畸過程中四種典型牙齒形態的計算結果表明：① 軟組織最大應變沿托槽尖銳邊緣分布，與臨床觀測的軟組織變形輪廓一致；② 隨著牙齒排齊，軟組織的最大應變也隨之減小，與矯治初期常見損傷和潰瘍以及矯治后期患者不適感減輕的臨床表現相符。本文的方法可為國內外口腔正畸醫療領域的相關量化分析研究提供參考，進一步有益于新型矯治裝置的產品研發分析。
引用本文： 華嘉皓, 吉利, 代清源, 梁振宇, 郭龍妹, 陳太聰. 矯牙過程托槽致唇頰軟組織非線性反應量化分析研究. 生物醫學工程學雜志, 2023, 40(2): 295-302. doi: 10.7507/1001-5515.202210016 復制
0 引言
近年來，生命科學領域針對人體組織生物力學機制的研究逐漸興起，相關醫療器械的應用研究也受到關注。常規的基于臨床案例統計研究方法屬于定性分析范疇，其需要樣本多、耗時長，不利于新產品和技術的研發。某些具有破壞性質的臨床應用試驗研究，出于人道主義要求，也不易進行。因此，基于有限元仿真的量化數值分析方法逐漸成為生物力學研究的重要工具，例如骨科中對骨折狀態下骨頭的力學分析以及假體植入后關節的應力分布分析[1-2]；口腔醫學領域中也有學者研究義齒修復的有限元模擬方法以及探討牙齒和牙槽骨的合理有限元模型[3-6]。
在口腔正畸方面，由于口內結構復雜且組織種類多樣，有限元應用略為滯后[7]。近十年來，相關有限元研究主要集中在微支抗、隱形矯治器、托槽等正畸器械對牙齒的力學作用分析[8-12]，缺少對于口內其它組織的影響研究，關于唇頰軟組織的生物力學影響更未見報道。而在正畸臨床中，除了牙齒矯正這一主要效應，唇頰軟組織損傷、潰瘍和患者佩戴不適感等不利反應也是臨床診療方案設計中的關注重點[13-14]，甚至會影響患者對正畸方案的選擇。因此，開展正畸器械對于唇頰軟組織的生物力學影響研究具有重要的臨床實踐意義，也是新型器械研發的重要依據。相關的有限元分析，需要處理已有研究中未涉及的多個新問題，包括唇頰軟組織的材料表征、復雜口腔活動的計算工況設計、器械與唇頰軟組織的大面積接觸等。
針對以上問題，本文研究矯牙過程中托槽致唇頰軟組織非線性反應的高精度有限元建模和分析方法，具體包括：① 根據唇頰生物組成特點，優化選取二階奧格登（Ogden）模型，對類脂肪材料的唇頰軟組織進行表征；② 根據口腔活動特點，建立托槽介入和正交滑動的兩階段仿真模型，并對關鍵接觸參數進行優化設置；③ 采用整體模型—子模型的兩層次分析方法，基于整體模型計算得到的位移邊界，實現子模型高精度應變的高效求解。最終對正畸過程中的四種典型牙齒形態開展計算，并與臨床觀測進行對比驗證。綜上所述，本文研究通過探明常用矯牙托槽與唇頰軟組織之間的作用機制以及造成患者不適的原因，以期為新型矯牙產品的研發提供有益參考。
1 問題背景
1.1 矯牙托槽的工作原理
托槽是用于牙齒正畸矯治的一種精密機械裝置，具有體積小和制造精度高的特點。如圖1所示為當前廣泛使用的方形托槽模型，基底長寬約為4 mm和3 mm，托槽高約2 mm。臨床使用時，托槽粘接在牙齒表面，弓絲穿過槽溝壓緊托槽，矯治力傳遞到所附著的牙齒，導致齒根處的牙槽骨受力發生骨細胞破壞和再生，牙槽骨重塑形態使得牙齒隨之移動，從而達到矯正牙齒位置的目的。
圖1
方形托槽結構
Figure1.
Structure of a square bracket
1.2 托槽矯治的臨床問題
在牙齒正畸的醫療實踐中，由于骨細胞重組過程緩慢，因此整個矯治過程耗時較長，通常為1～2年。在這一長周期內，患者會有較強的不適臨床反應。一方面，由于托槽介入唇頰與牙齒之間，導致唇頰軟組織變形張緊，刺激神經反射，引起患者長期的不舒適感。另一方面，在說話和咀嚼等大量的口腔活動中，托槽邊緣的刮蹭會造成唇頰軟組織的局部應力集中，嚴重時會使唇頰軟組織損傷，甚至發生潰瘍。
針對這些不利反應，現有的正畸研究采用案例統計的方法進行定性綜合分析，從臨床表征上評估矯治器的影響。為了更有效地對比不同矯治器，有時還需要在同一患者口內進行分口試驗，在不同牙齒上安裝不同的矯治器，然而這一試驗需要通過嚴格的醫療倫理審查。
因此，為了更有效地對矯治器進行分析評估，以及推進新型器械研發，有必要采用數值仿真手段，從人體生物力學角度開展定量分析。
本文將基于有限元分析軟件ANSYS Workbench（版本：19.2, ANSYS Inc., 美國），采用三維非線性有限元模擬方法，分析方形托槽對唇頰軟組織的力學影響，從而對正畸過程中的患者不舒適度進行定量評估。該方法已成功應用于本課題組研發的新型球形托槽的效果評估[15]，論證了其相對于方形托槽的優越性，該部分工作將另文論述。
2 唇頰軟組織材料的有限元表征
2.1 唇頰軟組織的組成和材料特點
根據人體唇頰部的局部生理解剖結構[16]，唇頰軟組織由皮膚、脂肪、肌肉和粘膜等組成，其中脂肪的體積占比達到一半以上。考慮到唇頰不同位置處的組織結構比例不同，以及脂肪的材料特性相對于其它軟組織得到了更為廣泛深入的實驗和理論研究[17]，因此為了建模和分析方便，本文中將唇頰軟組織統一近似為類脂肪材料進行處理。
Comley等[18-19]通過脂肪微觀結構研究，在37 ℃下將其理想化為不可壓縮的無粘性流體，低應變下彈性模量約為1 kPa，并通過實驗建立了脂肪應力—應變曲線。Alkhouli等[20]通過實驗研究，提出了脂肪材料特性的兩條漸近線理論，即初始彈性模量為1 kPa和0.3應變后達到最終彈性模量80 kPa，另外還發現應變高于0.3后，脂肪出現撕裂的風險很高，建議該值為脂肪組織的損傷閾值。因此，本文計算研究也將重點考察唇頰軟組織的應變分布結果。
2.2 有限元中的材料本構模型
針對脂肪彈性模量隨應變增加的特點，有限元中采用超彈性材料進行描述，常用的本構模型包括基于應變張量不變量定義的新胡克（Neo-Hookean）模型和門尼-里夫林（Mooney-Rivlin）模型，以及基于主伸長率定義的Ogden模型等[17]。
Sims等[21]的脂肪壓痕實驗顯示，Neo-Hookean模型更加適用于表達具有彈塑性變形行為的低中應變下的力學響應。Geerligs等[22]的脂肪剪切實驗表明，Mooney-Rivlin模型適用于卸載后能夠完全恢復形變的高應變情況。Comley等[19]通過壓縮實驗證明，Ogden模型不僅可以表達脂肪具有非線性行為和對稱的力學響應，還能夠體現其不可壓縮性。Mihai等[23]對比論證了Neo-Hookean模型和Mooney-Rivlin 模型不適用于在拉伸和剪切共同作用下的脂肪組織，而Ogden模型的模擬結果與實驗吻合良好，適用于三維有限元分析。
因此，本文參考材料實驗結果[18-20]，采用二階Ogden模型來描述唇頰軟組織的類脂肪材料。如圖2所示為最終擬合得到的材料單軸應力—應變曲線。由圖2可見，擬合得到的初始彈性模量約為1 kPa，且彈性模量隨著應變增加而變大，最終在0.32應變處達到最終彈性模量80 kPa，與材料特點相符。
圖2
擬合得到的單軸應力—應變曲線
Figure2.
Uniaxial stress-strain curve obtained by fitting
需要說明的是，Ogden模型有不同階次可供選擇，更高階次的模型對于實驗參考結果的擬合更準確，但有限元計算的收斂難度與耗時都會增加，尤其在與接觸、幾何等非線性耦合后更進一步加劇。因此，出于擬合精度與計算收斂性的綜合考慮，本文選用六參數的二階模型來進行材料模擬，相關的參數取值如下：材料常數MU1為37.524 Pa，材料常數A1為17.883，材料常數MU2為37.524 Pa，材料常數A2為17.883，不可壓縮參數D1為0 Pa，不可壓縮參數D2為0 Pa。
3 托槽介入與滑動計算分析
3.1 有限元計算方案
本研究采用擬動力分析方法，通過位移分步加載的方式，使用靜力分析模塊依次完成托槽介入和滑動兩個連續工況的有限元計算。
（1）幾何建模
本文選擇一名25歲女性患者作為研究對象，研究已獲得患者許可并簽署知情同意書，并已通過中山大學附屬第一醫院倫理委員會的倫理審查（批件號：倫審[2020]315-1號）。采用大視野錐形束投照計算機斷層掃描（cone beam computed tomography，CBCT）（i-CAT FLX, Kavo Inc.，美國）對其頜面部進行三維重建，采用口內掃描儀（Dentrix 60，Fussen Inc.，中國）掃描口內牙列。其中，CBCT能真實反映唇頰軟組織和牙根的形態和尺寸，口內掃描儀能精準獲取牙齒牙冠面的形態。后續通過醫學建模軟件Mimics（Materialise Inc.，比利時）和逆向工程軟件Geomagic Studio（Geomagic Inc.，美國）提取有限元分析所用的下唇組織（寬 × 高 × 厚 = 30 mm × 25 mm × 5 mm）和左下頜中切牙的三維幾何模型，導入有限元分析軟件ANSYS Workbench（版本：19.2，ANSYS Inc.，美國）完成幾何建模。然后，本文通過工業建模軟件SolidWorks（版本：2018，Dassault Systemes Inc.，美國）完整還原方形托槽的真實形態，以STP格式導入有限元分析軟件ANSYS Workbench（版本：19.2，ANSYS Inc.，美國）完成幾何建模。
最終組合完成的下唇—托槽—左下頜中切牙的有限元模型和網格剖分如圖3所示。其中，對于需要獲得精確應變分布的唇頰軟組織，選取二次六面體單元；對于牙齒及托槽，其剛度大、變形小，為適應幾何形狀選取線性四面體單元。
圖3
分析模型與網格剖分
Figure3.
Analysis model and mesh generation
（2）材料屬性
除了前面已定義的唇頰軟組織非線性材料模型外，還需要定義牙齒和托槽的材料模型。由于牙齒和托槽的剛度遠大于唇頰軟組織的剛度，實際變形小，因此設置牙齒和托槽為均質、各向同性的線彈性材料，相關彈性模量和泊松比分別取值如下[24]：左下頜中切牙彈性模量為20 GPa，泊松比為0.3；托槽彈性模量為200 Gpa，泊松比為0.3。
（3）邊界條件
首先，將牙齒根部設置為固定支座。
其次，考慮到無托槽情況下的說話、咀嚼等口腔活動的一般特點，即唇頰軟組織在對牙齒表面有一定壓迫的前提下，沿著牙齒表面發生滑動，因此在安裝托槽后的連續計算工況設計如下：
① 托槽介入工況，模擬托槽介入唇頰軟組織和牙齒之間的情況。實現如下，將唇頰軟組織內表面與托槽頂面發生臨界接觸時的狀態設置為模型計算的初始狀態，然后讓唇頰軟組織沿著牙齒表面的垂向壓入（約3 mm），使得在扣除托槽介入唇頰軟組織后（約2 mm），仍有牙齒介入唇頰軟組織（約1 mm），從而實現與無托槽情況下相近的唇頰軟組織對牙齒表面的壓迫作用。
② 托槽滑動工況，模擬托槽與唇頰軟組織間相對滑動的情況，出于最不利的受力情況考慮，采用正交位移的設置方法。實現如下，讓唇頰軟組織分別沿著牙齒表面的橫向和豎向移動2 mm，其中橫向移動的工況用于分析牙齒的橫向偏轉畸變情況，豎向移動的工況用于分析牙齒的豎向偏轉畸變情況。需要說明的是，本文計算中采用的2 mm滑動量要小于現實生活中的最大滑動量。其原因在于，本文作者在前期計算中發現，采用20 mm滑動量計算得到的唇頰應變結果與2 mm滑動量的計算結果幾乎一致，均表現為在相同范圍內的波動變化，但前者的計算耗時高達后者的10倍。因此，出于計算精度與計算效率的綜合考慮，本文建議采用2 mm的滑動量進行計算。
3.2 關鍵接觸參數設置
在口腔活動中，唇頰軟組織、托槽和牙齒之間存在接觸。其中，托槽與牙齒之間是粘接關系，設置為綁定接觸；唇頰軟組織與托槽、唇頰軟組織與牙齒之間為摩擦型接觸，口腔粘膜和唾液的影響導致摩擦效應微小，摩擦系數均取為0.05[25]。
由于唇頰軟組織的剛度遠小于托槽和牙齒，因此唇頰軟組織—托槽、唇頰軟組織—牙齒的接觸行為均設置為非對稱接觸，其中托槽和牙齒定義為目標面，唇頰軟組織定義為接觸面，以滿足接觸面不能穿透目標面的計算要求。并且由于工況位移與接觸單元尺寸相當，因此接觸間均設置為有限滑動，關閉小滑動開關。
接觸計算中的穿透容差參數不應大于唇頰軟組織厚度的1/100與有限元網格大小的1/10，考慮到計算模型中軟組織厚度為5 mm，以及有限元網格大小約為1 mm，因此穿透容差參數取為0.05 mm。
此外，由于唇頰軟組織為非線性材料，法向接觸剛度設置為在每次迭代后自動更新，以保障計算精度和收斂性。
4 優化計算方法
4.1 整體模型計算方案
在以插值位移場為基礎的有限元分析中，應變結果的收斂速度低于位移結果。因此，本文研究的唇頰軟組織應變反應，原則上需要劃分精細網格來保證計算精度。此外，由于托槽部件具有多處微小的邊緣倒角，其與唇頰軟組織的接觸會引起劇烈的應力集中，因此有限元模型更需要精細化的分析網格，才能有效捕捉唇頰軟組織的最大應變反應。
傳統的方法是在整體模型上直接加密分析網格，然而相應的計算自由度增加會導致電腦資源需求加大、計算速度變慢、收斂困難等問題。以本文受試者的近遠中扭轉30°牙齒形態為例，采用配置霄龍32核處理器（AMD EPYC 7452 32-Core）和256 GB內存的工作站電腦，針對整體模型開展了唇頰軟組織網格分別劃分為1.0、0.9、0.8、···、0.3 mm共8種情況的計算分析。非線性有限元的計算耗時由1.0 mm情況下的30 min，急劇升至0.3 mm情況下的1 000 min以上，超過17 h。這種細化方案不利于多種牙齒形態的重復分析，需要優化調整。
4.2 兩層次計算方案
本文根據位移和應變結果的精度區別，選擇兩種不同規模和精細度的網格模型分別計算位移和應變。其中將粗網格的整體模型計算得到的位移結果，作為細網格的子模型的位移邊界條件，進一步計算得到子模型的高精度應變結果。
如圖4所示，為本文采用的整體模型—子模型兩層次計算方案示意。其中，根據圣維南原理，接觸區域處的應力集中對遠離該位置處的位移和應變擾動影響微小。因此，本文以接觸區域為中心，將整體模型的唇頰軟組織切割出特定的區域（長 × 寬 × 厚 = 9 mm × 6 mm × 2.5 mm），并結合托槽的部分切割模型，組成第二層次分析用的子模型，然后對子模型進行更小尺寸的有限元網格剖分。
圖4
兩層次計算方案示意
Figure4.
Schematic diagram of two-level calculation scheme
4.3 方案對比
為了全面比對兩種計算方案的效果，以及確定合理的子模型網格精細度，本文開展了不同網格情況下的計算對比，最終得到的唇頰軟組織最大應變結果和總計算耗時如圖5所示。其中，黑色柱表示使用整體模型計算方案的結果，名稱“整1.0、整0.9、整0.8、···、整0.3”，分別對應唇頰軟組織網格尺寸為1.0、0.9、0.8、···、0.3 mm的情況；白色柱表示使用兩層次計算方案的結果，名稱“整1.0 +子0.3、整1.0 + 子0.2、整1.0 + 子0.1”，分別對應整體模型網格尺寸為1.0 mm和子模型網格尺寸為0.3 mm、整體模型網格尺寸為1.0 mm和子模型網格尺寸為0.2 mm、整體模型網格尺寸為1.0 mm和子模型網格尺寸為0.1 mm的情況。由圖5可見，采用整體模型1.0 mm + 子模型0.3 mm的兩層次計算方案，可在60 min的計算耗時內，得到收斂的應變結果，相對于整體模型計算方案，計算效率得到大幅提高。本文后續計算即采用該種網格設置進行。
圖5
不同網格尺寸下兩種計算方案的最大應變結果和計算耗時
Figure5.
Maximum strain result and computation time for the two calculation schemes with different meshes
5 典型牙齒形態分析結果
5.1 典型牙齒形態選擇
在本文受試者的近遠中扭轉30°牙齒形態之外，本文結合臨床上正畸患者切牙常見的幾種排列狀態[26]，另外補充了三種牙齒形態進行計算。最終，本文計算的四種典型牙齒形態為：無偏轉畸變（牙列整齊）、兩種橫向偏轉畸變（近遠中扭轉15°、近遠中扭轉30°）、一種豎向偏轉畸變（唇側傾斜15°）。如圖6所示，展示了四種典型牙齒形態分析模型的俯視圖和側視圖。四種形態下的計算結果，可用于量化評估牙齒的橫向正畸過程和豎向正畸過程的唇頰軟組織反應大小。
圖6
四種典型牙齒形態的計算模型
Figure6.
Calculation models for four typical tooth morphologies
5.2 分析結果
分別對這四種牙齒形態模型進行托槽介入和滑動兩個連續工況的計算，在有限元分析軟件中將作用全過程分解為40個分析步，其中前12步為托槽介入的過程，即每步對應0.25 mm的位移，后28步為托槽滑動的過程，即每步對應約0.071 mm的位移。得到唇頰軟組織在兩個工況下的最不利應變分布如圖7所示，其中最大應變數值隨各分析步變化的歷史如圖8所示。
圖7
四種典型牙齒形態下的最不利應變分布結果
Figure7.
Results of the most unfavourable strain distribution for four typical tooth morphologies
圖8
四種典型牙齒形態下的最大應變隨分析步變化曲線
Figure8.
Curves of maximum strain changing with loading step for four typical tooth morphologies
（1）在托槽介入階段，唇頰軟組織最大應變隨著壓入深度增加而逐步增大，高應變沿著托槽的尖銳邊緣分布，該處接觸存在明顯的應力集中現象。
（2）在托槽滑動階段，唇頰軟組織最大應變較介入階段有所加大，呈現波動變化，高應變同樣沿著托槽的尖銳邊緣分布，存在較強的應力集中。其中，對于唇側傾斜15°的牙齒形態，滑動會引起托槽介入深度的進一步加大，導致最大應變波動上揚。
（3）三種牙齒形態的變化，即從近遠中扭轉30°到近遠中扭轉15°，再到牙列整齊，可以代表牙齒的橫向正畸過程，唇頰軟組織最大應變由0.31減小到0.21，再到0.15，表示隨著橫向正畸的進行，唇頰軟組織的力學反應由大到小，患者不舒適感逐漸降低。其中，近遠中扭轉30°情況下應變為0.31，說明這種情況下有可能出現組織損傷[20]，存在潰瘍風險。
（4）兩種牙齒形態的變化，即唇側傾斜15°到牙列整齊，可以代表牙齒的豎向正畸過程，唇頰軟組織最大應變由0.22減小到0.15，表示隨著豎向正畸的進行，患者的不舒適感同樣會降低。
如圖9所示為該名患者佩戴托槽一周后的臨床照片，圖9中黃色箭頭標示出了托槽在嘴唇上留下的壓痕位置及對應關系，觀察可見：
圖9
矯治開始一周后的唇部牽張和未牽張照片
Figure9.
Retracted and unretracted photographs of the lip one week after the start of orthodontic treatment
（1）圖9右圖的自然未牽張狀態下，下唇內表面呈現了較明顯的托槽壓痕。其中，左下頜中切牙（近遠中扭轉30°）托槽對應的更深壓痕形狀與托槽的單側邊緣輪廓接近，左下頜側切牙（牙列整齊）托槽對應的等深度壓痕形狀與托槽的兩側邊緣輪廓接近，均與如圖7所示應變分布結果相符。
（2）圖9左圖的牽張狀態下，口腔粘膜下的唇頰軟組織健康狀態得以顯現。其中，左下頜中切牙托槽對應的下唇粘膜下出現了唇頰軟組織局部充血紅腫，顯示組織損傷，存在潰瘍風險，而左下頜側切牙托槽對應的下唇粘膜下的唇頰軟組織仍保持正常狀態。Leiva等[27]開展的雙盲隨機臨床試驗表明，在開始矯治的第一個月內，唇頰軟組織損傷和潰瘍高發，而后隨著矯治過程的進行，這一狀況逐步改善。前述有限元的計算結論與該臨床試驗結論及本案例觀察結果相符。
6 結語
本文采用耦合非線性有限元方法，首次對口腔活動中托槽導致的唇頰軟組織力學反應進行量化分析。其中，通過醫學掃描建立唇頰—托槽—牙齒的三維實體模型，引入Ogden模型進行唇頰軟組織的類脂肪材料表征，并采用整體模型—子模型的兩層次分析方法實現高效率計算，得到與臨床觀測相吻合的結果。該方法可為國內外口腔正畸醫療領域的相關量化分析研究提供有益的參考，可進一步應用于新型矯治裝置的產品研發分析。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：華嘉皓主要負責有限元模型建立、數據處理與分析、論文初稿撰寫；吉利主要負責研究技術路線輔助制定、人體試驗方案制定；代清源與梁振宇主要負責有限元模型計算；郭龍妹主要負責人體試驗開展和照片整理；陳太聰主要負責研究技術路線制定、有限元方案設計、論文審閱修訂。
倫理聲明：本研究通過了中山大學附屬第一醫院倫理委員會的審批（批文編號：倫審［2020］315-1號）。
0 引言
近年來，生命科學領域針對人體組織生物力學機制的研究逐漸興起，相關醫療器械的應用研究也受到關注。常規的基于臨床案例統計研究方法屬于定性分析范疇，其需要樣本多、耗時長，不利于新產品和技術的研發。某些具有破壞性質的臨床應用試驗研究，出于人道主義要求，也不易進行。因此，基于有限元仿真的量化數值分析方法逐漸成為生物力學研究的重要工具，例如骨科中對骨折狀態下骨頭的力學分析以及假體植入后關節的應力分布分析[1-2]；口腔醫學領域中也有學者研究義齒修復的有限元模擬方法以及探討牙齒和牙槽骨的合理有限元模型[3-6]。
在口腔正畸方面，由于口內結構復雜且組織種類多樣，有限元應用略為滯后[7]。近十年來，相關有限元研究主要集中在微支抗、隱形矯治器、托槽等正畸器械對牙齒的力學作用分析[8-12]，缺少對于口內其它組織的影響研究，關于唇頰軟組織的生物力學影響更未見報道。而在正畸臨床中，除了牙齒矯正這一主要效應，唇頰軟組織損傷、潰瘍和患者佩戴不適感等不利反應也是臨床診療方案設計中的關注重點[13-14]，甚至會影響患者對正畸方案的選擇。因此，開展正畸器械對于唇頰軟組織的生物力學影響研究具有重要的臨床實踐意義，也是新型器械研發的重要依據。相關的有限元分析，需要處理已有研究中未涉及的多個新問題，包括唇頰軟組織的材料表征、復雜口腔活動的計算工況設計、器械與唇頰軟組織的大面積接觸等。
針對以上問題，本文研究矯牙過程中托槽致唇頰軟組織非線性反應的高精度有限元建模和分析方法，具體包括：① 根據唇頰生物組成特點，優化選取二階奧格登（Ogden）模型，對類脂肪材料的唇頰軟組織進行表征；② 根據口腔活動特點，建立托槽介入和正交滑動的兩階段仿真模型，并對關鍵接觸參數進行優化設置；③ 采用整體模型—子模型的兩層次分析方法，基于整體模型計算得到的位移邊界，實現子模型高精度應變的高效求解。最終對正畸過程中的四種典型牙齒形態開展計算，并與臨床觀測進行對比驗證。綜上所述，本文研究通過探明常用矯牙托槽與唇頰軟組織之間的作用機制以及造成患者不適的原因，以期為新型矯牙產品的研發提供有益參考。
1 問題背景
1.1 矯牙托槽的工作原理
托槽是用于牙齒正畸矯治的一種精密機械裝置，具有體積小和制造精度高的特點。如圖1所示為當前廣泛使用的方形托槽模型，基底長寬約為4 mm和3 mm，托槽高約2 mm。臨床使用時，托槽粘接在牙齒表面，弓絲穿過槽溝壓緊托槽，矯治力傳遞到所附著的牙齒，導致齒根處的牙槽骨受力發生骨細胞破壞和再生，牙槽骨重塑形態使得牙齒隨之移動，從而達到矯正牙齒位置的目的。
圖1
方形托槽結構
Figure1.
Structure of a square bracket
1.2 托槽矯治的臨床問題
在牙齒正畸的醫療實踐中，由于骨細胞重組過程緩慢，因此整個矯治過程耗時較長，通常為1～2年。在這一長周期內，患者會有較強的不適臨床反應。一方面，由于托槽介入唇頰與牙齒之間，導致唇頰軟組織變形張緊，刺激神經反射，引起患者長期的不舒適感。另一方面，在說話和咀嚼等大量的口腔活動中，托槽邊緣的刮蹭會造成唇頰軟組織的局部應力集中，嚴重時會使唇頰軟組織損傷，甚至發生潰瘍。
針對這些不利反應，現有的正畸研究采用案例統計的方法進行定性綜合分析，從臨床表征上評估矯治器的影響。為了更有效地對比不同矯治器，有時還需要在同一患者口內進行分口試驗，在不同牙齒上安裝不同的矯治器，然而這一試驗需要通過嚴格的醫療倫理審查。
因此，為了更有效地對矯治器進行分析評估，以及推進新型器械研發，有必要采用數值仿真手段，從人體生物力學角度開展定量分析。
本文將基于有限元分析軟件ANSYS Workbench（版本：19.2, ANSYS Inc., 美國），采用三維非線性有限元模擬方法，分析方形托槽對唇頰軟組織的力學影響，從而對正畸過程中的患者不舒適度進行定量評估。該方法已成功應用于本課題組研發的新型球形托槽的效果評估[15]，論證了其相對于方形托槽的優越性，該部分工作將另文論述。
2 唇頰軟組織材料的有限元表征
2.1 唇頰軟組織的組成和材料特點
根據人體唇頰部的局部生理解剖結構[16]，唇頰軟組織由皮膚、脂肪、肌肉和粘膜等組成，其中脂肪的體積占比達到一半以上。考慮到唇頰不同位置處的組織結構比例不同，以及脂肪的材料特性相對于其它軟組織得到了更為廣泛深入的實驗和理論研究[17]，因此為了建模和分析方便，本文中將唇頰軟組織統一近似為類脂肪材料進行處理。
Comley等[18-19]通過脂肪微觀結構研究，在37 ℃下將其理想化為不可壓縮的無粘性流體，低應變下彈性模量約為1 kPa，并通過實驗建立了脂肪應力—應變曲線。Alkhouli等[20]通過實驗研究，提出了脂肪材料特性的兩條漸近線理論，即初始彈性模量為1 kPa和0.3應變后達到最終彈性模量80 kPa，另外還發現應變高于0.3后，脂肪出現撕裂的風險很高，建議該值為脂肪組織的損傷閾值。因此，本文計算研究也將重點考察唇頰軟組織的應變分布結果。
2.2 有限元中的材料本構模型
針對脂肪彈性模量隨應變增加的特點，有限元中采用超彈性材料進行描述，常用的本構模型包括基于應變張量不變量定義的新胡克（Neo-Hookean）模型和門尼-里夫林（Mooney-Rivlin）模型，以及基于主伸長率定義的Ogden模型等[17]。
Sims等[21]的脂肪壓痕實驗顯示，Neo-Hookean模型更加適用于表達具有彈塑性變形行為的低中應變下的力學響應。Geerligs等[22]的脂肪剪切實驗表明，Mooney-Rivlin模型適用于卸載后能夠完全恢復形變的高應變情況。Comley等[19]通過壓縮實驗證明，Ogden模型不僅可以表達脂肪具有非線性行為和對稱的力學響應，還能夠體現其不可壓縮性。Mihai等[23]對比論證了Neo-Hookean模型和Mooney-Rivlin 模型不適用于在拉伸和剪切共同作用下的脂肪組織，而Ogden模型的模擬結果與實驗吻合良好，適用于三維有限元分析。
因此，本文參考材料實驗結果[18-20]，采用二階Ogden模型來描述唇頰軟組織的類脂肪材料。如圖2所示為最終擬合得到的材料單軸應力—應變曲線。由圖2可見，擬合得到的初始彈性模量約為1 kPa，且彈性模量隨著應變增加而變大，最終在0.32應變處達到最終彈性模量80 kPa，與材料特點相符。
圖2
擬合得到的單軸應力—應變曲線
Figure2.
Uniaxial stress-strain curve obtained by fitting
需要說明的是，Ogden模型有不同階次可供選擇，更高階次的模型對于實驗參考結果的擬合更準確，但有限元計算的收斂難度與耗時都會增加，尤其在與接觸、幾何等非線性耦合后更進一步加劇。因此，出于擬合精度與計算收斂性的綜合考慮，本文選用六參數的二階模型來進行材料模擬，相關的參數取值如下：材料常數MU1為37.524 Pa，材料常數A1為17.883，材料常數MU2為37.524 Pa，材料常數A2為17.883，不可壓縮參數D1為0 Pa，不可壓縮參數D2為0 Pa。
3 托槽介入與滑動計算分析
3.1 有限元計算方案
本研究采用擬動力分析方法，通過位移分步加載的方式，使用靜力分析模塊依次完成托槽介入和滑動兩個連續工況的有限元計算。
（1）幾何建模
本文選擇一名25歲女性患者作為研究對象，研究已獲得患者許可并簽署知情同意書，并已通過中山大學附屬第一醫院倫理委員會的倫理審查（批件號：倫審[2020]315-1號）。采用大視野錐形束投照計算機斷層掃描（cone beam computed tomography，CBCT）（i-CAT FLX, Kavo Inc.，美國）對其頜面部進行三維重建，采用口內掃描儀（Dentrix 60，Fussen Inc.，中國）掃描口內牙列。其中，CBCT能真實反映唇頰軟組織和牙根的形態和尺寸，口內掃描儀能精準獲取牙齒牙冠面的形態。后續通過醫學建模軟件Mimics（Materialise Inc.，比利時）和逆向工程軟件Geomagic Studio（Geomagic Inc.，美國）提取有限元分析所用的下唇組織（寬 × 高 × 厚 = 30 mm × 25 mm × 5 mm）和左下頜中切牙的三維幾何模型，導入有限元分析軟件ANSYS Workbench（版本：19.2，ANSYS Inc.，美國）完成幾何建模。然后，本文通過工業建模軟件SolidWorks（版本：2018，Dassault Systemes Inc.，美國）完整還原方形托槽的真實形態，以STP格式導入有限元分析軟件ANSYS Workbench（版本：19.2，ANSYS Inc.，美國）完成幾何建模。
最終組合完成的下唇—托槽—左下頜中切牙的有限元模型和網格剖分如圖3所示。其中，對于需要獲得精確應變分布的唇頰軟組織，選取二次六面體單元；對于牙齒及托槽，其剛度大、變形小，為適應幾何形狀選取線性四面體單元。
圖3
分析模型與網格剖分
Figure3.
Analysis model and mesh generation
（2）材料屬性
除了前面已定義的唇頰軟組織非線性材料模型外，還需要定義牙齒和托槽的材料模型。由于牙齒和托槽的剛度遠大于唇頰軟組織的剛度，實際變形小，因此設置牙齒和托槽為均質、各向同性的線彈性材料，相關彈性模量和泊松比分別取值如下[24]：左下頜中切牙彈性模量為20 GPa，泊松比為0.3；托槽彈性模量為200 Gpa，泊松比為0.3。
（3）邊界條件
首先，將牙齒根部設置為固定支座。
其次，考慮到無托槽情況下的說話、咀嚼等口腔活動的一般特點，即唇頰軟組織在對牙齒表面有一定壓迫的前提下，沿著牙齒表面發生滑動，因此在安裝托槽后的連續計算工況設計如下：
① 托槽介入工況，模擬托槽介入唇頰軟組織和牙齒之間的情況。實現如下，將唇頰軟組織內表面與托槽頂面發生臨界接觸時的狀態設置為模型計算的初始狀態，然后讓唇頰軟組織沿著牙齒表面的垂向壓入（約3 mm），使得在扣除托槽介入唇頰軟組織后（約2 mm），仍有牙齒介入唇頰軟組織（約1 mm），從而實現與無托槽情況下相近的唇頰軟組織對牙齒表面的壓迫作用。
② 托槽滑動工況，模擬托槽與唇頰軟組織間相對滑動的情況，出于最不利的受力情況考慮，采用正交位移的設置方法。實現如下，讓唇頰軟組織分別沿著牙齒表面的橫向和豎向移動2 mm，其中橫向移動的工況用于分析牙齒的橫向偏轉畸變情況，豎向移動的工況用于分析牙齒的豎向偏轉畸變情況。需要說明的是，本文計算中采用的2 mm滑動量要小于現實生活中的最大滑動量。其原因在于，本文作者在前期計算中發現，采用20 mm滑動量計算得到的唇頰應變結果與2 mm滑動量的計算結果幾乎一致，均表現為在相同范圍內的波動變化，但前者的計算耗時高達后者的10倍。因此，出于計算精度與計算效率的綜合考慮，本文建議采用2 mm的滑動量進行計算。
3.2 關鍵接觸參數設置
在口腔活動中，唇頰軟組織、托槽和牙齒之間存在接觸。其中，托槽與牙齒之間是粘接關系，設置為綁定接觸；唇頰軟組織與托槽、唇頰軟組織與牙齒之間為摩擦型接觸，口腔粘膜和唾液的影響導致摩擦效應微小，摩擦系數均取為0.05[25]。
由于唇頰軟組織的剛度遠小于托槽和牙齒，因此唇頰軟組織—托槽、唇頰軟組織—牙齒的接觸行為均設置為非對稱接觸，其中托槽和牙齒定義為目標面，唇頰軟組織定義為接觸面，以滿足接觸面不能穿透目標面的計算要求。并且由于工況位移與接觸單元尺寸相當，因此接觸間均設置為有限滑動，關閉小滑動開關。
接觸計算中的穿透容差參數不應大于唇頰軟組織厚度的1/100與有限元網格大小的1/10，考慮到計算模型中軟組織厚度為5 mm，以及有限元網格大小約為1 mm，因此穿透容差參數取為0.05 mm。
此外，由于唇頰軟組織為非線性材料，法向接觸剛度設置為在每次迭代后自動更新，以保障計算精度和收斂性。
4 優化計算方法
4.1 整體模型計算方案
在以插值位移場為基礎的有限元分析中，應變結果的收斂速度低于位移結果。因此，本文研究的唇頰軟組織應變反應，原則上需要劃分精細網格來保證計算精度。此外，由于托槽部件具有多處微小的邊緣倒角，其與唇頰軟組織的接觸會引起劇烈的應力集中，因此有限元模型更需要精細化的分析網格，才能有效捕捉唇頰軟組織的最大應變反應。
傳統的方法是在整體模型上直接加密分析網格，然而相應的計算自由度增加會導致電腦資源需求加大、計算速度變慢、收斂困難等問題。以本文受試者的近遠中扭轉30°牙齒形態為例，采用配置霄龍32核處理器（AMD EPYC 7452 32-Core）和256 GB內存的工作站電腦，針對整體模型開展了唇頰軟組織網格分別劃分為1.0、0.9、0.8、···、0.3 mm共8種情況的計算分析。非線性有限元的計算耗時由1.0 mm情況下的30 min，急劇升至0.3 mm情況下的1 000 min以上，超過17 h。這種細化方案不利于多種牙齒形態的重復分析，需要優化調整。
4.2 兩層次計算方案
本文根據位移和應變結果的精度區別，選擇兩種不同規模和精細度的網格模型分別計算位移和應變。其中將粗網格的整體模型計算得到的位移結果，作為細網格的子模型的位移邊界條件，進一步計算得到子模型的高精度應變結果。
如圖4所示，為本文采用的整體模型—子模型兩層次計算方案示意。其中，根據圣維南原理，接觸區域處的應力集中對遠離該位置處的位移和應變擾動影響微小。因此，本文以接觸區域為中心，將整體模型的唇頰軟組織切割出特定的區域（長 × 寬 × 厚 = 9 mm × 6 mm × 2.5 mm），并結合托槽的部分切割模型，組成第二層次分析用的子模型，然后對子模型進行更小尺寸的有限元網格剖分。
圖4
兩層次計算方案示意
Figure4.
Schematic diagram of two-level calculation scheme
4.3 方案對比
為了全面比對兩種計算方案的效果，以及確定合理的子模型網格精細度，本文開展了不同網格情況下的計算對比，最終得到的唇頰軟組織最大應變結果和總計算耗時如圖5所示。其中，黑色柱表示使用整體模型計算方案的結果，名稱“整1.0、整0.9、整0.8、···、整0.3”，分別對應唇頰軟組織網格尺寸為1.0、0.9、0.8、···、0.3 mm的情況；白色柱表示使用兩層次計算方案的結果，名稱“整1.0 +子0.3、整1.0 + 子0.2、整1.0 + 子0.1”，分別對應整體模型網格尺寸為1.0 mm和子模型網格尺寸為0.3 mm、整體模型網格尺寸為1.0 mm和子模型網格尺寸為0.2 mm、整體模型網格尺寸為1.0 mm和子模型網格尺寸為0.1 mm的情況。由圖5可見，采用整體模型1.0 mm + 子模型0.3 mm的兩層次計算方案，可在60 min的計算耗時內，得到收斂的應變結果，相對于整體模型計算方案，計算效率得到大幅提高。本文后續計算即采用該種網格設置進行。
圖5
不同網格尺寸下兩種計算方案的最大應變結果和計算耗時
Figure5.
Maximum strain result and computation time for the two calculation schemes with different meshes
5 典型牙齒形態分析結果
5.1 典型牙齒形態選擇
在本文受試者的近遠中扭轉30°牙齒形態之外，本文結合臨床上正畸患者切牙常見的幾種排列狀態[26]，另外補充了三種牙齒形態進行計算。最終，本文計算的四種典型牙齒形態為：無偏轉畸變（牙列整齊）、兩種橫向偏轉畸變（近遠中扭轉15°、近遠中扭轉30°）、一種豎向偏轉畸變（唇側傾斜15°）。如圖6所示，展示了四種典型牙齒形態分析模型的俯視圖和側視圖。四種形態下的計算結果，可用于量化評估牙齒的橫向正畸過程和豎向正畸過程的唇頰軟組織反應大小。
圖6
四種典型牙齒形態的計算模型
Figure6.
Calculation models for four typical tooth morphologies
5.2 分析結果
分別對這四種牙齒形態模型進行托槽介入和滑動兩個連續工況的計算，在有限元分析軟件中將作用全過程分解為40個分析步，其中前12步為托槽介入的過程，即每步對應0.25 mm的位移，后28步為托槽滑動的過程，即每步對應約0.071 mm的位移。得到唇頰軟組織在兩個工況下的最不利應變分布如圖7所示，其中最大應變數值隨各分析步變化的歷史如圖8所示。
圖7
四種典型牙齒形態下的最不利應變分布結果
Figure7.
Results of the most unfavourable strain distribution for four typical tooth morphologies
圖8
四種典型牙齒形態下的最大應變隨分析步變化曲線
Figure8.
Curves of maximum strain changing with loading step for four typical tooth morphologies
（1）在托槽介入階段，唇頰軟組織最大應變隨著壓入深度增加而逐步增大，高應變沿著托槽的尖銳邊緣分布，該處接觸存在明顯的應力集中現象。
（2）在托槽滑動階段，唇頰軟組織最大應變較介入階段有所加大，呈現波動變化，高應變同樣沿著托槽的尖銳邊緣分布，存在較強的應力集中。其中，對于唇側傾斜15°的牙齒形態，滑動會引起托槽介入深度的進一步加大，導致最大應變波動上揚。
（3）三種牙齒形態的變化，即從近遠中扭轉30°到近遠中扭轉15°，再到牙列整齊，可以代表牙齒的橫向正畸過程，唇頰軟組織最大應變由0.31減小到0.21，再到0.15，表示隨著橫向正畸的進行，唇頰軟組織的力學反應由大到小，患者不舒適感逐漸降低。其中，近遠中扭轉30°情況下應變為0.31，說明這種情況下有可能出現組織損傷[20]，存在潰瘍風險。
（4）兩種牙齒形態的變化，即唇側傾斜15°到牙列整齊，可以代表牙齒的豎向正畸過程，唇頰軟組織最大應變由0.22減小到0.15，表示隨著豎向正畸的進行，患者的不舒適感同樣會降低。
如圖9所示為該名患者佩戴托槽一周后的臨床照片，圖9中黃色箭頭標示出了托槽在嘴唇上留下的壓痕位置及對應關系，觀察可見：
圖9
矯治開始一周后的唇部牽張和未牽張照片
Figure9.
Retracted and unretracted photographs of the lip one week after the start of orthodontic treatment
（1）圖9右圖的自然未牽張狀態下，下唇內表面呈現了較明顯的托槽壓痕。其中，左下頜中切牙（近遠中扭轉30°）托槽對應的更深壓痕形狀與托槽的單側邊緣輪廓接近，左下頜側切牙（牙列整齊）托槽對應的等深度壓痕形狀與托槽的兩側邊緣輪廓接近，均與如圖7所示應變分布結果相符。
（2）圖9左圖的牽張狀態下，口腔粘膜下的唇頰軟組織健康狀態得以顯現。其中，左下頜中切牙托槽對應的下唇粘膜下出現了唇頰軟組織局部充血紅腫，顯示組織損傷，存在潰瘍風險，而左下頜側切牙托槽對應的下唇粘膜下的唇頰軟組織仍保持正常狀態。Leiva等[27]開展的雙盲隨機臨床試驗表明，在開始矯治的第一個月內，唇頰軟組織損傷和潰瘍高發，而后隨著矯治過程的進行，這一狀況逐步改善。前述有限元的計算結論與該臨床試驗結論及本案例觀察結果相符。
6 結語
本文采用耦合非線性有限元方法，首次對口腔活動中托槽導致的唇頰軟組織力學反應進行量化分析。其中，通過醫學掃描建立唇頰—托槽—牙齒的三維實體模型，引入Ogden模型進行唇頰軟組織的類脂肪材料表征，并采用整體模型—子模型的兩層次分析方法實現高效率計算，得到與臨床觀測相吻合的結果。該方法可為國內外口腔正畸醫療領域的相關量化分析研究提供有益的參考，可進一步應用于新型矯治裝置的產品研發分析。
重要聲明
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
作者貢獻聲明：華嘉皓主要負責有限元模型建立、數據處理與分析、論文初稿撰寫；吉利主要負責研究技術路線輔助制定、人體試驗方案制定；代清源與梁振宇主要負責有限元模型計算；郭龍妹主要負責人體試驗開展和照片整理；陳太聰主要負責研究技術路線制定、有限元方案設計、論文審閱修訂。
倫理聲明：本研究通過了中山大學附屬第一醫院倫理委員會的審批（批文編號：倫審［2020］315-1號）。
