中國(guó)50th人體小腿有限元模型多工況優(yōu)化建模方法

張冠軍 王龍亮 鄧小朋 杜現(xiàn)平 官鳳嬌

1(湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410082)2(國(guó)防科學(xué)技術(shù)大學(xué)裝備綜合保障技術(shù)重點(diǎn)實(shí)驗(yàn)室, 長(zhǎng)沙 410073)

中國(guó)50th人體小腿有限元模型多工況優(yōu)化建模方法

張冠軍1王龍亮1鄧小朋1杜現(xiàn)平1官鳳嬌2*

1(湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410082)2(國(guó)防科學(xué)技術(shù)大學(xué)裝備綜合保障技術(shù)重點(diǎn)實(shí)驗(yàn)室, 長(zhǎng)沙 410073)

人體下肢是交通事故中最易受傷的部位之一，下肢有限元模型則成為研究下肢損傷機(jī)理和防護(hù)方法的重要工具，而要保證模型生物逼真度則需開(kāi)展多工況全面驗(yàn)證。人體模型驗(yàn)證所需的生物力學(xué)試驗(yàn)數(shù)據(jù)由于樣本尺寸、材料等多樣性造成力學(xué)響應(yīng)存在差異，難于用一個(gè)模型同時(shí)滿足多個(gè)試驗(yàn)數(shù)據(jù)。中國(guó)人體與歐美人體在幾何尺寸上的差異會(huì)造成生物力學(xué)響應(yīng)的不同，因此有必要開(kāi)發(fā)中國(guó)人體的有限元模型。通過(guò)CT、MRI掃描建立小腿幾何模型，利用由脛骨關(guān)鍵尺寸確定的縮放系數(shù)，將小腿幾何模型縮放至中國(guó)50百分位人體小腿。同樣采用縮放方法，將生物力學(xué)試驗(yàn)數(shù)據(jù)縮放至小腿尺寸所對(duì)應(yīng)的生物力學(xué)響應(yīng)數(shù)據(jù)，以解決試驗(yàn)樣本幾何差異的影響。以脛骨和腓骨的彈性模量、應(yīng)力-應(yīng)變曲線、失效應(yīng)變以及肉體體積模量為設(shè)計(jì)變量，采用優(yōu)化方法，同時(shí)擬合脛骨和腓骨在準(zhǔn)靜態(tài)、動(dòng)態(tài)載荷下不同加載位置的生物力學(xué)響應(yīng)，以及小腿動(dòng)態(tài)載荷下不同加載位置的生物力學(xué)響應(yīng)，解決試驗(yàn)樣本多樣性造成的模型驗(yàn)證難題。當(dāng)脛骨和腓骨皮質(zhì)骨彈性模量分別取18.43和18.23 GPa、失效應(yīng)變分別取1.2%和0.8%、小腿肉體體積模量取11.33 MPa時(shí)，所建小腿模型能夠較好地同時(shí)滿足多組試驗(yàn)中樣本的生物力學(xué)響應(yīng)，表明基于優(yōu)化的多工況人體模型驗(yàn)證方法能夠使有限元模型獲得較好的生物逼真度。

小腿；中國(guó)人體；優(yōu)化；有限元模型；驗(yàn)證

引言

交通安全問(wèn)題是一個(gè)世界性難題,全球每年超過(guò)117萬(wàn)人在交通事故中喪生，受傷者更是達(dá)到1 000萬(wàn)人[1]。下肢是交通事故中最易受到損傷的部位之一[2]，包括膝關(guān)節(jié)損傷、長(zhǎng)骨骨折等[3]。下肢損傷不僅給受害者造成巨大傷痛，也會(huì)給社會(huì)和家庭帶來(lái)沉重負(fù)擔(dān)。為更好地了解下肢損傷的機(jī)理，許多研究者采用尸體實(shí)驗(yàn)來(lái)研究下肢損傷，但由于實(shí)驗(yàn)成本高、樣本難以獲得且重復(fù)性差等原因，人體有限元模型獲得了廣泛應(yīng)用。

在國(guó)外，下肢有限元模型發(fā)展較早。1993年，Bermond等就開(kāi)發(fā)了包括股骨、脛骨和主要韌帶的行人下肢模型，但幾何外形不準(zhǔn)確[4]。1998年，Wykowski等開(kāi)發(fā)了乘員下肢有限元模型，只有骨骼并且將骨骼定義為剛體[5]。2005年，Untaroiu等基于LS-DYNA求解器，開(kāi)發(fā)了行人下肢有限元模型，定義了肌腱、肌肉和皮膚，并進(jìn)行了較詳細(xì)的驗(yàn)證[6]。2013年，Untaroiu等基于先前的模型，更詳細(xì)地提取骨骼幾何外形，建立了乘員下肢有限元模型[7]。在國(guó)內(nèi)，下肢有限元模型研究起步比較晚，2000年后一些學(xué)者根據(jù)國(guó)外模型進(jìn)行了材料改進(jìn)和驗(yàn)證[8-11]，并沒(méi)有根據(jù)中國(guó)人體的幾何外形建立中國(guó)人體下肢有限元模型。目前的下肢有限元模型大多都是根據(jù)歐美人體建立的，而中國(guó)人體在尺寸上與歐美人體存在較大差異，其生物力學(xué)響應(yīng)與中國(guó)人體存在多大差異目前并不明確。

