基于优化的中国50th人体大腿有限元模型验证方法的研究*

张冠军，王龙亮，胡跃群，杜现平，曹立波

(1.湖南大学，汽车车身先进设计制造国家重点实验室，长沙 410082； 2.中南大学湘雅三医院放射科，长沙 410013)

2016105

基于优化的中国50th人体大腿有限元模型验证方法的研究*

张冠军1，王龙亮1，胡跃群2，杜现平1，曹立波1

(1.湖南大学，汽车车身先进设计制造国家重点实验室，长沙 410082； 2.中南大学湘雅三医院放射科，长沙 410013)

目前的假人和人体有限元模型大多是根据欧美人体建立的，故由人体身材的差别引起生物力学的响应的差异值得探讨。建立中国50th人体有限元模型有助于提高中国人体的损伤防护水平。通过CT扫描数据获得大腿的几何模型,并将其缩放到中国50th的人体股骨尺寸。根据股骨解剖学结构将股骨头、股骨颈、股骨体和内外侧髁等的皮质骨和松质骨赋予不同的材料参数，并利用LS-OPT对大腿模型的材料参数进行优化，使仿真结果与缩放至中国50th人体股骨尺寸的实验数据吻合，以满足不同的加载部位、加载方向和加载速率等载荷工况的验证要求。仿真结果表明优化后的模型具有较高的生物逼真度，并能适应多种载荷工况。

大腿；股骨；有限元模型；验证

前言

世界卫生组织在2013年的报告中指出，全世界每年有接近124万的道路使用者在交通事故中丧生，并且受伤人数更是高达2 000～5 000万人[1]。下肢损伤在人体最容易受到损伤的8个部位中高居第2位[2]。下肢损伤不仅给受害者带来长期生活不便，而且也给社会和家庭带来沉重的负担。为更好地了解下肢损伤的机理，国外许多研究者采用尸体撞击实验进行研究，由于实验成本高、样本难以获得和重复性差，所以具有一定的局限性。随着计算机技术的发展，大量的数学模型应用于汽车安全性研究。人体下肢有限元模型由于能够详细地计算出骨骼内部的应力应变分布，并能获知人体内部组织与结构的力学响应，从而被广泛应用于乘员和行人的下肢损伤研究。

国外下肢有限元模型发展较早，文献[3]中就开发了包括股骨、胫骨和主要韧带的行人下肢模型，但几何外形不够准确。文献[4]中开发了乘员下肢有限元模型，但仅包含骨骼，且将骨骼定义为刚体。文献[5]中基于LS-DYNA求解器开发了行人下肢有限元模型，该模型定义了肌腱、肌肉和皮肤，并进行了较详细的验证。文献[6]中基于先前的模型更详细地提取骨骼几何外形，建立了乘员下肢有限元模型。

国内建立下肢有限元模型起步较晚，文献[7]～文献[10]中根据国外模型进行了材料改进和验证,但并没有根据中国人体的几何外形建立中国50百分位人体下肢有限元模型。

目前的下肢有限元模型大多是根据欧美人体建立的，由于中国人体在尺寸上与欧美人体有较大差别，模型的生物力学响应与中国人体存在多大差异目前并不明确。

为研究中国50th人体的生物力学响应，本文中利用中国人体下肢的医学数据建立中国50th人体大腿有限元模型。通过下肢的CT扫描数据，获得精确的大腿几何模型。根据大腿的解剖学结构和各部位的材料特性赋予其不同的材料参数，并利用优化方法对大腿材料参数进行优化，以使大腿模型在加载的不同部位、方向和速率下的生物力学响应与实验数据较好地吻合。

1 大腿有限元模型的建立

1.1 几何模型

基于计算机断层成像(CT)并综合运用医学图像处理软件Mimics和三维处理软件Geomagic Studio提取大腿几何模型，保证大腿具有详细的解剖学结构和准确几何形状。CT数据来源于一位骨骼正常的接近中国50th的男性血管疾病患者(身高173.1cm，体质量69.7kg)。CT图像和修复后的股骨和肌肉三维图如图1所示。

图1 大腿CT图像和修复后的股骨、肌肉三维图

1.2 有限元模型

