咀嚼機器人在義齒材料測試和下頜康復訓練領域具有廣闊的應用前景，而顳下頜關節的機構型式是影響咀嚼機器人性能表現的重要因素。鑒于目前彈性元件已在仿生機器人領域得到廣泛應用，因此本文在點接觸高副的基礎上采用彈性元件模擬顳下頜關節內關節盤緩沖振蕩的生物力學特點，形成彈性顳下頜關節機構型式；然后，本文探討了彈性顳下頜關節對咀嚼機器人在自由度、運動學與動力學等方面的影響，根據關節面幾何約束對顳下頜關節的位置與速度進行運動學分析，并進行基于拉格朗日方程的動力學分析；最后，使用響應面法對彈性元件預載荷以及剛度取值的影響進行分析。本文結果表明，彈性顳下頜關節能夠有效保證關節的靈活運動與穩定受力。綜上，本文所提出的彈性顳下頜關節機構方案進一步提升了咀嚼機器人仿生性，也為粘彈性關節盤的仿生設計提供了新思路。
引用本文： 秦文龍, 叢明, 任翔, 溫海營, 劉冬. 仿生咀嚼機器人彈性顳下頜關節設計與性能分析. 生物醫學工程學雜志, 2020, 37(3): 512-518, 526. doi: 10.7507/1001-5515.201812051 復制
引言
近年來，具有多自由度的機器人正越來越多地出現在醫療領域，如外科手術機器人、腕功能康復機器人和牙科機器人等[1-3]。其中，咀嚼機器人作為一種可以模擬人類咀嚼軌跡與咬合力的機械裝置受到牙科學研究人員的廣泛關注[4]。通過增設或改進不同的機器人專用模塊的方式，咀嚼機器人可以應用到牙科義齒測試、下頜康復訓練和學科教學等多種不同的場景[5-7]。
為保證咀嚼機器人運動與受力的仿生性，咀嚼機器人的設計通常是參照口頜內的骨—肌系統展開[8]。值得注意的是，口頜系統中顳下頜關節呈現旋轉與平移的復合運動特點，發揮著穩定下頜運動與受力的重要作用[9]。為了使得咀嚼機器人能夠復現這一生物學特點，如何構建符合人體顳下頜關節的仿生設計成為研究熱點[8]。早在咀嚼機器人引入仿生顳下頜關節之前，牙科學中用于模擬下頜運動的咬合架已經將顳下頜關節以鉸鏈型式、上窩下球型式和上球下窩型式加入到機構設計中，但均無法實現仿生髁突的復合運動[10]。吉林大學孫鐘雷等[11]在設計用于食品物性參數檢測的仿生咀嚼裝置（僅實現開閉口運動）時將下頜支、髁突與關節窩分別簡化為圓柱、球體與活動球槽，但是并未對顳下頜關節作進一步的機構學分析。新西蘭大學 Cheng 等[12-13]和 Mostashiri 等[14]在并聯咀嚼機器人（拓撲結構符號[15]：6—SOC{—R—S—S—}）的基礎上同樣將髁突簡化為球體，但將顳骨關節面簡化為貫通式滑道，進行自由度計算與運動學分析。在該研究中，滑道上曲面或下曲面變為重復限制運動的虛約束，盡管能在不影響機構運動的前提下增加機構的剛性，改善其受力狀況，但是由于虛約束對機構的幾何條件要求較高，對機構的加工和裝配精度亦提出了較高的要求。大連理工大學溫海營等[16-18]則針對并聯咀嚼機器人（拓撲結構符號[15]：6—SOC{—P—U—S—}）提出一種仿生性更強的點接觸高副型式（high kinematic pair，HKP）顳下頜關節，該關節除球形髁突外，特別提取人體顳骨關節面解剖生理結構尺寸用于設計仿生顳骨關節面，充分利用了口頜系統的冗余特性，體現出仿生的運動學特征，并且實現了驅動器能量的高效利用，但是單面剛性約束導致髁突易于脫離顳骨關節面，造成關節受力的非規律性波動。
目前的仿生策略忽略彈性關節盤，僅使用點接觸高副型式對剛性顳骨關節面與髁突進行接觸面形貌仿生，存在關節自由度嚴格受限與受力波動的問題。本文創新性將髁桿與下頜分離為獨立構件，并且使用彈性元件模擬關節盤，形成彈性顳下頜關節；進一步，對具有新關節的咀嚼機器人進行自由度、運動學與動力學分析，并對彈性元件的預載荷與剛度參數的影響規律進行討論。通過上述研究，期望可以為咀嚼機器人開發人員提供新的顳下頜關節機構型式，實現咀嚼機器人仿生性的進一步提升，另一方面也可以為人體其他位置關節盤的仿生設計提供新思路。
1 引入彈性顳下頜關節的咀嚼機器人機構
1.1 顳下頜關節的仿生設計
顳下頜關節由顳骨關節面、關節盤、髁突及其他軟組織構成。其中，關節盤是人體顳下頜關節的重要組成部分，具有緩沖振蕩、分散載荷、穩定關節的生物力學作用[19-20]。由于粘彈性關節盤的體積較小并且具有各項異性的力學特點，若選擇直接物理重建的仿生方式將對三維立體打印技術和粘彈性材料制備技術提出較高的要求。因此，本文僅對關節盤緩沖振蕩的生物力學功能進行間接的仿生。
首先，根據生理髁突的形貌數據，設計末端帶有直徑為 5 mm 球體的桿件模擬髁突，稱為髁桿（其中，球體稱為髁球）。同時，使用分段三次樣條曲線對生理髁突標記點 O 的矢狀面運動曲線進行插值，作為髁球中心運動軌跡 z = f（x）。將髁球中心 E 運動軌跡沿髁球半徑在矢狀面內偏置，得到仿生顳骨關節面（假設顳骨關節面任意法向與矢狀面始終保持平行）。在髁桿與仿生下頜之間使用壓縮彈簧模擬關節盤，同時在仿生下頜中設置滑道約束髁桿僅作相對平移運動。GB/T 2089-2009《普通圓柱螺旋壓縮彈簧尺寸及參數（兩端圈并緊磨平或制扁）》[21]中提及了圓柱螺旋壓縮彈簧標準，本文根據其中不同彈簧直徑下的彈簧剛度與最大工作載荷數據表，可確定設計變量的變化范圍，取預載荷 Fpre 變化范圍為 10～100 N，初值為 55 N；彈簧剛度 kspring 變化范圍為 10～50 N/mm，初值為 30 N/mm。基于以上參數，由顳骨關節面、髁桿與彈性元件共同構成彈性的顳下頜關節結構，并形成兩種運動副：點接觸高副、移動副。在咀嚼運動時，通過彈性元件伸縮（彈性元件下端點 A 與髁球中心 E 下移，髁球中心 E 與生理髁突標記點 O 不再重合），保持髁球與顳骨關節面不發生脫離，同時避免單純點接觸高副帶來的接觸面碰撞問題，如圖 1 所示。
圖1
彈性顳下頜關節仿生結構與工作原理
Figure1.
Biomimetic structure of elastic temporomandibular joint and its working principle
1.2 咀嚼機器人機構分析
在彈性顳下頜關節基礎上，采用 6 條驅動支鏈模擬咀嚼肌群，包括：左側咬肌、右側咬肌、左側顳肌、右側顳肌、左側翼外肌和右側翼外肌，形成新型咀嚼機器人機構型式，其三維模型如圖 2 所示，機構型式如圖 3 所示。根據庫茲貝奇準則，計算咀嚼機器人機構自由度，如式（1）所示：
圖2
咀嚼機器人三維模型
Figure2.
3D model of the masticatory robot
圖3
咀嚼機器人機構型式
Figure3.
Mechanism of the masticatory robot
 |