不同組織層次的全面驗(yàn)證是保證人體有限元模型生物逼真度的有效手段。但由于試驗(yàn)樣本在幾何尺寸、組織材料等方面的差異，造成不同樣本的生物力學(xué)響應(yīng)存在較大差異，這給模型驗(yàn)證帶來(lái)了較大的挑戰(zhàn)，尤其是要同時(shí)滿足不同組織層次、不同載荷工況的多目標(biāo)模型驗(yàn)證。

因此，為建立具有較高生物逼真度且良好適用性的中國(guó)50th人體有限元模型，本研究利用醫(yī)學(xué)影像數(shù)據(jù)建立了中國(guó)50th人體小腿有限元模型，基于優(yōu)化方法構(gòu)建了模型驗(yàn)證方法，對(duì)不同組織層次、不同加載方向、不同加載位置等多種工況下的小腿模型進(jìn)行了驗(yàn)證，使一個(gè)模型同時(shí)與多種試驗(yàn)數(shù)據(jù)吻合。

1 小腿有限元模型建立

1.1 幾何模型

精確的幾何模型是獲得高質(zhì)量有限元模型的基礎(chǔ)。利用CT和MRI醫(yī)學(xué)影像分別構(gòu)建小腿長(zhǎng)骨(脛骨和腓骨)和肉體幾何，如圖1所示。醫(yī)學(xué)影像數(shù)據(jù)來(lái)源于一位骨骼正常的接近中國(guó)50th的男性血管疾病患者(身高173.1 cm，體重69.7 kg)。

圖1 小腿醫(yī)學(xué)影像與脛骨、腓骨和肌肉的幾何模型Fig.1 The medical imaging and geometry of the shank

1.2 有限元模型

綜合使用ANSYSY ICEM CFD和Hypermesh對(duì)小腿模型進(jìn)行網(wǎng)格劃分。脛骨體皮質(zhì)骨采用4層六面體網(wǎng)格模擬；脛骨兩端皮質(zhì)骨較薄，用殼單元模擬，松質(zhì)骨采用體單元模擬。脛骨體與脛骨兩端皮質(zhì)骨過(guò)渡區(qū)域選用階梯形狀逐漸過(guò)渡，以減小應(yīng)力集中，如圖2所示。腓骨骨干皮質(zhì)骨采用2層六面體單元模擬，兩端皮質(zhì)骨采用厚度為1.5 mm的殼單元模擬，骨干皮質(zhì)骨與兩端皮質(zhì)骨的過(guò)渡同樣采用階梯形狀。小腿的肉體采用實(shí)體單元模擬，用CONTACT_TIED_NODES_TO_ SURFACE將肌肉與脛骨、腓骨連接。皮膚采用殼單元模擬，使用共節(jié)點(diǎn)方式與肌肉連接。

圖2 小腿網(wǎng)格的劃分Fig.2 The mesh division of legs

脛骨、腓骨模型雅克比小于0.60的單元比例不超過(guò)1%，最小值為0.42；翹曲度大于20°的單元比例不超過(guò)2%，最大值為120°；長(zhǎng)寬比大于3.5的單元比例小于5%，最大值為5.2；最小單元尺寸4.69 mm，最大單元尺寸6.42 mm。

1.3 有限元模型尺寸縮放

由于本模型的人體尺寸(身高173.1 cm，體重69.7 kg)大于中國(guó)50百分位人體尺寸[12](身高 170.8 cm，體重65 kg)，需要將小腿模型縮放到中國(guó)50百分位人體小腿尺寸。在碰撞過(guò)程中，脛骨是決定小腿生物力學(xué)響應(yīng)的主要部分，因此本研究依據(jù)脛骨尺寸對(duì)小腿模型進(jìn)行縮放。

軸向縮放系數(shù)主要依據(jù)中國(guó)50百分位人體與本樣本脛骨長(zhǎng)的比例，橫向縮放系數(shù)主要依據(jù)中國(guó)50百分位人體與本樣本股骨上內(nèi)關(guān)節(jié)面矢徑和上外關(guān)節(jié)面矢徑比例的平均值，縱向縮放系數(shù)主要依據(jù)中國(guó)50百分位人體與本樣本的上段寬和下段寬比例的平均值[12]。最終確定脛骨的軸向縮放系數(shù)為0.965，橫向縮放系數(shù)為0.969，縱向縮放系數(shù)為0.964，如表1所示。將脛骨模型按照3個(gè)方向縮放系數(shù)平均值縮放至中國(guó)50百分位男性脛骨的尺寸，并依據(jù)此比例對(duì)小腿腓骨、肉體和皮膚模型進(jìn)行縮放，最終獲得中國(guó)50百分位男性小腿有限元模型。

1.4 小腿材料模型

脛骨、腓骨的皮質(zhì)骨采用可分別定義拉、壓應(yīng)力-應(yīng)變曲線的彈塑性材料，應(yīng)力-應(yīng)變曲線初始值如圖3所示[13]。采用黏彈塑性材料模擬松質(zhì)骨。皮膚采用彈性材料模擬，密度為1 000 kg/m3，彈性模量為1 MPa，泊松比為0.3，厚度為1 mm[14-15]。肉體使用黏彈性材料模擬，密度為1 000 kg/m3，體積模量為20 MPa，超彈性系數(shù)C1為1.2×10-4MPa，超彈性系數(shù)C2為2.5×10-4MPa，普羅尼算法系數(shù)S1為1.162，普羅尼算法系數(shù)S2為0.808，普羅尼算法系數(shù)T1為10.43 ms，普羅尼算法系數(shù)T2為84.1 ms[16]。采用失效應(yīng)變作為失效準(zhǔn)則來(lái)模擬損傷，在材料達(dá)到設(shè)定的失效應(yīng)變后自動(dòng)刪除失效單元。根據(jù)相關(guān)文獻(xiàn)設(shè)定模型優(yōu)化前各參數(shù)取值，如表2所示[14, 17-24]。

