李凡+粟思橙+胡偉+黃晶+楊濟(jì)匡

基金項(xiàng)目：國家自然科學(xué)基金資助項(xiàng)目（51205117）

作者簡介：李 凡（1981-），男，湖南邵陽人，湖南大學(xué)助理教授，博士，碩士研究生導(dǎo)師

通訊聯(lián)系人，Email：lifandudu@163.com

（湖南大學(xué) 汽車車身先進(jìn)設(shè)計(jì)制造國家重點(diǎn)實(shí)驗(yàn)室，湖南 長沙 410082） 摘 要：研究了肌肉主動力有限元模型的力學(xué)性能與穩(wěn)定性，使其適用于頭頸部動力學(xué)響應(yīng)的仿真研究.運(yùn)用準(zhǔn)線性黏彈性被動單元與Hill主動單元相耦合的方法，建立了兔子脛前肌的有限元仿真模型.對比Myers實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證了該模型在不同應(yīng)變率下離心收縮的應(yīng)力應(yīng)變特性.同時(shí)研究了等長收縮與軸向壓縮模式下的肌肉的力學(xué)特性.仿真結(jié)果表明，仿真模型與實(shí)驗(yàn)輸出應(yīng)力應(yīng)變曲線具有較好的一致性，在相同應(yīng)變情況下，最大應(yīng)力誤差僅為0.07 MPa.該模型具有較好的計(jì)算穩(wěn)定性與生物逼真度，能夠滿足人體頸部肌肉建模的需求.

關(guān)鍵詞：生物力學(xué)；肌肉主動力；有限元分析；準(zhǔn)線性黏彈性材料；本構(gòu)模型；耦合模型

中圖分類號：R318.01

人體頸部肌肉的主被動作用在頭頸部動力學(xué)響應(yīng)中具有重要的意義.低速碰撞環(huán)境中，駕駛員由于視覺或觸覺引發(fā)身體應(yīng)激反應(yīng)，使得肌肉緊張產(chǎn)生主動力.肌肉主動力對頭頸部運(yùn)動的穩(wěn)定性乃至減小碰撞中頭頸部損傷具有積極的作用[1-4].因而在研究車輛碰撞特別是后碰撞中的頭頸部損傷時(shí)，有必要考慮頸部肌肉主動力的影響.

有限元方法被廣泛應(yīng)用于被動安全與損傷生物力學(xué)領(lǐng)域，常見的模擬肌肉主動力的方法是采用單軸的Hill肌肉本構(gòu)模型[5-8].該方法最早將肌肉簡化為起點(diǎn)與止點(diǎn)之間的一維直線段，并不考慮幾何以及肌肉間交互作用的影響.Wittek[9]以串聯(lián)的Hill單元分段線性地?cái)M合肌肉走向與收縮力作用方向，然而串聯(lián)的Hill單元導(dǎo)致了不穩(wěn)定問題.基于連續(xù)介質(zhì)的肌肉本構(gòu)模型研究中，大多提出以復(fù)合材料的形式將肌肉的整體響應(yīng)分解為線性疊加的主動與被動本構(gòu)模型[10-13].主動本構(gòu)模型方面，Johansson等[14]最早將單軸的Hill模型修改成遞歸公式以應(yīng)用于復(fù)合本構(gòu)模型，隨后的多數(shù)研究都沿用了該方式[10-12].被動本構(gòu)模型往往以應(yīng)變能密度函數(shù)的形式，分別定義體積應(yīng)變響應(yīng)、等容超彈性響應(yīng)和等容黏彈性響應(yīng)并疊加得到被動響應(yīng).通過對各分量定義不同的應(yīng)變能函數(shù)，形成了多樣化的被動本構(gòu)模型.實(shí)際應(yīng)用中，Behr[15]以及Hedenstierna等[16]依據(jù)Wittek＼[9＼]的建議提出了在結(jié)構(gòu)而非本構(gòu)模型層面耦合肌肉主被動響應(yīng)的方法，并分別應(yīng)用于人體下肢與頸部有限元研究.該方法采用共節(jié)點(diǎn)的方式以實(shí)體單元抑制串聯(lián)Hill單元引起的不穩(wěn)定，是一種簡單、高效的建模方式.然而Behr等＼[15＼]并未考慮肌肉組織的應(yīng)變率效應(yīng)，且為了控制模型的不穩(wěn)定，采用了平緩的激活與加載曲線.這既降低了模型的逼真度也限制了其應(yīng)用范圍.Hedenstierna等＼[16＼]分別采用Ogden模型（超彈性橡膠材料）與QLV模型評估了該種模型在多種應(yīng)變率下的性能并維持了高應(yīng)變率下模型的穩(wěn)定性.結(jié)果表明QLV模型的精度更高，穩(wěn)定性卻較差.但是為了控制激活態(tài)下模型的穩(wěn)定性Hedenstierna調(diào)高了較大應(yīng)變（10%以上）時(shí)的應(yīng)力輸出.