其中，m 是自由度數目；n 是構件數目，共 16 個，分別為 6 個連桿、6 個滑塊、2 個髁桿、1 個下頜平臺和 1 個靜平臺；ji 是具有 i 個自由度的關節數目，包括：移動副 j1 = 8（髁桿與下頜之間 2 個移動副、滑塊與絲杠之間 6 個移動副）、萬向節 j2 = 6、球副 j3 = 6、4 自由度運動副 j4 = 0、點接觸高副 j5 = 2。由此，可計算得到機器人的自由度是 6。尤其是根據并聯彈性驅動器的定義，采用 2 條仿顳下頜關節彈性支鏈與 6 條仿肌肉驅動支鏈，共同構成并聯彈性驅動器[22]。
2 彈性顳下頜關節運動學分析
2.1 坐標系與咀嚼軌跡定義
為了便于描述接下來的運動學與動力學方程，首先建立咀嚼機器人動靜坐標系。左側髁突基坐標系為{G}，其原點 OG 位于左側生理髁突標記點，x、y 和 z 向分別為咀嚼機器人前方、左側方和上方，下文簡稱基坐標系。下頜連體坐標系為{L}，其原點 O 與下頜固結，隨下頜運動。在閉口時，下頜連體坐標系與基坐標系重合。
咀嚼時序為從正中頜位起始，通過右側咀嚼運動，回歸正中頜位。綜合考慮切牙“淚滴狀”的軌跡特征、非工作側髁突在矢狀面與水平面的投影特征以及工作側髁突的錐形面位移約束特征，進行右側咀嚼軌跡的定義，如圖 4 所示。
圖4
切牙咀嚼軌跡（相對于基坐標系{G}）
Figure4.
Chewing trajectory of the incisor (in the ground coordinate system {G})
2.2 運動學分析
2.2.1 位置分析
為了便于運動學分析，使用下頜連體坐標系相對于基坐標系的 6 個位姿參數進行下頜運動描述，即：x、y、z、α、β 和 γ。其中，α、β 和 γ 是歐拉角，依次為繞基坐標系 x 軸翻滾、y 軸俯仰和 z 軸偏航，x、y 和 z 是下頜連體坐標系原點 O 的平移量。
下頜任意一點 P 在基坐標系{G}的位置向量可以通過如下坐標變換矩陣求得，如式（2）所示：
 |