综合使用ANSYSY ICsEM CFD和Hypermesh软件对大腿模型进行网格划分。股骨体皮质骨使用六面体网格模拟，股骨两端皮质骨较薄，选用壳单元模拟。股骨两端的松质骨采用体单元模拟。肌肉使用体单元模拟,用CONTACT_TIED_NODES_TO_SURFACE将肌肉和股骨连接。皮肤采用壳单元模拟，使用共节点方式与肌肉连接。为减小应力集中，在股骨体与股骨两端皮质骨过渡区域选用阶梯形状逐渐过渡，如图2所示。

图2 股骨网格的划分

股骨模型雅克比小于0.6的单元比例不超过1%，最小值为0.42；翘曲度大于20°的单元比例不超过2%，最大值为120°；长宽比大于3.5的单元比例小于5%，最大值为5.2；最小单元尺寸4.69mm；最大单元尺寸6.424mm。大腿模型单元总数为59 634，其中体单元54 964个，壳单元4 670个。

1.3 有限元模型的尺寸缩放

由于本文中模型的人体尺寸大于中国50th人体尺寸[11](身高170.8cm，体质量65kg)，所以需要将大腿模型缩放到中国50th人体大腿尺寸。

轴向缩放系数主要依据中国50th与患者的股骨长的比例；横向缩放系数主要依据中国50th与患者的股骨颈横径、股骨体中部横径、股骨髁宽比例的平均值；径向缩放系数主要依据中国50th与患者的内侧髁长和外侧髁长比例的平均值[12]。最终确定的股骨轴向缩放系数为0.972，横向缩放系数为0.978，径向缩放系数为0.983，如表1所示。

表1 股骨参数缩放表

将股骨模型按照表1中的3个方向的缩放系数平均值放至中国50th男性股骨的尺寸，并依据此比例对大腿肌肉和皮肤模型进行缩放，最终获得了中国50th男性大腿有限元模型。

1.4 大腿材料模型

