个体化骨肌多体动力学和有限元联合建模的肩胛骨锁定板生物力学评估方法

范勋健,陈瑱贤,曹卓,夏志辉,陈军,靳忠民

(1.中国工程物理研究院机械制造工艺研究所,621900,四川绵阳;2.西安交通大学机械制造系统工程国家重点实验室,710049,西安;3.长安大学公路建设技术与装备教育部重点实验室,710064,西安;4.空军军医大学第二附属医院骨科,710038,西安;5.西南交通大学机械学院,610031,成都)

目前,取出物分析和数值仿真是植入物临床前失效评估的重要手段。取出物分析成本高、周期长,数值仿真能有效弥补实验方法的局限性,便于参数化研究[1-2]。数值仿真中的有限元方法已广泛用于骨折板的生物力学分析,但由于体内载荷环境复杂,以往采用的由文献或者简单实验确定的简化边界条件不符合体内的生理载荷状况[3-4]。Behren等发现了单一的静态加载不足以表示人体运动过程的受力情况[5]。Duda等证实了仿真结果高度依赖边界条件,考虑全部肌肉力与仅考虑主要肌肉力计算的骨应变差异高达26%[6]。

人体骨骼肌肉多体动力学建模系统允许考虑患者骨骼和肌肉的个体差异,可以同时计算体内关节力和力矩、肌肉力和韧带力,为预测人体日常活动下的体内载荷提供了平台[7-8]。最近,研究人员采用骨肌建模方法联合有限元方法,调查人体下肢在步行和骑自行车运动中股骨[9]以及植入物(骨折板[10]和人工髋关节[11])的受力。具体来讲,采用多体动力学方法预测人体日常活动下的生理载荷,并将其作为有限元分析的边界条件计算植入物的应力应变。这种联合建模方法能够在人体日常活动时体内生理载荷下,评估植入物的生物力学性能,相比以往采用简化边界条件更符合人体真实的力学环境,为改善植入物的临床效果提供了途径。但是,这些研究都没有考虑个体的解剖差异。Gerus等的研究表明,相比采用通用的骨肌模型,个体化的骨肌建模可以将膝关节的内侧接触力预测精度提高47%,将外侧接触力预测精度提高7%[12]。Pellikaan等开发和验证了骨变形缩放技术,允许更好的缩放通用骨肌模型来匹配个体的骨解剖特征[13]。

肩胛骨在肩关节的运动中有重要作用,肩胛骨如果固定不动,上臂只能主动抬高至60°,因此丧失肩胛骨活动时,肩部活动至少减少1/3[14-15]。同时,肩胛骨为扁平状不规则骨,术中复位固定困难,选择合适的内固定物对肩胛骨体部碎片的复位及骨折愈合具有重要意义[15-17]。目前,肩胛骨体部骨折的内固定物存在以下3点不足[14-18]:①钢板非整体结构,难以对抗剪力;②没有合适的外型,需要手术预弯,增加了手术难度,降低了钢板强度;③通常采用盆骨重建板,厚度为3～3.5 mm,过厚。针对上述问题,本文拟设计一种专用于肩胛骨体部复杂骨折的一体式锁定钢板,该钢板具有以下优点:①采用整体结构抵抗剪力;②外型符合肩胛骨体部的解剖特点;③厚度满足肩胛骨生物力学要求。

针对临床上肩胛骨骨折板的现存问题,本文采用个体化骨肌多体动力学和有限元联合建模方法,建立患者个体化的骨肌多体动力学模型和肩胛骨骨折固定有限元模型,在体内生理载荷下调查肩胛骨体部骨折板的力学性能,从而为肩胛骨体部骨折板的设计提供理论依据,也为医生手术提供一定的参考。

图1 骨肌多体动力学和有限元方法联合建模技术路线图

1 骨肌多体动力学与有限元建模方法

本文采用骨肌多体动力学联合有限元方法调查肩胛骨骨折板的生物力学性能,技术路线如图1所示,具体步骤包括:①个体化骨肌建模,提取个体CT数据建立肩关节的三维实体模型,采用基于骨解剖点的骨变形缩放技术和肌肉募集准则,将Anybody通用上肢模型精确缩放成与中国人个体相匹配的上肢骨骼肌肉模型;②个体化骨肌动力学建模,采集并处理个体的三维运动学数据,并将其输入骨肌模型,通过多体逆向动力学方法计算上臂抬高运动中肩膀的肌肉力、关节力和力矩;③肌肉力和关节力输出转换,编程提取骨肌动力学软件计算的附着于肩胛骨上的全部肌肉力和关节力,输出转换为能被有限元软件采用的边界条件;④肩胛骨骨折固定有限元建模,建立与个体相对应的肩胛骨骨折固定有限元模型,将所有附着于肩胛骨上的肌肉力和关节力加载至有限元模型,在体内生理载荷下对比一体式骨折板与传统重建板的力学性能。相比以往的植入物失效预测方法,本文方法采用患者体内真实的生理载荷评估骨折板的生物力学性能,结果更加可靠。本文建立的个体化骨肌多体动力学模型和肩胛骨骨折固定有限元载荷模型见图2,图2b中各点为载荷加载点标识,有限元模型中所有肩胛骨受力的加载点标识、最大幅值和肌肉的生理载荷面积见表1。

