变色龙舌头弹射过程的有限元模拟与仿真

张 一，桑建兵，刘帅龙，王 颖

（河北工业大学 机械工程学院，天津 300401）

0 引言

生物软组织的性质已成为生物医学工程领域研究的一个重要课题，对于开发出性能卓越的生物材料具有十分重要的工程意义，所以研究生物软组织的力学性能就显得尤为重要。在对生物软组织的研究中，发现变色龙有一种独特的能力，它能以很快的速度将舌头伸长到一倍体长。1977年，Harkness 等[1]研究发现变色龙通过调整其眼睛晶状体焦点可从其聚焦机构获得单眼距离信息，来作为调节信号以判断距离，与此同时将舌头缓慢伸到嘴边。1991年，Wainwright等[2]通过高速摄影技术研究了变色龙舌头弹射过程，研究发现弹射有加速弹出过程，最大加速度可达500 m/s2。2000年，Herrel等[3]研究了变色龙捕食过程，发现变色龙舌头加速后有一段时间保持较为恒定的速度，之后减速并且猎物被舌头吸附，最后将猎物收回口中。至此研究完成了变色龙捕食过程的4 个阶段：估计猎物距离舌头慢慢伸出（阶段I），舌头加速前进（阶段II），在保持一段恒定速度（阶段III），舌垫减速（阶段IV）。2003年Groot 等[4]通过实验得到了变色龙捕食过程中的速度，本文对比该曲线进行研究。变色龙舌的物质是不均匀的、各向异性的，它由具有超弹性性质的几乎或完全不可压缩材料组成。大多数生物软组织由胶原蛋白、弹性体和平滑肌组成，这3种组织成分在生物软组织中的比例及其相互交织形成的三维结构对其整体力学性能有很大影响。为了了解变色龙舌的力学特性，考虑了均匀各向同性材料，从而得到了简单可靠的变形模拟结果[5-6]。传统的力学试验研究需要校准好的实验设备和优质的材料资源，使成本大大增加。随着计算机的发展，生物力学实验的数值模拟可以得到近似解，其中有限元分析方法在生物力学中得到了广泛的应用，尤其是显式求解方法[7]。因此，本文对变色龙舌喷射过程进行了高应变速率的数值模拟，主要采用有限元分析方法来模拟短时间内发生的弹射等动态过程。

1 超弹性材料的基本理论

随着生物力学的发展，研究人员对生物软组织的力学特性有了更深入的了解[8-9]。建立符合生物力学特征的软组织模型涉及到复杂的材料参数。此外，由于难以实现复杂的纤维结构，大多数研究将生物软组织材料简化为各向同性不可压缩超弹性材料。应力-应变关系的非线性可以用应变能函数表示，各本构模型都是应变能函数的一种特殊形式。因此，超弹性模型的本构关系可以用应变能[10]表示为

式中，I1、I2、I3为Cauchy-Green变形张量不变量。

式中：λ1，λ2，λ3为主伸长比；γi为主应变；对于不可压缩材料

Cauchy应力张量σ的表达式可以由应变能函数W给出：

式中：p为静水压力，B是左Cauchy-Green应变张量。

本文选择缩减多项式本构模型（reduced polynomial），缩减多项式的应变能函数为

对于不可压缩材料D=0，采用二阶缩减多项式模型即Yeoh模型如式（6）所示：

2 变色龙舌头的动态仿真

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 材料参数的确定

采用Abaqus软件中的显式方法对变色龙舌头的弹出过程进行了非线性动态仿真，虽然隐式方法是求解非线性问题最有效的方法，但显式方法因其在求解大变形、接触和材料非线性方面具有明显的优势，因此本论文采用显式方法进行求解。

变色龙舌骨主要起3种作用：连接舌头与口腔并且支撑舌部组织，控制舌头伸出的方向以及与人类一样在吞咽时起作用。2015年，张冠军等[12]建立了包含骨骼等详细解剖学结构的行人下肢的有限元模型，并且利用行人下肢弯曲的生物力学实验，对有限元模型进行了验证，给出了几种下肢骨骼的材料。参考了该文献中下肢骨骼的材料参数作为本有限元模型中的舌骨的材料参数。经过有限元仿真分析，初步确定了在时间增量步收敛情况下加速肌的材料参数。舌骨和加速肌部分的材料参数如表1所示。

