基于本构模型参数的肌肉主动力响应有限元分析

2014-11-20李凡粟思橙胡伟黄晶杨济匡

湖南大学学报·自然科学版 2014年10期

李凡+粟思橙+胡伟+黄晶+杨济匡

基金项目：国家自然科学基金资助项目（51205117）

作者简介：李凡（1981-），男，湖南邵阳人，湖南大学助理教授，博士，硕士研究生导师

通讯联系人，Email：lifandudu@163.com

（湖南大学汽车车身先进设计制造国家重点实验室，湖南长沙 410082）摘要：研究了肌肉主动力有限元模型的力学性能与稳定性，使其适用于头颈部动力学响应的仿真研究.运用准线性黏弹性被动单元与Hill主动单元相耦合的方法，建立了兔子胫前肌的有限元仿真模型.对比Myers实验数据验证了该模型在不同应变率下离心收缩的应力应变特性.同时研究了等长收缩与轴向压缩模式下的肌肉的力学特性.仿真结果表明，仿真模型与实验输出应力应变曲线具有较好的一致性，在相同应变情况下，最大应力误差仅为0.07 MPa.该模型具有较好的计算稳定性与生物逼真度，能够满足人体颈部肌肉建模的需求.

关键词：生物力学；肌肉主动力；有限元分析；准线性黏弹性材料；本构模型；耦合模型

中图分类号：R318.01

人体颈部肌肉的主被动作用在头颈部动力学响应中具有重要的意义.低速碰撞环境中，驾驶员由于视觉或触觉引发身体应激反应，使得肌肉紧张产生主动力.肌肉主动力对头颈部运动的稳定性乃至减小碰撞中头颈部损伤具有积极的作用[1-4].因而在研究车辆碰撞特别是后碰撞中的头颈部损伤时，有必要考虑颈部肌肉主动力的影响.

有限元方法被广泛应用于被动安全与损伤生物力学领域，常见的模拟肌肉主动力的方法是采用单轴的Hill肌肉本构模型[5-8].该方法最早将肌肉简化为起点与止点之间的一维直线段，并不考虑几何以及肌肉间交互作用的影响.Wittek[9]以串联的Hill单元分段线性地拟合肌肉走向与收缩力作用方向，然而串联的Hill单元导致了不稳定问题.基于连续介质的肌肉本构模型研究中，大多提出以复合材料的形式将肌肉的整体响应分解为线性叠加的主动与被动本构模型[10-13].主动本构模型方面，Johansson等[14]最早将单轴的Hill模型修改成递归公式以应用于复合本构模型，随后的多数研究都沿用了该方式[10-12].被动本构模型往往以应变能密度函数的形式，分别定义体积应变响应、等容超弹性响应和等容黏弹性响应并叠加得到被动响应.通过对各分量定义不同的应变能函数，形成了多样化的被动本构模型.实际应用中，Behr[15]以及Hedenstierna等[16]依据Wittek＼[9＼]的建议提出了在结构而非本构模型层面耦合肌肉主被动响应的方法，并分别应用于人体下肢与颈部有限元研究.该方法采用共节点的方式以实体单元抑制串联Hill单元引起的不稳定，是一种简单、高效的建模方式.然而Behr等＼[15＼]并未考虑肌肉组织的应变率效应，且为了控制模型的不稳定，采用了平缓的激活与加载曲线.这既降低了模型的逼真度也限制了其应用范围.Hedenstierna等＼[16＼]分别采用Ogden模型（超弹性橡胶材料）与QLV模型评估了该种模型在多种应变率下的性能并维持了高应变率下模型的稳定性.结果表明QLV模型的精度更高，稳定性却较差.但是为了控制激活态下模型的稳定性Hedenstierna调高了较大应变（10%以上）时的应力输出.

国内关于肌肉的有限元研究较少.多数生物力学有限元模型中，肌肉作用直接以离散载荷施加于模型中[17-18].在动态研究中，肌肉同样被化简为一维单元，且往往仅保留了被动响应[19].部分模型采用实体单元模拟肌肉，但仍未添加主动响应[20].本构模型研究方面，龚亚琦等[21]将Johansson等[14]与Tang等[12]的研究相结合提出了包含疲劳因素的连续介质本构模型.

