牙頜模型作為口頜系統與咀嚼機器人系統的代表性部件，其建模方式是影響此類多體系統動力學模型準確度的重要因素。本文以咀嚼機器人右側第一磨牙為研究對象，提出一種由 V 型面與球面構成的點接觸高副等效建模方式。首先，利用有限單元法對原始模型的 3 種靜態接觸工況（咬入接觸、正中接觸、咬出接觸）與動態咀嚼工況進行咬合動力學分析，得到期望咬合力；其次，運用赫茲接觸模型建立點接觸高副等效模型靜態接觸工況的咬合力解析式，根據期望咬合力對解析式內法向量與咬合剛度等參數進行等效設計；最后，對等效模型在靜態接觸與動態咀嚼工況中的咬合力表現進行評價。本文結果表明，等效模型能夠實現靜態接觸工況 8 項期望咬合力的等效，同時對動態咀嚼工況初期與末期咬合力變化吻合良好，中期引入一定程度的沖擊，但可通過后續軌跡規劃進行弱化。本文所提出的牙頜模型等效建模方案，進一步提升了多體系統動力學模型準確性，也為人體其他復雜接觸的動力學建模提供了新思路。
引用本文： 秦文龍, 叢明, 任翔, 劉冬. 基于咬合動力學的磨牙點接觸高副等效建模與評價. 生物醫學工程學雜志, 2020, 37(4): 614-621. doi: 10.7507/1001-5515.201906021 復制
版權信息： ?四川大學華西醫院華西期刊社《生物醫學工程學雜志》版權所有，未經授權不得轉載、改編
引言
研究人體口頜系統的肌力模式、關節受力與咬合性能問題，離不開準確的多體系統動力學模型。同樣地，研究咀嚼機器人系統（仿生口頜系統）的咀嚼運動與咬合力控制問題亦然。近年來，牙頜模型作為口腔代表性部件正越來越多地作為子結構方式引入到多體系統動力學建模中，以進一步提高系統整體動力學模型質量，而不是僅僅依靠咬合力數值的方式[1-2]。目前，可采用的動力學建模策略包括多體系統仿真、有限元仿真以及聯合仿真。其中，有限元理論通過設置合理的網格密度，可解決原始咬合接觸的難收斂問題，最大限度地還原咬合面與食品形貌，模擬咬合面碰撞與食品變形行為，不足之處在于增加了上下頜、咀嚼肌肉以及顳下頜關節等非必要組件的網格建模工作量和后續動力學模型的解算難度[3-6]。為了盡可能地降低多體系統建模工作量與解算難度，同時也為了保證牙頜模型建模精度，廣泛采用的是多體系統仿真與聯合仿真的方案[7]。
多體系統仿真是指在多體系統仿真平臺下對牙頜模型進行建模，僅采取分析力學解算，具有建模簡單與求解速度快的優勢。可執行仿真的軟件平臺較多，例如：機械系統動力學仿真軟件 Adams（MSC Inc., 美國）、下頜運動評估軟件 Kinematics 3D（Zebris Medical Inc., 德國）、人體肌肉骨骼仿真軟件 OpenSim（The National Center for Simulation in Rehabilitation Research, 美國）和人體運動仿真軟件 Anybody Modeling System（AnyBody Technology Inc., 丹麥）[7-9]。但是，多體系統仿真針對接觸區域僅進行有限度的離散化，當面對磨牙所存在的尖、窩、脊等復雜接觸面形貌特征時，接觸判斷準則將難以正確判斷接觸是否發生，繼而發生接觸穿透現象，并致使模型求解失敗[10]。少部分研究人員對磨牙咬合動力學模型進行了簡化等效，例如：2018 年，王貴飛[11]在機械系統動力學仿真軟件 Adams（MSC Inc., 美國）內使用簡化的平面—球面—平面型式構造兩個點接觸高副（higher kinematic pair, HKP）描述上頜磨牙—食品—下頜磨牙三體接觸關系。然而，該點接觸高副并未體現出磨牙正中接觸的三向受力規律，也無法適用于動態咀嚼軌跡的加載，使得咀嚼機器人咬合力控制并不全面。針對復雜接觸問題，報道的另一種建模方案是將復雜接觸面進一步離散化，例如：Tuijt 等[12]與 Skipper 等[13]進行口頜多體系統建模時，將顳下頜關節接觸面離散為致密點云，用作后續接觸分析，實現在單一仿真環境下復雜接觸的多體動力學分析。可見，盡管多體系統仿真方案支持平臺較多，但存在牙頜模型建模精度或離散化程度不足的問題。
聯合仿真是指在多體系統仿真基礎上按需對牙頜模型進行有限元仿真，同時采用分析力學與彈性力學解算。醫用多體系統仿真軟件 ArtiSynth（University of British Columbia, 加拿大）能夠實現剛體與軟組織的聯合動力學仿真，即剛體動力學分析與有限元分析的聯合，現已應用于開閉口運動下肌力模式分析與下頜旋轉中心分析[14-16]。2008 年，Hannam 等[17]在運用該軟件對口頜骨—肌系統進行動力學分析時，預先根據仿真食品被搗碎過程中力學特性設計咬合力，最終口頜系統仿真結果與臨床數據吻合。同樣地，該仿真策略也已成功應用于顳下頜關節等類似復雜接觸問題。2019 年，Sagl 等[18]提出一種綜合動態咀嚼模型，由肌肉驅動的下頜剛體、顳下頜關節盤有限元模型與彈性基軟骨模型構成，結果顯示這種聯合模型能夠在保證仿真效果的同時，有效降低計算時間。然而，能夠實現聯合仿真建模的類似軟件平臺并不多。另一方面，聯合仿真中牙頜模型有限元分析部分更多是獨立地在有限元軟件平臺上實施。Benazzi 等[19]與 Qin 等[20]分別使用有限元軟件 Abaqus（SIMULIA Inc., 美國）對口頜系統中牙頜模型的咬合接觸問題進行建模，輸出咀嚼軌跡與咬合力之間的數學關系。Pascale 等[21]使用 C 語言結合可視化軟件 Vega（MultiGen Paradigm Inc., 美國）庫函數編寫牙頜有限元模型，同時采用分區的碰撞檢測算法計算咬合接觸受力。后續，該團隊使用擴展型拉普拉斯光滑算法對咬合面形貌進行簡化，進一步實現簡化的咬合接觸應力分析[22]，但并未給出將以上獨立的仿真數據應用（聯合）至多體系統模型中的方法。
基于以上研究，現有多體系統仿真平臺除非在分析力學基礎上增強模型離散能力，或補充彈性力學算法轉為聯合仿真平臺，否則將無法實現牙頜模型子結構模型的準確建模。為了解決這一問題，本文提出一種新的仿真策略，即設計一種適應現有多體系統仿真平臺離散程度的簡化牙頜模型，該模型由 V 型面與球面構成點接觸高副，使其分析力學結果與原始模型的彈性力學結果能夠保持等效。本研究期望可為多體系統的牙頜模型提供新的動力學建模策略，實現多體系統的動力學建模精度的進一步提升，另一方面通過本文的嘗試，也希望可以為人體其他復雜接觸的動力學建模提供新思路。
1 原始模型咬合動力學分析
1.1 三維與有限元模型的建立
第一磨牙（以下簡稱磨牙）作為牙頜模型中承受極大咬合力的牙位而常被用于各類咬合問題的研究，本文選取咀嚼機器人系統右側磨牙作為研究對象，進行咬合動力學簡化建模與評價[23]。面向義齒斷裂或磨損測試的仿生型咀嚼機器人系統一般選用標準牙頜模型作為基牙安裝于仿生口腔環境內。使用光柵掃描設備 AutoScan DS100（杭州先臨三維科技股份有限公司，中國）對日進牙科教學模型（Nissin Education Products，中國）的磨牙部件進行掃描，掃描精度 0.015 mm。其次，進行咬合面三維重建，并根據預備工藝進行預備體曲面造型，最終得到基牙與義齒的三維模型，如圖1 所示。定義靜坐標系{G}以上頜磨牙中央窩為原點，以磨牙后方、側方與上方為x、y與z軸正向。
圖1
咀嚼機器人磨牙三維與有限元模型
Figure1.
Geometrical and finite element models of molars assembled in chewing robot
使用網格劃分軟件 HyperMesh 14.0（Altair Engineering Inc., 美國）對磨牙全瓷冠以 0.1 mm 為全局單元尺寸進行網格劃分，形成網格模型。設置基牙材料為復合樹脂（密度 1 200 kg/m3，彈性模量 11.5 GPa，泊松比 0.24）；金屬冠材料為鈷鉻合金（密度 8 000 kg/m3，彈性模量 220 GPa，泊松比 0.30）；義齒類型為全瓷冠（密度 6 100 kg/m3，彈性模量 220 GPa，泊松比 0.30）[24-25]。將其導入有限元軟件 Abaqus 11.0（SIMULIA Inc., 美國）中，設置咬合面接觸類型為面—面接觸類型，解算機制為“運動學模式”。鑒于唾液的潤滑作用，本文不考慮切向咬合力。
1.2 靜態接觸與動態咀嚼工況設計
定義下頜磨牙動坐標系{L}隨下頜磨牙運動而運動，在正中接觸時，與靜坐標系{G}重合。下頜運動軌跡使用動坐標系在靜坐標系中的位姿進行描述。設計 3 種典型接觸工況，分別表征咬入接觸（intrusive contact）、正中接觸（centric occlusion）與咬出接觸（extrusive contact）。完整的靜態接觸工況定義應包含該接觸狀態對應咀嚼軌跡（僅咬合相）過程點在水平投影上的位置、初始頜位高度、垂直加載速率與仿真時間等數據。假設咀嚼軌跡水平投影滿足咬入角θ1 與咬出角θ2 為 10° 的角度特征，咬入接觸、正中接觸與咬出接觸對應咀嚼軌跡過程點在水平方向上的投影上分別為 A 點、O 點與 B 點，其中 OA 表示O點與 A 點間的距離，為 1.5 mm；OB表示O點與 B 點間的距離，為 0.5 mm；并定義A點、O 點與 B 點初始頜位高度分別為 ? 0.6、0、? 0.6 mm，如圖2 所示。統一設置垂直加載速率為 2 mm/s，仿真時間為 0.2 s，歷史輸出量為法向咬合力（以下簡稱咬合力）的三向分力。
圖2
靜態接觸與動態咀嚼工況設計
Figure2.
Cases design for static contact and dynamic chewing
本文選取多個靜態接觸工況間接表征動態咀嚼工況。為了盡可能表現咬合力的變化細節，新增 4 個咀嚼軌跡中間點，其水平投影分別C～F點，彼此之間的距離分別為：OC = 0.90 mm；OD = 0.30 mm；OE = 0.10 mm；OF = 0.30 mm，如圖2 所示，綠色帶箭頭實線表示磨牙咀嚼軌跡；紅色帶箭頭虛線表示中間點垂直軌跡；黑色帶箭頭實線表示起止點垂直軌跡。采用與文獻[19-20]相一致的咀嚼軌跡規劃原則，要求中間點（O 點，對應投影為C～F點）的z向侵入位移保持一致（取值為 0.004 mm），并定義起始點與終止點（對應投影為A點與B點）的z向侵入位移為 0 mm。
1.3 原始模型動力學仿真結果討論
以［? 1 000 N, 1 000 N］為 z 向咬合力區間，在此區間范圍內觀察靜態接觸工況的咬合力變化規律。刪除接觸間隙后，咬合力相對下頜磨牙垂直侵入位移的變化規律，如圖 3 所示。整體上來看，咬合力呈現近似線性增加，后期梯度逐漸變大。從咬合力分量上來看，相較于 y 和 z 向，x 向受力不明顯，即頰舌向受力明顯，近遠中向受力不明顯。特別地，正中接觸的 y 向分力小于 0，x 向分力大于 0，表明上頜磨牙受到頰側指向舌側、近中指向遠中的側向力，呈現典型的正中接觸咬合力特征[26]。
圖3
原始模型靜態接觸與動態咀嚼工況咬合力
Figure3.
Bite force of the original model in static contact and dynamic chewing cases
基于系列靜態接觸工況的咬合力數據，通過線性插值得到動態咀嚼過程中咬合力相對下頜磨牙水平滑移距離的變化規律，如圖3 所示。整體來看，呈現緩升—平穩—陡降—緩降的波動特征，沒有出現明顯咬合力沖擊，體現了咬合面形貌生物力學優越性。在咬入階段，咬合力呈現近似線性增加。觀察y向分力可知，在下頜水平滑移 0.6 mm 后，下頜磨牙舌斜面與上頜磨牙頰斜面逐漸開始接觸，實現了咬入接觸與正中接觸之間的漸進過渡。在正中接觸附近，咬合力保持穩定，表明圍繞正中接觸存在短距離的等剛度階段。在咬出階段，咬合力則在經歷速降與兩階段緩降后回歸零值，其速降原因在于下頜磨牙頰斜面僅在較短的滑移距離下就脫離受力。前述等剛度與速降的變化特征屬于典型的咀嚼過程中咬合力變化規律[19-20]。
2 等效模型幾何建模與受力分析
2.1 幾何建模
在曲面造型軟件 CREO 3.0（PTC Inc., 美國）中進行點接觸高副等效模型的 V 型面與球面造型設計，如圖4 所示。具體步驟如下：
圖4
等效模型幾何模型與受力分析
Figure4.
Geometric model and contact force analysis of the equivalent model
（1）以上下頜磨牙牙頸邊緣線拉伸創建等效模型的外周輪廓面。
（2）采用與原始模型相一致的靜坐標系，通過咬入斜面單位法向量 nα = ［xα, yα, zα］與咬出斜面單位法向量 nβ = ［xβ, yβ, zβ］，創建過原點的咬入接觸面α與咬出接觸面β。
（3）根據下頜磨牙牙尖形貌確定球體半徑為R，這里取 1.8 mm。圍繞下頜球心坐標 (xS, yS, zS) 使用旋轉命令創建磨牙球面。
（4）通過圓角命令創建各曲面間過渡面。
2.2 受力分析
2.2.1 咬合力定義
基于赫茲接觸模型定義接觸點Pi處咬合力與侵入位移的冪指數關系，動力學方程如式（1）所示：
 |