股骨皮质骨在受拉和受压时表现出不同的力学特性[13]。因此，选择可分别定义拉、压应力-应变曲线的弹塑性材料(#124)模拟皮质骨[13]。同时，使用弹塑性材料模拟松质骨，使用黏弹性材料模拟皮肤和肌肉组织。采用失效应变作为失效准则来模拟损伤，在材料达到设定的失效应变后自动删除失效单元。股骨近心端材料的弹性模量并不完全一致[14]，故将股骨近心端皮质骨和松质骨分别分成3个部分，如图2所示。根据相关文献设定模型优化前各参数取值如表2所示[6,15-19]。

皮肤密度为1 000kg/m3，弹性模量为1MPa，泊松比为0.3，厚度为1mm[15,20]；肌肉选用黏弹性材料(#92)，密度为1 000kg/m3，体积模量为20MPa，C1为0.12kPa，C2为0.25kPa，S1为1.162，S2为0.808，T1为10.43ms，T2为84.1ms[21]。

2 大腿有限元模型的优化验证

验证分为准静态验证和动态验证，如表3所示。

表2 股骨模型材料参数设置

表3 有限元模型验证

注：A-P表示载荷由前向后方向；L-M表示载荷由外侧向内侧方向。

传统的验证方式采用试错法人工调节材料参数对模型进行验证，效率低且准确度较差。本文中采用优化方法自动获得最佳材料参数值，能同时保证模型的生物力学响应与多个实验结果吻合。

为更准确地获得材料参数，先对股骨近心端、中部、远心端3个加载位置的动态仿真进行股骨体皮质骨的弹性模量和拉、压应力-应变曲线进行优化(此时未定义失效应变)，将得到的最佳材料参数代入到大腿模型。然后对大腿中部、远心端两个加载位置的动态仿真进行肌肉材料参数优化，使用优化结果更新大腿模型。最后，对股骨近心端、中部、远心端和大腿中部、远心端5个加载位置的动态仿真进行股骨体皮质骨的失效应变优化。最终获取一组最佳材料参数，对股骨的准静态实验进行校核，如图3所示。

图3 大腿模型验证流程图

材料参数优化分3步进行：(1)在股骨体皮质骨的弹性模量和拉、压应力-应变曲线优化中,设计变量分别为弹性模量和拉、压应力-应变曲线X、Y轴的缩放系数；(2)在肌肉材料参数优化中，设计变量为体积模量；(3)在股骨体皮质骨的失效应变优化中，设计变量为失效应变。

各设计变量的取值范围如表4所示。

表4 设计变量的取值范围

前两步优化的目标是仿真与实验的力-位移曲线的均方差f(X)[29]的最小化，即

(1)

式中：X为设计变量，k=1,2,…,K，K为实验曲线的条数；m=1,2,…,Pk，Pk为第k条曲线中计算点的个数；Wm为权重系数，本文中各个工况全为1；fm(X)为响应面近似模型的计算值；Gm为实验测试点的值。

股骨体皮质骨失效应变的优化目标是各工况仿真与实验的力-位移曲线的力极大值点对应的位移之差绝对值f(X)的最小化，即

(2)

在优化中，如果某个工况下有多个目标实验曲线，则采用这些实验曲线的平均曲线作为该工况优化的目标曲线。平均曲线的计算方法参照文献[30]。

2.1 实验数据缩放

股骨动态仿真采用文献[25]中的股骨长度为467mm的数据，大腿动态仿真采用文献[27]中缩放至美国50th人体的数据。为了准确验证中国50th人体的有限元模型，还需将实验数据缩放到与中国50th人体几何相对应的数值，以消除尺寸差异对生物力学响应的影响。

文献[2]中为了准确获得美国50th人体的股骨和大腿的损伤耐受限度和弯曲响应，采用股骨长度缩放系数λL计算位移缩放系数λD和力缩放系数λF，然后对实验数据进行缩放：

λL=L50th/L样本

(3)

λD=λL

(4)

(5)

式中：L50th为美国50th人体的股骨长度；L样本为实验样本的股骨长度。

实验曲线中的位移、力、弯矩乘以相应的缩放系数即可得到与美国50th人体相对应的曲线。

由于股骨准静态实验所用股骨长度尺寸不能确定，本文中只对股骨和大腿动态仿真的实验数据进行缩放。股骨和大腿动态验证的实验数据来源于缩放后的数据，所以在股骨样本的近心端1/3处、中间、远心端1/3处的缩放系数相同，大腿样本的中间、远心端1/3处的缩放系数相同，如表5所示。实验样本的位移、力、弯矩乘以相应缩放系数即可得到与中国50th人体相对应的缩放实验曲线。

表5 实验样本的缩放系数和股骨长度

2.2 股骨和大腿模型动态仿真

鉴于行人侧面遭受撞击的几率远大于其他方向[27]，本文中仅进行L-M方向的动态仿真。文献[32]和文献[27]中对股骨的近心端1/3处、中部、远心端1/3处和大腿中部、远心端1/3处进行了三点弯曲实验，文献[28]和文献[25]中则使用其数据开展了模型验证。参照上述实验和仿真建立股骨和大腿的动态三点弯曲验证模型，如图4和图5所示。股骨和大腿的加载速度分别为1.2和1.5m/s。

图4 股骨近心端1/3、中部、远心端1/3处动态三点弯曲验证

图5 大腿中部、远心端1/3处动态三点弯曲验证

2.3 股骨模型准静态仿真

股骨准静态三点弯曲仿真可以验证大腿模型材料参数的设置。依据文献[28]和文献[22]中的相关描述建立仿真模型。用直径为25mm的刚性圆筒冲击器，以0.01m/s的速度对股骨中部进行加载。根据载荷加载方向的不同，股骨准静态三点弯曲仿真分为A-P和L-M两个方向的验证。仿真设置如图6所示。

图6 股骨准静态三点弯曲验证

2.4 股骨头准静态压溃仿真

在汽车正面碰撞中，乘员的股骨近心端是易受伤害部位[6]。为确保大腿有限元模型的生物逼真度，有必要对股骨头进行准静态压溃验证。文献[26]中进行了18组股骨头准静态压溃实验见图7。股骨轴向与竖直方向成20°，骨干区域完全约束，用直径为30mm的圆柱形刚性冲击器以0.5mm/s的速度对股骨头加载，直至断裂。仿真设置如图7(b)所示。

图7 股骨头准静态压溃验证

3 模型仿真结果

3.1 大腿材料优化结果

大腿材料优化参数收敛过程如图8所示。由图可见，随着迭代次数的增加，股骨体皮质骨的弹性模量、应力-应变曲线的X、Y轴缩放系数以及肌肉体积模量的兴趣域空间逐渐缩小，各参数取值逐渐稳定。

图8 优化材料参数收敛过程

目标函数优化历程如图9所示。由图可见，随着迭代次数的增加，股骨材料优化以及肌肉材料优化的目标函数逐渐稳定并趋于最小化，当迭代次数达到15次后，目标函数值达到要求。

图9 目标函数优化历程

最终股骨体皮质骨弹性模量为14.83GPa，拉、压应力-应变曲线X、Y轴的缩放系数分别为0.601 3和0.464 1，肌肉体积模量为11.33MPa，失效应变为1.156%。

3.2 动态仿真优化结果

股骨和大腿材料参数优化结果如图10所示。由于股骨中部和远心端1/3处实验曲线有3条，本文在优化中将其拟合成一条均值线。股骨动态三点弯曲仿真的力-位移曲线与目标曲线吻合程度较高；大腿模型三点弯曲仿真的力-位移曲线在上半段与目标曲线拟合较好，虽然在下半段存在一些偏差，但曲线的走势还是一致的。以上结果表明：通过优化，股骨近心端、中部、远心端和大腿中部、远心端5个加载位置的动态仿真都很好地与目标曲线吻合。

图10 股骨和大腿优化仿真结果

3.3 准静态仿真结果

使用优化后的材料参数进行股骨A-P、L-M方向的准静态三点弯曲仿真，由于实验涉及的股骨长度尺寸不能确定，在准静态验证时未对实验曲线进行缩放，其仿真曲线和实验曲线如图11所示。仿真结果在实验曲线范围内，并与实验曲线保持了较好的一致性。模型发生损伤时的变形量和碰撞力也与实验曲线吻合较好。

图11 股骨在A-P、L-M加载方向准静态三点弯曲仿真与实验对比

3.4 股骨头仿真结果

仿真结果如图12所示。由图可见，股骨头在碰撞力达到7.3kN时发生断裂，文献[26]中的18组实验结果显示，碰撞力在3.1～15.0kN时股骨颈发生断裂，平均值为8.4±3.0kN。虽然仿真结果比实验均值小，但仍在实验范围之内。

图12 股骨头断裂时刻的碰撞力

4 结论

由CT扫描数据获取几何模型，划分网格后根据中国50百分位人体股骨的尺寸将模型缩放至中国50th人体尺寸，并将实验数据缩放到与中国50th人体几何相对应的数值。使用优化方法对股骨体皮质骨材料的弹性模量、拉压应力-应变曲线、失效应变和肌肉的体积模量进行优化，使股骨近心端、中部、远心端3个加载位置的动态仿真，大腿中部、远心端两个加载位置的动态仿真，股骨中部A-P和L-M两个加载方向的准静态仿真和股骨头准静态压溃仿真中模型的生物力学响应与实验结果吻合良好，所建立的中国50th人体大腿有限元模型具有较好的生物逼真度。因此，基于优化的模型验证方法可用于多实验工况的模型验证，验证效率高且所获得的模型生物逼真度高。

[1] BURTON A, HARVEY A, BLAKEMAN D, et al. Global status report on road safety 2013: supporting a decade of action[R]. Geneva, Switzerland: World Health Organization (WHO),2013.

[2] KIM Y S, CHOI H H, CHO Y N, 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] 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,1993:117-129.

[4] 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,1998:101-116.

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

[6] UNTAROIU C D, YUE N, SHIN J. A finite element model of the lower limb for simulating automotive impacts[J]. Ann Biomed Eng,2013,41(3):513-526.

[7] 张冠军.行人下肢的碰撞损伤特性及相关参数研究[D].长沙:湖南大学,2009.

[8] 杨济匡,方海峰.人体下肢有限元动力学分析模型的建立和验证[J].湖南大学学报(自然科学版),2005,32(5):31-36.

[9] 李正东,刘宁国,黄平,等.下肢有限元模型的建立及损伤机制重建[J].中国司法鉴定,2012(6):37-42.

[10] 蒋小晴,杨济匡,王丙雨,等.乘员股骨在轴向压力—弯矩下的损伤生物力学机理研究[J].力学学报,2014(3):465-474.

[11] 刘俊先,张兴和.中国正常人体测量值[M].北京:中国医药科技出版社,1994.

[12] CRANDALL J R, PORTIER L, PETIT P, et al. Biomechanical response and physical properties of the leg, foot, and ankle[C]. Stapp Car Crash Conference,1996.

[13] NOVITSKAYA E, CHEN P Y, 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] LOTZ J C, GERHART T N, HAYES W C. Mechanical properties of metaphyseal bone in the proximal femur[J]. Journal of Biomechanics,1991,24(5):317-329.

