放射式冲击波在软组织中传播特性的有限元分析

2021-08-11刘慧力王照智

科学技术与工程 2021年20期

赵彻，刘慧力，王照智，冯春

(1.常州工学院航空与机械工程学院，常州 213032； 2.吉林大学工程仿生教育部重点实验室，长春 130025; 3.沈阳工业大学机械工程学院，沈阳 110870)

随着人口老龄化的深入，关节炎和肌腱炎等软组织慢性损伤性疾病呈高发态势，已成为康复医学研究的重点课题[1]。与传统的药物和按摩治疗相比，放射式体外冲击波治疗(radial shock wave therapy，rESWT)具有易于操作、安全和治疗费用低等特点，在软组织损伤性疾病的治疗方面已得到广泛认可[2]。大量的临床研究已经证明，rESWT在足底筋膜炎[3]、肩钙化肌腱炎[4]和慢性跟腱炎[5]等多种软组织损伤性疾病的治疗中疗效显著。

然而，rESWT修复软组织损伤、缓解疼痛的具体机制尚不明确，治疗方案的制定主要依靠临床医生的经验和患者的主观感受。Walewicz等[6]采用随机和单盲试验评估了rESWT对下腰痛患者的影响，表明与传统的运动改善计划相比，使用rESWT对减轻疼痛和改善总体功能状态具有长期的重大影响。Zhang等[7]将rESWT与药物在治疗慢性盆腔疼痛综合征的效果进行了对比，表明rESWT的治疗效果更佳，症状复发率也更低。然而，要想制定更合理的治疗方案，开发更先进的治疗仪，仍需要对病患目标部位的生物学反应与特定治疗剂量(压强和能流密度等)间的关系进行量化分析[8-9]。

为解决上述问题，首先要了解放射式冲击波的传播特性。近年来，中外科研人员围绕这一问题开展了研究工作。Perez等[10]使用光纤探头水听器对ESWT体外冲击波治疗(extracorporeal shock wave therapy，ESWT)设备在水中产生的声场进行表征，得出焦点处的峰值正压范围为2～45 MPa。Payne等[11]使用聚二甲基硅氧烷(polydimethylsiloxane，PDMS)胶垫模仿人体软组织，对冲击波在其中的传播数据进行了测量，结果表明介质的力学性能对冲击波的传播特性具有显著影响。Liu等[12]使用柔性膜传感器分别对放射式冲击波在金属薄片、水和动物组织中的传播数据进行了测量，结果表明不同的测量介质对冲击波的传播数据影响较大。

综上所述，为降低测量难度和规避医学伦理风险，科研人员主要应用水或硅胶等标准介质代替软组织测量放射式冲击波的传播数据，其结果与体内的真实情况差距较大。此外，目前仍缺少量化冲击波在非均匀介质中三维传播数据的技术手段。因此，亟需新的理论与方法来解决这一难题。

有限元分析法已广泛应用于软组织的生物力学分析，尤其适合研究常规手段无法观测到的微观力学响应。李飞等[13]通过数值模拟接触分析研究了髋臼周围截骨术中动步态的软骨层接触压力和接触面积，表明髋臼周围截骨术后能明显的改善软骨层接触压力分布，为髋臼周围截骨术提供生物力学参考数据。杨得富等[14]基于有限元分析方法，建立耦合电磁场和生物传热过程的微波消融模型，研究了单缝同轴微波消融天线在多种肿瘤中的SAR值、消融区域温度分布及消融灶形态，为提高肿瘤消融治疗效果提供了理论支持。王方等[15]采用人体胸部有限元模型，基于胸腔骨骼和软组织的多个材料参数，进行正面冲击载荷下的胸部碰撞模拟，分别输出主要的胸部响应参数值，并采用标准多项式拟合方法详细对比分析了胸部材料参数发生变化时胸腔的不同响应。

