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

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

0 引言

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

1 超彈性材料的基本理論

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

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

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

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

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

本文選擇縮減多項(xiàng)式本構(gòu)模型（reduced polynomial），縮減多項(xiàng)式的應(yīng)變能函數(shù)為

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

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

2.1 有限元模型的建立

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

圖1 高速攝像機(jī)下的變色龍F(tuán)ig.1 Chameleon imaged using a high-speed camera

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

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

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

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

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

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

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

圖3 實(shí)驗(yàn)測定的變色龍舌頭的速度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使舌頭達(dá)到一定伸長，之后瞬間去掉外界加載的速度。通過分析舌尖的速度下降情況得到與實(shí)驗(yàn)數(shù)據(jù)較為適當(dāng)?shù)牟牧蠀?shù)。根據(jù)人體軟組織的材參數(shù)，設(shè)置復(fù)合體的密度為970 kg/m3。為了減少各種因素對有限元模擬的影響，忽略了舌骨和復(fù)合體接觸之中的摩擦，并且采用了材料參數(shù)中C10=C20的數(shù)據(jù)。

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

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

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

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

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

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

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

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

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

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

2.3 加載力的確定

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

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

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

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

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

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

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

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

圖12 估計(jì)與實(shí)驗(yàn)速度對比Fig.12 The estimated speed is compared with the experimental speed

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

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

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

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

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

表4 預(yù)測的加載結(jié)果Tab.4 Predicted loading results

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

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

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

圖14 速度預(yù)測曲線與實(shí)驗(yàn)曲線的對比Fig.14 The comparison between the velocity prediction curve and the experimental curve

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

圖15 位移預(yù)測曲線與實(shí)驗(yàn)曲線的對比Fig.15 The comparison between the position prediction curve and the experimental curve

4 結(jié)論

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