本文目的在于研究一种以实体单元与梁单元共节点方式耦合主被动响应的肌肉有限元建模方法，并分析讨论模型的力学特性与稳定性特征.通过调整主动与被动本构模型的参数改善模型的生物逼真度并提高其在被动与激活态时的计算稳定性，以满足不同应变率下肌肉生物力学响应的模拟仿真要求.

1 方法与材料

建立了兔子胫前肌（Tibialis Anterior， TA）几何模型.将准线性黏弹性（QLV）本构模型与Hill主动本构模型耦合，建立了TA的肌肉主动力模型.根据Myers实验验证了该模型在不同应变率下肌肉离心收缩的应力应变特性，同时研究了等长收缩与轴向压缩模式下的肌肉的力学特性.

1.1 肌肉几何模型

胫前肌的几何模型根据Myers等[22-23]实验样本的横截面轮廓重建而得.胫前肌的总长度与等效生理横截面积（Physiologic Cross Sectional Area， PCSA）分别为90 mm与45 mm2.两端肌腱长14 mm＼[11＼]，起点与止点的横截面积分别为21 mm2与12.5 mm2.兔子胫前肌有限元模型如图1所示.

1.2 肌肉材料模型

材料实验研究表明[24-25]，肌肉的被动响应具有高度非线性的黏弹性.因而实体单元材料采用了LSDYNA中MAT_176 （QLV Model）模型[26].QLV模型的黏弹性响应是根据应变历史对非线性的瞬时弹性响应进行遗传积分而得，其瞬时弹性响应与黏弹性应力为：

1.3 模型验证与工况模拟

Myers等人[22-23]为了研究骨骼肌对人体颈部在冲击环境中动力学响应的影响，以兔子胫前肌为样本进行了多种应变率（1/s，10/s，25/s）的活体单轴拉伸实验.Myers量化了兔子胫前肌在不同应变率下的主、被动响应，并认为其可用于人体颈部动力学的数值研究.该实验因其丰富的数据与适合的实验条件设置而经常为人体颈部数值研究所引用[1-2，16].

建立的有限元模型进行离心收缩、等长收缩和轴向压缩共5组9次仿真.仿真边界条件设置见表2.除D组以外，所有仿真都分别在肌肉未激活与完全激活2种状态下进行加载（未激活与完全激活状态分别以0和1表示，A0即1/s应变率下未激活时的离心收缩仿真，其余以此类推）.A，B，C 3组仿真分别以3种应变率加载，并与Myers的实验结果进行对比以验证模型的有效性；D，E 2组仿真仅研究模型在不同载荷下的力学特征与稳定性.

2 结果

仿真结果与实验数据的对比如图3所示.图中3类不同走向的曲线分别代表了TA在未激活、完全激活状态下的工程应力曲线（被动响应与整体响应）以及两者的差值（主动响应）.

图3中分别给出了Myers实验数据、TA模型仿真结果以及Hedenstierna的研究结果（分别以Myers、仿真曲线以及Hedenstierna标记），为简化起见，仅给出了Hedenstierna的主动响应以及25/s应变率下的被动响应曲线.由图可知，未激活时TA模型的计算结果在加载初期存在一定偏差.在B0与C0组仿真中该偏差的峰值分别为0.07 MPa与0.06 MPa，A0组的误差较小.C0组仿真结果对实验数据拟合较好，未出现Hedenstierna研究中应力偏低以及不稳定的状况.当完全激活时，B1与C1组仿真的整体响应最大误差皆约为0.07 MPa，对比B0与C0组的结果可知该误差主要源自于被动响应.A1组仿真中，较大应变下的主动响应出现了下滑趋势，25%应变时误差为0.05 MPa.等长收缩仿真中，在初始激活时刻模型内产生了高频振动的速度场，然而其幅值非常小，轴向最大振动速度约为0.3 m/s，且在10 ms内迅速衰减至0.05 m/s.肌腹处的截面输出力也很稳定.当轴向压缩时，模型激活与否的响应差距较大，如图4所示.被动的模型在压缩加载的早期便屈曲失稳，呈C形状，而激活的模型在工程应变达到20%时仍能保持均匀稳定的压缩状态.

3 讨论