为此，将有限元分析法应用于放射式冲击波在软组织中传播特性的表征与分析中：建立包含皮肤、脂肪，肌肉和骨骼的数值计算模型，揭示放射式冲击波在多层软组织中的传播特性，获取不同目标位置的压强和能流密度等关键信息，驱动气压、治疗探头与皮肤的耦合方式，以及治疗探头几何形状等因素对冲击波作用效果的影响。相关成果将为解决放射式冲击波在软组织中传播特性难以量化分析的难题提供新的解决方案，也将为患者提供针对性更强、更加安全有效的治疗方案，以及新型rESWT设备的研发，提供了重要参考。

1 建模与仿真

1.1 仿真模型的建立

如图1所示，rESWT的原理是子弹在压缩空气的推进下在弹道中以高速撞击治疗探头并产生放射式冲击波，进而传入人体软组织。根据这一原理，仅对子弹高速碰撞治疗探头的瞬间进行建模仿真。在CATIA V5中完成子弹、治疗探头(D20型)、保护壳，密封圈和阻尼器的建模，并组成装配体，保持子弹和治疗探头同心。其中，子弹直径为6 mm，高度为16 mm，治疗探头的最大直径为20 mm，凸面高度为0.7 mm。之后，建立包含皮肤层、脂肪层，肌肉层和皮质骨的多层生物组织模型，作为放射式冲击波治疗仪的作用对象。由于模型具有对称性，只取其1/2部分导入有限元分析软件ABAQUS中。

图1 放射式体外冲击波治疗仪的基本原理

为提高计算效率，降低沙漏效应和体积自锁效应的影响，采用8节点六面体网格、线性缩减积分、沙漏控制单元(C3D8R)对模型进行网格划分。经过网格的收敛性检验，并综合考虑计算效率和精度，除治疗探头与软组织接触的部位采用网格加密处理外，其余网格单元尺寸均设定为0.8 mm。有限元模型整体如图2所示，共包含349 628个网格单元。

图2 有限元模型

1.2 赋予材料属性

在放射式冲击波治疗仪中，子弹和治疗探头均为304不锈钢材料(杨氏模量2.1×1011Pa，泊松比0.3)。密封圈和阻尼器由硅橡胶材料制成，其力学性能可用超弹性材料模型Mooney-Rivlin中的3个常量C10(1.933×106Pa)、C01(0.483×106Pa)和D1(1.0×10-9Pa-1)表征。为节约计算资源，将保护壳的材料属性设置为刚体。根据生物软组织力学性能的相关研究结果[16-17]，皮肤、脂肪和肌肉可用一阶Ogden橡胶材料模型表征，特征常量为μ和α，具体参数如表1所示。相对刚性的皮质骨简化为线弹性模型(杨氏模量1.2×1010Pa，泊松比0.3)。

表1 软组织的力学性能参数

1.3 边界条件设置与数值计算方法

在ABAQUS软件中对有限元模型施加载荷和边界条件。根据流体动力学试验结果[18]，在预定义场中设定不同驱动气压作用下子弹的初速度。子弹与治疗探头，治疗探头与密封圈和阻尼器之间均采用表面-表面法向硬接触。为研究耦合方式对放射式冲击波的影响，采用两种方式设置治疗探头与皮肤层间的接触类型；一种是绑定约束，用来模拟皮肤表面涂有超声波耦合剂，与治疗探头完全贴合的情况。一种是表面-表面法向硬接触，用来模拟治疗探头直接作用与干燥的皮肤表面。用绑定约束模拟皮肤、脂肪、肌肉和皮质骨间的粘连状态。最后，在ABAQUS的Dynamic/Explicit求解器中进行数值计算，最小分析步设定为1.0×10-8s。

2 仿真结果与分析

2.1 放射式冲击波在软组织中的传播规律

图3为驱动气压为0.3 MPa时D20型治疗探头产生的冲击波在多层组织模型中的传播信息，t为碰撞开始后的时间。当t=2 μs时，冲击波仍在治疗探头内传播。当t=10 μs时，随着子弹与治疗探头的分离，驻波开始传入皮肤层。当t=20 μs时，首个包含正压(红色)和负压(蓝色)的冲击波在软组织中产生，波长为18 mm，且波前已到达脂肪层。当t=40 μs时，冲击波呈明显的放射式传播，且影响范围逐渐扩大，波前到达肌肉层。