[15] BEILLAS P, BEGEMAN P C, YANG K H, et al. Lower limb: advanced FE model and new experimental data[J]. Stapp Car Crash J,2001,45:469-494.

[16] BURSTEIN A H, REILLY D T, MARTENS M. Aging of bone tissue: mechanical properties[J]. J Bone Joint Surg Am,1976,58(1):82-96.

[17] GALBUSERA F, FREUTEL M, DURSELEN L, et al. Material models and properties in the finite element analysis of knee ligaments: a literature review[J]. Front Bioeng Biotechnol,2014,2:54.

[18] MARTENS M, Van AUDEKERCKE R, DELPORT P, et al. The mechanical characteristics of cancellous bone at the upper femoral region[J]. J Biomech,1983,16(12):971-983.

[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] KARIMI A, NAVIDBAKHSH M. Material properties in unconfined compression of gelatin hydrogel for skin tissue engineering applications[J]. Biomed Tech (Berl),2014,59(6):479-486.

[21] SNEDEKER J G, MUSER M H, WALZ F H. 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-457.

[22] YAMADA H, EVANS F G. Strength of biological materials[M]. Baltimore: Williams & Wilkins,1970:297.

[23] MATHER B S. Variation with age and sex in strength of the femur[J]. Med Biol Eng,1968,6(2):129-132.