其中，SP 為任意一點 P 在基坐標系{G}中的坐標（列向量）；SO 為下頜連體坐標系{L}中原點 O 在基坐標系{G}中的坐標（列向量）；LSOSP 是在下頜連體坐標系{L}中任意一點 P 相對于原點 O 的位置向量（列向量）；GRL 為下頜連體坐標系{L}在基坐標系{G}中的方位，如式（3）所示：
 |
其中，c 代表 cos（·），s 代表 sin（·）。
以左側髁桿為例進行位置分析。彈簧下端點 A 在基坐標系下的坐標如式（4）所示：
 |
推導出髁桿方向向量 n，如式（5）所示：
 |
髁球中心 E 具有兩個幾何約束：球心始終處于髁球中心運動軌跡曲線上，滿足 z = f（x）；球心受髁桿方向向量約束。據此聯立方程組，求解髁球中心坐標，如式（6）所示。
 |
其中，e1 與 e3 分別為 x 向與 z 向單位向量（行向量）。根據 E 與 A 的坐標可換算出左側彈簧長度變化量 H，同理可求出右側彈簧長度變化量，如圖 5 所示。
圖5
彈簧變形量及變形速度
Figure5.
Deformation and changing rate of the spring
2.2.2 速度分析
同樣地，以左側髁桿為例進行速度分析。已知髁球中心的瞬時速度矢量與髁球中心軌跡平面相切。將待求的髁球中心速度 VE 沿切平面進行矢量分解，并代入由切平面產生的運動速度約束，如式（7）所示：
 '/> |
進一步地，髁球中心在不同坐標系中的相對速度關系如式（8）所示：
 |
其中，v 為彈簧變形速度；LVE 為髁球中心在下頜連體坐標系中的速度分量，可根據不同坐標系下速度變換關系進行求解，如式（9）所示：
 |
其中，ωmandi 為仿生下頜在基坐標系下角速度。聯立式（7）～（9），推導得到左側彈簧變形速度 v，同理可求出右側彈簧變形速度，如圖 5 所示。
結果表明，工作側顳下頜關節彈性元件的變化量未超過 1.7 mm，近似呈現正弦波動。非工作側變化量超過工作側髁突變化量，波動程度較小。相比于單純的點接觸高副型式顳下頜關節，彈性元件的引入允許髁突在顳下頜關節面下實現更靈活的空間運動，而不發生干涉。該變化量亦與人體關節盤較小的結構尺寸相適應。最大變形速度發生在非工作側髁突的快速閉口階段，接近 40 mm/s。較大的變化速率體現出人體關節盤具備優良的緩沖髁突速度沖擊的能力。同樣地，非工作側髁突變化速率超過工作側髁突，且兩者均呈現波動。
3 彈性顳下頜關節動力學分析
為了進一步研究顳下頜關節對咀嚼機器人受力的影響，本文采用第二類拉格朗日方程建立動力學方程，如式（10）所示：
 |
其中，U 為拉格朗日函數，定義為動能之和 T 減去勢能之和 W；qk 為咀嚼機器人第 k 項廣義坐標（即下頜 6 個位姿參數）；Qk 為咀嚼機器人非保守力在第 k 項廣義坐標下的廣義力之和。
3.1 動能和勢能的計算
3.1.1 仿生下頜動能和勢能
仿生下頜的動能 Kmandi 如式（11）所示：
 |
其中，mmandi 為仿生下頜的質量；Vmandi 為仿生下頜質心在基坐標系下的線速度；Imandi 為仿生下頜在基坐標系中關于質心的慣性矩陣；ωmandi 為仿生下頜在基坐標系中的角速度。
仿生下頜的勢能 Pmandi 如式（12）所示：
 |
其中，Zmandi 為仿生下頜質心在基坐標系中 z 方向的坐標；g 為重力加速度。
3.1.2 彈性元件勢能
左側彈性元件的勢能 Pspring 如式（13）所示：
 |
同理，計算出右側彈性元件的勢能。
3.1.3 6 條驅動支鏈動能和勢能
第 i 條驅動支鏈的動能 Kchain，i 如式（14）所示：
 |