2 个体上肢骨骼肌肉多体动力学分析

2.1 实验数据

与空军军医大学第二附属医院合作采集了一例中国成年男性的相关实验数据进行建模和研究。该男子年龄26岁、身高1.73 m、体重65 kg、胸腔宽0.278 m、上臂长0.282 m、前臂长0.276 m。采集的实验数据包括个体右侧肩胛骨和肱骨的CT数据(层厚0.6 mm、矩阵512×512,视野385 mm×385 mm)、个体右臂匀速外展和前屈的三维运动学数据。以上数据由西安交通大学机械制造系统工程国家重点实验室VICON三维运动捕捉系统采集。

2.2 个体化上肢骨肌多体动力学建模

2.2.1 骨肌模型的个体化精确缩放 在Anybody 6.0.7软件中,通过修改软件提供的通用模型(模型库版本1.6.3,代表欧洲人的平均身体尺寸)来建立个体化上肢骨肌多体动力学模型[19]。采用的Anybody上肢骨肌动力学模型已经过文献的广泛验证,并用于预测日常活动下肩膀的肌肉力和关节力[20-23]。利用Mimics 16.0和Geomagic 13.0软件反求CT数据中个体肩胛骨和肱骨的三维实体几何模型,结果如图3所示。

表1 肩胛骨受力的加载点标识、最大受力幅值和肌肉生理载荷面积

注:Fx、Fy、Fz分别为x、y、z方向的最大受力。

(a)上肢骨肌多体动力学模型 (b)肩胛骨骨折固定有限元载荷模型图2 个体化骨肌多体动力学与有限元分析联合的建模方法示意图

(a)CT数据 (b)三维实体几何模型图3 个体肩胛骨和肱骨三维实体建模

采用Meshlab 1.3.2软件在个体肩胛骨和肱骨几何表面确定骨解剖点,在Anybody通用模型中的肩胛骨和肱骨表面选取对应的骨解剖点。根据选取的骨解剖点信息和骨几何精确缩放技术,将通用模型右侧的肩胛骨和肱骨进行个体化缩放。首先,通过仿射变换将通用模型的骨骼(待缩放骨)粗略缩放成个体的骨骼几何(目标骨);然后,基于非线性径向基插值函数对骨骼几何进行精确缩放,同时定义边界函数优化缩放效果;最后,基于骨解剖点和反向刚性转换,将缩放后骨的参考坐标系由个体CT数据的坐标系转换为Anybody骨肌建模系统中的坐标系。缩放中使用的非线性径向基插值函数缩放方程[24-26]为

(1)

式中:f表示缩放后的模型解剖点;φ表示非线性径向基插值函数;cj是φ的系数,根据目标骨和待缩放骨的骨解剖点计算得到;‖y-xj‖表示目标点与中心点之间的距离;q是多项式p的阶数。个体化精确缩放后的肩关节骨骼模型如图4所示。

图4 个体化肩关节精确缩放

由于肌肉通过起始点和插入点依附在骨骼上,进行骨骼缩放的同时,骨骼上的肌肉附着点和肌肉形状也需要根据骨骼几何进行调整,从而更符合个体的真实解剖。肌肉的缩放遵循线性几何缩放定理[27],公式为

s=Sp+t

(2)

式中:s是被缩放后的点的位置向量;S是缩放矩阵;p表示初始点的位置;t是平移向量。

对模型中未缩放的骨骼体节按照长度-质量-脂肪缩放准则进行缩放,依据个体身高、体重和脂肪比例调整骨骼长度和肌肉强度。身高体重相似的两个人,脂肪比例高的人肌肉强度低,人体的脂肪比例F的具体计算公式[28]为

(3)

式中:B表示身体质量指数;M表示人体体重;H表示人体身高。在Anybody骨肌建模系统中,计算每束肌肉的强度时都会乘以表1中对应的生理横截面积。

2.2.2 运动学数据采集及处理 Anybody通用模型符合欧洲人的运动习惯,利用Vicon三维运动数据采集系统获得个体右臂外展90°和前屈45°的运动学数据,采集过程见图5。将采集的个体运动数据以C3D格式导入精确缩放后的个体骨肌模型,进行运动学参数优化,并计算上肢各关节的运动角和手位移等数据,作为逆向动力学分析的输入。