本文目的在于研究以共节点方式耦合肌肉主被动响应的相关问题.所建模型从应力预测的精确性与模型稳定性两方面研究.TA模型的离心收缩仿真结果表明，QLV模型能够较好地模拟肌肉组织的拉伸特性.这与Hedenstierna＼[16＼]的结论一致.在未激活状态的仿真中，QLV模型精确地预测了3种应变率下肌肉组织的应力曲线，且C0组仿真并未出现Hedenstierna＼[16＼]的研究中明显的应力低估＼[图3（c）中ε＼]以及不稳定现象.这是由于Hedenstierna＼[16＼]确定的黏性参数是以Best[33]的松弛实验调整而得，本文是依据模型本构方程拟合获得.此外Hedenstierna＼[16＼]以Davis＼[27＼]测定的应力曲线为目标拟合了瞬时弹性响应参数，这既限定了QLV模型黏性响应的变动空间，也与Myers的实验数据存在不兼容.B0与C0仿真中，加载初期较大的偏差是源自于以阶跃函数设置的等应变率速度边界条件.TA模型在由静止转变为高速拉伸的运动过程中产生了附加动载荷，计算应力时并未对该波动进行处理.完全激活时，Hill模型也能精确地描述胫前肌的主动收缩特性.Hedenstierna＼[16＼]在研究中为了控制模型稳定性，人为调高在较大应变时主动响应（图3中δ），导致整体响应偏高约15%.本文避免了该问题.B1与C1仿真在加载初期的应力误差除了继承于被动响应的部分，也受到了主动响应的影响.该部分误差是由于主动单元在高速伸长时主动收缩力受到Hill模型速度曲线的放大.A1仿真中应变超过10%时主动响应的下滑主要是由于TA模型中肌肉与肌腱接合部分的刚度过渡较大，近似于准静态的载荷条件使得接合处的主动单元未被充分拉伸，而中央肌腹部分的主动单元被过度拉伸，在长度曲线的影响下，造成了主动收缩力的下滑.

与Ogden模型在永久的弹性响应上叠加线性黏性应力不同，QLV模型是对非线性的瞬时弹性响应遗传积分得到黏性应力，所以通过改变瞬时弹性响应，QLV模型能够更广泛地拟合非线性黏性特性.然而，由于缺乏稳定的弹性响应，基于QLV模型建立的肌肉有限元模型存在较大的不稳定风险.这些不稳定一方面是源自于Hill模型的特性，另一方面又与QLV模型迅速衰减的松弛特性相关.为了拟合肌肉组织的高度非线性拉伸特性，QLV模型的应力水平在松弛阶段最初的几十毫秒内将迅速衰减到峰值的40%，并且继续随应变的变动而变动，这种迅速松弛的特征很容易破坏模型的平衡状态.同时，串联的Hill单元会导致不稳定.按照不稳定性的来源可以将其分为长度与速度不稳定.根据长度曲线，串联的Hill单元在均匀伸长时，节点位置的扰动将导致相邻Hill单元产生不同的收缩力，节点力的不平衡将对扰动形成正反馈促使其发散.长度曲线的斜率越大，模型越不稳定.同理，模型长度不变时，节点的速度扰动将导致节点力不平衡，而该不平衡将对扰动形成负反馈导致高频振动.两类不稳定严重限制了模型的应用.等长收缩时观察到TA模型的高频振动即速度不稳定.表现为位移场与截面力的高频振动.在离心收缩以及较大应变的等长收缩时，可观察到沙漏状不均匀网格变形，即长度不稳定.该类不稳定往往从某些易发区域产生并迅速扩展到整个模型.不稳定的扩散过程与肌肉的不可压缩性相关：某Hill单元的发散性伸长将强制其依附的实体单元棱边的相邻或相对棱边收缩，形成连锁效应.模型的速度不稳定可以通过适当提高Vmax，增加并联以及减少串联的Hill单元个数3种方法共同控制.而长度不稳定则需要调整Hill单元特性分布，如图2所示，沿轴向由两端往中央逐步增加长度特性的斜率进行控制.此外，将QLV的部分刚度转移至Hill单元，由Hill单元的非线性弹簧提供弹性应力以及避免过高的Vmax等都有助于提高串联Hill单元的稳定性.由图3可知以上稳定性控制措施对模型的力学响应影响较小.经调整后的TA模型在离心收缩、等长收缩以及轴向压缩等工况下都能避免不稳定问题，输出稳定的结果.

本文中兔子胫前肌的几何模型是通过文献中的横截面积与长度重建而成，存在一定的几何误差.此外，由于目的是研究在冲击环境下颈部肌肉的力学响应，所以在拟合模型参数时，未考虑肌肉组织在长时间范围（100 ms以上）内的松弛特性以及压缩特性.