表1 舌骨和加速肌的材料参数Tab.1 Material parameters of the hyoid and the accelerator muscle

2003年，Groot 等[4]通过实验获得了变色龙捕食的速度如图3所示。

从图3可以看出，当速度达到0.6 m/s 后，舌头有一段匀速的过程。根据这个基本过程，为了估计材料参数，首先对变色龙舌头的单轴拉伸过程进行了有限元模拟，确定的速度加载方式及边界条件如图4所示。软胶原组织的力学行为与抓取方法密切相关，约束的应用影响着载荷在整个材料中的传递，从而影响组织的力学性能[13]。经过多次加载尝试，最终确定的边界条件为：舌骨左端为固定端约束，复合体左端限制了y方向的位移以及x和z方向的转动，加速肌和复合体之间采用绑定约束。

图3 实验测定的变色龙舌头的速度Fig.3 The velocity of the chameleon's tongue was measured experimentally

图4 确定材料参数的加载方式Fig.4 Material parameters are determined by means of loading

将加速肌加载恒定速度0.6 m/s持续0.05 s使舌头达到一定伸长，之后瞬间去掉外界加载的速度。通过分析舌尖的速度下降情况得到与实验数据较为适当的材料参数。根据人体软组织的材参数，设置复合体的密度为970 kg/m3。为了减少各种因素对有限元模拟的影响，忽略了舌骨和复合体接触之中的摩擦，并且采用了材料参数中C10=C20的数据。

选择相同的节点输出速度变换曲线，输出速度的节点选择如图5所示，红色的点即为选择的节点。

图5 输出速度的节点Fig.5 Output speed node

为了确定材料参数的范围，首先设置了材料参数分别为10 Pa、100 Pa、200 Pa、300 Pa、400 Pa。通过以上的速度加载方式，输出了0.05 s之后几种材料参数对应的速度，速度曲线如图6所示。

由图6 可以看出，时间0.05 s 之后，各种材料参数的速度曲线开始下降。当材料参数C10=10 Pa，到0.07 s时速度由0.6 m/s降至0.42 m/s，速度下降的速度较小。当材料参数扩大到100 Pa，0.05 s到0.06 s有较小的降速，但是0.06 s时有明显的突变，并且0.07 s时速度降到负值，说明模型开始收缩，不满足需要的条件。当材料参数大于100 Pa 到400 Pa，速度均于0.055 s 左右降至0，之后为负值，说明模型开始收缩，并且材料参数越大模型的收缩越明显。可以从图7总体看出，材料参数最好设置在10～100 Pa之间。

图6 速度随时间变化的曲线Fig.6 Curves of velocity versus time

根据之前的结论，设置了材料参数分别为10 Pa、20 Pa、30 Pa、40 Pa、50 Pa。通过以上的速度加载方式，输出了0.05 s 之后几种材料参数对应的速度，速度曲线如图7所示。

由图7 可以看出，材料参数越大，速度下降越明显，尤其是C10=50 Pa时，速度能够在0.02 s之内降至0.32 m/s。每次增大同样大小的材料参数后，速度下降幅度变小。当C10=50 Pa 时，伸长后的模型的应力云图如图10所示。

图7 速度随时间变化的曲线Fig.7 Curves of velocity versus time

由图8 可以看出，伸长后的模型应力均匀，并且整个模型均匀伸长，没有出现局部的断裂与褶皱现象，体现了在高应变率下对超弹性材料分析的可行性。通过加载的速度与实验速度的比较，根据伸长情况，考虑到实验速度与加载速度的区别以及其他因素的影响，最后确定材料参数为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

分阶段加载力得到舌尖的速度曲线，将曲线和实验曲线相比对得到一组预估的加载数据。根据实验速度曲线，将力的加载分成了13个阶段。通过反演的方式，结合有限元模拟与仿真确定出了每一时间段的加载力，最后得到了整个弹出过程的加载力数值如表2所示。

表2 变色龙舌头弹出过程的加载力数值Tab.2 The loading force of the chameleon tongue during ejection

将表2中不同时间段的压力数值输入Abaqus 进行计算，得到变色龙舌头弹射过程的Mises 应力云图如图11 所示。可以看出，伸长后的模型应力均匀，并且整个模型均匀伸长，没有出现局部的断裂与褶皱现象，体现了压力驱动模型的可行性。得到了一组速度曲线，将速度曲线与实验得到的速度曲线对比，如图12所示。