图5 个体右臂抬高时的三维运动数据采集过程

2.2.3 个体右臂抬高运动逆向动力学分析 将运动学分析计算的关节角等数据结合牛顿第二定律与肌肉募集准则进行逆向动力学分析,计算得到的肩关节力和力矩如图6所示,可以看出:上臂外展运动中,盂肱关节力和力矩在外展90°时达到最大值,分别为428 N和10.3 N·m,盂肱关节力最大值大约等于人的体重;上臂前屈45°运动中,盂肱关节力和力矩的最大值分别为379 N和7.4 N·m。预测的盂肱关节力和力矩与实验测量结果[20-23]一致。

(a)外展和前屈运动的盂肱关节接触力

(b)外展和前屈运动的盂肱关节力矩图6 右臂外展90°和前屈45°的盂肱关节接触力和力矩

3 个体肩胛骨骨折固定生物力学评估

3.1 肩胛骨骨折固定实体建模

肩胛骨骨折固定的三维装配模型在SolidWorks 2013中完成,断裂的肩胛骨分别由传统重建板固定和新设计的一体式骨折板固定,建立的三维模型如图7所示。依据肩胛骨解剖和骨折手术治疗原则,新设计骨折板的外形近似三角形,该设计由资深的骨科临床医生确定。钢板的设计参数包括:螺钉直径5 mm、螺孔间距13 mm、钢板宽度10 mm、钢板厚度3 mm。重建板的设计参数依据文献[29]进行设定:螺钉直径5 mm、螺孔间距13 mm、钢板宽度10 mm、钢板厚度3 mm。骨折形式由骨科医生确定,模拟肩胛骨体部粉碎性骨折,钢板和螺钉的植入位置及数量由医生确定。

(a)传统重建板 (b)一体式骨折板图7 采用重建板和一体式骨折板的肩胛骨骨折固定三维模型

(a)肩胛骨骨折固定 (b)载荷模型 (c)皮质骨和松质骨图8 肩胛骨骨折固定有限元模型

3.2 肩胛骨骨折固定有限元建模

采用Abaqus 6.14建立肩胛骨骨折固定有限元模型,结果如图8所示。肩胛骨模型考虑了皮质骨和松质骨,根据文献[3],肩胛骨皮质骨的厚度为0.6 mm,本文采用从骨表层向内偏置0.6 mm模拟皮质骨,剩余部分为松质骨,骨质参数见表2,E为弹性模量。相比以往采用向外偏置模拟皮质骨,内加皮质骨的方法确保了骨头真实的几何维度,保证了仿真结果的准确性。骨折板和螺钉采用钛合金(TiAl6V4)材料,该材料具有良好生物相容性,弹性模量为115 GPa,泊松比为0.3[30]。

表2 骨质参数

通过HyperMesh 13.0建立肩胛骨骨折固定网格模型。骨折板和螺钉的网格大小由网格敏感性分析确定,分别选择0.5 mm、1 mm和1.5 mm的网格尺寸。网格敏感性分析表明:与0.5 mm网格相比,采用1 mm网格的最大范氏应力相对偏差小于5%,而1.5 mm网格的偏差大于5%;采用1 mm网格的计算时间是0.5 mm网格的一半。综合考虑计算精度和效率,确定骨折板和螺钉的网格大小为1 mm,对应单元数分别为15 765和15 356。肩胛骨的网格敏感性分析显示:网格大小从1 mm增加到2 mm,最大范氏应力的变化小于5%。为了缩短计算时间,采用2 mm网格,皮质骨和松质骨的网格单元数分别是17 432和112 345。

有限元载荷采用个体骨肌模型计算的右臂外展和前屈运动中体内的生理载荷,包括盂肱关节力和力矩以及所有附着于肩胛骨表面上的225个肌肉力。肌肉力和关节力及力矩分别作用在肌肉附着点和关节点上,为了避免应力集中,每个肌肉附着点和关节点(见图2)与肩胛骨表面临近区域相耦合。根据文献[1],肩膀运动时肩峰相对固定,故在肩峰处设定全约束。

根据文献[3]、文献[30-31],断裂骨段间和骨板与骨界面间设置为面面接触,摩擦因数分别为0.15和0.4。因为锁定钢板的螺钉能够紧固连接,螺钉与螺孔和骨头间采用绑定连接。

3.3 肩胛骨骨折板的生物力学分析