表2 小腿模型材料參數(shù)設(shè)置Tab.2 Material parameters of leg model

2 小腿有限元模型的優(yōu)化驗(yàn)證

充分全面的驗(yàn)證是保證有限元模型適應(yīng)性和逼真度的基礎(chǔ)，本研究進(jìn)行的模型驗(yàn)證包括準(zhǔn)靜態(tài)驗(yàn)證和動(dòng)態(tài)驗(yàn)證，如表3所示。傳統(tǒng)的驗(yàn)證方式采用試錯(cuò)法手動(dòng)調(diào)節(jié)材料參數(shù)對(duì)模型進(jìn)行驗(yàn)證，效率低且準(zhǔn)確度有限。本研究采用優(yōu)化方法自動(dòng)獲得最佳材料參數(shù)值，能同時(shí)保證模型的生物力學(xué)響應(yīng)與多個(gè)試驗(yàn)結(jié)果吻合。

圖4 小腿模型驗(yàn)證流程Fig.4 Flow chart of the shank model validation

先對(duì)脛骨近心端、中部、遠(yuǎn)心端3個(gè)加載位置的動(dòng)態(tài)仿真進(jìn)行脛骨體皮質(zhì)骨的彈性模量和拉、壓應(yīng)力-應(yīng)變曲線優(yōu)化(此時(shí)未定義失效應(yīng)變)，將得到的最佳材料參數(shù)更新到脛骨模型，然后對(duì)脛骨近心端、中部、遠(yuǎn)心端3個(gè)加載位置的動(dòng)態(tài)仿真進(jìn)行脛骨失效應(yīng)變優(yōu)化。同樣,對(duì)腓骨進(jìn)行類似的優(yōu)化驗(yàn)證。接著將得到的最佳材料參數(shù)更新到小腿模型，對(duì)小腿近心端、中部、遠(yuǎn)心端3個(gè)加載位置的動(dòng)態(tài)仿真進(jìn)行肉體材料的體積模量?jī)?yōu)化。最終獲取一組最佳材料參數(shù)，再對(duì)脛骨和腓骨的準(zhǔn)靜態(tài)試驗(yàn)進(jìn)行校核。優(yōu)化流程如圖4所示，設(shè)計(jì)變量的取值范圍如表4所示[7, 14, 19, 30-32]。

表3 小腿有限元模型驗(yàn)證Tab.3 Finite element model validation of the shank

注：A-P為載荷由前向后方向，L-M為載荷由外側(cè)向內(nèi)側(cè)。

Note：A-P: loading direction from anterior to posterior；L-M: loading direction from lateral to medial.

表4 設(shè)計(jì)變量的取值范圍Tab.4 The range of design variables

在脛骨和腓骨體皮質(zhì)骨的彈性模量和拉、壓應(yīng)力-應(yīng)變曲線優(yōu)化中,設(shè)計(jì)變量分別彈性模量和拉、壓應(yīng)力-應(yīng)變曲線橫(應(yīng)變)、縱(應(yīng)力)坐標(biāo)的縮放系數(shù)。在肌肉材料參數(shù)優(yōu)化中，設(shè)計(jì)變量為體積模量。目標(biāo)函數(shù)為仿真與實(shí)驗(yàn)的力-位移曲線的均方差f(X)的最小值[33]，即

(1)

式中：X為設(shè)計(jì)變量，k=1,2, …，K，K為實(shí)驗(yàn)曲線的條數(shù)；m=1,2, …，Pk，Pk為第k條曲線中計(jì)算點(diǎn)的個(gè)數(shù)；Wm為權(quán)重系數(shù)，本研究中各工況權(quán)重為1；fm(X)為響應(yīng)面近似模型的計(jì)算值；Gm為實(shí)驗(yàn)測(cè)試點(diǎn)的值。

在脛骨和腓骨體皮質(zhì)骨的失效應(yīng)變優(yōu)化中，設(shè)計(jì)變量為失效應(yīng)變。目標(biāo)函數(shù)為各工況仿真與實(shí)驗(yàn)的力-位移曲線的力極大值點(diǎn)對(duì)應(yīng)的位移之差絕對(duì)值f(X)的最小值，即

(2)

式中：n=1,2,…，N，N為工況的個(gè)數(shù)；Xn為第n個(gè)工況中實(shí)驗(yàn)的力-位移曲線的位移最大值；Xsn為第n個(gè)工況中仿真的力-位移曲線的力極大值對(duì)應(yīng)的位移。

在優(yōu)化中，如果某個(gè)工況下的目標(biāo)實(shí)驗(yàn)曲線有多個(gè)，則采用這些實(shí)驗(yàn)曲線的平均曲線作為該工況優(yōu)化的目標(biāo)曲線。平均曲線的計(jì)算方法參照Lessley等的計(jì)算(2004)[34]。