图3 放射式冲击波在软组织模型中的传播信息

应用有限元分析法可展示多层软组织模型中任意单元格内冲击波的波形数据。图4为冲击波在皮肤层与治疗探头接触处网格单元中的波形，由于皮肤的声阻抗要远小于不锈钢材料，所以冲击波的振幅值缩小了一个数量级，但皮肤层中的最大正压强(P+)仍达到6.22 MPa。皮肤中密布大量的感觉神经末梢，所以在使用放射式冲击波体外治疗仪进行治疗时，患者普遍会感到刺痛。

图4 放射式冲击波在皮肤层中的波形图

图5为放射式冲击波传播至肌肉层时网格单元中的波形图。从图5可以看出，由于皮肤、脂肪和肌肉组织的阻抗效应，在深度为50 mm时，P+急剧下降，仅为1.26 MPa。由图5中的数据还可计算出放射式冲击波在软组织中传播的平均速度为1 500 m/s。将放射式冲击波的压强随传播深度变化如图6所示，可以看出，冲击波的P+和P-(最大负压强)受软组织阻抗效应的影响随传播深度的增加呈现梯度减小的趋势。

图5 放射式冲击波在深度为50 mm的肌肉层中的波形图

图6 放射式冲击波的压强随传播深度变化的规律

能流密度(EFD)是衡量放射式冲击波治疗效果的重要量化指标之一，其计算公式为

(1)

式(1)中：ρ为传播介质的密度；CL为冲击波在介质中的传播速度(CL=1 500 m/s)；P为冲击波的压强。根据有限元仿真结果并结合式(1)计算得到的能流密度如图7所示。从图7可以看出，与压强相比，放射式冲击波的能流密度受生物软组织阻抗效应的影响更加显著，从皮肤表层传播至深度50 mm处时，EFD从最高时的0.47 mJ/mm2逐渐降低至0.024 mJ/mm2，且受多层软组织内材料密度变化的影响，曲线的斜率也呈现逐渐降低的趋势。

图7 放射式冲击波的能流密度随传播深度变化的规律

2.2 耦合方式对冲击波传播特性的影响

在临床实践中，需要先在患处皮肤表面涂抹超声波耦合剂，使治疗探头与皮肤完全贴合后再启动仪器。然而，一些患者在家庭中使用治疗仪时往往省略这一步骤。因此，有必要研究治疗探头与皮肤间的耦合方式对治疗效果的具体影响。

图8为治疗探头在两种耦合方式下皮肤层产生的冲击波波形。如图8中红色曲线所示，子弹的高速碰撞会使治疗探头产生振动(周期为500 μs)，在缺少耦合剂黏合的情况下会与干燥的皮肤表面产生短暂分离，导致后续的冲击波无法连续传入组织。虽然最大正压强P+仍可达6.22 MPa，但由式(1)计算得到的能流密度值仅为0.016 mJ/mm2，几乎不会产生治疗效果。

图8 耦合方式对放射式冲击波波形的影响

结果表明，在治疗探头与皮肤表面接触不良的情况下，虽然患者会感受到刺痛，但放射式冲击波产生的能量无法传入组织内部，导致治疗失败。因此，在皮肤表面涂抹耦合剂，是放射式体外冲击波治疗仪发挥治疗效果的必要条件。

2.3 驱动气压对冲击波传播特性的影响

通过改变压缩空气的气压来改变子弹的初速度，进而对冲击波进行调制，可以更好地针对不同部位的软组织损伤进行治疗。然而，目前仍没有有效的方法来评估这种调制方式的有效性和准确性。因此，本文利用有限元仿真法针对这一问题进行研究。参照流体动力学的试验结果[18]，在ABAQUS软件的LOAD模块中将子弹的初速度设置为9、12.3、14.9、17.81 m/s,对应的驱动气压分别为0.1、0.2、0.3、0.4 MPa。