图11 输入估计力的舌头拉伸前后模型Fig.11 Input estimated force after tongue stretch before and after model

由图12可以看出，有限元分析所得到的速度曲线和文献中的实验曲线基本吻合，0～0.04 s时舌头加速过程与实验曲线吻合较好，之后的最大速度相差大约为20%，存在较大的误差，所以引入神经网络优化程序对不同阶段的加载力进一步的优化。

图12 估计与实验速度对比Fig.12 The estimated speed is compared with the experimental speed

3 基于神经网络加载力的优化

近年来，神经网络以其强大的非线性和自适应特性以及学习能力成为解决复杂工程问题的一种常用工具。对变色龙舌头弹出过程进行动态仿真的过程中，通过人工干预对每一个阶段的加载力进行准确预测是一项耗时巨大的工作，很难得到最优解。为了得到较为精确的结果，准确预测变色龙舌头弹射时的力学性能，引入了神经网络，通过神经网络和有限元相结合对加载力进行优化。

神经网络的过程分为两阶段，一是训练学习、二是操作预测。训练神经网络需要能够体现整体特点的样本，所以将13个阶段的估计力分别取到等差数列得到一个13因素12水平的数据集如表3所示，如果将所有数据的组合都考虑的到，则有1213组数据，样本量巨大。

表3 加载力的等差数据集Tab.3 An arithmetic data set of loading forces

因此利用正交实验法建立一套13因素12水平的正交实验方案L144(1213)，该方案共有144组数据。将这144组仿真数据分别输入到Abaqus中进行仿真计算，得到了每一组仿真数据对应的舌尖的速度-时间数据。将144组仿真数据为训练样本，时间和速度数据作为输入参数，压力作为输出，对神经网络进行训练。最后得到加载力的预测结果如表4所示。

表4 预测的加载结果Tab.4 Predicted loading results

将预测结果输入Abaqus 进行仿真计算。拉伸前后舌头模型对比如图13所示。通过对舌头弹射拉伸前后的图像进行对比，发现弹射后舌头的长度是原来长度的两倍多。应力云图表明，弹射后的复合体的应力高度均匀，反映了有限元建模在高应变率下进行超弹性材料单轴拉伸分析的可行性。

图13 舌头弹射前后对比Fig.13 Contrast before and after tongue ejection

并将有限元模拟的时间-速度曲线与实验得到的速度进行对比，如图14所示。从图中可以看出，预测结果和实验数据具有良好的一致性。在第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前的仿真结果与实验数据吻合较好。由于神经网络中ReLU激活函数的单侧抑制特性，导致了0.089 s后的预测值与实验结果之间出现了较为明显的偏差。

图14 速度预测曲线与实验曲线的对比Fig.14 The comparison between the velocity prediction curve and the experimental curve

将有限元模拟的时间-位移曲线与实验得到的位移进行对比，如图15所示。从图15中可以看出，位移预测仿真结果与实验结果基本一致。这些结果表明，基于有限元的力学性能分析和神经网络相结合可以可靠地预测软组织的力学性能。只要网络经过适当的训练，在很短的时间内就可以预测到所需要的速度所对应的压力。因此，与仅使用有限元法进行仿真相比，该组合方法在可靠、准确预测力学性能的同时，大大减少了计算的时间成本和工作量。

图15 位移预测曲线与实验曲线的对比Fig.15 The comparison between the position prediction curve and the experimental curve

4 结论

根据高速摄像机下的变色龙捕食过程中舌头的比例建立了舌头的有限元模型，分为3个部分：加速肌、复合体以及舌骨。对舌头弹射过程进行了有限元仿真，通过在加速肌部分加载速度的方法，通过对比初步确定了复合体材料的本构模型为Yeoh二阶模型以及材料的本构参数为C10=C20=50 Pa。考虑到变色龙舌头的弹出过程主要是由舌尖（加速肌）表面受到的压力驱动而产生，建立了变色龙舌头的“压力驱动模型”，将力的加载分成了13个阶段，通过反演的方式，结合有限元模拟与仿真确定出了每一时间段的加载力，得到了整个弹射过程的速度曲线。将有限元模拟与神经网络预测相结合，将所得结果与文献中实验所得到的速度曲线和位移曲线进行对比，结果具有较好的一致性，验证了有限元分析和神经网络相结合解决高应变率下非线性问题的正确性。