變色龍舌頭彈射過程的有限元模擬與仿真

張 一，桑建兵，劉帥龍，王 穎

（河北工業(yè)大學 機械工程學院，天津 300401）

0 引言

生物軟組織的性質已成為生物醫(yī)學工程領域研究的一個重要課題，對于開發(fā)出性能卓越的生物材料具有十分重要的工程意義，所以研究生物軟組織的力學性能就顯得尤為重要。在對生物軟組織的研究中，發(fā)現(xiàn)變色龍有一種獨特的能力，它能以很快的速度將舌頭伸長到一倍體長。1977年，Harkness 等[1]研究發(fā)現(xiàn)變色龍通過調整其眼睛晶狀體焦點可從其聚焦機構獲得單眼距離信息，來作為調節(jié)信號以判斷距離，與此同時將舌頭緩慢伸到嘴邊。1991年，Wainwright等[2]通過高速攝影技術研究了變色龍舌頭彈射過程，研究發(fā)現(xiàn)彈射有加速彈出過程，最大加速度可達500 m/s2。2000年，Herrel等[3]研究了變色龍捕食過程，發(fā)現(xiàn)變色龍舌頭加速后有一段時間保持較為恒定的速度，之后減速并且獵物被舌頭吸附，最后將獵物收回口中。至此研究完成了變色龍捕食過程的4 個階段：估計獵物距離舌頭慢慢伸出（階段I），舌頭加速前進（階段II），在保持一段恒定速度（階段III），舌墊減速（階段IV）。2003年Groot 等[4]通過實驗得到了變色龍捕食過程中的速度，本文對比該曲線進行研究。變色龍舌的物質是不均勻的、各向異性的，它由具有超彈性性質的幾乎或完全不可壓縮材料組成。大多數(shù)生物軟組織由膠原蛋白、彈性體和平滑肌組成，這3種組織成分在生物軟組織中的比例及其相互交織形成的三維結構對其整體力學性能有很大影響。為了了解變色龍舌的力學特性，考慮了均勻各向同性材料，從而得到了簡單可靠的變形模擬結果[5-6]。傳統(tǒng)的力學試驗研究需要校準好的實驗設備和優(yōu)質的材料資源，使成本大大增加。隨著計算機的發(fā)展，生物力學實驗的數(shù)值模擬可以得到近似解，其中有限元分析方法在生物力學中得到了廣泛的應用，尤其是顯式求解方法[7]。因此，本文對變色龍舌噴射過程進行了高應變速率的數(shù)值模擬，主要采用有限元分析方法來模擬短時間內發(fā)生的彈射等動態(tài)過程。

1 超彈性材料的基本理論

隨著生物力學的發(fā)展，研究人員對生物軟組織的力學特性有了更深入的了解[8-9]。建立符合生物力學特征的軟組織模型涉及到復雜的材料參數(shù)。此外，由于難以實現(xiàn)復雜的纖維結構，大多數(shù)研究將生物軟組織材料簡化為各向同性不可壓縮超彈性材料。應力-應變關系的非線性可以用應變能函數(shù)表示，各本構模型都是應變能函數(shù)的一種特殊形式。因此，超彈性模型的本構關系可以用應變能[10]表示為

式中，I1、I2、I3為Cauchy-Green變形張量不變量。

式中：λ1，λ2，λ3為主伸長比；γi為主應變；對于不可壓縮材料

Cauchy應力張量σ的表達式可以由應變能函數(shù)W給出：

式中：p為靜水壓力，B是左Cauchy-Green應變張量。

本文選擇縮減多項式本構模型（reduced polynomial），縮減多項式的應變能函數(shù)為

對于不可壓縮材料D=0，采用二階縮減多項式模型即Yeoh模型如式（6）所示：

2 變色龍舌頭的動態(tài)仿真

2.1 有限元模型的建立

變色龍的舌頭主要由3部分組成：加速肌（AM）、復合體（RC）和舌骨（HY）。2016年Asher等[11]用高速攝像機拍攝了變色龍吐舌過程，如圖1所示，本文利用圖像中變色龍舌頭在彈出過程的尺寸信息，構建了變色龍舌頭的有限元模型如圖2 所示。其中舌骨在彈出過程中其變形相對于加速肌和復合體來說很小，故將舌骨考慮為線彈性材料，而加速肌和復合體考慮為超彈性材料。

