行人膝关节韧带的建模及验证研究

2012 年(第 34 卷)第 1 期 汽 车 工 程 Automotive Engineering 2012(Vol． 34)No． 1 2012012 行人膝关节韧带的建模及验证研究* * 中央高校基本科研业务费(531107040162)、湖南大学汽车车身先进设计制造国家重点实验室自主研究课题(61075004) 和上海市数字化汽车车身重点实验室开放课题基金(201004)资助。 原稿收到日期为 2010 年 11 月 26 日，修改稿收到日期为 2011 年 4 月 8 日。 张冠军1，曹立波1，官凤娇1，张 勇2，Yang King H． 3 (1. 湖南大学，汽车车身先进设计制造国家重点实验室，长沙 410082; 2. 华侨大学机电及自动化学院，厦门 361021; 3. Bioengineering Center，Wayne State University，Detroit，US 48201) ［摘要］ 基于人体模型 THUMS中的膝韧带模型，改进了材料参数和单元属性。采用近年较新的生物力学实 验数据，逐个验证韧带模型在不同应变率下的生物逼真度。结果明，膝韧带 ACL、PCL、MCL和 LCL在不同应变率 下的拉力与位移关系曲线与试验结果吻合较好，说明模型具有较好的生物逼真度，能较准确地模拟韧带的损伤和拉 伸响应。 关键词:汽车与行人碰撞;膝关节;韧带;有限元模型;验证 A Study on the Modeling and Its Validation for Pedestrian's Knee Ligaments Zhang Guanjun1，Cao Libo1，Guan Fengjiao1，Zhang Yong2 ＆ Yang King H． 3 1. Hunan University，State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body，Changsha 410082; 2. College of Mechanical Engineering and Automation，Huaqiao University，Xiamen 361021; 3. Bioengineering Center，Wayne State University，Detroit，US 48201 ［Abstract］ Based on the knee ligament model in total human model for safety (THUMS) ，with material pa- rameters and element properties modified，and using the newer experimental data of biomechanics in recent years， the biofidelities of ligament models under different strain rates are validated one by one． The results indicate that the relation curves of tension force vs． displacement of ACL，PCL，MCL and LCL knee ligaments in different strain dates agree well with experiment data，meaning that the models have good biofidelity and can accurately simulate the injury and tension response of ligaments． Keywords:vehicle-to-pedestrian crash;knee joint;ligament;finite element model;validation 前言 行人是道路交通参与者中的弱者，在交通事故 中遭受的伤害比较严重。统计表明，全世界每年约 有 117 万人死于交通事故，其中行人约占 65%。我 国每年交通事故造成约 2. 5 万行人死亡，约占交通 事故死亡总数的 25%，行人安全形势非常严峻［1］。 在汽车与行人碰撞过程中，通常是下肢最先与汽车 接触，下肢是最易发生损伤的部位，而膝关节韧带损 伤则是行人下肢的主要损伤形式且容易导致长期 残疾。 有限元分析法在膝关节韧带的损伤研究中获得 了较广泛的应用。用铰链模拟膝关节的膝关节韧带 有限元模型，如文献［2］和文献［3］中开发的用于乘 员安全研究的下肢模型，简单且计算速度快，但无法 预测韧带损伤。早期的韧带有限元模型，如文献 ［4］和文献［5］中开发的用于行人安全研究的下肢 模型，采用一维单元模拟，通常是弹簧阻尼单元，材 料参数简单，但存在不能定义韧带与人体其它组织 的接触等局限性。近期的韧带模型如文献［6］～文 献［8］中开发的行人下肢有限元模型，通常采用膜 (或壳)单元或体单元模拟，能够模拟韧带与其周围 组织的接触，生物逼真度较好，但这类模型需要较复 · 58 · 汽 车 工 程 2012 年(第 34 卷)第 1 期 杂的韧带组织的材料参数，需要以深入的韧带生物 力学实验为基础。韧带属于各向异性黏弹性材料， 主要承受拉力，且其力学性能受应变率的影响较大。 由于膝关节韧带的生物力学实验数据相对匮乏，不 同时期的模型很难进行较全面的验证。 全人体安全模型(total human model for safety， THUMS)或简称人体模型，主要针对乘员载荷情况 开展了一些验证，未对膝关节韧带进行独立的验证。 为获得良好的生物逼真度，本文中基于 THUMS膝韧 带有限元模型，对韧带材料的弹性模量、屈服应力和 单元属性等进行了改进，并综合较新的生物力学实 验数据，对膝关节主要韧带———ACL(前交叉韧带)、 PCL(后交叉韧带)、MCL(胫侧副韧带)和 LCL(腓侧 副韧带)进行了不同应变率下的验证，为后续行人下 肢有限元模型的开发奠定基础，同时也为行人下肢 损伤机理和安全防护研究提供准确高效的研究手段。 1 韧带的材料特性 韧带属于各向异性黏弹性材料，主要承受拉力。 文献［9］中通过对 MCL 的研究表明，胶原纤维横向 上的力学特性只为纵向的 1 /20。文献［10］中对日 本雌性兔子“股骨-MCL-胫骨”的动态拉伸实验研究 则表明，应变率对 MCL 的应力-应变曲线有较大 影响。 文献［9］、文献［11］～ 文献［14］中对 MCL 和 LCL材料的特性做了研究，获得了不同应变率下韧 带的弹性模量、极限应力和延伸率等参数，见表 1。 文献［15］中则在文献［13］实验的基础上对膝关节 主要韧带做了进一步的研究，获得了不同应变率下 膝韧带的最大拉力，见表 2。这些生物力学实验数据 为有限元模型的开发验证提供了必要的数据基础。 表 1 膝韧带MCL和 LCL的生物力学特性 文献 ［11］ ［13］ ［14］ ［9］ ［12］ 实验条件 年龄 20 ～ 75 55 ～ 90 21 ～ 30 62 ± 18 40 ～ 66 载荷 2mm /s(8mm /s) 1. 7mm /s 应变率 1 /s 应变率 1% /s 1. 6mm /s(1 600mm /s) MCL 弹性模量 /MPa － － － 332 ± 58 94. 0 ± 68. 2(182. 5 ± 78. 2) 极限应力 /MPa － － － 39 ± 5 15. 1 ± 9. 5(16. 9 ± 2. 8) 最大拉力 /N 468 ± 106(665 ± 236) 465 ± 190 － － 1 305 ± 278(1 215 ± 209) 延伸率 /% 23 ± 8(24 ± 4) 13 ± 6 － 17 ± 2 20. 3 ± 3. 9(11. 5 ± 5. 3) LCL 弹性模量 /MPa － － 345 ± 22. 4 － 367 ± 146(626 ± 131) 极限应力 /MPa － － 36 ± 12 － 23. 7 ± 9. 3(60. 2 ± 25. 7) 最大拉力 /N － 425 ± 95 － － 387 ± 92(571 ± 148) 延伸率 /% － 19 ± 8 15 ± 4 － 7. 1 ± 0. 9(10. 5 ± 4. 6) 表 2 膝韧带的耐受限度 韧带 载荷 /(mm/s) 最大拉力 /kN MCL 0. 016 1. 0 ± 0. 34 1. 6 1. 3 ± 0. 28 1 600 1. 2 ± 0. 21 LCL 0. 016 0. 24 ± 0. 02 1. 6 0. 39 ± 0. 09 1 600 0. 57 ± 0. 15 aACL 1 600 0. 99 ± 0. 56 pACL 1 600 1. 0 ± 0. 19 aPCL 1 600 0. 65 ± 0. 29 pPCL 1 600 0. 29 ± 0. 094 2 膝关节韧带有限元模型 为提高人体模型的生物逼真度，进行单个组织 的验证十分必要。因此，在 THUMS膝关节韧带模型 的基础上，对材料参数和单元属性等进行改进，并进 行不同应变率下的全面验证。模型使用 LS-DYNA 有限元分析软件进行求解分析。 韧带的厚度一般较小，通常采用壳单元模拟。 同时，壳单元相对于体单元，计算效率较高。由于膜 单元在大模型计算中的稳定性不易把握，为尽量保 证韧带在膝关节模型中具有较好的稳定性，采用 LS- DYNA中的 16 号全积分壳单元模拟膝关节 ACL、 PCL、MCL和 LCL 4 条主要韧带，并在韧带厚度方向 上设置 5 个积分点，同时使用 8 号沙漏模式。韧带 与骨骼间通过节点重合进行连接，如图 1 所示。 根据近期发表的相关文献提供的材料参数，结 合我国人体韧带的厚度，重新设定各韧带的材料参 2012(Vol． 34)No． 1 张冠军，等:行人膝关节韧带的建模及验证研究 · 59 · 图 1 膝关节韧带有限元模型 数，并在韧带厚度偏差允许的范围内对各韧带的厚 度进行适当的调整，使仿真结果与实验结果尽可能 吻合。其中，ACL 的厚度为 8. 0mm、PCL 的厚度为 6. 0mm、MCL 的厚度为 3. 3mm、LCL 的厚度为 2. 9mm。鉴于膝横韧带、板股后韧带和髌骨与髌韧 带对行人膝关节动力学响应的影响较小，不再对其 进行单独的验证。 由于应变率对韧带材料的力学性能影响较大，4 条主要韧带均采用#19 塑性材料(* MAT_STRAIN_ RATE_DEPENDENT_PLASTICITY) ，可分别定义材 料在不同应变率下的弹性模量、屈服应力和极限应 力(应变)等。通过定义有效塑性应变模拟韧带损 伤。根据生物力学实验数据［9，11 － 13］，设定密度均为 1 100kg /m3，泊松比均为 0. 49，4 条主要韧带的基本 材料参数如表 3 所示。 表 3 膝韧带的材料参数 韧带 参数 应变率 / s － 1 0 0. 3 3 1010 MCL 弹性模量 /MPa 155 160 200 230 屈服应力 /MPa 21 30 37 38 切线模量 /MPa 100 110 117 150 极限应变 0. 02 0. 02 0. 02 0. 02 LCL 弹性模量 /MPa 80 110 120 170 屈服应力 /MPa 15 26 30 32 切线模量 /MPa 40 41 44 120 极限应变 0. 02 0. 02 0. 02 0. 02 ACL 弹性模量 /MPa 175 180 210 260 屈服应力 /MPa 33 45 50 55 切线模量 /MPa 100 100 110 140 极限应变 0. 04 0. 04 0. 04 0. 04 PCL 弹性模量 /MPa 155 160 200 230 屈服应力 /MPa 21 30 37 38 切线模量 /MPa 100 110 117 150 极限应变 0. 04 0. 04 0. 04 0. 04 3 膝关节韧带有限元模型验证 3. 1 膝关节韧带的生物力学实验 文献［13］和文献［15］中先后进行了人体膝关 节韧带的拉伸实验，对人体韧带样本以一定的速度 (0. 016、1. 6 和 1 600mm /s)进行拉伸，研究应变率对 人体韧带力学响应的影响。本文中将主要依据这些 较新的生物力学实验，分别对膝关节的主要韧带 ACL、PCL、MCL和 LCL进行不同应变率下的生物逼 真度验证。 3. 2 韧带模型验证 根据文献［13］和文献［15］中的生物力学实验 建立膝关节各韧带的验证模型，如图 2 所示。由于 长骨皮质骨的弹性模量比松质骨大很多，而韧带均 附着在皮质骨上，因此验证模型中仅保留了骨骼的 皮质骨，而去除了松质骨，以提高计算效率。约束胫 骨与腓骨下部节点的 6 个自由度和股骨上部节点的 5 个自由度，而仅释放其向上的移动自由度。拉伸 速度根据实验设置，分为准静态(0. 016mm /s)、低速 (1. 6mm /s)和高速(1 600mm /s)3 种。仿真中输出 韧带拉力和股骨的位移，绘制拉力与位移的关系 曲线。 图 2 膝关节主要韧带验证模型 4 验证结果及分析 4. 1 MCL验证结果及分析 MCL在高速(1 600mm /s)下拉伸的仿真过程如 图 3 所示。约在 8. 9ms 时，MCL 在股骨连接点附近 发生撕裂，此时股骨位移约为 13. 2mm，拉力约为 1. 30kN。韧带拉力与其位移的关系曲线如图 4 所 示。由图可知，MCL有限元模型的仿真结果能很好 地与实验曲线吻合。 · 60 · 汽 车 工 程 2012 年(第 34 卷)第 1 期 图 3 MCL高速拉伸仿真过程 图 4 MCL高速拉伸力-位移关系曲线 MCL在低速(1. 6mm /s)下受拉，在股骨连接点 附近发生撕裂，此时股骨位移约为 12. 0mm，拉力约 为 1. 01kN。 MCL在准静态下受拉，在股骨连接点附近发生 撕裂，此时股骨位移约为 10. 7mm，拉力约为 0. 72kN。韧带拉力与其位移的关系曲线如图 5 所 示，仿真结果与实验结果基本吻合。 图 5 MCL准静态拉伸力-位移关系曲线 4. 2 LCL验证结果及分析 LCL在高速(1 600mm /s)下拉伸的仿真过程如 图 6 所示。约在 6. 0ms 时，LCL 在股骨连接点附近 发生撕裂，此时股骨位移约为 9. 1mm，拉力约为 0. 45kN。 LCL在低速(1. 6mm /s)下受拉，在股骨连接点 附近发生撕裂，此时股骨位移约为 9. 0mm，拉力约为 图 6 LCL高速拉伸仿真过程 0. 31kN。 LCL在准静态(0. 016mm /s)下受拉，在股骨连 接点附近发生撕裂，此时股骨位移约为 7. 7mm，拉力 约为 0. 21kN。 高速、低速和准静态载荷下韧带拉力与位移的 关系曲线分别如图 7 ～图 9 所示。结果显示，准静态 载荷下，LCL模型的响应能与实验结果较好地吻合。 随着应变率的提高，LCL 模型的极限拉力比实验结 果小，但仍在实验区间范围内，如表 3 所示。 图 7 LCL高速拉伸力-位移关系曲线 图 8 LCL中速拉伸力-位移关系曲线 4. 3 ACL和 PCL验证结果及分析 目前，交叉韧带仅有高速条件下的力学特性数 据［13，15］，如表 4 所示。其中 aACL 和 pACL 分别为 ACL的前束和后束，aPCL 和 pPCL 分别为 PCL 的前 束和后束。由于前后交叉韧带仅在功能上存在分 束，在解剖学结构上并无明显的分束，本文中仍将 2012(Vol． 34)No． 1 张冠军，等:行人膝关节韧带的建模及验证研究 · 61 · 图 9 LCL准静态拉伸力-位移关系曲线 前、后交叉韧带作为一个整体建模。所以，仿真中获 得的 ACL(或 PCL)的最大拉力应不小于其任一束的 最大拉力，且不大于其中两束韧带的合力。验证中 所得到的 ACL和 PCL的最大拉力分别为 1. 27kN和 1. 10kN，均在实验数据范围内。 表 4 膝韧带仿真与实验数据对比 韧带 载荷 / (mm/s) 最大拉力 /kN 实验 仿真 MCL 0. 016 1. 0 ± 0. 34 0. 72 1. 6 1. 3 ± 0. 28 1. 01 1 600 1. 2 ± 0. 21 1. 30 LCL 0. 016 0. 24 ± 0. 02 0. 21 1. 6 0. 39 ± 0. 09 0. 31 1 600 0. 57 ± 0. 15 0. 45 ACL aACL pACL 1 600 0. 99 ± 0. 56 1. 00 ± 0. 19 1. 27 PCL aPCL pPCL 1 600 0. 65 ± 0. 29 0. 29 ± 0. 09 1. 10 将上述各韧带(MCL、LCL、ACL、PCL)的验证结 果与实验结果列于表 4 以进行比较。结果显示，各 韧带的生物力学特性与实验结果基本吻合，模型具 有较好的生物逼真度。 虽然仿真结果与实验结果吻合较好，但依然存 在一些差异，主要的原因如下。 (1)通过扫描等手段准确获取韧带的外观几何 参数和连接位置有一定的难度，因此在有限元模型 中较准确地模拟韧带的大小和位置仍是一个难题。 (2)生物力学实验所使用的样本数量有限，且 样本间存在一定差异，致使实验所获得的材料参数 差异较大。在数学模型中设置相关参数须首先确定 参数的变化区间，然后凭经验选择合适的参数值。 因此，很难通过一套特定的参数模拟不同人体韧带 的力学特性，仿真与实验结果仍存在一些差异。 由于行人在碰撞过程中，韧带拉伸时的应变率 处于动态的变化过程，而不同应变率下的韧带材料 特性不同。因此，研究不同应变率下韧带的损伤具 有重要意义。而基于良好的韧带生物逼真度所建立 的膝关节模型，能有效地预测行人韧带的损伤。 5 结论 根据较新的韧带材料参数，采用壳单元建立了 膝关节韧带有限元模型，利用较新的生物力学实验， 分别对胫侧副韧带 MCL、腓侧副韧带 LCL、前交叉韧 带 ACL和后交叉韧带 PCL 进行了拉伸验证。采用 基于应变率的弹塑性材料模型和基于有效塑性应变 的损伤准则，能够有效模拟膝关节韧带的力学特性 和损伤。通过对膝关节主要韧带的验证，能够为开 发较高生物逼真度的膝关节模型奠定基础。 参考文献 ［1］ Zhang G，Cao L，Hu J，et al． A Field Data Analysis of Risk Fac- tors Affecting the Injury Risks in Vehicle-to-Pedestrian Crashes ［C］． 52nd AAAM Annual Conference，San Diego，2008． ［2］ Bedewi P G． The Biomechanics of Human Lower Extremity Injury in the Automotive Crash Environment［D］． Washington D． C: George Washington Universtiy，1998． ［3］ Wykowski E，Sinnhuber R，Appel H． Finite Element Model of Hu- man Lower Extremities in a Frontal Impact［C］． 1998 International IRCOBI Conference on the Biomechanics of Impact． Goteborg， Sweden，1998． ［4］ Schuster P J，Chou C C，Prasad P，et al． Development and Vali- dation of a Pedestrian Lower Limb Non-linear 3-d Finite Element Model［J］． Stapp Car Crash J，2000，44:315 － 334． ［5］ Yang J K，Wittek A，Kajzer J． Finite Element Model of the Hu- man Lower Extremity Skeleton System in a Lateral Impact［C］． 1996 International IRCOBI Conference on the Biomechanics of Im- pact，1996． ［6］ Beillas P，Begeman P C，Yang K H，et al． Lower Limb:Ad- vanced FE Model and New Experimental Data［J］． Stapp Car Crash J，2001，45:469 － 94． ［7］ Takahashi Y，Kikuchi Y，Mori F，et al． Advanced FE Lower Limb Model for Pedestrians［C］． The 18th International Technical Con- ference on the Enhanced Safety of Vehicles (ESV)． Nagoya，Ja- pan:US Department of Transportation，National Highway Traffic Safety Administration，2003． ［8］ 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． ［9］ Quapp K M，Weiss J A． Material Characterization of Human Medi- al Collateral Ligament［J］． Journal of Biomechanical Engineering， 1998，120(6) :757 － 763． (下转第 75 页) 2012(Vol． 34)No． 1 刚宪约，等:基于位移法的多轴汽车轮胎载荷计算的研究 · 75 · 5 数值算例 参考文献［4］中的数据，某四轴汽车满载总质 量为 27 397kg，质心高度为 1 350mm。根据本文中 定义的坐标系，各轮胎参数见表 1。计算如下工况 的轮胎载荷:(1)弯曲工况;(2)扭转工况，轮胎 1 抬 高 100mm;(3)制动工况，ax = － 7m /s 2; (4)转弯工 况，ay = 4m /s 2。结果见表 2。 表 1 四轴汽车轮胎参数 轮胎序号 等效刚度 /(N /mm) x /mm y /mm 1 440. 5 3 175 1 065 2 440. 5 3 175 － 1 035 3 440. 5 1 651 1 065 4 440. 5 1 651 － 1 035 5 366. 5 － 2 159 1 065 6 366. 5 － 2 159 － 1 035 7 366. 5 － 3 683 1 065 8 366. 5 － 3 683 － 1 035 6 结论 基于位移法建立的弯曲、扭转、起动 /制动和转 弯等典型工况的轮胎载荷计算模型，无须建立繁杂 的变形协调方程，格式统一，仅包含简单的矩阵乘法 表 2 各工况下的轮胎载荷 N 轮胎序号 弯曲工况 扭转工况 制动工况 转弯工况 转弯工况* 1 36 260 60 302 50 998 14 793 17 032 2 37 307 29 321 52 045 58 774 56 534 3 36 191 20 005 43 875 14 724 16 963 4 37 237 33 074 44 922 58 705 56 465 5 29 967 24 449 21 690 12 107 13 970 6 30 838 35 323 22 561 48 699 46 836 7 29 910 27 571 15 764 12 049 13 912 8 30 781 38 445 16 635 48 642 46 778 注:表 2 中“转弯工况* ”为不考虑整车质心侧向位移引起的重 力侧倾力矩的影响计算得到的轮胎载荷数据。 和三阶线性方程组求解运算，特别方便计算机求解。 对于类似于转弯制动的复合工况或者汽车在颠簸路 面上行驶时的动态工况，仅须将本文中相应类型的 整车载荷向量进行叠加，即可方便地求解相应的轮 胎载荷。 参考文献 ［1］ 纪晓鹏．轿车制动力分配与轴荷分配的匹配关系［J］． 上海汽 车，2002(9) :4 － 5． ［2］ 尹俊峰． 多轴专用车底盘轴荷估算方法［J］． 专用汽车，2004 (3) :27 － 28． ［3］ 柴新伟，杨世文．三轴汽车轴荷计算及轴距选择［J］．机械管理 与开发，2010，25(1) :53 － 54，56． ［4］ 王国军，陈欣，徐安桃，等． 多轴汽车轴荷的合理分配与确定 ［J］．军事交通学院学报，2009，11(3) :  53 － 56． (上接第 61 页) ［10］ Yamamoto S，Saito A，Nagasaka K，et al． The Strain-Rate De- pendence of Mechanical Properties of Rabbit Knee Ligaments ［C］． The 19th International Technical Conference on the En- hanced Safety of Vehicles (ESV)． Nagoya，Japan:US Depart- ment of Transportation，National Highway Traffic Safety Adminis- tration，2003． ［11］ Kennedy J C，Hawkins R J，Willis R B，et al． Tension Studies of Human Knee Ligaments:Yield Point，Ultimate Failure，and Dis- ruption of the Cruciate and Tibial Collateral Ligaments［J］． Jour- nal of Bone and Joint Surgery，1976，58(3) :350 － 355． ［12］ Kerrigan J R，Ivarsson B J，Bose D，et al． Rate-Sensitive Consti- tutive and Failure Properties of Human Collateral Knee Ligaments ［C］． 2003 International IRCOBI Conference on the Biomechan- ics of Impact． Lisbon，Portugal，2003． ［13］ Marinozzi G，Pappalardo S，Steindler R． Human Knee Ligaments: Mechanicla Tests and Ultrastructural Observations［J］． Italian Journal of Orthopaedics and Traumatology，1983，9(2) :231 － 240． ［14］ Butler D L，Kay M D，Stouffer D C． Comparison of Material Properties in Fascicle-Bone Units from Human Patellar Tendon and Knee Ligaments［J］． Journal of Biomechanics，1986，19(6) : 425 － 432． ［15］ van Dommelen J A W，Ivarsson B J，Jolandan M M，et al． Char- acterization of the Rate-Dependent Mechanical Properties and Failure of Human Knee Ligaments［C］． SAE Paper 2005 － 01 － 0293．