2.1 生物力學(xué)實(shí)驗(yàn)數(shù)據(jù)縮放

脛骨和腓骨動(dòng)態(tài)仿真采用Takahashi等實(shí)驗(yàn)(2003)中脛骨長(zhǎng)度為467 mm的數(shù)據(jù)[35]，小腿動(dòng)態(tài)仿真采用Kerrigan(2008)縮放至美國(guó)50百分位人體的數(shù)據(jù)[36]。為了準(zhǔn)確驗(yàn)證中國(guó)50百分位人體的有限元模型，還需將實(shí)驗(yàn)數(shù)據(jù)縮放到與中國(guó)50百分位人體幾何相對(duì)應(yīng)的數(shù)值，以消除尺寸差異對(duì)生物力學(xué)響應(yīng)的影響。

Kerrigan(2008)為準(zhǔn)確獲得美國(guó)50百分位人體小腿的損傷耐受限度和彎曲響應(yīng)，采用脛骨長(zhǎng)度縮放系數(shù)λL計(jì)算位移縮放系數(shù)λD和力縮放系數(shù)λF，然后對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行縮放[2]，具體計(jì)算如下：

式中，L50th為美國(guó)50百分位人體的脛骨長(zhǎng)度，L樣本為試驗(yàn)樣本的脛骨長(zhǎng)度。

實(shí)驗(yàn)曲線中的位移、力、彎矩乘以相應(yīng)的縮放系數(shù)，即可得到與美國(guó)50百分位相對(duì)應(yīng)的曲線。按照上述實(shí)驗(yàn)數(shù)據(jù)縮放方法，只需將美國(guó)50百分位人體的脛骨長(zhǎng)度縮放至中國(guó)50百分位人體的脛骨長(zhǎng)度。

由于脛骨準(zhǔn)靜態(tài)試驗(yàn)所用脛骨長(zhǎng)度尺寸不能確定，本研究只對(duì)脛骨、腓骨和小腿動(dòng)態(tài)仿真的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行縮放。文獻(xiàn)中脛骨、腓骨和小腿動(dòng)態(tài)實(shí)驗(yàn)數(shù)據(jù)已經(jīng)過(guò)縮放，所以本研究在脛骨、腓骨樣本的近心端1/3處、中間、遠(yuǎn)心端1/3處的縮放系數(shù)相同，小腿樣本的近心端1/3處、中間、遠(yuǎn)心端1/3處的縮放系數(shù)相同，如表5所示。實(shí)驗(yàn)樣本的位移、力乘以相應(yīng)縮放系數(shù)即可得到與中國(guó)50百分位人體尺寸相對(duì)應(yīng)的縮放實(shí)驗(yàn)曲線。

表5 實(shí)驗(yàn)樣本的縮放系數(shù)和脛骨長(zhǎng)度

Tab.5 The scaling coefficients of the experimental samples and the length of the tibia

縮放參數(shù)脛骨L樣本小腿L樣本Prox.Mid.Dist.Prox.Mid.Dist.中國(guó)50thL50th脛骨長(zhǎng)/mm391.00391.00391.00378.70378.70378.70354.02λL0.9050.9050.9050.9350.9350.935—λD0.9050.9050.9050.9350.9350.935—λF0.8190.8190.8190.8740.8740.874—

2.2 脛骨、腓骨和小腿模型動(dòng)態(tài)模型

鑒于人體側(cè)面遭受撞擊的幾率遠(yuǎn)大于其他方向[36]，本研究?jī)H進(jìn)行L-M方向的動(dòng)態(tài)仿真。Kerrigan等(2003)和Kerrigan(2008)對(duì)脛骨、腓骨和小腿的近心端1/3處、中部、遠(yuǎn)心端1/3處進(jìn)行了三點(diǎn)彎曲實(shí)驗(yàn)[28, 36]。在脛骨和腓骨的實(shí)驗(yàn)中，沖擊塊被一層厚25 mm的ConforTM 泡沫包裹，沖擊速度為1.45 m/s，在小腿實(shí)驗(yàn)中，加載速度為1.5 m/s。Untaroiu等(2005)[2]和Takahashi等(2003)[16]使用其數(shù)據(jù)開(kāi)展了模型驗(yàn)證。本研究參照上述實(shí)驗(yàn)和仿真建立脛骨、腓骨和小腿的動(dòng)態(tài)三點(diǎn)彎曲驗(yàn)證模型，分別如圖5～7所示。金屬方盒、弧形板、支撐平板和沖擊塊均定義為剛體，泡沫定義為柔性體(Material #57)。長(zhǎng)骨與金屬盒的連接定義為剛性連接(*CONSTRAINED_EXTRA_NODES _SET)[2,16]。

圖5 脛骨動(dòng)態(tài)三點(diǎn)彎曲仿真模型Fig.5 Simulation of tibia dynamic 3-point bending test

圖7 小腿動(dòng)態(tài)三點(diǎn)彎曲仿真模型Fig.7 Simulation of shank dynamic 3-point bending test

圖6 腓骨動(dòng)態(tài)三點(diǎn)彎曲仿真模型
Fig.6 Simulation of fibula dynamic 3-point bending test

2.3 脛骨和腓骨模型準(zhǔn)靜態(tài)仿真模型

依據(jù)Yamada、Asang、Ehler相關(guān)文獻(xiàn)[25-27]建立仿真模型。用直徑為25 mm的剛性圓筒沖擊器，以0.01 m/s的速度對(duì)長(zhǎng)骨中部進(jìn)行加載。根據(jù)載荷加載方向的不同，脛骨和腓骨準(zhǔn)靜態(tài)三點(diǎn)彎曲仿真分為A-P和L-M兩個(gè)方向的驗(yàn)證。仿真設(shè)置如圖8、9所示。

圖8 脛骨準(zhǔn)靜態(tài)三點(diǎn)彎曲仿真模型。(a) A-P；(b) L-MFig.8 Quasi-static 3-point bending simulation of tibia. (a) A-P; (b) L-M

圖9 脛骨準(zhǔn)靜態(tài)三點(diǎn)彎曲仿真模型。(a) A-P；(b) L-MFig.9 Quasi-static 3-point bending simulation of fibula. (a) A-P; (b) L-M

3 模型仿真結(jié)果

3.1 優(yōu)化結(jié)果

圖10列出了典型的設(shè)計(jì)變量在優(yōu)化過(guò)程中的收斂歷程。隨著迭代次數(shù)的增加，各設(shè)計(jì)變量興趣域空間逐漸縮小，各參數(shù)取值逐漸穩(wěn)定。最終脛骨體皮質(zhì)骨彈性模量為18.43 GPa，拉、壓應(yīng)力-應(yīng)變曲線橫、縱坐標(biāo)的縮放系數(shù)分別為2.560和0.684，失效應(yīng)變?yōu)?.2%；腓骨體皮質(zhì)骨彈性模量為18.23 GPa，拉、壓應(yīng)力-應(yīng)變曲線橫、縱坐標(biāo)的縮放系數(shù)分別為2.888和0.725，失效應(yīng)變?yōu)?.8%；小腿肉體積模量為11.33 MPa。