为了评估骨折板的力学性能,采用最大范氏应力作为量化指标,选取传统重建板作为优化的参照标准,通过与重建板进行对比分析来考察设计的有效性,完成一体式骨折板的厚度和宽度设计。一体式骨折板和重建板的初始厚度为3 mm、宽度为10 mm,本文先后调查了两种骨折板在右臂外展90°和前屈45°运动下的生物力学性能。在骨折板的设计优化过程中,遵循控制变量原则。

在右臂外展90°运动下进行骨折板的设计,骨折板和螺钉的范式应力云图如图9所示。对于骨折板的厚度设计,可以分析得出:3 mm厚一体式骨折板的最大范氏应力(106 MPa)相比3 mm厚重建板的最大范氏应力(128 MPa)减小了17.2%,表明一体式骨折板的初始设计有效;最大范氏应力的位置出现在肩胛骨关节盂附近,说明外展运动中骨折板在关节盂处容易发生断裂,也证实了之前的临床观察[32]和数值研究[33-34]结果;2.6 mm厚一体式骨折板的最大范氏应力(118 MPa)相比3 mm厚重建板的最大范氏应力(128 MPa)下降了7.8%,表明一体式骨折板的厚度优化有效;最大范氏应力位置仍出现在肩胛骨关节盂附近。

(a)3 mm厚10 mm宽 (b)3 mm厚10 mm宽重建板 一体式骨折板

(c)2.6 mm厚10 mm宽 (d)3 mm厚6 mm宽一体式骨折板 一体式骨折板图9 外展90°时重建板和一体式骨折板的范氏应力云图

对于骨折板的宽度设计,分析可以得出:10 mm宽一体式骨折板的最大范氏应力(106 MPa)小于10 mm宽重建板的的最大范氏应力(128 MPa),表明一体式骨折板的初始设计有效;6 mm宽的一体式骨折板的最大范氏应力(245 MPa)大于10 mm宽重建板的最大范氏应(128 MPa),表明6 mm宽度一体式骨板的力学性能相比重建板差;同时,最大范氏应力位置发生了变化,大约出现在肩胛骨体部居中的位置。由此可知,右臂外展90°运动时,厚度2.6 mm、宽度10 mm的一体式骨折板为目前最优设计。

在右臂前屈45°运动下,将外展运动中获得的一体式骨折板的有效设计方案与重建板进行对比分析,骨折板和螺钉的范氏应力云图如图10所示,可以分析得出:3 mm和2.6 mm厚的一体式骨折板的最大范氏应力分别为96 MPa和113 MPa,与传统重建板的最大范氏应力(142 MPa)相比分别下降32.4%和20.4%。综上,右臂外展90°和前屈45°运动时,厚度2.6 mm和宽度10 mm为一体式骨折板的理想的参数组合。

(a)3 mm厚 (b)3 mm厚 (c)2.6 mm厚10 mm宽 10 mm宽一 10 mm宽一重建板 体式骨折板 体式骨折板图10 前屈45°时重建板和一体骨折板的范氏应力云图

同时,本文考察了重建板和一体式骨折板(厚2.6 mm、宽10 mm)的骨折固定模型在外展90°和前屈45°时的肩胛骨综合位移,对比分析了两种固定方式的稳定性,并与正常肩胛骨的位移进行比较,结果发现:在外展90°时,重建板、一体式骨折板固定和正常肩胛骨的最大综合位移分别是1.14 mm、1.07 mm和0.86 mm;在前屈45°时,重建板、一体式骨折板固定和正常肩胛骨的最大综合位移分别是2.17 mm、2.11 mm和1.89 mm。分析可知,一体式骨折板相比重建板在固定复杂肩胛骨骨折时,其一体式结构具有更好的抗弯和抗扭性能,能够对骨折断端提供更加稳定的把持力和更牢靠的固定[35-36]。

4 结 论

在人体生物力学数值研究中,真实的边界条件是确保仿真结果准确可靠的关键因素[5-6],但是之前的研究多采用简化的边界条件[3-4]。本文建立的个体化骨肌多体动力学和有限元分析联合建模方法,允许在个体患者体内生理载荷状况下评估植入物的力学性能,确保仿真结果更加真实合理。本文得出的主要结论如下。

(1)建立了个体化上肢骨肌多体动力学和肩胛骨骨折固定有限元联合模型,应用于调查上臂不同抬高运动体内载荷状况下一体式肩胛骨骨折板和传统重建板的生物力学性能差异。

(2)2.6 mm厚10 mm宽一体式骨折板相比3 mm厚10 mm宽重建板具有更好的生物力学性能,外展和前屈运动中骨折板的最大范氏应力分别下降了7.8%和20.4%。

(3)与传统植入物失效评估方法相比,该方法允许考虑个体的解剖学差异和体内生理载荷状况,能够在个体日常活动体内载荷状况下调查植入物的生物力学性能,为假体设计、手术技术和术后康复提供了理论依据。