图9为在同一时刻(50 μs)4种驱动气压下的rESWT产生的放射式冲击波的仿真结果。从中可以发现，虽然rESWT产生的冲击波在三维形态上未发生显著变化，但随着驱动气压的增加，其影响范围逐渐扩大，同一位置处的强度也随之增大。图10为4种驱动气压下放射式冲击波在皮肤层中的波形，虽然波形在整体形态上具有相似性，但其压强随驱动气压的增加而明显增强，P+从驱动气压为0.1 MPa时的3.02 MPa提高至驱动气压为0.4 MPa时的8.74 MPa，即驱动气压越高，患者感受到的刺痛感越强。

图9 驱动气压为0.1～0.4 MPa的治疗仪输出的放射式冲击波

图10 不同驱动气压作用下冲击波的波形图

图11为4种驱动气压下下放射式冲击波的能流密度(EFD)随传播深度变化的规律。可见，驱动气压的提高可显著增加放射式冲击波的传播深度。Rompe等[19]将体外冲击波治疗仪产生的能流密度划分为3级：高能量区大于0.60 mJ/mm2、中等能量区在0.28～0.60 mJ/mm2范围、低能量区在0.08～

图11 不同驱动气压作用下冲击波的能流密度随传播深度变化的规律

0.28 mJ/mm2范围。如能流密度小于0.08 mJ/mm2则几乎不会产生治疗效果。根据这一标准，当驱动气压较低时(0.1 MPa)，即使是皮肤层产生的冲击波仍属于低能量区。而将驱动气压提高至0.4 MPa时，软组织模型中11～25 mm深的部位仍为中等能量区，低能量区更是可以深入41 mm处。

仿真结果在rESWT的临床实践中具有一定的指导意义。对于跟腱炎和网球肘等浅层的软组织损伤，以及病患为儿童或老人的情况，可以采用低驱动气压(0.1～0.2 MPa)进行治疗，在达到治疗效果的同时还能尽量减少病患的痛苦。而对于肩膀和大腿肌肉深层的软组织损伤，以及病患为肌肉发达的运动员或肥胖人群的情况，应采用高驱动气压(0.3～0.4 MPa),使放射式冲击波产生的能量能够传达患处，以保证治疗效果。

2.4 治疗探头几何结构对冲击波传播特性的影响

除调节驱动气压外，改变rESWT治疗探头的几何结构也可以对放射式冲击波的传播特性起到调制作用，从而提高针对不同部位软组织损伤的治疗效果。然而，由于缺少对冲击波在人体内传播数据的直接测量方法，新型治疗探头的研发仍主要依靠患者的主观反馈和医生的临床经验，其效率和有效性都难以保证。

所提出的有限元分析法可以为解决这一难题提供新思路，即通过对具有不同几何结构的治疗探头产生的放射式冲击波进行模拟仿真，研究其传播特性。作为对比，将R15即(即头部直径为15 mm)型和D35型(即头部直径为35 mm)治疗探头的几何信息导入ABAQUS软件中进行有限元仿真，当驱动气压为0.3 MPa时，其结果如图12所示。

从图12(a)可以看出，与头部直径为20 mm的D20型治疗探头[图12(b)]相比，在相同条件下，R15型治疗探头产生的放射式冲击波形态更加规整，能量更为集中，但影响范围有所减少。由式(1)的计算结果可知，R15型产生的能流密度(EFD)最大值为0.38 mJ/mm2，有效传播深度为28 mm，均小于D20型(0.48 mJ/mm2，38 mm)。因此，R15型探头更适合治疗肌腱炎和跟骨骨刺等深度较浅且范围较小的软组织损伤。

从图12(b)可以看出，与D20型和R15型相比，在相同条件下，D35型治疗探头产生的放射式冲击波的影响范围更广，且由于结构中存在较大的几何突变，冲击波在探头内部反射、叠加，造成其在组织中的传播形态更为复杂。由式(1)的计算结果可知，D35型产生的能流密度(EFD)最大值为0.46 mJ/mm2，有效传播深度为37 mm，与D20型十分接近。因此，D35型探头更适用于面积较大的深层肌肉痛点。