圖10 材料參數(shù)收斂過(guò)程。(a)脛骨彈性模量；(b) 脛骨應(yīng)力-應(yīng)變曲線橫坐標(biāo)(應(yīng)變)縮放系數(shù)；(c) 脛骨應(yīng)力-應(yīng)變曲線縱坐標(biāo)(應(yīng)力)縮放系數(shù)；(d) 小腿體積模量Fig.10 Convergence process of material parameters. (a) Elastic modulus of the tibia; (b) Abscissa scaling factor of the tibial stress-strain curve; (c) Ordinate scaling factor of the tibial stress strain curve; (d) Bulk modulus of the shank

3.2 動(dòng)態(tài)仿真優(yōu)化結(jié)果

脛骨、腓骨和小腿材料參數(shù)優(yōu)化后，仿真結(jié)果與縮放后試驗(yàn)曲線的對(duì)比分別如圖11～13所示。

圖11 脛骨動(dòng)態(tài)三點(diǎn)彎曲優(yōu)化驗(yàn)證結(jié)果。(a)近心端1/3處；(b)中部；(c)遠(yuǎn)心端1/3處Fig.11 Dynamic 3-point bending validation results for tibia against test results with optimization method. (a) Proximal 1/3; (b) Middle-shaft; (c) Distal 1/3

圖12 腓骨動(dòng)態(tài)三點(diǎn)彎曲優(yōu)化驗(yàn)證結(jié)果。(a)近心端1/3處；(b)中部；(c)遠(yuǎn)心端1/3處Fig.12 Dynamic 3-point bending validation results for fibula against test results with optimization method. (a) Proximal 1/3; (b) Middle-shaft; (c) Distal 1/3

圖13 小腿動(dòng)態(tài)優(yōu)化驗(yàn)證與縮放后的試驗(yàn)曲線。(a)近心端1/3處；(b)中部；(c)遠(yuǎn)心端1/3處Fig.13 Dynamic 3-point bending validation results for shank against test results with optimization method. (a) Proximal 1/3; (b) Middle-shaft; (c) Distal 1/3

由圖11～圖13可知，脛骨動(dòng)態(tài)三點(diǎn)彎曲仿真的近心端1/3處和遠(yuǎn)心端1/3處力-位移曲線與目標(biāo)曲線吻合程度較好；脛骨中部三點(diǎn)彎曲仿真的力-位移曲線在上半段與目標(biāo)曲線擬合較好，雖然在下半段存在一些偏差，但曲線的走勢(shì)基本一致。腓骨近心端、中部、遠(yuǎn)心端3處加載位置的動(dòng)態(tài)仿真都與目標(biāo)曲線吻合良好。小腿近心端和中部加載位置的動(dòng)態(tài)仿真與目標(biāo)曲線吻合較好，但在遠(yuǎn)心端1/3加載位置的仿真結(jié)果不夠理想。

3.3 準(zhǔn)靜態(tài)仿真結(jié)果

使用優(yōu)化后的材料參數(shù)進(jìn)行脛骨和腓骨A-P、L-M兩個(gè)加載方向的準(zhǔn)靜態(tài)三點(diǎn)彎曲仿真。由于實(shí)驗(yàn)涉及的脛骨長(zhǎng)度無(wú)法確定，在準(zhǔn)靜態(tài)仿真時(shí)未對(duì)實(shí)驗(yàn)曲線進(jìn)行縮放。脛骨和腓骨的驗(yàn)證結(jié)果分別如圖14、15所示。仿真結(jié)果在實(shí)驗(yàn)曲線范圍內(nèi)，并與實(shí)驗(yàn)曲線保持了較好的一致性。模型發(fā)生損傷時(shí)的變形量和碰撞力也與實(shí)驗(yàn)曲線吻合較好。

