骑行时不同屈曲角度膝关节软骨受力分析

马雪梅 张春秋* 高丽兰 叶金铎 张西正

1(天津市先进机电系统设计与智能控制重点实验室，天津理工大学机械工程学院 天津 300384)2(军事医学科学院卫生装备研究所，天津 300161)

骑行时不同屈曲角度膝关节软骨受力分析

马雪梅1张春秋1*高丽兰1叶金铎1张西正2

1(天津市先进机电系统设计与智能控制重点实验室，天津理工大学机械工程学院 天津 300384)2(军事医学科学院卫生装备研究所，天津 300161)

在对骑行人车系统动力学研究的基础上，分析在骑行周期内膝关节软骨的受力状态，以期获得关节软骨受力的规律，增加对骑行时膝关节生物力学特性的理解。基于有限元分析法，建立包括股骨、胫骨、腓骨、髌骨、关节软骨、半月板及韧带组织的人体全膝关节有限元模型。将骨组织刚体化，并对其施加屈曲位移边界条件，包括胫骨相对股骨的内旋、内收、前移和外移以及髌骨相对股骨的屈曲、内旋、内倾以及外移。通过显式动力学分析计算，获得膝关节屈曲60°、80°和100°相位，同时得到骑行状态上述屈曲位处膝关节软骨的应力分布。结果 通过有限元分析，获得骑行姿态下膝关节相关相位的力学分布规律。结果表明，相同载荷下，最大von-Mises应力出现在屈曲100°位置，股骨软骨应力增幅达71.25%，髌骨软骨增大29.36%；随着骑行屈曲角度的增加，胫股关节高应力区逐渐向膝关节后部转移，髌股关节软骨受力逐渐上移。骑行时高应力发生在膝关节屈曲角度较大位置，胫骨平台软骨后侧、髌骨软骨上侧承受更大应力。

膝关节；软骨；骑行；有限元模型；屈曲

引言

骑行运动在日常生活中非常普遍。但是目前对骑行时人体膝关节生物力学特性的研究较少，而骑行过程中膝关节的受力规律以及膝关节软骨损伤与骑行的关系却是在崇尚健康骑行的当下值得关注的。

膝关节作为人体的一个重要负重关节，一直是骨关节生物力学领域的研究重点，人们利用仿真与实验方法对膝关节运动及承载状态进行探索。已有实验对膝关节做深度屈曲运动时产生的曲转矩和关节力进行研究，并与步行和上楼梯等日常运动进行对比，证明深度屈曲比其他步态产生更大的关节力，但该实验未考察软骨的受力状态[1]。也有专业人士对汽车驾驶、上楼梯以及深度屈曲时肌肉活动和关节力状态进行了研究，关注重点主要在关节活动而非关节软骨的接触力[2-4]。软骨产生过大应力是诱发关节炎的主要原因，所以关注关节软骨的应力状态将更有利于判断关节损伤与膝关节运动的关系。目前，应用实验方法测量膝关节内软骨接触应力依然困难，有限元技术的应用弥补了上述不足。建立人体下肢的肌骨模型，利用有限元技术分析步行时下肢肌肉力、膝关节韧带张力以及不同步态阶段关节软骨上的应力；通过改变材料属性，模拟患关节炎疾病的膝关节在相同条件下的受力状态[5]。Fan等建立完整的膝关节模型，研究处于跪姿的膝关节软骨及半月板的应力应变，并与相同条件下的站姿进行对比，发现膝关节在跪时将产生更大的应力，更易诱发膝关节损伤[6]。

现有研究表明，骑行周期内膝关节关节力最大值为831 N,出现在骑行周期的200°～230°范围内[7]。本研究在此基础上，选取骑行最低点(屈曲60°)以及骑行时膝关节力最大区域两个相位(屈曲80°和屈曲100°)作为考察点，利用ABAQUS显式动力学分析骑行时膝关节的软骨应力值及分布规律。

1 材料与方法

1.1 几何建模