[24] STROMSOE K, HOISETH A, ALHO A, et al. Bending strength of the femur in relation to non-invasive bone mineral assessment[J]. J Biomech,1995,28(7):857-861.

[25] TAKAHASHI Y, KIKUCHI Y, MORI F, et al. Advanced FE lower limb model for pedestrians[C]. 18th International ESV Conference,2003.

[26] KEYAK J H, ROSSI S A, JONES K A, et al. Prediction of femoral fracture load using automated finite element modeling[J]. J Biomech,1998,31(2):125-133.

[27] KERRIGAN J R. A computationally efficient mathematical model of the pedestrian lower extremity[D]. University of Virginia,2008.

[28] UNTAROIU C D. 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]. University of Virginia,2005.

[29] 官凤娇.冲击载荷下的生物组织材料参数反求及损伤研究[D].长沙:湖南大学,2011.

[30] LESSLEY D, CRANDALL J, SHAW G, et al. A normalization technique for developing corridors from individual subject responses[C]. SAE Paper 2004-01-0288.

[31] EHLER E, LÖSCHE H. Die menschliche tibia unter biegebelastung[J]. Beitr. Orthop,1970,17(5):291-304.

[32] KERRIGAN J R, BHALLA K S, MADELEY N J, et al. Experiments for Establishing Pedestrian-Impact Lower Limb Injury Criteria[C]. SAE 2003 World Congress & Exhibition,2003.

A Study on the Validation Method of the 50th Percentile ChineseThigh Finite Element Model Based on Optimization

Zhang Guanjun1, Wang Longliang1, Hu Yuequn2, Du Xianping1& Cao Libo1

1.HunanUniversity,StateKeyLaboratoryofAdvancedDesignandManufacturingforVehicleBody,Changsha410082; 2.DepartmentofRadiology,TheThirdXiangyaHospitalofCentralSouthUniversity,Changsha410013

At present, crash dummies and finite element human models are mainly created based on occidental human statistical data, so the differences in biomechanical response caused by the differences in human stature are worth exploring. Establishing the 50th percentile Chinese model is conducive to improving the injury protection level of Chinese people. Based on CT scan data, thigh geometry model is obtained and scaled to the 50th percentile Chinese human femur size. Different material parameters are given to the cortical and cancellous bones of femoral head, femoral neck, femoral body and lateral condyles according to femoral anatomic structure. The material parameters of thigh model are optimized by using LS-OPT to make simulation results well agree with experimental data of the 50th percentile Chinese thigh size and meet the validation requirements for loading conditions with different locations, directions and rates of loading. The results of simulation show that the model optimized has higher biological fidelity and can be adapted to a variety of loading conditions.

thigh; femur; finite element model; validation

*国家自然科学基金(51205118)、湖南大学汽车车身先进设计制造国家重点实验室自主研究课题和中央高校基本科研业务费资助。

原稿收到日期为2015年5月14日，修改稿收到日期为2015年7月8日。