國內(nèi)關(guān)于肌肉的有限元研究較少.多數(shù)生物力學(xué)有限元模型中，肌肉作用直接以離散載荷施加于模型中[17-18].在動態(tài)研究中，肌肉同樣被化簡為一維單元，且往往僅保留了被動響應(yīng)[19].部分模型采用實(shí)體單元模擬肌肉，但仍未添加主動響應(yīng)[20].本構(gòu)模型研究方面，龔亞琦等[21]將Johansson等[14]與Tang等[12]的研究相結(jié)合提出了包含疲勞因素的連續(xù)介質(zhì)本構(gòu)模型.

本文目的在于研究一種以實(shí)體單元與梁單元共節(jié)點(diǎn)方式耦合主被動響應(yīng)的肌肉有限元建模方法，并分析討論模型的力學(xué)特性與穩(wěn)定性特征.通過調(diào)整主動與被動本構(gòu)模型的參數(shù)改善模型的生物逼真度并提高其在被動與激活態(tài)時(shí)的計(jì)算穩(wěn)定性，以滿足不同應(yīng)變率下肌肉生物力學(xué)響應(yīng)的模擬仿真要求.

1 方法與材料

建立了兔子脛前?。═ibialis Anterior， TA）幾何模型.將準(zhǔn)線性黏彈性（QLV）本構(gòu)模型與Hill主動本構(gòu)模型耦合，建立了TA的肌肉主動力模型.根據(jù)Myers實(shí)驗(yàn)驗(yàn)證了該模型在不同應(yīng)變率下肌肉離心收縮的應(yīng)力應(yīng)變特性，同時(shí)研究了等長收縮與軸向壓縮模式下的肌肉的力學(xué)特性.

1.1 肌肉幾何模型

脛前肌的幾何模型根據(jù)Myers等[22-23]實(shí)驗(yàn)樣本的橫截面輪廓重建而得.脛前肌的總長度與等效生理橫截面積（Physiologic Cross Sectional Area， PCSA）分別為90 mm與45 mm2.兩端肌腱長14 mm＼[11＼]，起點(diǎn)與止點(diǎn)的橫截面積分別為21 mm2與12.5 mm2.兔子脛前肌有限元模型如圖1所示.

1.2 肌肉材料模型

材料實(shí)驗(yàn)研究表明[24-25]，肌肉的被動響應(yīng)具有高度非線性的黏彈性.因而實(shí)體單元材料采用了LSDYNA中MAT_176 （QLV Model）模型[26].QLV模型的黏彈性響應(yīng)是根據(jù)應(yīng)變歷史對非線性的瞬時(shí)彈性響應(yīng)進(jìn)行遺傳積分而得，其瞬時(shí)彈性響應(yīng)與黏彈性應(yīng)力為：

1.3 模型驗(yàn)證與工況模擬