对受试者临床检测，确认其无膝关节疾病。对该受试者膝关节进行螺旋CT扫描，扫描过程中膝关节处于伸直位置，告知受试者保持放松，以确保膝关节内部无预应力。扫描层厚0.8 mm，层数为1 450层。受试者知情同意，并签署知情同意书。

将以DICOM格式保存的扫描数据导入到软件MIMICS(Materialise Inc., Belgium)中。应用阈值分割法对组织进行预分割，使图像显示效果最佳。在此基础上，依据膝关节解剖结构[8-9]，对图像进行编辑，将编辑好的数据在MIMICS中经三维计算生成全膝关节的几何模型。该几何模型包括股骨、胫骨、腓骨、髌骨以及它们的软骨层，半月板，前后交叉韧带，内外侧副韧带和髌腱。将建好的几何模型以STL格式导入逆向工程软件中进行表面优化，最终拟合成NURBS曲面(见图1)。将拟合好的曲面以IGES格式保存后导入Pro/Engineer软件中将模型实体化。

图1 全膝关节几何模型Fig.1 The geometric model of knee joint

1.2 有限元模型建立

将实体化的全膝关节几何模型导入到有限元分析软件ABAQUS(Simulia Inc.,USA)中进行网格划分，由于此模型形状不规则，并且含有接触，所以选择修正二次四面体单元C3D10M。 软骨层的网格密度为2.5 mm，模型节点数为80 042，单元数为44 628。为了验证有限元模型网格密度的合理性，对相同加载条件下软骨层网格密度为1.5 mm时的模型进行计算，两种网格密度下的计算结果接近，选取2.5 mm合理。

对骨及其软组织的材料属性进行简化，根据相关的研究结果[10-11]赋值如表1所示。将韧带组织定义为超弹性材料，材料的可压缩性使用默认值，应变势能选用Ogden模型。根据现有的结果[12]给出该材料单轴测试的名义应力和名义应变。定义韧带的初始应变[13]如表2所示。

表1 膝关节组织的材料属性[10-11]Tab.1 The material properties of knee joint[10-11]

表2 不同部位韧带初始应变Tab.2 Initial strain of the ligaments

1.3 接触和边界条件设定

模型中所有软骨及韧带均被绑定在相应的骨组织上，将内、外侧半月板的前后角绑定在胫骨平台上，并定义如下接触：股骨软骨与胫骨软骨接触，股骨软骨与半月板接触，半月板与胫骨软骨接触，股骨软骨与髌骨软骨接触。设定软骨与软骨之间接触摩擦系数为0.002 5[14]。

获得螺旋CT扫描数据时，膝关节处于伸直位置，需要观察骑行屈曲60°、80°以及100°三相位膝关节的力学状态，为了获得上述3个大角度屈曲位，选择ABAQUS显示动力学求解。屈曲模型中，对所有骨组织施加刚体约束。在股骨上建立屈曲轴线(FEA)，如图2所示。

将屈曲轴线竖直向下平移至胫骨平台，建立胫骨运动坐标系参考轴。为了对膝关节的运动进行完整的描述，依据Grood的结论[15]，以股骨屈曲轴中点、胫骨运动坐标系参考轴中点为坐标原点建立运动坐标系，将股骨坐标系原点移动到髌骨中心点，建立髌骨坐标系(见图2)。

将处于力线位置的股骨坐标原点作为载荷施加点，施加骑行时的最大载荷831 N[7]，由于屈曲运动，载荷方向也将随之变化，如图3。屈曲过程中膝关节胫骨髌骨相对股骨运动的定义参考相关文

献[4,16-17]。胫骨相对股骨发生的旋转运动(内旋、内收)及平移运动(前移、外移)的定义如表3所示，髌骨相对股骨的旋转(屈曲、外倾和内旋)及外移的定义如表4所示。

图2 膝关节坐标系Fig.2 Coordinate system of knee joint

图3 膝关节不同骑行相位模型。 (a)屈曲60°；(b) 屈曲80°；(c) 屈曲100°Fig.3 3D model of knee under different phases of riding. (a)60° flexion; (b)80° flexion; (c)100° flexion