其中，Fi 為咬合力；di 為垂向位移（方向沿 z 軸正向）；ai 與 ei 分別為與接觸點相關的咬合剛度和冪指數，當 ei 取 1 或 1.5 時，表征剛性圓柱體或球體與彈性半空間的接觸；η 為接觸點 Pi 所在接觸面；nη 為接觸面 η 的單位法向量。值得注意的是，赫茲接觸模型是多體系統仿真軟件的常見內置接觸算法，如機械系統動力學仿真軟件 Adams（MSC Inc., 美國）中的“沖擊模式”算法，因此基于此公式所設計的等效模型可直接轉至仿真平臺使用。
2.2.2 接觸咬合力計算
對于咬入接觸狀態，球面與 V 型面僅存在單個接觸點P1（位于接觸面α）；對于正中接觸狀態，球面與 V 型面同時存在 2 個接觸點P3 與 P4；對于咬出接觸狀態，球面與 V 型面又恢復為僅存在單個接觸點 P2（位于接觸面 β）如圖 4 所示。根據式（1），i 分別取值 1～4，η對應取 α 與 β，此時咬合力表達式如式（2）～（4）所示：
 |
 |
 |
其中，Fin，Fco 和 Fex，分別為咬入接觸、正中咬合和咬出接觸的咬合力；din，dco 和 dex 分別為咬入接觸、正中咬合和咬出接觸的垂向位移。
3 基于靜態接觸工況的等效模型參數設計
3.1 期望咬合力擬合
使用最小二乘法對原始模型咬合力進行擬合，得到期望咬合剛度與冪指數，并選用均方根誤差（root mean square error, RMSE）為擬合質量評價指標，擬合結果如表1 所示（ain, x、ain, y 和 ain, z 分別為咬入接觸的 x 向、y 向和 z 向咬合力期望咬合剛度，aco, x、aco, y 和 aco, z 分別為正中咬合的 x 向、y 向和 z 向咬合力期望咬合剛度，aex, x、aex, y 和 aex, z 分別為咬出接觸的 x 向、y 向和 z 向咬合力期望咬合剛度，期望冪指數符號同理）。擬合結果表明，x 向冪指數高于 z 向，z 向冪指數高于 y 向，表明三個方向上非線性程度順序為 x>z>y。冪指數基本介于 1～1.5 之間，表明咬合接觸在各個方向上的接觸近似介于圓柱面與球面對彈性半空間的擠壓，體現出咬合接觸的復雜性。為了便于后續參數設計，同時使用固定冪指數 1 對各向分力進行擬合。
表1
基于赫茲接觸模型期望咬合力擬合結果
Table1.
Fitness of expected bite force based on Hertz contact model
3.2 斜面法向量與參數設計
3.2.1 法向量與咬合剛度設計
為簡化計算，這里采用統一的冪指數。根據表 1 中固定冪指數的 9 項期望咬合剛度，建立關于單位法向量與咬合剛度等 10 個未知數的非線性方程組，如式（5）所示：
 |