圖14 脛骨準(zhǔn)靜態(tài)三點(diǎn)彎曲驗(yàn)證結(jié)果。(a) A-P；(b) L-MFig.14 Quasi-static 3-point bending validation results for tibia against test results with optimization method

圖15 腓骨準(zhǔn)靜態(tài)三點(diǎn)彎曲驗(yàn)證結(jié)果。(a) A-P；(b) L-MFig.15 Quasi-static 3-point bending validation results for fibula against test results with optimization method. (a) A-P; (b) L-M

4 討論

本研究采用優(yōu)化方法對(duì)所建立的中國(guó)50th男性人體小腿有限元模型進(jìn)行了全面的驗(yàn)證，股骨和脛骨的彈性模量、應(yīng)力-應(yīng)變曲線縮放系數(shù)、失效應(yīng)變和小腿肉的體積模量等設(shè)計(jì)變量均能夠有效收斂于最優(yōu)值。脛骨和腓骨皮質(zhì)骨的彈性模量分別為18.43和18.23 GPa，這與多個(gè)文獻(xiàn)給出的14.3～21.1GPa[19, 37-38]的范圍一致。脛骨和腓骨皮質(zhì)骨的失效應(yīng)變分別為1.2%和0.8%，這與文獻(xiàn)[19]基本一致。小腿肉的體積模量為11.33 MPa，也在文獻(xiàn)給出的1.33～29.2 MPa[16,39]的范圍內(nèi)。

腓骨近心端、中部和遠(yuǎn)心端3處不同加載位置的動(dòng)態(tài)仿真都與試驗(yàn)曲線吻合良好，A-P、L-M兩個(gè)加載方向的準(zhǔn)靜態(tài)仿真也與試驗(yàn)曲線吻合較好。脛骨近心端和遠(yuǎn)心端加載位置的動(dòng)態(tài)仿真以及A-P、L-M兩個(gè)加載方向的準(zhǔn)靜態(tài)仿真均與試驗(yàn)結(jié)果吻合良好，但中部加載位置的動(dòng)態(tài)仿真與試驗(yàn)結(jié)果有一定差異。小腿中部加載的動(dòng)態(tài)仿真與試驗(yàn)結(jié)果吻合較好，但近心端和遠(yuǎn)心端加載位置的動(dòng)態(tài)仿真與試驗(yàn)結(jié)果有一定差異。上述驗(yàn)證所涉及的生物力學(xué)試驗(yàn)對(duì)應(yīng)的試樣樣本較多，樣本年齡分布較廣，個(gè)體差異帶來(lái)較大的生物力學(xué)響應(yīng)的差異。因此，用一套材料參數(shù)同時(shí)滿足基于不同樣本的多組試驗(yàn)曲線存在較大的現(xiàn)實(shí)難度。通過(guò)多工況同時(shí)優(yōu)化并設(shè)置相同的權(quán)重，能夠最大限度地保證仿真結(jié)果的總體擬合效果。

由于缺少中國(guó)人體的生物力學(xué)試驗(yàn)數(shù)據(jù)，本研究開(kāi)展驗(yàn)證所用的試驗(yàn)數(shù)據(jù)均來(lái)源于歐美人體。但中國(guó)人體與歐美人體在解剖學(xué)結(jié)構(gòu)方面存在一定的差異，而基于縮放方法的試驗(yàn)數(shù)據(jù)縮放，能夠在一定程度上減小解剖學(xué)結(jié)構(gòu)差異對(duì)模型驗(yàn)證的影響。不管中國(guó)人體還是歐美人體均屬人類，兩者在力學(xué)特性及材料參數(shù)方面雖存在差異，但可能并不大。但是，將來(lái)仍需要針對(duì)中國(guó)人體開(kāi)展相應(yīng)的生物力學(xué)試驗(yàn)，以揭示與歐美人體在力學(xué)特性方面的差異，并為中國(guó)人體數(shù)學(xué)模型的研究提供基礎(chǔ)性數(shù)據(jù)支持。

5 結(jié)論

本研究通過(guò)CT和MRI醫(yī)學(xué)影像數(shù)據(jù)建立了小腿骨骼和肉體幾何模型，并根據(jù)人體統(tǒng)計(jì)學(xué)數(shù)據(jù)將小腿模型縮放到中國(guó)50百分位人體尺寸。使用優(yōu)化方法，對(duì)脛骨、腓骨和小腿不同加載位置的仿真進(jìn)行優(yōu)化，獲得了理想的材料參數(shù)，并能夠滿足準(zhǔn)靜態(tài)下脛骨和腓骨不同加載方向的生物力學(xué)逼真度。結(jié)果表明，所建立的中國(guó)人體50百分位小腿模型具有較高的生物逼真度，基于優(yōu)化和縮放方法的人體有限元模型驗(yàn)證方法能夠有效解決試驗(yàn)樣本幾何、材料多樣性帶來(lái)的多工況試驗(yàn)數(shù)據(jù)擬合問(wèn)題，保證模型具有較好的生物逼真度和較廣泛的適用性。

[1] Zhang Guanjun, Cao Libo, Hu Jingwen, et al. A field data analysis of risk factors affecting the injury risks in vehicle-to-pedestrian crashes [J]. Ann Adv Automot Med, 2008, 52: 199-214.

