侧向撞击时肩部有限元分析

张 晶，郭 媛，张绪树，张玉文

(太原理工大学 应用力学与生物医学工程研究所，太原 030024)

肩关节是人体结构中最复杂、活动度最大的关节之一[1]。肩部损伤常发生在各种碰撞中，无论是汽车行驶安全方面还是骨科领域均将其视为一种严重的关节损伤。由于肩部经常是乘客和车厢或气囊等约束系统作用时的一个主要载荷路径，肩部的响应会对身体其他部分造成影响，传向胸部的载荷直接决定于肩部的响应[2-4]。研究肩关节在侧向碰撞工况下的响应对预测碰撞对乘客的影响有重要意义。

近年来，采用有限元等数值方法来模拟肩关节的力学行为的研究受到普遍关注[3]。HAUG et al[5]构建了完整的运动人体肩关节刚体有限元模型，该模型采用非线性弹簧模拟肩部关节的连接，可用来研究上肢受力和关节扭矩。INGRAM D[6]用有限元方法对肩部肌骨系统进行了研究，分析了肩关节形态学结构的功能，所建模型可再现肩部的动态行为。DUPREY et al[7]根据COMPIGNE et al[8]尸体滑车实验建立模型，研究了肩部在不同冲击力作用下的力学响应并与实验进行了对比。虽然这些模型和尸体实验进行了验证，但存在大量简化，如均未考虑软骨和肌肉的影响，同时模型网格较为粗糙，未能精确反应肩关节几何特征，影响仿真结果的准确性。国内杨晓霞[9]构建了肩胛骨、锁骨的三维有限元模型，研究了盂肱关节在4种外展角度工况下肩锁关节的应力分布情况。杨志勇等[10]建立了锁骨的三维有限元模型并进行轴向加载，发现锁骨应力主要集中在锁骨干中约1/3区域，且下表面应力集中较上表面明显。禹铭杨[11]建立了不包括肌肉的肩部模型，研究了肩部形态学参数对肩峰撞击症的影响。国内研究主要集中于用静态加载来模拟肩部各种受力工况，且没有关于肩部侧向撞击的报道。

本文旨在建立肩部的有限元模型，模拟肩部在侧向冲击作用下的力学响应，为预测汽车碰撞或其他肩部侧向冲击情形下的肩部损伤提供参考。

1 材料和方法

1.1 几何模型的建立

对一位肩部无畸形、病变成年男性志愿者的左侧肩部进行CT扫描，扫描层厚0.5 mm，得到CT断层扫描图像，以DICOM格式保存。

在Mimics中新建工程，导入包含左侧肩部的CT扫描数据，设定分割阈值生成蒙板。对蒙板进行裁剪、编辑、缩小范围等操作，使锁骨、肩胛骨、肱骨近端、肩部软组织能明显区分，在完成光滑、开闭等形态学操作之后进行三维重建。对三维模型进行光滑、三角形缩减、包裹等操作，以提高其质量[12-14]。

根据韧带解剖位置起止点，在骨模型表面选取相应位置重建韧带，考虑到韧带对骨施加面载荷，使用一组弹簧单元模拟韧带[14]。最终在Abaqus中重建出肩部主要韧带：喙锁韧带、肩锁韧带、喙肩韧带、盂肱韧带[15]。

基于相邻骨之间的空隙均由软骨填充，且软骨厚度以及空隙最小间距的一半作为假设，建立软骨几何模型[16]。将上述三维模型导入3-Matic中，根据软骨解剖位置在相应骨表面标记相应面，将标记面从骨面上分离后，沿外法线方向偏移相应距离生成软骨，对实体进行倒圆角、光滑等操作以提高其几何质量。

侧向撞击钢板在Abaqus中创建为边长200 mm的正方形壳。

1.2 三维有限元模型的建立

将所建立的实体模型载入3-Matic中，生成体网格。因网格质量对计算结果有重要影响，所以在生成体网格之前进行了过滤锋锐三角形、减少三角面片数量、网格质量检测、自动面网格生成等操作，使网格质量满足计算要求[12]。肩部有限元模型参数见表1，完整的模型如图1所示。

表1 肩部有限元模型参数Table 1 Parameters of the shoulder finite element model

图1 肩部侧向碰撞的有限元模型Fig.1 Finite element model of shoulder lateral impact

1.3 材料赋值

由于软件和实验数据的局限性，本文涉及的生物材料属性均假设为连续、均质、线弹性。将划分好网格的有限元模型以INP文件格式导入到Mimics中，基于CT图像的灰度值，使用经验公式以10种材料属性对骨性部分进行赋值[17]。根据参考文献[18-19]，对锁骨、肩胛依据经验式(1)和(2)赋值：

D=1.067×VG-961.608 .

(1)

E=0.01×D1.86.

(2)

对肱骨依据经验式(3)和(4)赋值：

D=0.624×VG-465.976 .

(3)

E=0.32×D1.41.

(4)

式中：D为材料密度，kg/m3；VG为CT图像的灰度值；E为材料弹性模量，MPa.

肌肉、软骨、钢板的材料属性见表2.韧带设置为只能承受拉力的非线性弹簧单元，肩部主要韧带参数见表3.