屈曲/(°)内旋/(°)内收/(°)前移/mm外移/mm604.00.58.03.5804.51.012.03.81005.02.514.04.0

表4 膝关节髌骨相对股骨的旋转和平移Tab.4 Relative rotation and translation of patella vs femur

2 结果

得到不同相位的屈曲模型如图3所示。处于不同骑行位置的膝关节，其软骨上的von-Mises应力分布(见图4)及大小(见图5)也随屈曲角度不同而不同，在加载831 N时，屈曲60°位置软骨最大应力为4.57 MPa，发生在股骨软骨上(见图4(a))；屈曲80°时，膝关节软骨最大应力为6.042 MPa，发生在胫骨平台上(见图4(b))；屈曲100°时，膝关节软骨的最大应力为7.83 MPa，发生在股骨软骨上(见图4(c))。

在所关注的骑行阶段，软骨的von-Mises应力因屈曲角度不同而不同，屈曲角度增大时，软骨的应力值上升，股胫关节的应力分布由膝关节前侧向后侧转移，髌股关节应力分布向髌骨软骨上方转移。同时，根据计算结果发现，屈曲60°到屈曲100°过程中，膝关节内侧应力值逐渐增大，这可能是由于屈曲时胫骨相对于股骨发生的内旋与内倾运动造成的。

图4 不同骑行相位关节软骨应力状态(由左至右分别为股骨软骨、胫骨平台软骨和髌骨软骨)。(a)屈曲60°；(b)屈曲80°；(c)屈曲100°Fig.4 The distribution of von Mises in different postures (From the left to the right are femur cartilage, tibia cartilage and patella cartilage). (a)60°flexion; (b)80°flexion; (c)100°flexion

图5 骑行时不同相位软骨von-Mises应力值Fig.5 von Mises stress on cartilage in different phases

3 讨论和结论

在本研究中，建立完整膝关节几何模型，软骨厚度分布符合生理结构特征。建模前查阅膝关节的相关解剖资料，对其解剖学特点有全面的了解。建模后，对重要考察的部位如关节软骨、半月板的几何尺寸与相关文献的研究结果[8-9]做对比，软骨厚度及分布，半月板矢径横径及前后角高度均在标准尺寸范围内。

建立有限元模型，分析在骑行最低点以及骑行周期中膝关节合力最大区域的两个相位处关节软骨的应力状态。对有限元分析结果进行直接验证难以实现，本研究用所建立的膝关节模型模拟步行运动，得到步态周期中50%和75%两相位处膝关节胫骨软骨上的应力分别为5.0和7.9 MPa，这与文献[5]的6.1 和7.8 MPa非常接近。

较步行和爬楼梯等日常运动，骑行时膝关节负重小，但屈曲范围大。大角度屈曲时，胫股关节接触面积减小，软骨受力向后转移；同时，内侧软骨应力增大，外侧软骨应力减小。可见骑行时膝关节软骨的应力分布及大小不断变化，变化幅度和范围更大，频率更快。由此可知通过优化自行车的结构尺寸而改变膝关节屈曲范围，进而控制软骨应力，可以达到健康骑行的目的。另一方面，如果骑速过快，软骨易发生疲劳损伤，损伤位置分布在膝关节股骨软骨内侧后部的可能性更大。

随着健康产业的发展，更应重视骑行的舒适感，故本工作仍需进一步探索。首先，本课题只关注骑行过程中屈曲60°、80°和100°这3个具体位置膝关节软骨的应力分布，是一种静态研究，尚需对整个骑行过程做动力学模拟。其次，考虑不同骑行速度对膝关节生物力学特性的影响，休闲骑行与竞技骑行有较大差别。最后，文中对相关生物组织的材料属性进行了简化，今后的研究工作可以采用更实际的材料特性。虽然存在上述不足，对于了解骑行时膝关节软骨的应力分布以及关节软骨损伤与骑行之间的关系，本研究的结果仍具有参考作用。

[1] Takeo N, Chris OD, Eugene JA, et al. Mechanical loads at the knee joint during deep flexion [J]. Journal of Orthopaedic Research, 2002, 20(4): 881-886.