圖1 高速攝像機下的變色龍Fig.1 Chameleon imaged using a high-speed camera

圖2 變色龍舌頭模型Fig.2 Model of chameleon tongue

2.2 材料參數(shù)的確定

采用Abaqus軟件中的顯式方法對變色龍舌頭的彈出過程進行了非線性動態(tài)仿真，雖然隱式方法是求解非線性問題最有效的方法，但顯式方法因其在求解大變形、接觸和材料非線性方面具有明顯的優(yōu)勢，因此本論文采用顯式方法進行求解。

變色龍舌骨主要起3種作用：連接舌頭與口腔并且支撐舌部組織，控制舌頭伸出的方向以及與人類一樣在吞咽時起作用。2015年，張冠軍等[12]建立了包含骨骼等詳細解剖學結構的行人下肢的有限元模型，并且利用行人下肢彎曲的生物力學實驗，對有限元模型進行了驗證，給出了幾種下肢骨骼的材料。參考了該文獻中下肢骨骼的材料參數(shù)作為本有限元模型中的舌骨的材料參數(shù)。經(jīng)過有限元仿真分析，初步確定了在時間增量步收斂情況下加速肌的材料參數(shù)。舌骨和加速肌部分的材料參數(shù)如表1所示。

表1 舌骨和加速肌的材料參數(shù)Tab.1 Material parameters of the hyoid and the accelerator muscle

2003年，Groot 等[4]通過實驗獲得了變色龍捕食的速度如圖3所示。

從圖3可以看出，當速度達到0.6 m/s 后，舌頭有一段勻速的過程。根據(jù)這個基本過程，為了估計材料參數(shù)，首先對變色龍舌頭的單軸拉伸過程進行了有限元模擬，確定的速度加載方式及邊界條件如圖4所示。軟膠原組織的力學行為與抓取方法密切相關，約束的應用影響著載荷在整個材料中的傳遞，從而影響組織的力學性能[13]。經(jīng)過多次加載嘗試，最終確定的邊界條件為：舌骨左端為固定端約束，復合體左端限制了y方向的位移以及x和z方向的轉動，加速肌和復合體之間采用綁定約束。

圖3 實驗測定的變色龍舌頭的速度Fig.3 The velocity of the chameleon's tongue was measured experimentally

圖4 確定材料參數(shù)的加載方式Fig.4 Material parameters are determined by means of loading

將加速肌加載恒定速度0.6 m/s持續(xù)0.05 s使舌頭達到一定伸長，之后瞬間去掉外界加載的速度。通過分析舌尖的速度下降情況得到與實驗數(shù)據(jù)較為適當?shù)牟牧蠀?shù)。根據(jù)人體軟組織的材參數(shù)，設置復合體的密度為970 kg/m3。為了減少各種因素對有限元模擬的影響，忽略了舌骨和復合體接觸之中的摩擦，并且采用了材料參數(shù)中C10=C20的數(shù)據(jù)。

選擇相同的節(jié)點輸出速度變換曲線，輸出速度的節(jié)點選擇如圖5所示，紅色的點即為選擇的節(jié)點。

圖5 輸出速度的節(jié)點Fig.5 Output speed node

為了確定材料參數(shù)的范圍，首先設置了材料參數(shù)分別為10 Pa、100 Pa、200 Pa、300 Pa、400 Pa。通過以上的速度加載方式，輸出了0.05 s之后幾種材料參數(shù)對應的速度，速度曲線如圖6所示。

由圖6 可以看出，時間0.05 s 之后，各種材料參數(shù)的速度曲線開始下降。當材料參數(shù)C10=10 Pa，到0.07 s時速度由0.6 m/s降至0.42 m/s，速度下降的速度較小。當材料參數(shù)擴大到100 Pa，0.05 s到0.06 s有較小的降速，但是0.06 s時有明顯的突變，并且0.07 s時速度降到負值，說明模型開始收縮，不滿足需要的條件。當材料參數(shù)大于100 Pa 到400 Pa，速度均于0.055 s 左右降至0，之后為負值，說明模型開始收縮，并且材料參數(shù)越大模型的收縮越明顯。可以從圖7總體看出，材料參數(shù)最好設置在10～100 Pa之間。

圖6 速度隨時間變化的曲線Fig.6 Curves of velocity versus time