表2 材料属性[20]Table 2 Mechanical properties of materials

1.4 接触、载荷和边界条件设置

软骨、软组织与所附着骨面之间的接触均设为固定约束，以保证它们之间没有相对位移。使用Abaqus中的面-面接触算法模拟关节与软骨之间的接触，由于人体中关节面之间存在润滑，因此将关节面之间设为无摩擦的接触关系。同样地使用面-面接触算法模拟肩部软组织与钢板的接触，其中法向采用硬接触，切向接触设置为摩擦系数为0.6的库伦摩擦[21]。

将钢板置于yoz平面，并将其6个自由度完全约束。调整肩与钢板的相对位置，使肩部软组织外侧最大程度接近钢板，以缩短计算时间。肩部施加沿x轴方向的速度场，大小为3.5 m/s，分析时间设置为0.02 s.

根据运动生物力学知识可知，肩胛骨活动受到肩胛胸壁关节的限制，使肩胛骨在运动中和胸壁紧密切合，防止了翼状飘起，近似可认为其沿肩胛胸壁方向无线位移，故本文限制了肩胛骨前骨面上部分节点沿y轴方向的线位移。锁骨由胸锁关节固定于第1肋骨上，锁骨近端受到胸锁关节的限制使其旋转运动较多而线位移较少，模拟中限制锁骨近端z轴方向的线位移来代表胸锁关节对其的约束[9,15,22]。整个动态碰撞求解过程在Abaqus显示求解器中完成，耗时29 h.求解完成之后进入Abaqus可视化结果处理页面，对结果进行分析。

模型验证参照杨晓霞[9]肩锁关节力学测试实验的固定和加载方式(如图2)，将肩胛骨与胸壁接触一侧的部分节点进行固定约束，分别用70,80,90,100 N的拉力沿x轴负向作用于锁骨近端，使用Abaqus求解器中的隐式方法完成模拟。选取如下测量点：点1为锁骨外端中点，点2为锁骨外端中点在肩峰上对应点(内外方向)，点3为喙突尖顶点。测量标志点的位移，将标志点位移曲线与实验结果对比。

图2 沿x轴横向拉力加载的实验和模拟Fig.2 Experiment and simulation of x axial tension

2 实验结果

2.1 静态加载位移情况

施加横向拉力下目标点的位移结果见表4.从模拟和实验的对比图(图3)可见，随着拉力的增加，目标点位移也在增加，喙突尖顶点位移最大，而锁骨外端中点在肩峰上对应点位移较小，位移整体趋势与实验结果一致。因本模型材料属性均设定为线弹性，故位移结果表现出一定的线性特征。

表4 沿x轴横向拉力加载方式下目标点的位移Table 4 Displacement of target point under x axial tension load

图3 沿x轴拉力加载下目标点位移Fig.3 Displacement of target point under x axial tension force

2.2 侧向冲击后肩部整体受力情况

本文研究的碰撞过程历时20 ms：0～12 ms时肩部与墙面发生接触；13～20 ms时肩部进入无外力作用，即变形恢复阶段。模拟展现了碰撞的整个过程。

碰撞中von Mises应力最大值发生在11 ms时，大小为70.12 MPa，位置在肩胛颈内侧，这与临床观察肩胛颈易发生骨折位置相符[23]。整个碰撞过程中最大von Mises应力变化如图4所示。结合整个碰撞的动态过程可知，最大von Mises应力起初出现于肩部肌肉，肩部肌肉的变形引起肱骨的受力变化，最大应力由肱骨向肩胛骨转移；3～11 ms最大应力主要发生在肩胛骨上，这反映出在碰撞中肩胛骨起着主要承担载荷的作用。碰撞作用的后期，肩部反弹与刚性面脱离接触，肩部储存的能量释放，肩部各组成部分之间相互作用，最大应力发生在肱骨和锁骨远端软骨上。

图4 最大von Mises应力-时间曲线Fig.4 Max von Mises stress versus time curve

从图5最大法向接触力-时间曲线可看出，肩部与刚性墙接触后，法向接触力开始上升，此时肩部速度在下降，肩部肌肉变形压缩，在3 ms时法向接触力达到峰值12.77 kN，这表示软组织被最大程度压缩，肩部各组织大部分区域速度降为0.随后肩部反弹，法向接触力开始下降；在7 ms时法向接触力下降到低谷后又开始上升，从碰撞运动过程推测这可能是因为肱骨和肌肉沿x轴正向恢复变形速度大于肩部沿x轴负向反弹速度，肩部和刚性墙再次发生了接触压缩，从而使接触力上升。在13 ms接触力下降为0，意味着肩部肌肉最下层和钢板脱离，冲击结束。

图5 最大法向接触力-时间曲线Fig.5 Max normal contact force versus time curve

2.3 侧向冲击后肩部各骨骼受力情况

2.3.1 肱骨骨应力分布