同時考慮單位向量的模為 1，聯合構成 11 個方程。增加剛度較小的 aex, x 作為補充未知數，保證方程組有解，這意味著方程組的解能夠實現靜態接觸工況 8 項咬合剛度與相應咬合力的復現。
3.2.2 下頜磨牙球心坐標優化
約束正中接觸時，球心距離咬入與咬出斜面的距離均等于 R，如式（6）所示。為使上下頜的動態接觸范圍盡可能達到最大，這里以 xS 與 yS 平方和取最小值為優化目標，使用中等規模算法（medium-scale algorithm）對球心坐標進行求解。
 |
3.2.3 上頜磨牙剛度分布
定義輔助點 Pref 為通過球心坐標 (x, y, z) 的軌跡平行線與斜面交線的交點 (xref, yref, zref)，如圖 5 所示，綠色有箭頭實線表示咀嚼軌跡；綠色虛線表示軌跡平行線；黑色實線表示斜面交線。假設咬入斜面 P1 點與 P3 點之間以及咬出斜面 P2 點與 P4 點之間剛度線性分布，沿斜線交線平行方向為均勻分布，得到斜面剛度與球心坐標之間的計算式，如式（7）所示：
圖5
等效模型咬合剛度定義和求解結果
Figure5.
Occlusal stiffness definition and distribution of the equivalent model
 |