根據(jù)之前的結論，設置了材料參數(shù)分別為10 Pa、20 Pa、30 Pa、40 Pa、50 Pa。通過以上的速度加載方式，輸出了0.05 s 之后幾種材料參數(shù)對應的速度，速度曲線如圖7所示。

由圖7 可以看出，材料參數(shù)越大，速度下降越明顯，尤其是C10=50 Pa時，速度能夠在0.02 s之內降至0.32 m/s。每次增大同樣大小的材料參數(shù)后，速度下降幅度變小。當C10=50 Pa 時，伸長后的模型的應力云圖如圖10所示。

圖7 速度隨時間變化的曲線Fig.7 Curves of velocity versus time

由圖8 可以看出，伸長后的模型應力均勻，并且整個模型均勻伸長，沒有出現(xiàn)局部的斷裂與褶皺現(xiàn)象，體現(xiàn)了在高應變率下對超彈性材料分析的可行性。通過加載的速度與實驗速度的比較，根據(jù)伸長情況，考慮到實驗速度與加載速度的區(qū)別以及其他因素的影響，最后確定材料參數(shù)為C10=C20=50 Pa。

圖8 模型伸長后的應力云圖Fig.8 Stress nephogram of the model after elongation

2.3 加載力的確定

高速攝像機拍攝的變色龍吐舌過程中舌頭伸出階段和舌頭彈出階段如圖9所示[14]，從圖9可以看出，加速肌的形狀接近凸臺型，當變色龍把舌頭從嘴里伸出來時，舌頭彈出。所以提出了一種壓力驅動模型，即在加速器的斜面施加一個壓力，依靠壓力對加速肌進行彈射。這種壓力驅動模型的加載方式如圖10所示。

圖9 高速攝像機下變色龍舌頭彈射過程Fig.9 Chameleon tongue ejection process under high-speed camera

圖10 變色龍舌頭的加載模型Fig.10 Load model of chameleon tongue

分階段加載力得到舌尖的速度曲線，將曲線和實驗曲線相比對得到一組預估的加載數(shù)據(jù)。根據(jù)實驗速度曲線，將力的加載分成了13個階段。通過反演的方式，結合有限元模擬與仿真確定出了每一時間段的加載力，最后得到了整個彈出過程的加載力數(shù)值如表2所示。

表2 變色龍舌頭彈出過程的加載力數(shù)值Tab.2 The loading force of the chameleon tongue during ejection

將表2中不同時間段的壓力數(shù)值輸入Abaqus 進行計算，得到變色龍舌頭彈射過程的Mises 應力云圖如圖11 所示。可以看出，伸長后的模型應力均勻，并且整個模型均勻伸長，沒有出現(xiàn)局部的斷裂與褶皺現(xiàn)象，體現(xiàn)了壓力驅動模型的可行性。得到了一組速度曲線，將速度曲線與實驗得到的速度曲線對比，如圖12所示。

圖11 輸入估計力的舌頭拉伸前后模型Fig.11 Input estimated force after tongue stretch before and after model

由圖12可以看出，有限元分析所得到的速度曲線和文獻中的實驗曲線基本吻合，0～0.04 s時舌頭加速過程與實驗曲線吻合較好，之后的最大速度相差大約為20%，存在較大的誤差，所以引入神經(jīng)網(wǎng)絡優(yōu)化程序對不同階段的加載力進一步的優(yōu)化。

圖12 估計與實驗速度對比Fig.12 The estimated speed is compared with the experimental speed

3 基于神經(jīng)網(wǎng)絡加載力的優(yōu)化

近年來，神經(jīng)網(wǎng)絡以其強大的非線性和自適應特性以及學習能力成為解決復雜工程問題的一種常用工具。對變色龍舌頭彈出過程進行動態(tài)仿真的過程中，通過人工干預對每一個階段的加載力進行準確預測是一項耗時巨大的工作，很難得到最優(yōu)解。為了得到較為精確的結果，準確預測變色龍舌頭彈射時的力學性能，引入了神經(jīng)網(wǎng)絡，通過神經(jīng)網(wǎng)絡和有限元相結合對加載力進行優(yōu)化。

神經(jīng)網(wǎng)絡的過程分為兩階段，一是訓練學習、二是操作預測。訓練神經(jīng)網(wǎng)絡需要能夠體現(xiàn)整體特點的樣本，所以將13個階段的估計力分別取到等差數(shù)列得到一個13因素12水平的數(shù)據(jù)集如表3所示，如果將所有數(shù)據(jù)的組合都考慮的到，則有1213組數(shù)據(jù)，樣本量巨大。