整个碰撞中肱骨所受最大等效应力发生在15 ms，大小为47.15 MPa.从图6肱骨截面图可见最大应力位于内表面，并且在整个碰撞中肱骨应力集中主要发生在其内部。碰撞中应力逐渐由肱骨干移向肱骨近端，随着力的传递，应力集中在颈干交界部位和肱骨干上部，在其内侧产生聚积。同时，与关节盂相接触的肱骨关节面部分，应力也在逐渐增加，但肱骨头上的应力明显大于肱骨干，据此推测肱骨头和肱骨颈干交界部位更易发生骨折。

图6 最大应力发生时肱骨应力云图Fig.6 Von Mises stress of humerus at the time of maximum stress

2.3.2 锁骨应力分布

锁骨上的最大等效应力发生在14 ms，大小为13.96 MPa，整个碰撞过程中锁骨上应力分布较为均匀。从14 ms时的锁骨应力云图(图7)可见锁骨下表面应力集中更明显，最大应力发生在锁骨干中部约1/3处，且近端较远端更易发生应力集中。同时观察应力云图发现，下表面比上表面更易发生应力集中，说明撞击中锁骨下表面断裂风险更高。

图7 最大应力发生时锁骨应力云图Fig.7 Von Mises stress of clavicle at the time of maximum stress

2.3.3 肩胛骨的应力分布

整个碰撞过程中最大应力发生在肩胛骨上，在最大应力发生时刻，肩胛骨上的应力分布如图8所示。肩胛颈处有明显的应力集中，且整个过程中肩胛颈相对于周围的冈下窝、肩胛冈均有较大的应力值和更明显的应力集中现象；表明力经肩锁和盂肱关节传导时，肩胛颈易发生骨折，与WILLIAMS et al[24]统计的大多数肩胛骨骨折为外科颈骨骨折的结论一致。因此在肩部骨折固定时，应首先考虑肩胛颈的完整性。

图8 最大应力发生时肩胛骨应力云图Fig.8 Von Mises stress of scapula at the time of maximum stress

3 讨论

近年来，有限元分析法在生物力学研究中的应用日益广泛和深入，并且在处理复杂生物力学问题时它的优势得到体现，已成为分析和预测骨应力、应变的重要工具。采用有限元分析法，肩关节的结构形状、组织性能、载荷边界条件等均可用数学形式概括，可通过改变其中任一参数以观察其对整体结构的影响，从而解释损伤中的力学变化[1,6,20]。现在人体生物力学有限元模型的建立大多依靠CT等影像学图像和Abaqus等有限元软件[1]，这样一方面可以增加模型模拟的准确性，使模型更趋于真实；另一方面简化了建模过程，节省人力物力。

本文基于逆向工程原理建立了人体左肩部复合体，在Mimics中将CT扫描所得连续断层图像进行重建、分割、优化处理，建模过程中始终保证软件对数据的自动识别优化。模型采用了骨广泛的力学参数，不再将骨简单区分为松质骨和密质骨。根据张国栋等[17]的研究，将骨有限元模型赋予10种材料属性后可准确反映其真实应力状态，因此本文骨的材质同样设为10种。本文所建模型包括了从CT图像中提取的软组织，更符合实际的侧向撞击情况，体现了软组织对冲击能量的吸收作用。肩部实体模型均划分为4节点四面体单元，韧带用非线性弹簧单元代替，整个模型共135 677个节点，516 392个单元，高度精细划分。模型的建立较好地体现了肩部的复杂几何形状和骨密度的不均匀性，保证了模型从数据源上的准确性和有效性。

静态加载模拟所得结果与杨晓霞[9]所测结果吻合，证明所建模型合理有效。动态碰撞所得结果表明碰撞中肩胛颈骨折风险最大，肱骨头、肱骨颈、锁骨下表面均为损伤易发区，由此为肩部碰撞救治损伤位置的确定提供了指导。同时最大等效应力-时间曲线和接触力-时间曲线也可为肩关节护具的开发提供参考，护具设计可围绕降低Von Mises峰值应力和接触力来进行设计。

需要指出的是，因缺少肩部侧向撞击的实验数据，所以所建模型只能通过对比部分结构静态加载位移进行验证，这显然未能完整显示出模型与真实肩部的接近程度，今后需要进行力、位移以及其他参数的综合比对。侧向撞击速度选择COMPIGNE et al[8]统计的肩部侧向撞击平均速度3.5 m/s，设置为纯侧向碰撞，实际碰撞中可能有更高的速度且方向多变，今后需要进行更高速度不同方向碰撞结果的对比。本文的结果对于了解肩部侧向撞击过程中各组织的生物力学响应具有一定的学术价值，对于肩部撞击的预防保护以及撞击后的临床治疗具有一定的参考价值。

4 结束语

本文采用组织逆向重建方法，利用Mimics、3-Matic、Abaqus等软件构建了人体左肩部各组织三维有限元模型，并验证了模型的有效性。纯侧向撞击结果表明，肩部骨折固定时，应首先考虑肩胛颈的完整性，肱骨头和肱骨颈干交界为肱骨骨折高风险区，锁骨干中部约1/3处应力集中解释了该处易发骨折的原因。今后还需完善模型，尝试以实体单元来构建韧带和使用超弹性材料属性，进一步研究不同速度侧向碰撞对肩部以及上躯干的影响。