Myers等人[22-23]為了研究骨骼肌對人體頸部在沖擊環(huán)境中動力學(xué)響應(yīng)的影響，以兔子脛前肌為樣本進(jìn)行了多種應(yīng)變率（1/s，10/s，25/s）的活體單軸拉伸實(shí)驗(yàn).Myers量化了兔子脛前肌在不同應(yīng)變率下的主、被動響應(yīng)，并認(rèn)為其可用于人體頸部動力學(xué)的數(shù)值研究.該實(shí)驗(yàn)因其豐富的數(shù)據(jù)與適合的實(shí)驗(yàn)條件設(shè)置而經(jīng)常為人體頸部數(shù)值研究所引用[1-2，16].

建立的有限元模型進(jìn)行離心收縮、等長收縮和軸向壓縮共5組9次仿真.仿真邊界條件設(shè)置見表2.除D組以外，所有仿真都分別在肌肉未激活與完全激活2種狀態(tài)下進(jìn)行加載（未激活與完全激活狀態(tài)分別以0和1表示，A0即1/s應(yīng)變率下未激活時(shí)的離心收縮仿真，其余以此類推）.A，B，C 3組仿真分別以3種應(yīng)變率加載，并與Myers的實(shí)驗(yàn)結(jié)果進(jìn)行對比以驗(yàn)證模型的有效性；D，E 2組仿真僅研究模型在不同載荷下的力學(xué)特征與穩(wěn)定性.

2 結(jié) 果

仿真結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的對比如圖3所示.圖中3類不同走向的曲線分別代表了TA在未激活、完全激活狀態(tài)下的工程應(yīng)力曲線（被動響應(yīng)與整體響應(yīng)）以及兩者的差值（主動響應(yīng)）.

圖3中分別給出了Myers實(shí)驗(yàn)數(shù)據(jù)、TA模型仿真結(jié)果以及Hedenstierna的研究結(jié)果（分別以Myers、仿真曲線以及Hedenstierna標(biāo)記），為簡化起見，僅給出了Hedenstierna的主動響應(yīng)以及25/s應(yīng)變率下的被動響應(yīng)曲線.由圖可知，未激活時(shí)TA模型的計(jì)算結(jié)果在加載初期存在一定偏差.在B0與C0組仿真中該偏差的峰值分別為0.07 MPa與0.06 MPa，A0組的誤差較小.C0組仿真結(jié)果對實(shí)驗(yàn)數(shù)據(jù)擬合較好，未出現(xiàn)Hedenstierna研究中應(yīng)力偏低以及不穩(wěn)定的狀況.當(dāng)完全激活時(shí)，B1與C1組仿真的整體響應(yīng)最大誤差皆約為0.07 MPa，對比B0與C0組的結(jié)果可知該誤差主要源自于被動響應(yīng).A1組仿真中，較大應(yīng)變下的主動響應(yīng)出現(xiàn)了下滑趨勢，25%應(yīng)變時(shí)誤差為0.05 MPa.等長收縮仿真中，在初始激活時(shí)刻模型內(nèi)產(chǎn)生了高頻振動的速度場，然而其幅值非常小，軸向最大振動速度約為0.3 m/s，且在10 ms內(nèi)迅速衰減至0.05 m/s.肌腹處的截面輸出力也很穩(wěn)定.當(dāng)軸向壓縮時(shí)，模型激活與否的響應(yīng)差距較大，如圖4所示.被動的模型在壓縮加載的早期便屈曲失穩(wěn)，呈C形狀，而激活的模型在工程應(yīng)變達(dá)到20%時(shí)仍能保持均勻穩(wěn)定的壓縮狀態(tài).

3 討 論

