基于本構模型參數的肌肉主動力響應有限元分析

2014-11-20 17:20:45李凡粟思橙胡偉黃晶楊濟匡

湖南大學學報·自然科學版 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.