表3 加載力的等差數(shù)據(jù)集Tab.3 An arithmetic data set of loading forces

因此利用正交實驗法建立一套13因素12水平的正交實驗方案L144(1213)，該方案共有144組數(shù)據(jù)。將這144組仿真數(shù)據(jù)分別輸入到Abaqus中進行仿真計算，得到了每一組仿真數(shù)據(jù)對應的舌尖的速度-時間數(shù)據(jù)。將144組仿真數(shù)據(jù)為訓練樣本，時間和速度數(shù)據(jù)作為輸入?yún)?shù)，壓力作為輸出，對神經(jīng)網(wǎng)絡進行訓練。最后得到加載力的預測結果如表4所示。

表4 預測的加載結果Tab.4 Predicted loading results

將預測結果輸入Abaqus 進行仿真計算。拉伸前后舌頭模型對比如圖13所示。通過對舌頭彈射拉伸前后的圖像進行對比，發(fā)現(xiàn)彈射后舌頭的長度是原來長度的兩倍多。應力云圖表明，彈射后的復合體的應力高度均勻，反映了有限元建模在高應變率下進行超彈性材料單軸拉伸分析的可行性。

圖13 舌頭彈射前后對比Fig.13 Contrast before and after tongue ejection

并將有限元模擬的時間-速度曲線與實驗得到的速度進行對比，如圖14所示。從圖中可以看出，預測結果和實驗數(shù)據(jù)具有良好的一致性。在第1階段，從0～0.02 s，舌頭從口中緩慢地低速伸出。在第2階段，從0.02～0.04 s，舌頭開始有較快的加速，直到第3階段，速度達到0.6 m/s 左右。通過對該過程中2條曲線的比較，預測曲線與實驗曲線吻合較好，最大相對誤差小于8%。從而證明了模型的正確性。第3階段，舌頭以相對均勻的速度向前噴射約0.05 s。在第4階段，舌頭的速度下降，直到達到0 m/s。由圖可知，0.089 s前的仿真結果與實驗數(shù)據(jù)吻合較好。由于神經(jīng)網(wǎng)絡中ReLU激活函數(shù)的單側抑制特性，導致了0.089 s后的預測值與實驗結果之間出現(xiàn)了較為明顯的偏差。

圖14 速度預測曲線與實驗曲線的對比Fig.14 The comparison between the velocity prediction curve and the experimental curve

將有限元模擬的時間-位移曲線與實驗得到的位移進行對比，如圖15所示。從圖15中可以看出，位移預測仿真結果與實驗結果基本一致。這些結果表明，基于有限元的力學性能分析和神經(jīng)網(wǎng)絡相結合可以可靠地預測軟組織的力學性能。只要網(wǎng)絡經(jīng)過適當?shù)挠柧殻诤芏痰臅r間內就可以預測到所需要的速度所對應的壓力。因此，與僅使用有限元法進行仿真相比，該組合方法在可靠、準確預測力學性能的同時，大大減少了計算的時間成本和工作量。

圖15 位移預測曲線與實驗曲線的對比Fig.15 The comparison between the position prediction curve and the experimental curve

4 結論

根據(jù)高速攝像機下的變色龍捕食過程中舌頭的比例建立了舌頭的有限元模型，分為3個部分：加速肌、復合體以及舌骨。對舌頭彈射過程進行了有限元仿真，通過在加速肌部分加載速度的方法，通過對比初步確定了復合體材料的本構模型為Yeoh二階模型以及材料的本構參數(shù)為C10=C20=50 Pa。考慮到變色龍舌頭的彈出過程主要是由舌尖（加速肌）表面受到的壓力驅動而產生，建立了變色龍舌頭的“壓力驅動模型”，將力的加載分成了13個階段，通過反演的方式，結合有限元模擬與仿真確定出了每一時間段的加載力，得到了整個彈射過程的速度曲線。將有限元模擬與神經(jīng)網(wǎng)絡預測相結合，將所得結果與文獻中實驗所得到的速度曲線和位移曲線進行對比，結果具有較好的一致性，驗證了有限元分析和神經(jīng)網(wǎng)絡相結合解決高應變率下非線性問題的正確性。