4 结论

QLV本构模型模拟肌肉组织的被动拉伸特性具有较好的生物逼真性.通过分别调整串联与并联的Hill单元个数以及调整Hill单元长度特性曲线斜率的分布可以在不影响仿真精度的情况下提高激活态下主被动耦合模型的稳定性.本文建立的肌肉有限元主被动耦合模型可以应用于汽车碰撞中的头颈部生物力学研究.

参考文献

[1] MERTZ H J，PATRICK L M. Strength and response of the human neck [C]//Proceedings of Stapp Car Crash Conference. Coronado， California：Society of Automotive Engineers， 1971：207-255.

[2] BRAULT J R， SIEGMUND G P， WHEELER J B. Cervical muscle response during whiplash： evidence of a lengthening muscle contraction[J]. Clinical Biomechanics， 2000， 15（6）： 426-435.

[3 ] SIEGMUND G P， SANDERSON D J， MYERS B S， et al. Awareness affects the response of human subjects exposed to a single whiplashlike perturbation[J]. Spine， 2003， 28（7）： 671-679.

[4] KUMAR S， FERRARI R， NARAYAN Y. Kinematic and electromyographic response to whiplash loading in lowvelocity whiplash impacts—a review[J]. Clinical Biomechanics， 2005， 20（4）： 343-356.

[5] DE JAGER M K J. Mathematical headneck models for acceleration impacts[D]. Eindhoven Dutch： Eindhoven University of Technology， 1996：88-92.

[6] CHANCEY V C， NIGHTINGALE R W， VAN EE C A， et al. Improved estimation of human neck tensile tolerance： reducing the range of reported tolerance using anthropometrically correct muscles and optimized physiologic initial conditions[J]. Stapp Car Crash Journal， 2003， 47： 135-153.

[7] BROLIN K， HALLDIN P， LEIJONHUFVUD I. The effect of muscle activation on neck response[J]. Traffic Injury Prevention， 2005， 6（1）： 67-76.

[8] FICE J B， CRONIN D S. Investigation of whiplash injuries in the upper cervical spine using a detailed neck model[J]. Journal of Biomechanics， 2012， 45（6）： 1098-1102.

[9] WITTEK A. Mathematical modeling of the muscle effects on the human body responses under transient loads example of the headneck complex [D]. Chalmers： Department of Applied Mechanics， Chalmers University of Technology， 2000：1-20.

[10] TSUI C P， TANG C Y， LEUNG C P， et al. Active finite element analysis of skeletal muscletendon complex during isometric， shortening and lengthening contraction [J]. BioMedical Materials and Engineering， 2004， 14（3）： 271-279.

[11] LU Y T， ZHU H X， RICHMOND S， et al. A viscohyperelastic model for skeletal muscle tissue under high strain rates[J]. Journal of Biomechanics， 2010， 43（13）： 2629-2632.

[12] TANG C Y， TSUI C P， STOJANOVIC B， et al. Finite element modelling of skeletal muscles coupled with fatigue[J]. International Journal of Mechanical Sciences， 2007， 49（10）： 1179-1191.

[13] EHRET A E， BL M， ITSKOV M. A continuum constitutive model for the active behaviour of skeletal muscle[J]. Journal of the Mechanics and Physics of Solids， 2011， 59（3）： 625-636.

[14] JOHANSSON T， MEIER P， BLICKHAN R. A finiteelement model for the mechanical analysis of skeletal muscles [J]. Journal of Theoretical Biology， 2000， 206（1）： 131-149.

[15] BEHR M， ARNOUX P J， SERRE T， et al. Tonic finite element model of the lower limb[J]. Journal of Biomechanical Engineering， 2006， 128（2）： 223.

[16] HEDENSTIERNA S， HALLDIN P， BROLIN K. Evaluation of a combination of continuum and truss finite elements in a model of passive and active muscle tissue[J]. Computer Methods in Biomechanics and Biomedical Engineering， 2008， 11（6）： 627-639.

[17] 魏高峰. 人体骨肌系统的整体生物力学建模与仿真分析研究[D]. 上海交通大学机械与动力工程学院， 2010：56-75.

WEI Gaofeng. Research on the biomechanical modeling and simulation of human musculoskeletal system [D]. Shanghai： School of Mechanical and Power Engineering， Shanghai Jiaotong University， 2010：56-75. （In Chinese）