其中，mslid，i 為滑塊的質量；Vslid，i 為滑塊的速度；Vlink，i 為連桿質心在基坐標系下的線速度；mlink，i 為連桿的質量；Ilink，i 為連桿在基坐標系下對質心的慣量矩陣；ωlink，i 為仿生下頜在基坐標系中的角速度。
第 i 條驅動支鏈的勢能 Pchain，i 如式（15）所示：
 |
其中，Zslid，i 與 Zlink，i 為滑塊和連桿質心在基坐標系中 z 方向的坐標。
3.2 廣義力的計算
咀嚼機器人力學體系中非保守力包括 6 組驅動力 Fslid，i 與咬合力 Fmolar，對應的是廣義力為 Qslid，k 與 Qmolar，k；保守力包括仿生下頜所受到的彈簧力，可以根據彈簧預載荷與變形量直接算出。
3.2.1 驅動力廣義力
根據驅動力 τi 與移動副位移為 Sslid，i，換算出 6 組驅動力對應的廣義力，如式（16）所示：
 |
3.2.2 咬合力廣義力
牙科模型咬合力 Fmolar 作用在右側第一磨牙上，咬合力在正中頜時達到最大值 150 N，在咬合相階段呈現正弦變化規律，方向保持垂直向上；作用點 z 軸位移為 Zmolar，換算出咬合力對應的廣義力，如式（17）所示：
 |
將式（11）～（17）代入拉格朗日方程式（10），得到完整的咀嚼機器人動力學模型，并解算得到驅動力 Fslid，i 與驅動功率 Pslid，i，如圖 6 所示。進一步分析能量利用率，如式（18）所示：
圖6
電機驅動力與驅動功率
Figure6.
Driving force and driving power
 |
其中，η 為能量利用率；Wall 為消耗的總功；Wpos 為正功；Wneg 為負功。
結果表明，在咬入相與咬出相階段由于咬合力的原因，使得該時間段驅動力較大。顳肌相對于咬肌與翼外肌具有不同的方向，同時具有較大的數值。最大驅動力發生在左側顳肌，達到 768.6 N。盡管右側咀嚼運動存在不對稱性，但是雙側驅動力之間的變化規律與變動范圍存在較強的相關性。由于彈性元件的引入，使得關節內力僅與彈性元件參數有關，圍繞預載荷穩定地上下浮動，最大內力為 72.9 N。電機驅動功率的變化范圍在 ? 25～20 W。由于負功率的存在，此時能量利用率為 49.82%。
4 彈性顳下頜關節參數分析
本節采用中心復合設計準則（central composite design，CCD）對彈性元件的預載荷與剛度參數進行分析，分析能量利用率、最大驅動力、最大關節內力的受控情況。從參數取值范圍中選擇 9 個設計點進行動力學仿真實驗，設計點選取如圖 7 所示，max（·）表示取最大值，min（·）表示取最小值。為了更好地分析樣本空間，使用含交叉項的二次型函數法建立參數與動力學響應之間的函數關系，擬合形成響應曲面，如式（19）所示：
圖7
實驗設計與響應面擬合
Figure7.
Design of experiments and fit of response surfaces
 |