[2] Kim YS, Choi HH, Cho YN, et al. Numerical investigations of interactions between the knee-thigh-hip complex with vehicle interior structures [J]. Stapp Car Crash J, 2005, 49: 85-115.

[3] Yang Jikuang. Injury biomechanics in car-pedestrian collisions: development, validation and application of human-body mathematical models [D]. Gothenburg: Chalmers University of Technology, 1997.

[4] Bermond F, Ramet M, Bouquet R, et al. A finite element model of the pedestrian knee joint in lateral impact [C] // Proceedings of the International Research Council on the Biomechanics of Injury Conference. Eindhoven: INRETS-Bron, 1993: 117-129.

[5] Wykowski E, Sinnhuber R, Appel H. Finite element model of human lower extremities in a frontal impact [C]//Proceedings of the International Research Council on the Biomechanics of Injury Conference. Goteborg: IRCOBI Secretariat, 1998: 101-116.

[6] Untaroiu CD. Development and validation of a finite element model of human lower limb: including detailed geometry, physical material properties, and component validations for pedestrian injuries [D]. Charlottesville: University of Virginia, 2005.

[7] Untaroiu CD, Yue N, Shin J. A finite element model of the lower limb for simulating automotive impacts [J]. Annals of Biomedical Engineering, 2013, 41(3): 513-526.

[8] 楊濟(jì)匡, 方海峰. 人體下肢有限元?jiǎng)恿W(xué)分析模型的建立和驗(yàn)證 [J]. 湖南大學(xué)學(xué)報(bào)(自然科學(xué)版), 2005, 32(5): 31-36.

[9] 張冠軍. 行人下肢的碰撞損傷特性及相關(guān)參數(shù)研究 [D].長(zhǎng)沙: 湖南大學(xué), 2009.

[10] 李正東, 劉寧國(guó), 黃平, 等. 下肢有限元模型的建立及損傷機(jī)制重建 [J]. 中國(guó)司法鑒定, 2012(6): 37-42.

[11] 蔣小晴, 楊濟(jì)匡, 王丙雨, 等. 乘員股骨在軸向壓力-彎矩下的損傷生物力學(xué)機(jī)理研究 [J]. 力學(xué)學(xué)報(bào), 2013, 46(3): 465-474.

[12] 劉俊先, 張興和. 中國(guó)正常人體測(cè)量值 [M]. 北京: 中國(guó)醫(yī)藥科技出版社, 1994: 62-63.

[13] Novitskaya E, Chen PY, Hamed E, et al. Recent advances on the measurement and calculation of the elastic moduli of cortical and trabecular bone: a review [J]. Theoretical and Applied Mechanics, 2011, 38(3): 209-297.

[14] Beillas P, Begeman PC, Yang KH, et al. Lower limb: advanced FE model and new experimental data [J]. Stapp Car Crash J, 2001, 45: 469-494.

[15] Karimi A, Navidbakhsh M. Material properties in unconfined compression of gelatin hydrogel for skin tissue engineering applications [J]. Biomedical Engineering/ Biomedizinische Technik, 2014, 59(6): 479-486.

[16] Snedeker JG, Muser MH, Walz FH. Assessment of pelvis and upper leg injury risk in car-pedestrian collisions: comparison of accident statistics, impactor tests and a human body finite element model [J]. Stapp Car Crash J, 2003, 47: 437-458.

[17] Iwamoto M, Tamura A, Furusu K, et al. Development of a finite element model of the human lower extremity for analyses of automotive crash injuries [C] //SAE 2000 World Congress. Detroit: SAE International, 2000: 2000-01-0621.

[18] Schuster PJ, Chou CC, Prasad P, et al. Development and validation of a pedestrian lower limb non-linear 3-d finite element model [J]. Stapp Car Crash J, 2000, 44: 315-334.

[19] Takahashi Y, Kikuchi Y, Konosu A, et al. Development and validation of the finite element model for the human lower limb of pedestrians [J]. Stapp Car Crash J, 2000, 44: 335-355.

[20] Chawls A, Mukherjee S, Mohan D, et al. Validation of lower extremity model in THUMS [C] //Proceedings of the International Research Council on the Biomechanics of Injury conference. Graz: IRCOBI Secretariat, 2004: 155-166.

[21] Arnoux PJ, Cesari D, Behr M, et al. Pedestrian lower limb injury criteria evaluation: a finite element approach [J]. Traffic Injury Prevention, 2005, 6(3): 288-297.

[22] Untaroiu C, Darvish K, Crandall J, et al. A finite element model of the lower limb for simulating pedestrian impacts [J]. Stapp Car Crash J, 2005, 49: 157-181.

[23] Silvestri C, Ray M. Development of a finite element model of the knee-thigh-hip of a 50th percentile male including ligaments and muscles [J]. International Journal of Crashworthiness, 2009, 14(2): 215-229.

[24] Neale M, Thomas R, Bateman H, et al. A finite element modelling investigation of lower leg injuries [C] //Proceedings of the 20th International Technical Conference on the Enhanced Safety of Vehicles (ESV). Lyon: NHTSA, 2007: 07-0077.

[25] Yamada H. Strength of biological materials [M]. Baltimore: The Williams & Wilkins Company, 1970: 49-73.