其中，aα 為咬入斜面剛度；aβ 為咬出斜面剛度；zi 為接觸點 Pi 的 z 向坐標。
實例計算結果顯示，反向求解得到 aex, x 的設計值為 ? 0.226 × 103 N/mm，與原始剛度相比變小，但對整體咬合力影響不大。咬入斜面法向量［0.277, 0.486, 0.829］，單側咬入剛度 2.138×104 N/mm，正中時咬入剛度變為 2.606 × 104 N/mm；咬出斜面法向量［? 0.016, ? 0.656, 0.754］，單側咬出剛度 1.861 × 104 N/mm，正中時咬出剛度變為 5.341 × 104 N/mm。優化得到下頜球心坐標 (0.034, 0.139, ? 2.264)，達到最小距離 0.143 mm。
4 等效模型動力學仿真與評價
采用與前述靜態接觸與動態咀嚼工況一致的軌跡，對等效模型進行動力學分析。假設下頜球心沿某軌跡進行運動，首先對該軌跡進行離散，其次根據軌跡離散點坐標計算該軌跡點與各斜面的垂向距離，建立等效模型咬合力動力學模型，如式（8）所示：
 |
其中，dα 和 dβ 分別為球心軌跡點距離斜面 α 和 β 的距離；max(a, b) 表示取 a 與 b 兩者的較大值。
等效模型分析結果如圖6 所示。前述可知，等效模型的參數設計是基于靜態接觸工況的 8 項期望咬合力數據，因此等效模型很好地實現了相應咬合力的復現，而具體吻合程度是由擬合質量決定。
圖6
等效模型靜態接觸與動態咀嚼工況咬合力
Figure6.
Bite force of the equivalent model in static contact and dynamic chewing cases
對于動態咀嚼工況而言，咬合力呈現緩升—陡升—陡降—緩降的波動特征。等效模型在滑移前期（0～0.6 mm）與后期（1.8～2.0 mm）咬合力與原始模型吻合良好。然而，滑移中期由咬入接觸進入正中接觸、正中接觸進入咬出接觸時，咬合力引入一定程度的沖擊。原因在于 V 型面的平面形貌，使得單雙側接觸過渡距離較短。可通過定義合理的滑移位移—咀嚼時間曲線，即延長過渡時間，實現對于沖擊的弱化。
5 結論
本文為了實現多體系統動力學模型對于咬合接觸的直接模擬，基于靜態接觸工況提出一種由 V 型面與球面構成的點接觸高副，最終得到以下結論：原始模型靜態接觸工況咬合力數據呈現近似線性增加；動態咀嚼工程咬合力呈現緩升—平穩—陡降—緩降的波動特征；等效模型能夠實現 8 項靜態咬合力的復現；等效模型在動態咀嚼工況呈現緩升—陡升—陡降—緩降的波動特征，一定程度上實現了動態咬合力變化趨勢的復現。整體來看，等效模型利用相對簡單的形貌盡可能復現了原始模型的動力學特征。
若需要對進一步緩解或消除等效模型在咀嚼工況中期的動力學沖擊，可考慮使用弧形曲面代替 V 型面，這將在接下來的研究中進一步完善。此外，磨牙的點接觸高副等效方法可為其他人體關節的等效建模提供參考。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
研究人體口頜系統的肌力模式、關節受力與咬合性能問題，離不開準確的多體系統動力學模型。同樣地，研究咀嚼機器人系統（仿生口頜系統）的咀嚼運動與咬合力控制問題亦然。近年來，牙頜模型作為口腔代表性部件正越來越多地作為子結構方式引入到多體系統動力學建模中，以進一步提高系統整體動力學模型質量，而不是僅僅依靠咬合力數值的方式[1-2]。目前，可采用的動力學建模策略包括多體系統仿真、有限元仿真以及聯合仿真。其中，有限元理論通過設置合理的網格密度，可解決原始咬合接觸的難收斂問題，最大限度地還原咬合面與食品形貌，模擬咬合面碰撞與食品變形行為，不足之處在于增加了上下頜、咀嚼肌肉以及顳下頜關節等非必要組件的網格建模工作量和后續動力學模型的解算難度[3-6]。為了盡可能地降低多體系統建模工作量與解算難度，同時也為了保證牙頜模型建模精度，廣泛采用的是多體系統仿真與聯合仿真的方案[7]。
多體系統仿真是指在多體系統仿真平臺下對牙頜模型進行建模，僅采取分析力學解算，具有建模簡單與求解速度快的優勢。可執行仿真的軟件平臺較多，例如：機械系統動力學仿真軟件 Adams（MSC Inc., 美國）、下頜運動評估軟件 Kinematics 3D（Zebris Medical Inc., 德國）、人體肌肉骨骼仿真軟件 OpenSim（The National Center for Simulation in Rehabilitation Research, 美國）和人體運動仿真軟件 Anybody Modeling System（AnyBody Technology Inc., 丹麥）[7-9]。但是，多體系統仿真針對接觸區域僅進行有限度的離散化，當面對磨牙所存在的尖、窩、脊等復雜接觸面形貌特征時，接觸判斷準則將難以正確判斷接觸是否發生，繼而發生接觸穿透現象，并致使模型求解失敗[10]。少部分研究人員對磨牙咬合動力學模型進行了簡化等效，例如：2018 年，王貴飛[11]在機械系統動力學仿真軟件 Adams（MSC Inc., 美國）內使用簡化的平面—球面—平面型式構造兩個點接觸高副（higher kinematic pair, HKP）描述上頜磨牙—食品—下頜磨牙三體接觸關系。然而，該點接觸高副并未體現出磨牙正中接觸的三向受力規律，也無法適用于動態咀嚼軌跡的加載，使得咀嚼機器人咬合力控制并不全面。針對復雜接觸問題，報道的另一種建模方案是將復雜接觸面進一步離散化，例如：Tuijt 等[12]與 Skipper 等[13]進行口頜多體系統建模時，將顳下頜關節接觸面離散為致密點云，用作后續接觸分析，實現在單一仿真環境下復雜接觸的多體動力學分析。可見，盡管多體系統仿真方案支持平臺較多，但存在牙頜模型建模精度或離散化程度不足的問題。
聯合仿真是指在多體系統仿真基礎上按需對牙頜模型進行有限元仿真，同時采用分析力學與彈性力學解算。醫用多體系統仿真軟件 ArtiSynth（University of British Columbia, 加拿大）能夠實現剛體與軟組織的聯合動力學仿真，即剛體動力學分析與有限元分析的聯合，現已應用于開閉口運動下肌力模式分析與下頜旋轉中心分析[14-16]。2008 年，Hannam 等[17]在運用該軟件對口頜骨—肌系統進行動力學分析時，預先根據仿真食品被搗碎過程中力學特性設計咬合力，最終口頜系統仿真結果與臨床數據吻合。同樣地，該仿真策略也已成功應用于顳下頜關節等類似復雜接觸問題。2019 年，Sagl 等[18]提出一種綜合動態咀嚼模型，由肌肉驅動的下頜剛體、顳下頜關節盤有限元模型與彈性基軟骨模型構成，結果顯示這種聯合模型能夠在保證仿真效果的同時，有效降低計算時間。然而，能夠實現聯合仿真建模的類似軟件平臺并不多。另一方面，聯合仿真中牙頜模型有限元分析部分更多是獨立地在有限元軟件平臺上實施。Benazzi 等[19]與 Qin 等[20]分別使用有限元軟件 Abaqus（SIMULIA Inc., 美國）對口頜系統中牙頜模型的咬合接觸問題進行建模，輸出咀嚼軌跡與咬合力之間的數學關系。Pascale 等[21]使用 C 語言結合可視化軟件 Vega（MultiGen Paradigm Inc., 美國）庫函數編寫牙頜有限元模型，同時采用分區的碰撞檢測算法計算咬合接觸受力。后續，該團隊使用擴展型拉普拉斯光滑算法對咬合面形貌進行簡化，進一步實現簡化的咬合接觸應力分析[22]，但并未給出將以上獨立的仿真數據應用（聯合）至多體系統模型中的方法。
基于以上研究，現有多體系統仿真平臺除非在分析力學基礎上增強模型離散能力，或補充彈性力學算法轉為聯合仿真平臺，否則將無法實現牙頜模型子結構模型的準確建模。為了解決這一問題，本文提出一種新的仿真策略，即設計一種適應現有多體系統仿真平臺離散程度的簡化牙頜模型，該模型由 V 型面與球面構成點接觸高副，使其分析力學結果與原始模型的彈性力學結果能夠保持等效。本研究期望可為多體系統的牙頜模型提供新的動力學建模策略，實現多體系統的動力學建模精度的進一步提升，另一方面通過本文的嘗試，也希望可以為人體其他復雜接觸的動力學建模提供新思路。
1 原始模型咬合動力學分析
1.1 三維與有限元模型的建立
第一磨牙（以下簡稱磨牙）作為牙頜模型中承受極大咬合力的牙位而常被用于各類咬合問題的研究，本文選取咀嚼機器人系統右側磨牙作為研究對象，進行咬合動力學簡化建模與評價[23]。面向義齒斷裂或磨損測試的仿生型咀嚼機器人系統一般選用標準牙頜模型作為基牙安裝于仿生口腔環境內。使用光柵掃描設備 AutoScan DS100（杭州先臨三維科技股份有限公司，中國）對日進牙科教學模型（Nissin Education Products，中國）的磨牙部件進行掃描，掃描精度 0.015 mm。其次，進行咬合面三維重建，并根據預備工藝進行預備體曲面造型，最終得到基牙與義齒的三維模型，如圖1 所示。定義靜坐標系{G}以上頜磨牙中央窩為原點，以磨牙后方、側方與上方為x、y與z軸正向。
圖1
咀嚼機器人磨牙三維與有限元模型
Figure1.
Geometrical and finite element models of molars assembled in chewing robot
使用網格劃分軟件 HyperMesh 14.0（Altair Engineering Inc., 美國）對磨牙全瓷冠以 0.1 mm 為全局單元尺寸進行網格劃分，形成網格模型。設置基牙材料為復合樹脂（密度 1 200 kg/m3，彈性模量 11.5 GPa，泊松比 0.24）；金屬冠材料為鈷鉻合金（密度 8 000 kg/m3，彈性模量 220 GPa，泊松比 0.30）；義齒類型為全瓷冠（密度 6 100 kg/m3，彈性模量 220 GPa，泊松比 0.30）[24-25]。將其導入有限元軟件 Abaqus 11.0（SIMULIA Inc., 美國）中，設置咬合面接觸類型為面—面接觸類型，解算機制為“運動學模式”。鑒于唾液的潤滑作用，本文不考慮切向咬合力。
1.2 靜態接觸與動態咀嚼工況設計
定義下頜磨牙動坐標系{L}隨下頜磨牙運動而運動，在正中接觸時，與靜坐標系{G}重合。下頜運動軌跡使用動坐標系在靜坐標系中的位姿進行描述。設計 3 種典型接觸工況，分別表征咬入接觸（intrusive contact）、正中接觸（centric occlusion）與咬出接觸（extrusive contact）。完整的靜態接觸工況定義應包含該接觸狀態對應咀嚼軌跡（僅咬合相）過程點在水平投影上的位置、初始頜位高度、垂直加載速率與仿真時間等數據。假設咀嚼軌跡水平投影滿足咬入角θ1 與咬出角θ2 為 10° 的角度特征，咬入接觸、正中接觸與咬出接觸對應咀嚼軌跡過程點在水平方向上的投影上分別為 A 點、O 點與 B 點，其中 OA 表示O點與 A 點間的距離，為 1.5 mm；OB表示O點與 B 點間的距離，為 0.5 mm；并定義A點、O 點與 B 點初始頜位高度分別為 ? 0.6、0、? 0.6 mm，如圖2 所示。統一設置垂直加載速率為 2 mm/s，仿真時間為 0.2 s，歷史輸出量為法向咬合力（以下簡稱咬合力）的三向分力。
圖2
靜態接觸與動態咀嚼工況設計
Figure2.
Cases design for static contact and dynamic chewing
本文選取多個靜態接觸工況間接表征動態咀嚼工況。為了盡可能表現咬合力的變化細節，新增 4 個咀嚼軌跡中間點，其水平投影分別C～F點，彼此之間的距離分別為：OC = 0.90 mm；OD = 0.30 mm；OE = 0.10 mm；OF = 0.30 mm，如圖2 所示，綠色帶箭頭實線表示磨牙咀嚼軌跡；紅色帶箭頭虛線表示中間點垂直軌跡；黑色帶箭頭實線表示起止點垂直軌跡。采用與文獻[19-20]相一致的咀嚼軌跡規劃原則，要求中間點（O 點，對應投影為C～F點）的z向侵入位移保持一致（取值為 0.004 mm），并定義起始點與終止點（對應投影為A點與B點）的z向侵入位移為 0 mm。
1.3 原始模型動力學仿真結果討論
以［? 1 000 N, 1 000 N］為 z 向咬合力區間，在此區間范圍內觀察靜態接觸工況的咬合力變化規律。刪除接觸間隙后，咬合力相對下頜磨牙垂直侵入位移的變化規律，如圖 3 所示。整體上來看，咬合力呈現近似線性增加，后期梯度逐漸變大。從咬合力分量上來看，相較于 y 和 z 向，x 向受力不明顯，即頰舌向受力明顯，近遠中向受力不明顯。特別地，正中接觸的 y 向分力小于 0，x 向分力大于 0，表明上頜磨牙受到頰側指向舌側、近中指向遠中的側向力，呈現典型的正中接觸咬合力特征[26]。
圖3
原始模型靜態接觸與動態咀嚼工況咬合力
Figure3.
Bite force of the original model in static contact and dynamic chewing cases
基于系列靜態接觸工況的咬合力數據，通過線性插值得到動態咀嚼過程中咬合力相對下頜磨牙水平滑移距離的變化規律，如圖3 所示。整體來看，呈現緩升—平穩—陡降—緩降的波動特征，沒有出現明顯咬合力沖擊，體現了咬合面形貌生物力學優越性。在咬入階段，咬合力呈現近似線性增加。觀察y向分力可知，在下頜水平滑移 0.6 mm 后，下頜磨牙舌斜面與上頜磨牙頰斜面逐漸開始接觸，實現了咬入接觸與正中接觸之間的漸進過渡。在正中接觸附近，咬合力保持穩定，表明圍繞正中接觸存在短距離的等剛度階段。在咬出階段，咬合力則在經歷速降與兩階段緩降后回歸零值，其速降原因在于下頜磨牙頰斜面僅在較短的滑移距離下就脫離受力。前述等剛度與速降的變化特征屬于典型的咀嚼過程中咬合力變化規律[19-20]。
2 等效模型幾何建模與受力分析
2.1 幾何建模
在曲面造型軟件 CREO 3.0（PTC Inc., 美國）中進行點接觸高副等效模型的 V 型面與球面造型設計，如圖4 所示。具體步驟如下：
圖4
等效模型幾何模型與受力分析
Figure4.
Geometric model and contact force analysis of the equivalent model
（1）以上下頜磨牙牙頸邊緣線拉伸創建等效模型的外周輪廓面。
（2）采用與原始模型相一致的靜坐標系，通過咬入斜面單位法向量 nα = ［xα, yα, zα］與咬出斜面單位法向量 nβ = ［xβ, yβ, zβ］，創建過原點的咬入接觸面α與咬出接觸面β。
（3）根據下頜磨牙牙尖形貌確定球體半徑為R，這里取 1.8 mm。圍繞下頜球心坐標 (xS, yS, zS) 使用旋轉命令創建磨牙球面。
（4）通過圓角命令創建各曲面間過渡面。
2.2 受力分析
2.2.1 咬合力定義
基于赫茲接觸模型定義接觸點Pi處咬合力與侵入位移的冪指數關系，動力學方程如式（1）所示：
|
其中，Fi 為咬合力；di 為垂向位移（方向沿 z 軸正向）；ai 與 ei 分別為與接觸點相關的咬合剛度和冪指數，當 ei 取 1 或 1.5 時，表征剛性圓柱體或球體與彈性半空間的接觸；η 為接觸點 Pi 所在接觸面；nη 為接觸面 η 的單位法向量。值得注意的是，赫茲接觸模型是多體系統仿真軟件的常見內置接觸算法，如機械系統動力學仿真軟件 Adams（MSC Inc., 美國）中的“沖擊模式”算法，因此基于此公式所設計的等效模型可直接轉至仿真平臺使用。
2.2.2 接觸咬合力計算
對于咬入接觸狀態，球面與 V 型面僅存在單個接觸點P1（位于接觸面α）；對于正中接觸狀態，球面與 V 型面同時存在 2 個接觸點P3 與 P4；對于咬出接觸狀態，球面與 V 型面又恢復為僅存在單個接觸點 P2（位于接觸面 β）如圖 4 所示。根據式（1），i 分別取值 1～4，η對應取 α 與 β，此時咬合力表達式如式（2）～（4）所示：
|
|
|
其中，Fin，Fco 和 Fex，分別為咬入接觸、正中咬合和咬出接觸的咬合力；din，dco 和 dex 分別為咬入接觸、正中咬合和咬出接觸的垂向位移。
3 基于靜態接觸工況的等效模型參數設計
3.1 期望咬合力擬合
使用最小二乘法對原始模型咬合力進行擬合，得到期望咬合剛度與冪指數，并選用均方根誤差（root mean square error, RMSE）為擬合質量評價指標，擬合結果如表1 所示（ain, x、ain, y 和 ain, z 分別為咬入接觸的 x 向、y 向和 z 向咬合力期望咬合剛度，aco, x、aco, y 和 aco, z 分別為正中咬合的 x 向、y 向和 z 向咬合力期望咬合剛度，aex, x、aex, y 和 aex, z 分別為咬出接觸的 x 向、y 向和 z 向咬合力期望咬合剛度，期望冪指數符號同理）。擬合結果表明，x 向冪指數高于 z 向，z 向冪指數高于 y 向，表明三個方向上非線性程度順序為 x>z>y。冪指數基本介于 1～1.5 之間，表明咬合接觸在各個方向上的接觸近似介于圓柱面與球面對彈性半空間的擠壓，體現出咬合接觸的復雜性。為了便于后續參數設計，同時使用固定冪指數 1 對各向分力進行擬合。
表1
基于赫茲接觸模型期望咬合力擬合結果
Table1.
Fitness of expected bite force based on Hertz contact model
3.2 斜面法向量與參數設計
3.2.1 法向量與咬合剛度設計
為簡化計算，這里采用統一的冪指數。根據表 1 中固定冪指數的 9 項期望咬合剛度，建立關于單位法向量與咬合剛度等 10 個未知數的非線性方程組，如式（5）所示：
|
同時考慮單位向量的模為 1，聯合構成 11 個方程。增加剛度較小的 aex, x 作為補充未知數，保證方程組有解，這意味著方程組的解能夠實現靜態接觸工況 8 項咬合剛度與相應咬合力的復現。
3.2.2 下頜磨牙球心坐標優化
約束正中接觸時，球心距離咬入與咬出斜面的距離均等于 R，如式（6）所示。為使上下頜的動態接觸范圍盡可能達到最大，這里以 xS 與 yS 平方和取最小值為優化目標，使用中等規模算法（medium-scale algorithm）對球心坐標進行求解。
|
3.2.3 上頜磨牙剛度分布
定義輔助點 Pref 為通過球心坐標 (x, y, z) 的軌跡平行線與斜面交線的交點 (xref, yref, zref)，如圖 5 所示，綠色有箭頭實線表示咀嚼軌跡；綠色虛線表示軌跡平行線；黑色實線表示斜面交線。假設咬入斜面 P1 點與 P3 點之間以及咬出斜面 P2 點與 P4 點之間剛度線性分布，沿斜線交線平行方向為均勻分布，得到斜面剛度與球心坐標之間的計算式，如式（7）所示：
圖5
等效模型咬合剛度定義和求解結果
Figure5.
Occlusal stiffness definition and distribution of the equivalent model
|
其中，aα 為咬入斜面剛度；aβ 為咬出斜面剛度；zi 為接觸點 Pi 的 z 向坐標。
實例計算結果顯示，反向求解得到 aex, x 的設計值為 ? 0.226 × 103 N/mm，與原始剛度相比變小，但對整體咬合力影響不大。咬入斜面法向量［0.277, 0.486, 0.829］，單側咬入剛度 2.138×104 N/mm，正中時咬入剛度變為 2.606 × 104 N/mm；咬出斜面法向量［? 0.016, ? 0.656, 0.754］，單側咬出剛度 1.861 × 104 N/mm，正中時咬出剛度變為 5.341 × 104 N/mm。優化得到下頜球心坐標 (0.034, 0.139, ? 2.264)，達到最小距離 0.143 mm。
4 等效模型動力學仿真與評價
采用與前述靜態接觸與動態咀嚼工況一致的軌跡，對等效模型進行動力學分析。假設下頜球心沿某軌跡進行運動，首先對該軌跡進行離散，其次根據軌跡離散點坐標計算該軌跡點與各斜面的垂向距離，建立等效模型咬合力動力學模型，如式（8）所示：
|
其中，dα 和 dβ 分別為球心軌跡點距離斜面 α 和 β 的距離；max(a, b) 表示取 a 與 b 兩者的較大值。
等效模型分析結果如圖6 所示。前述可知，等效模型的參數設計是基于靜態接觸工況的 8 項期望咬合力數據，因此等效模型很好地實現了相應咬合力的復現，而具體吻合程度是由擬合質量決定。
圖6
等效模型靜態接觸與動態咀嚼工況咬合力
Figure6.
Bite force of the equivalent model in static contact and dynamic chewing cases
對于動態咀嚼工況而言，咬合力呈現緩升—陡升—陡降—緩降的波動特征。等效模型在滑移前期（0～0.6 mm）與后期（1.8～2.0 mm）咬合力與原始模型吻合良好。然而，滑移中期由咬入接觸進入正中接觸、正中接觸進入咬出接觸時，咬合力引入一定程度的沖擊。原因在于 V 型面的平面形貌，使得單雙側接觸過渡距離較短。可通過定義合理的滑移位移—咀嚼時間曲線，即延長過渡時間，實現對于沖擊的弱化。
5 結論
本文為了實現多體系統動力學模型對于咬合接觸的直接模擬，基于靜態接觸工況提出一種由 V 型面與球面構成的點接觸高副，最終得到以下結論：原始模型靜態接觸工況咬合力數據呈現近似線性增加；動態咀嚼工程咬合力呈現緩升—平穩—陡降—緩降的波動特征；等效模型能夠實現 8 項靜態咬合力的復現；等效模型在動態咀嚼工況呈現緩升—陡升—陡降—緩降的波動特征，一定程度上實現了動態咬合力變化趨勢的復現。整體來看，等效模型利用相對簡單的形貌盡可能復現了原始模型的動力學特征。
若需要對進一步緩解或消除等效模型在咀嚼工況中期的動力學沖擊，可考慮使用弧形曲面代替 V 型面，這將在接下來的研究中進一步完善。此外，磨牙的點接觸高副等效方法可為其他人體關節的等效建模提供參考。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