[18] 仲照明，宫赫，肖智韬，等. 髋臼力和肌肉力对股骨有限元分析结果的影响[J]. 哈尔滨工业大学学报， 2011， 43（S1）： 303-308.

ZHONG Zhaoming， GONG He， XIAO Zhitao， et al. Effects of hip joint force and muscle forces on the results of finite element analysis of human femur [J]. Journal of Harbin Institute of Technology， 2011，43（S1）： 303-308. （In Chinese）

[19] 杨济匡，姚剑峰. 人体颈部动力学响应分析有限元模型的建立和验证[J]. 湖南大学学报：自然科学版， 2003，30（4）： 40-46.

YANG Jikuang， YAO Jianfeng. Development and validation of a human neck FE model in impact loading condition [J]. Journal of Hunan University：Natural Sciences， 2003，30（4）：40-46. （In Chinese）

[20] 张治纲，李曙光，范立冬，等. 人体胸部有限元模型研究[J]. 第三军医大学学报， 2009， 31（6）： 535-537.

ZHANG Zhigang， LI Shuguang， FAN Lidong， et al. Establishment of finite element model for human thorax [J]. Journal of Third Military Medical University， 2009，31（6）： 535-537. （In Chinese）

[21] 龚亚琦，樊建平. 基于纤维增强超弹性基体模型的骨骼肌有限元模拟[J]. 中国生物医学工程学报， 2009， 28（5）： 792-796.

GONG Yaqi， FAN Jianping. Numerical analysis of skeletal muscle using model of hyperelastic matrix and active contraction fiber [J]. Chinese Journal of Biomedical Engineering， 2009， 28（5）：792-796. （In Chinese）

[22] MYERS B S， VANEE C A， CAMACHO D L， et al. On the structural and material properties of mammalian skeletal muscle and its relevance to human cervical impact dynamics [C]//Proceedings of the 39th Stapp Car Crash Conference. San Diego， California：Society of Automotive Engineers， 1995：No.952723.

[23] MYERS B S， WOOLLEY C T， SLOTTER T L， et al. The influence of strain rate on the passive and stimulated engineering stresslarge strain behavior of the rabbit tibialis anterior muscle[J]. Journal of Biomechanical Engineering， 1998， 120（1）： 126-132.

[24] ABRAHAM A C， KAUFMAN K R， HAUT DONAHUE T L. Phenomenological consequences of sectioning and bathing on passive muscle mechanics of the New Zealand white rabbit tibialis anterior[J]. Journal of the Mechanical Behavior of Biomedical Materials， 2013，17：290-295.

[25] MEYER G A， MCCULLOCH A D， LIEBER R L. A nonlinear model of passive muscle viscosity[J]. Journal of Biomechanical Engineering， 2011， 133（9）： 091007.

[26] HALLQUIST J O. LSDYNA theory manual[M]. Livermore， California： Livermore Software Technology Corporation，2006.

[27] DAVIS J， KAUFMAN K R， LIEBER R L. Correlation between active and passive isometric force and intramuscular pressure in the isolated rabbit tibialis anterior muscle[J]. Journal of Biomechanics， 2003， 36（4）： 505-512.

[28] MAGANARIS C N， PAUL J P. In vivo human tendon mechanical properties[J]. The Journal of Physiology， 1999， 521（1）： 307-313.

[29] GRAS L L， MITTON D， VIOT P， et al. Hyperelastic properties of the human sternocleidomastoideus muscle in tension[J]. Journal of the Mechanical Behavior of Biomedical Materials， 2012，15：131-140.

[30] STH J， BROLIN K， HAPPEE R. Active muscle response using feedback control of a finite element human arm model [J]. Computer Methods in Biomechanics and Biomedical Engineering， 2012， 15（4）： 347-361.

[31] WITTEK A， KAJZER J. Modeling the muscle influence on the kinematics of the headneck complex in impacts[J].Memoirs of the School of Engineering， Nagoya Univ，1997，49：155-205.

[32] ZAJAC F E. Muscle and tendon： properties， models， scaling， and application to biomechanics and motor control[J]. Critical Reviews in Biomedical Engineering， 1988， 17（4）： 359-411.

[33] BEST T M， MCELHANEY J， GARRETT W E， et al. Characterization of the passive responses of live skeletal muscle using the quasilinear theory of viscoelasticity[J]. Journal of biomechanics， 1994， 27（4）： 413-419.