其中，x 為參數組成的向量；n 為向量 x 的長度；xi 或 xj 為向量 x 的第 i 或 j 個元素；y 為某動力學響應；a0 代表常數項待定系數，aj 代表一次項待定系數，aij 代表二次項待定系數。使用最小二乘法求解待定系數，得到響應面擬合結果，如圖 7 所示。
結果表明，隨著彈性元件參數的變化，能量利用率、最大驅動力與最大關節內力呈現出不同強弱程度的可控性。隨著預載荷和剛度的增加，驅動力與關節力呈近似線性同步增加，顯示出可控性較強，且前者的增加梯度大于后者。同時，兩者受預載荷的影響明顯大于剛度參數。預載荷（或剛度）較小時，剛度（或預載荷）的增加使得能量利用率上升明顯，但隨著預載荷（或剛度）的逐漸變大，剛度（或預載荷）對能量利用率的影響越小，靈敏度下降。整體來看，盡管增加關節盤的預載荷與剛度有利于能量利用率的上升，但基本維持在 49%，顯示出該指標可控性弱。一方面，較弱的可控性代表該指標主動優化空間不足，但另一方面，未發生顯著增減，又反映出指標具備穩定性。
5 結論
本文考慮關節盤的生物力學作用，提出彈性顳下頜關節結構，并進行了自由度、運動學與動力學分析，最后進行了彈性元件參數的討論，得到以下結論：① 在點接觸高副型式顳下頜關節基礎上，將髁桿與下頜拆分為獨立構件，此時咀嚼機器人為 6 自由度，且允許髁桿沿移動副發生小量的自由位移；② 在動力學上有效地穩定了顳下頜關節的受力，使其始終圍繞預載荷變化，同時驅動電機存在負功率現象；③ 彈性元件的預載荷與剛度的持續增加，驅動力與關節力呈近似線性的上升趨勢，而能量利用率幾乎維持不變。
鑒于本文主要使用彈性元件進行關節盤彈性特點的仿生，相較于人體真實的關節盤，忽略了內部阻尼所發揮的能量耗散作用，這將在接下來的研究中進一步完善。此外，彈性顳下頜關節結構的分析可為現有的人工顳下頜關節的改善提供參考。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
近年來，具有多自由度的機器人正越來越多地出現在醫療領域，如外科手術機器人、腕功能康復機器人和牙科機器人等[1-3]。其中，咀嚼機器人作為一種可以模擬人類咀嚼軌跡與咬合力的機械裝置受到牙科學研究人員的廣泛關注[4]。通過增設或改進不同的機器人專用模塊的方式，咀嚼機器人可以應用到牙科義齒測試、下頜康復訓練和學科教學等多種不同的場景[5-7]。
為保證咀嚼機器人運動與受力的仿生性，咀嚼機器人的設計通常是參照口頜內的骨—肌系統展開[8]。值得注意的是，口頜系統中顳下頜關節呈現旋轉與平移的復合運動特點，發揮著穩定下頜運動與受力的重要作用[9]。為了使得咀嚼機器人能夠復現這一生物學特點，如何構建符合人體顳下頜關節的仿生設計成為研究熱點[8]。早在咀嚼機器人引入仿生顳下頜關節之前，牙科學中用于模擬下頜運動的咬合架已經將顳下頜關節以鉸鏈型式、上窩下球型式和上球下窩型式加入到機構設計中，但均無法實現仿生髁突的復合運動[10]。吉林大學孫鐘雷等[11]在設計用于食品物性參數檢測的仿生咀嚼裝置（僅實現開閉口運動）時將下頜支、髁突與關節窩分別簡化為圓柱、球體與活動球槽，但是并未對顳下頜關節作進一步的機構學分析。新西蘭大學 Cheng 等[12-13]和 Mostashiri 等[14]在并聯咀嚼機器人（拓撲結構符號[15]：6—SOC{—R—S—S—}）的基礎上同樣將髁突簡化為球體，但將顳骨關節面簡化為貫通式滑道，進行自由度計算與運動學分析。在該研究中，滑道上曲面或下曲面變為重復限制運動的虛約束，盡管能在不影響機構運動的前提下增加機構的剛性，改善其受力狀況，但是由于虛約束對機構的幾何條件要求較高，對機構的加工和裝配精度亦提出了較高的要求。大連理工大學溫海營等[16-18]則針對并聯咀嚼機器人（拓撲結構符號[15]：6—SOC{—P—U—S—}）提出一種仿生性更強的點接觸高副型式（high kinematic pair，HKP）顳下頜關節，該關節除球形髁突外，特別提取人體顳骨關節面解剖生理結構尺寸用于設計仿生顳骨關節面，充分利用了口頜系統的冗余特性，體現出仿生的運動學特征，并且實現了驅動器能量的高效利用，但是單面剛性約束導致髁突易于脫離顳骨關節面，造成關節受力的非規律性波動。
目前的仿生策略忽略彈性關節盤，僅使用點接觸高副型式對剛性顳骨關節面與髁突進行接觸面形貌仿生，存在關節自由度嚴格受限與受力波動的問題。本文創新性將髁桿與下頜分離為獨立構件，并且使用彈性元件模擬關節盤，形成彈性顳下頜關節；進一步，對具有新關節的咀嚼機器人進行自由度、運動學與動力學分析，并對彈性元件的預載荷與剛度參數的影響規律進行討論。通過上述研究，期望可以為咀嚼機器人開發人員提供新的顳下頜關節機構型式，實現咀嚼機器人仿生性的進一步提升，另一方面也可以為人體其他位置關節盤的仿生設計提供新思路。
1 引入彈性顳下頜關節的咀嚼機器人機構
1.1 顳下頜關節的仿生設計
顳下頜關節由顳骨關節面、關節盤、髁突及其他軟組織構成。其中，關節盤是人體顳下頜關節的重要組成部分，具有緩沖振蕩、分散載荷、穩定關節的生物力學作用[19-20]。由于粘彈性關節盤的體積較小并且具有各項異性的力學特點，若選擇直接物理重建的仿生方式將對三維立體打印技術和粘彈性材料制備技術提出較高的要求。因此，本文僅對關節盤緩沖振蕩的生物力學功能進行間接的仿生。
首先，根據生理髁突的形貌數據，設計末端帶有直徑為 5 mm 球體的桿件模擬髁突，稱為髁桿（其中，球體稱為髁球）。同時，使用分段三次樣條曲線對生理髁突標記點 O 的矢狀面運動曲線進行插值，作為髁球中心運動軌跡 z = f（x）。將髁球中心 E 運動軌跡沿髁球半徑在矢狀面內偏置，得到仿生顳骨關節面（假設顳骨關節面任意法向與矢狀面始終保持平行）。在髁桿與仿生下頜之間使用壓縮彈簧模擬關節盤，同時在仿生下頜中設置滑道約束髁桿僅作相對平移運動。GB/T 2089-2009《普通圓柱螺旋壓縮彈簧尺寸及參數（兩端圈并緊磨平或制扁）》[21]中提及了圓柱螺旋壓縮彈簧標準，本文根據其中不同彈簧直徑下的彈簧剛度與最大工作載荷數據表，可確定設計變量的變化范圍，取預載荷 Fpre 變化范圍為 10～100 N，初值為 55 N；彈簧剛度 kspring 變化范圍為 10～50 N/mm，初值為 30 N/mm。基于以上參數，由顳骨關節面、髁桿與彈性元件共同構成彈性的顳下頜關節結構，并形成兩種運動副：點接觸高副、移動副。在咀嚼運動時，通過彈性元件伸縮（彈性元件下端點 A 與髁球中心 E 下移，髁球中心 E 與生理髁突標記點 O 不再重合），保持髁球與顳骨關節面不發生脫離，同時避免單純點接觸高副帶來的接觸面碰撞問題，如圖 1 所示。
圖1
彈性顳下頜關節仿生結構與工作原理
Figure1.
Biomimetic structure of elastic temporomandibular joint and its working principle
1.2 咀嚼機器人機構分析
在彈性顳下頜關節基礎上，采用 6 條驅動支鏈模擬咀嚼肌群，包括：左側咬肌、右側咬肌、左側顳肌、右側顳肌、左側翼外肌和右側翼外肌，形成新型咀嚼機器人機構型式，其三維模型如圖 2 所示，機構型式如圖 3 所示。根據庫茲貝奇準則，計算咀嚼機器人機構自由度，如式（1）所示：
圖2
咀嚼機器人三維模型
Figure2.
3D model of the masticatory robot
圖3
咀嚼機器人機構型式
Figure3.
Mechanism of the masticatory robot
|
其中，m 是自由度數目；n 是構件數目，共 16 個，分別為 6 個連桿、6 個滑塊、2 個髁桿、1 個下頜平臺和 1 個靜平臺；ji 是具有 i 個自由度的關節數目，包括：移動副 j1 = 8（髁桿與下頜之間 2 個移動副、滑塊與絲杠之間 6 個移動副）、萬向節 j2 = 6、球副 j3 = 6、4 自由度運動副 j4 = 0、點接觸高副 j5 = 2。由此，可計算得到機器人的自由度是 6。尤其是根據并聯彈性驅動器的定義，采用 2 條仿顳下頜關節彈性支鏈與 6 條仿肌肉驅動支鏈，共同構成并聯彈性驅動器[22]。
2 彈性顳下頜關節運動學分析
2.1 坐標系與咀嚼軌跡定義
為了便于描述接下來的運動學與動力學方程，首先建立咀嚼機器人動靜坐標系。左側髁突基坐標系為{G}，其原點 OG 位于左側生理髁突標記點，x、y 和 z 向分別為咀嚼機器人前方、左側方和上方，下文簡稱基坐標系。下頜連體坐標系為{L}，其原點 O 與下頜固結，隨下頜運動。在閉口時，下頜連體坐標系與基坐標系重合。
咀嚼時序為從正中頜位起始，通過右側咀嚼運動，回歸正中頜位。綜合考慮切牙“淚滴狀”的軌跡特征、非工作側髁突在矢狀面與水平面的投影特征以及工作側髁突的錐形面位移約束特征，進行右側咀嚼軌跡的定義，如圖 4 所示。
圖4
切牙咀嚼軌跡（相對于基坐標系{G}）
Figure4.
Chewing trajectory of the incisor (in the ground coordinate system {G})
2.2 運動學分析
2.2.1 位置分析
為了便于運動學分析，使用下頜連體坐標系相對于基坐標系的 6 個位姿參數進行下頜運動描述，即：x、y、z、α、β 和 γ。其中，α、β 和 γ 是歐拉角，依次為繞基坐標系 x 軸翻滾、y 軸俯仰和 z 軸偏航，x、y 和 z 是下頜連體坐標系原點 O 的平移量。
下頜任意一點 P 在基坐標系{G}的位置向量可以通過如下坐標變換矩陣求得，如式（2）所示：
|
其中，SP 為任意一點 P 在基坐標系{G}中的坐標（列向量）；SO 為下頜連體坐標系{L}中原點 O 在基坐標系{G}中的坐標（列向量）；LSOSP 是在下頜連體坐標系{L}中任意一點 P 相對于原點 O 的位置向量（列向量）；GRL 為下頜連體坐標系{L}在基坐標系{G}中的方位，如式（3）所示：
|
其中，c 代表 cos（·），s 代表 sin（·）。
以左側髁桿為例進行位置分析。彈簧下端點 A 在基坐標系下的坐標如式（4）所示：
|
推導出髁桿方向向量 n，如式（5）所示：
|
髁球中心 E 具有兩個幾何約束：球心始終處于髁球中心運動軌跡曲線上，滿足 z = f（x）；球心受髁桿方向向量約束。據此聯立方程組，求解髁球中心坐標，如式（6）所示。
|
其中，e1 與 e3 分別為 x 向與 z 向單位向量（行向量）。根據 E 與 A 的坐標可換算出左側彈簧長度變化量 H，同理可求出右側彈簧長度變化量，如圖 5 所示。
圖5
彈簧變形量及變形速度
Figure5.
Deformation and changing rate of the spring
2.2.2 速度分析
同樣地，以左側髁桿為例進行速度分析。已知髁球中心的瞬時速度矢量與髁球中心軌跡平面相切。將待求的髁球中心速度 VE 沿切平面進行矢量分解，并代入由切平面產生的運動速度約束，如式（7）所示：
| '/> |
進一步地，髁球中心在不同坐標系中的相對速度關系如式（8）所示：
|
其中，v 為彈簧變形速度；LVE 為髁球中心在下頜連體坐標系中的速度分量，可根據不同坐標系下速度變換關系進行求解，如式（9）所示：
|
其中，ωmandi 為仿生下頜在基坐標系下角速度。聯立式（7）～（9），推導得到左側彈簧變形速度 v，同理可求出右側彈簧變形速度，如圖 5 所示。
結果表明，工作側顳下頜關節彈性元件的變化量未超過 1.7 mm，近似呈現正弦波動。非工作側變化量超過工作側髁突變化量，波動程度較小。相比于單純的點接觸高副型式顳下頜關節，彈性元件的引入允許髁突在顳下頜關節面下實現更靈活的空間運動，而不發生干涉。該變化量亦與人體關節盤較小的結構尺寸相適應。最大變形速度發生在非工作側髁突的快速閉口階段，接近 40 mm/s。較大的變化速率體現出人體關節盤具備優良的緩沖髁突速度沖擊的能力。同樣地，非工作側髁突變化速率超過工作側髁突，且兩者均呈現波動。
3 彈性顳下頜關節動力學分析
為了進一步研究顳下頜關節對咀嚼機器人受力的影響，本文采用第二類拉格朗日方程建立動力學方程，如式（10）所示：
|
其中，U 為拉格朗日函數，定義為動能之和 T 減去勢能之和 W；qk 為咀嚼機器人第 k 項廣義坐標（即下頜 6 個位姿參數）；Qk 為咀嚼機器人非保守力在第 k 項廣義坐標下的廣義力之和。
3.1 動能和勢能的計算
3.1.1 仿生下頜動能和勢能
仿生下頜的動能 Kmandi 如式（11）所示：
|
其中，mmandi 為仿生下頜的質量；Vmandi 為仿生下頜質心在基坐標系下的線速度；Imandi 為仿生下頜在基坐標系中關于質心的慣性矩陣；ωmandi 為仿生下頜在基坐標系中的角速度。
仿生下頜的勢能 Pmandi 如式（12）所示：
|
其中，Zmandi 為仿生下頜質心在基坐標系中 z 方向的坐標；g 為重力加速度。
3.1.2 彈性元件勢能
左側彈性元件的勢能 Pspring 如式（13）所示：
|
同理，計算出右側彈性元件的勢能。
3.1.3 6 條驅動支鏈動能和勢能
第 i 條驅動支鏈的動能 Kchain，i 如式（14）所示：
|
其中，mslid，i 為滑塊的質量；Vslid，i 為滑塊的速度；Vlink，i 為連桿質心在基坐標系下的線速度；mlink，i 為連桿的質量；Ilink，i 為連桿在基坐標系下對質心的慣量矩陣；ωlink，i 為仿生下頜在基坐標系中的角速度。
第 i 條驅動支鏈的勢能 Pchain，i 如式（15）所示：
|
其中，Zslid，i 與 Zlink，i 為滑塊和連桿質心在基坐標系中 z 方向的坐標。
3.2 廣義力的計算
咀嚼機器人力學體系中非保守力包括 6 組驅動力 Fslid，i 與咬合力 Fmolar，對應的是廣義力為 Qslid，k 與 Qmolar，k；保守力包括仿生下頜所受到的彈簧力，可以根據彈簧預載荷與變形量直接算出。
3.2.1 驅動力廣義力
根據驅動力 τi 與移動副位移為 Sslid，i，換算出 6 組驅動力對應的廣義力，如式（16）所示：
|
3.2.2 咬合力廣義力
牙科模型咬合力 Fmolar 作用在右側第一磨牙上，咬合力在正中頜時達到最大值 150 N，在咬合相階段呈現正弦變化規律，方向保持垂直向上；作用點 z 軸位移為 Zmolar，換算出咬合力對應的廣義力，如式（17）所示：
|
將式（11）～（17）代入拉格朗日方程式（10），得到完整的咀嚼機器人動力學模型，并解算得到驅動力 Fslid，i 與驅動功率 Pslid，i，如圖 6 所示。進一步分析能量利用率，如式（18）所示：
圖6
電機驅動力與驅動功率
Figure6.
Driving force and driving power
|
其中，η 為能量利用率；Wall 為消耗的總功；Wpos 為正功；Wneg 為負功。
結果表明，在咬入相與咬出相階段由于咬合力的原因，使得該時間段驅動力較大。顳肌相對于咬肌與翼外肌具有不同的方向，同時具有較大的數值。最大驅動力發生在左側顳肌，達到 768.6 N。盡管右側咀嚼運動存在不對稱性，但是雙側驅動力之間的變化規律與變動范圍存在較強的相關性。由于彈性元件的引入，使得關節內力僅與彈性元件參數有關，圍繞預載荷穩定地上下浮動，最大內力為 72.9 N。電機驅動功率的變化范圍在 ? 25～20 W。由于負功率的存在，此時能量利用率為 49.82%。
4 彈性顳下頜關節參數分析
本節采用中心復合設計準則（central composite design，CCD）對彈性元件的預載荷與剛度參數進行分析，分析能量利用率、最大驅動力、最大關節內力的受控情況。從參數取值范圍中選擇 9 個設計點進行動力學仿真實驗，設計點選取如圖 7 所示，max（·）表示取最大值，min（·）表示取最小值。為了更好地分析樣本空間，使用含交叉項的二次型函數法建立參數與動力學響應之間的函數關系，擬合形成響應曲面，如式（19）所示：
圖7
實驗設計與響應面擬合
Figure7.
Design of experiments and fit of response surfaces
|
其中，x 為參數組成的向量；n 為向量 x 的長度；xi 或 xj 為向量 x 的第 i 或 j 個元素；y 為某動力學響應；a0 代表常數項待定系數，aj 代表一次項待定系數，aij 代表二次項待定系數。使用最小二乘法求解待定系數，得到響應面擬合結果，如圖 7 所示。
結果表明，隨著彈性元件參數的變化，能量利用率、最大驅動力與最大關節內力呈現出不同強弱程度的可控性。隨著預載荷和剛度的增加，驅動力與關節力呈近似線性同步增加，顯示出可控性較強，且前者的增加梯度大于后者。同時，兩者受預載荷的影響明顯大于剛度參數。預載荷（或剛度）較小時，剛度（或預載荷）的增加使得能量利用率上升明顯，但隨著預載荷（或剛度）的逐漸變大，剛度（或預載荷）對能量利用率的影響越小，靈敏度下降。整體來看，盡管增加關節盤的預載荷與剛度有利于能量利用率的上升，但基本維持在 49%，顯示出該指標可控性弱。一方面，較弱的可控性代表該指標主動優化空間不足，但另一方面，未發生顯著增減，又反映出指標具備穩定性。
5 結論
本文考慮關節盤的生物力學作用，提出彈性顳下頜關節結構，并進行了自由度、運動學與動力學分析，最后進行了彈性元件參數的討論，得到以下結論：① 在點接觸高副型式顳下頜關節基礎上，將髁桿與下頜拆分為獨立構件，此時咀嚼機器人為 6 自由度，且允許髁桿沿移動副發生小量的自由位移；② 在動力學上有效地穩定了顳下頜關節的受力，使其始終圍繞預載荷變化，同時驅動電機存在負功率現象；③ 彈性元件的預載荷與剛度的持續增加，驅動力與關節力呈近似線性的上升趨勢，而能量利用率幾乎維持不變。
鑒于本文主要使用彈性元件進行關節盤彈性特點的仿生，相較于人體真實的關節盤，忽略了內部阻尼所發揮的能量耗散作用，這將在接下來的研究中進一步完善。此外，彈性顳下頜關節結構的分析可為現有的人工顳下頜關節的改善提供參考。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