本文目的在于研究以共節(jié)點(diǎn)方式耦合肌肉主被動響應(yīng)的相關(guān)問題.所建模型從應(yīng)力預(yù)測的精確性與模型穩(wěn)定性兩方面研究.TA模型的離心收縮仿真結(jié)果表明，QLV模型能夠較好地模擬肌肉組織的拉伸特性.這與Hedenstierna＼[16＼]的結(jié)論一致.在未激活狀態(tài)的仿真中，QLV模型精確地預(yù)測了3種應(yīng)變率下肌肉組織的應(yīng)力曲線，且C0組仿真并未出現(xiàn)Hedenstierna＼[16＼]的研究中明顯的應(yīng)力低估＼[圖3（c）中ε＼]以及不穩(wěn)定現(xiàn)象.這是由于Hedenstierna＼[16＼]確定的黏性參數(shù)是以Best[33]的松弛實(shí)驗(yàn)調(diào)整而得，本文是依據(jù)模型本構(gòu)方程擬合獲得.此外Hedenstierna＼[16＼]以Davis＼[27＼]測定的應(yīng)力曲線為目標(biāo)擬合了瞬時(shí)彈性響應(yīng)參數(shù)，這既限定了QLV模型黏性響應(yīng)的變動空間，也與Myers的實(shí)驗(yàn)數(shù)據(jù)存在不兼容.B0與C0仿真中，加載初期較大的偏差是源自于以階躍函數(shù)設(shè)置的等應(yīng)變率速度邊界條件.TA模型在由靜止轉(zhuǎn)變?yōu)楦咚倮斓倪\(yùn)動過程中產(chǎn)生了附加動載荷，計(jì)算應(yīng)力時(shí)并未對該波動進(jìn)行處理.完全激活時(shí)，Hill模型也能精確地描述脛前肌的主動收縮特性.Hedenstierna＼[16＼]在研究中為了控制模型穩(wěn)定性，人為調(diào)高在較大應(yīng)變時(shí)主動響應(yīng)（圖3中δ），導(dǎo)致整體響應(yīng)偏高約15%.本文避免了該問題.B1與C1仿真在加載初期的應(yīng)力誤差除了繼承于被動響應(yīng)的部分，也受到了主動響應(yīng)的影響.該部分誤差是由于主動單元在高速伸長時(shí)主動收縮力受到Hill模型速度曲線的放大.A1仿真中應(yīng)變超過10%時(shí)主動響應(yīng)的下滑主要是由于TA模型中肌肉與肌腱接合部分的剛度過渡較大，近似于準(zhǔn)靜態(tài)的載荷條件使得接合處的主動單元未被充分拉伸，而中央肌腹部分的主動單元被過度拉伸，在長度曲線的影響下，造成了主動收縮力的下滑.