[26] Ehler E, L?sche H. Die menschliche tibia unter biegebelastung [J]. Beitr. Orthop, 1970, 17(5): 291-304.

[27] Asang E. Experimental biomechanics of the human leg: a basis for interpreting typical skiing injury mechanisms [J]. The Orthopedic Clinics of North America, 1976, 7(1): 63-73.

[28] Kerrigan J, Ivansson B, Bose D, et al. Rate-sensitive constitutive and failure properties of human collateral knee ligaments [C]//Proceedings of the International Research Council on the Biomechanics of Injury conference. Lisbon: IRCOBI Secretariat, 2003: 177-190.

[29] Begeman P, Paravasthu N. Static and dynamic compression loading of the lower leg [C]//Proceedings of the Seventh Injury Prevention Through Biomechanics Symposium. Detroit: Wayne State University, 1997.

[30] Burstein AH, Reilly DT, Martens M. Aging of bone tissue: mechanical properties [J]. J Bone Joint Surg Am, 1976, 58(1): 82-86.

[31] Martens M, van Audekercke R, Delport P, et al. The mechanical characteristics of cancellous bone at the upper femoral region [J]. Journal of Biomechanics, 1983, 16(12): 971-983.

[32] Galbusera F, Freutel M, Dürselen L, et al. Material models and properties in the finite element analysis of knee ligaments: a literature review[J]. Frontiers in Bioengineering and Biotechnology, 2014, 2: 1-11.

[33] 官鳳嬌. 沖擊載荷下的生物組織材料參數(shù)反求及損傷研究 [D]. 長(zhǎng)沙: 湖南大學(xué), 2011.

[34] Lessley D, Crandall J, Shaw G, et al. A normalization technique for developing corridors from individual subject responses [C] //SAE 2000 World Congress. Detroit: SAE International, 2004: 2004-01-0288.

[35] Takahashi Y, Kikuchi Y, Mori F, et al. Advanced finite element lower limb model for pedestrians [C] //Proceedings of the 18th International Technical Conference on the Enhanced Safety of Vehicles (ESV). Lyon: NHTSA, 2003: No 218.

[36] Kerrigan JR. A computationally efficient mathematical model of the pedestrian lower extremity [D]. Charlottesville: University of Virginia, 2008.

[37] Rho JY, Ashman RB, Turner CH. Young's modulus of trabecular and cortical bone material: ultrasonic and microtensile measurements [J]. J Biomech, 1993, 26(2): 111-119.

[38] Snyder SM, Schneider E. Estimation of mechanical properties of cortical bone by computed tomography [J]. J Orthop Res, 1991, 9(3): 422-431.

[39] Ruan J, El-Jawahri R, Chai L, et al. Prediction and analysis of human thoracic impact response and injuries in cadaver impacts using a full human body finite element model [J]. Stapp Car Crash J, 2003, 47: 299-321.

Study on FE Modeling Method of the 50th Percentile Chinese Shank with Optimization Methods under Multi-Loading Condition

Zhang Guanjun1Wang Longliang1Deng Xiaopeng1Du Xianping1Guan Fengjiao2*

1(StateKeyLaboratoryofAdvancedDesignandManufacturingforVehicleBody,HunanUniversity,Changsha410082，China)2(CollegeofMechatronicEngineeringandAutomation,NationalUniversityofDefenseTechnology,Changsha410073,China)

The lower extremity of the human body is one of the most vulnerable parts in traffic accidents, and the finite element (FE) model of the lower limb has become an important tool to study mechanisms and protective methods of lower limb injuries. In order to ensure the fidelity of the model, it is necessary to carry out a comprehensive validation under multi-loading conditions. The biomechanical test data are different in the mechanical response due to the sample size and material diversity, so it is difficult to use one model to meet the multiple test data. The differences in geometric dimensions between Chinese and occidental human body would result in differences of biomechanical response, so it is necessary to develop the FE model of Chinese human body. In this paper, Geometry model based on the CT and MRI data was scaled to the 50 percentile male according to scaling factor derived from key dimensions of the tibia. The scaling method was also used to scale biomechanical test data to the biomechanical response data corresponding to the size of the leg. The elastic modulus, stress-strain curve, failure strain of the tibia and fibula and physical bulk modulus of the muscle were selected as design variables, the optimization method was used to fit tibia, fibula and calf biomechanical responses under different load positions in the quasi-static and dynamic loading conditions, which solved the problem with diversity caused by the test sample. When the elastic modulus of tibia and fibula were 18.43 and 18.23GPa, respectively, the failure strains were 1.156% and 0.8%, respectively, and physical bulk modulus of the calf was 11.33MPa, the calf model could fit the biomechanical responses well in many groups of tests. The FE model verification methods based on the optimization of multi-loading conditionsmade the FE model acquire better biofidelity effectively.

shank; Chinese human body; optimization; finite element model; validation

10.3969/j.issn.0258-8021. 2017. 02.010

2016-06-29， 錄用日期:2016-09-07

國(guó)家自然科學(xué)基金(11402296)；中央高校基本科研業(yè)務(wù)費(fèi)(531107040162)；湖南大學(xué)汽車車身先進(jìn)設(shè)計(jì)制造國(guó)家重點(diǎn)實(shí)驗(yàn)室自主研究課題(51475002)

R318

A

0258-8021(2017) 02-0195-010

*通信作者(Corresponding author)，E-mail: gfj0921@163.com