[2] 丁晨，王君泽，高瞻. 基于AnyBody的汽车驾驶中人体腿部逆向动力学仿真 [J]. 中国生物医学工程学报, 2013, 32(1): 124-128.

[3] Li Jingsheng, Ali H, Lucile C, et al. Kinematic characteristics of the tibiofemoral joint during a step-up activity [J].Gait Posture, 2013, 38(4): 712-716.

[4] Wei Qi, Ali H, Tsung YT, et al. In vivo kinematics of the knee during weight bearing in high flexion [J]. Journal of Biomech, 2013, 46(9): 1576-1582.

[5] Adouni M, Shirazi-Adl A. Evaluation of knee joint muscle forces and tissue stresses-strains during gait in severe OA versus normal subjects [J]. J Orthop Res, 2014,32(1):69-78.

[6] Wang Yuxing, Fan Yubo, Zhang Ming. Comparison of stress on knee cartilage during kneeling and standing using finite element models [J]. Medical Engineering&Physics, 2014, 36(4): 439-447.

[7] 郭艺华. 骑行时人车系统的动力学及相关有限元分析 [D]. 天津: 天津理工大学, 2015.

[8] 孙彩英，崔建岭，李石玲，等. MRI测量正常人膝关节软骨厚度 [J]. 实用放射学杂志, 2004, 20(11): 1007-1010.

[9] 徐伟，周庭永，钱学华，等. 膝关节半月板的三维断面解剖学研究 [J]. 中国临床解剖学杂志, 2010, 28(3): 272-276.

[10] Li Guoan, Lopez O, Rubash H. Variability of a three-dimensional finite element model constructed using magnetic resonance images of a knee for joint contact stress analysis [J]. Journal of Biomechanical Engineering, 2001, 123(4): 341-346.

[11] Yao Jie, Snibbe J, Maloney M, et al. Stresses and strains in the medial meniscus of an acl deficient knee under anterior loading: a finite element analysis with image-based experimental validation [J]. Journal of Biomechanical Engineering, 2006, 128(1): 135-141.

[12] Mesfar W, Shirazi-Adl A. Biomechanics of the knee joint in flexion under various quadriceps forces [J]. Knee, 2005,12(6): 424-434.

[13] Butler DL, Kay MD, Stouffer DC. Comparison of material properties in fascicle-bone units from human patellar tendon and knee ligaments [J]. Journal of Biomech, 1986, 19(6): 425-432.

[14] 王亦璁. 膝关节外科的基础和临床 [M]. 北京: 人民卫生出版社, 1999:20-22.

[15] Grood ES, Suntay WJ. A joint coordinate system for the clinical description of three-dimensional motions: application to the knee [J]. J Biomech Eng, 1983,105(2):136-144.

[16] Lu Tungwu, Tsai TY, Kuo Meiying, et al.Invivothree-dimensional kinematics of the normal knee during active extension under unloaded and loaded conditions using single-plane fluoroscopy [J]. Medical Engineering&Physics, 2008, 30(8): 1004-1012.

[17] Amis AA, Senavongse W, Bull AM. Patellofemoral kinematics during knee flexion-extension: an in vitro study [J]. Journal of Orthopaedic Research, 2006, 24(12): 2201-2211.

The Analysis of Stress on Knee Cartilage in Different Flexion Angles During Riding

Ma Xuemei1Zhang Chunqiu1*Gao Lilan1Ye Jinduo1Zhang Xizheng2

1(TianjinKeyLaboratoryforAdvancedMechatronicSystemDesignandIntelligentControl,SchoolofMechanicalEngineering,TianjinUniversityofTechnology,Tianjin300384,China)2(InstituteofMedicalEquipment,AcademyofMilitaryMedicalSciences,Tianjin300161,China)

knee joint; cartilage; riding; finite element model; flexion

10.3969/j.issn.0258-8021. 2016. 01.015

2015-08-03， 录用日期:2015-09-28

国家自然科学基金重点项目(11432016)；国家自然科学基金(11172208，11402171)

R318

D

0258-8021(2016) 01-0119-05

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