與Ogden模型在永久的彈性響應(yīng)上疊加線性黏性應(yīng)力不同，QLV模型是對非線性的瞬時(shí)彈性響應(yīng)遺傳積分得到黏性應(yīng)力，所以通過改變瞬時(shí)彈性響應(yīng)，QLV模型能夠更廣泛地?cái)M合非線性黏性特性.然而，由于缺乏穩(wěn)定的彈性響應(yīng)，基于QLV模型建立的肌肉有限元模型存在較大的不穩(wěn)定風(fēng)險(xiǎn).這些不穩(wěn)定一方面是源自于Hill模型的特性，另一方面又與QLV模型迅速衰減的松弛特性相關(guān).為了擬合肌肉組織的高度非線性拉伸特性，QLV模型的應(yīng)力水平在松弛階段最初的幾十毫秒內(nèi)將迅速衰減到峰值的40%，并且繼續(xù)隨應(yīng)變的變動而變動，這種迅速松弛的特征很容易破壞模型的平衡狀態(tài).同時(shí)，串聯(lián)的Hill單元會導(dǎo)致不穩(wěn)定.按照不穩(wěn)定性的來源可以將其分為長度與速度不穩(wěn)定.根據(jù)長度曲線，串聯(lián)的Hill單元在均勻伸長時(shí)，節(jié)點(diǎn)位置的擾動將導(dǎo)致相鄰Hill單元產(chǎn)生不同的收縮力，節(jié)點(diǎn)力的不平衡將對擾動形成正反饋促使其發(fā)散.長度曲線的斜率越大，模型越不穩(wěn)定.同理，模型長度不變時(shí)，節(jié)點(diǎn)的速度擾動將導(dǎo)致節(jié)點(diǎn)力不平衡，而該不平衡將對擾動形成負(fù)反饋導(dǎo)致高頻振動.兩類不穩(wěn)定嚴(yán)重限制了模型的應(yīng)用.等長收縮時(shí)觀察到TA模型的高頻振動即速度不穩(wěn)定.表現(xiàn)為位移場與截面力的高頻振動.在離心收縮以及較大應(yīng)變的等長收縮時(shí)，可觀察到沙漏狀不均勻網(wǎng)格變形，即長度不穩(wěn)定.該類不穩(wěn)定往往從某些易發(fā)區(qū)域產(chǎn)生并迅速擴(kuò)展到整個(gè)模型.不穩(wěn)定的擴(kuò)散過程與肌肉的不可壓縮性相關(guān)：某Hill單元的發(fā)散性伸長將強(qiáng)制其依附的實(shí)體單元棱邊的相鄰或相對棱邊收縮，形成連鎖效應(yīng).模型的速度不穩(wěn)定可以通過適當(dāng)提高Vmax，增加并聯(lián)以及減少串聯(lián)的Hill單元個(gè)數(shù)3種方法共同控制.而長度不穩(wěn)定則需要調(diào)整Hill單元特性分布，如圖2所示，沿軸向由兩端往中央逐步增加長度特性的斜率進(jìn)行控制.此外，將QLV的部分剛度轉(zhuǎn)移至Hill單元，由Hill單元的非線性彈簧提供彈性應(yīng)力以及避免過高的Vmax等都有助于提高串聯(lián)Hill單元的穩(wěn)定性.由圖3可知以上穩(wěn)定性控制措施對模型的力學(xué)響應(yīng)影響較小.經(jīng)調(diào)整后的TA模型在離心收縮、等長收縮以及軸向壓縮等工況下都能避免不穩(wěn)定問題，輸出穩(wěn)定的結(jié)果.

本文中兔子脛前肌的幾何模型是通過文獻(xiàn)中的橫截面積與長度重建而成，存在一定的幾何誤差.此外，由于目的是研究在沖擊環(huán)境下頸部肌肉的力學(xué)響應(yīng)，所以在擬合模型參數(shù)時(shí)，未考慮肌肉組織在長時(shí)間范圍（100 ms以上）內(nèi)的松弛特性以及壓縮特性.

4 結(jié) 論

QLV本構(gòu)模型模擬肌肉組織的被動拉伸特性具有較好的生物逼真性.通過分別調(diào)整串聯(lián)與并聯(lián)的Hill單元個(gè)數(shù)以及調(diào)整Hill單元長度特性曲線斜率的分布可以在不影響仿真精度的情況下提高激活態(tài)下主被動耦合模型的穩(wěn)定性.本文建立的肌肉有限元主被動耦合模型可以應(yīng)用于汽車碰撞中的頭頸部生物力學(xué)研究.

參考文獻(xiàn)

[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] 魏高峰. 人體骨肌系統(tǒng)的整體生物力學(xué)建模與仿真分析研究[D]. 上海交通大學(xué)機(jī)械與動力工程學(xué)院， 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] 仲照明， 宮赫， 肖智韜，等. 髖臼力和肌肉力對股骨有限元分析結(jié)果的影響[J]. 哈爾濱工業(yè)大學(xué)學(xué)報(bào)， 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ì)匡， 姚劍峰. 人體頸部動力學(xué)響應(yīng)分析有限元模型的建立和驗(yàn)證[J]. 湖南大學(xué)學(xué)報(bào)：自然科學(xué)版， 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]. 第三軍醫(yī)大學(xué)學(xué)報(bào)， 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] 龔亞琦， 樊建平. 基于纖維增強(qiáng)超彈性基體模型的骨骼肌有限元模擬[J]. 中國生物醫(yī)學(xué)工程學(xué)報(bào)， 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.

[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.

[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.