爆炸冲击波对肺损伤的数值模拟*

周 杰，陶 钢，王 健

(南京理工大学能源与动力工程学院，江苏 南京210094)

爆炸冲击波对人体的损伤主要集中在含有空气的生物器官:耳膜、肺及肠胃［1］。G.J.Argyros［2］指出，冲击波对肺的损伤是使人致命的关键因素。冲击波与人体胸部作用时，诱发生物组织产生应力波，应力波在肺组织中传播时，会导致肺泡破裂、撕裂、内爆裂［3-4］。肺在受到应力作用后的损伤位置、损伤程度和损伤机理一直受到人们的关注，有限元方法已成为研究损伤问题最主要的方法。

A.I.D＇yachenko 等［1］建立了一维有限元肺模型，并选用粘弹性材料模型研究了冲击波传播时造成的损伤。A.D.Geer［3］建立了人体胸部二维有限元模型，根据软组织和硬组织的特性，选择合理的材料模型，研究了肺在冲击波作用下的损伤过程。

本文中，基于CT 图片，建立人体胸部三维有限元模型。模型简化成由肌肉、骨骼及肺3 部分组成。采用LS-DYNA 有限元程序中流固耦合方法，对人体胸部在自由空间爆炸场中受冲击波作用的力学过程进行数值模拟。利用计算获得的冲击波入射超压峰值和正压持续时间，参照Bowen 曲线［5］，验证模型的损伤状态。通过观察肺部应力变化过程，分析肺部前表面应力的变化规律，得到损伤最严重的区域，从而了解损伤的力学机制。

1 几何模型的建立

本文中CT 图像源自网络。CT 图像为一套健康男性(45 岁、身高170 cm、体重68 kg)的DICOM 格式扫描影像，扫描层厚0.6 mm，分辨率512×512 像素。选择其中胸部的686 张CT 图像作为处理对象。由于人体胸部肌肉、骨骼和肺等组织材料属性的差异，使得反映在CT 图像上的灰度值不同，如图1 所示，肺的灰度值在－1 024 ～－750 之间，肌肉的灰度值在－750 ～1 250 之间，而骨骼的灰度值在1 250 ～3 056 之间。图2 是建立人体胸部三维模型的流程图。首先利用Mimics 软件，根据CT 图像灰度值的不同，分别提取人体胸部的肌肉、骨骼和肺，获得ASII 点云数据。采用CATA 软件对点云数据分别进行去噪、滤波和光滑等处理，去除了不存在的奇点和多余的拐点，使点云的曲率变化变小和二阶几何连续。通过Geomagic 软件，将各区域的点云数据按照给定的精度拟合成曲面，再根据胸部肌肉、骨骼和肺各部分的曲面曲率变化情况划分成多个区域，最后将各区域曲面拟合生成CAD 几何模型。由于对原始数据进行了去噪、滤波和光滑处理，因此生成的人体胸部几何模型中肌肉、骨骼和肺的几何关系也发生了一定的变化，三者之间会出现间隙或重叠等不匹配的情况，再利用Solidworks 软件进行几何模型修复，得到封闭的CAD 几何模型。

图1 CT 图像Fig.1 CT image

图2 建模流程图Fig.2 Process of establishing model

2 有限元模型的建立

基于人体CT 图像建立人体三维有限元模型。为了提高工作的可行性和计算效率，建立的有限元模型只保留了人体胸部的基本特征。因为心脏和肌肉材料的静态力学特征大致相同，因此用肌肉代替了心脏。采用映射方法，对人体胸部模型进行六面体结构化网格划分，由于骨骼的外形十分复杂，在CAE 软件处理过程中困难较大，因此进行了大量简化。人体的胸部网格模型是由肌肉、骨骼和肺3 个模型组成的，如图3 所示。由于肌肉体积较大，选用的最大单元尺寸为5 mm，而骨骼体积较小且形状不规则。这里采用的最大单元尺寸为2 mm，肺的结构比较规则，选用的最大单元尺寸为5 mm，整个模型网格数为43 477。胸部与颈部的接触平面、胸部与腹部的接触平面设定为透射边界，即应力波可透过接触平面而不发生反射。

准确模拟爆炸冲击波对人体胸部的作用过程，关键在于选取合理的人体胸部生物组织的材料模型和参数。对于模型中的胸部软组织部分:肌肉组织采用粘弹性模型［6］，肺组织采用连续介质的Fung-Vawter 模型［7］。硬组织骨骼则采用线弹性模型［7］。

3 计算模型和方法

模拟自由空间爆炸场中冲击波对人体胸部损伤效应的计算模型如图4 所示:人体胸部正面朝向爆炸源，为了获得一定的冲击波压力，模型距离爆炸源0.7 m;边界的几何中心处加入一小块TNT 炸药，运用点起爆方式获得球形爆炸冲击波，炸药的网格数为344，药量为70 g;设定长方体空气域为1.2 m×1.4 m×0.45 m，网格数量为515 656，为模拟无限大空间效应，6 个边界均设为透射边界，这些边界允许流体介质流出空气域。

图4 胸部损伤效应的计算模型Fig.4 Calculation model of chest damage effect

采用流固耦合计算方法，对计算区域进行单元离散后，利用LS-DYNA 显式积分法求解。计算域内所有网格采用八结点映射网格单元，空气和炸药为多物质流体ALE 网格，人体胸部生物组织为固体Lagrangian 网格，流固耦合采用罚函数约束方式追踪结构和流体位置间的相对位移，计算界面力并分布到流体结点上实现耦合。

图5 Bowen 损伤曲线Fig.5 Bowen damage curves

4 结果与分析

4.1 肺的损伤评估

Bowen 损伤曲线是基于爆炸冲击波对动物的损伤实验数据、根据相应的准则把动物的损伤状态折合成人在自由空间爆炸场中的损伤状态［5］，如图5(人体胸部正面朝向爆炸源)所示。通过冲击波入射的超压峰值和正压持续时间，参照Bowen 损伤曲线就可以评估肺临界损伤(Dc)到99%致命损伤(D99)范围内的损伤状态，D1、D99表示在受爆炸冲击波作用后的24 h 内的死亡率分别为1%和99%。

根据计算获得的冲击波入射超压峰值pmax=1 160 kPa、正压持续时间t=0.3 ms，如图6 所示。参照图5 的Bowen 损伤曲线获得Dc、D1、D10和D50损伤状态下的超压值。正压时间0.3 ms 时，肺的临界损伤状态的超压值为900 kPa。1%致命的损伤状态超压值为2 500 kPa，10%致命的损伤状态超压值为2 900 kPa，50%致命的损伤状态的超压值为3 600 kPa。由此可以得出本文计算模型中肺处于Dc与D1损伤状态之间。

4.2 力学响应与损伤效应分析

选择紧靠胸廓前的空气单元中的冲击波的波形，由于冲击波是从低阻抗材料空气向高阻抗的胸廓传播的，冲击波在胸廓表面会产生反射，所以波形会呈现曲线上升的趋势。参考理想冲击波的波形［3］，可以获得入射冲击波的最大超压峰值为1 160 kPa，正压持续时间为0.3 ms。如图6 所示，在0.46 ms 之前肺表面的压力为环境压力，应力波还没有作用到肺组织;当t=0.46 ms 时应力波开始作用到肺表面，随后压力迅速增大;当t=1.0 ms时肺表面的超压峰值达到最大值1 302 kPa，随后压力迅速下降到环境压力以下;在随后的时间内，压力振动衰减到环境压力。由此可知，在球形冲击波作用下，肺组织中的最大正应力峰值大于入射冲击波的超压峰值。

图6 冲击波和肺组织中最大受压单元压力变化曲线Fig.6 Incident shock wave and stress wave of the lung tissue maximum pressure element

图7 典型时刻的肺的应力云图Fig.7 Stress contours of the lung in certain time

图7 是典型时刻肺的应力云图。t=0.45 ms时，应力波还没有作用到肺组织上，肺的应力大小与环境压力相同(图7(a));当t=0.61 ms 时，应力波最初作用到肺组织表面区域的正应力峰值达到1 230 kPa(图7(b));当t=1.0 ms 时，应力波的正应力峰值集中分布在肺的前表面，此时正应力峰值为1 302 kPa 达到最大，正应力峰值集中的部分也可能是肺组织损伤最为严重的区域(图7(c))。为了更好地观察肺组织内部应力变化的情况，图7(d)～(e)中对肺部做了切片处理。当t=2.59 ms 时，肺的正应力峰值下降到575 kPa，是因为应力波在人体胸部传播受到人体生物组织的阻尼作用产生了衰减(图7(d));当t=4.14 ms 时，肺的正应力峰值为534 kPa，虽然肺组织对应力波传播有阻尼作用，但是由于胸后肋骨对应力波有反射作用，所以在肺后半段与肋骨接触的部分也会出现比较高的正应力值(图7(e));当t=5.0 ms 时，应力峰值下降到184 kPa，且大部分区域的值都比较低，在－200 ～120 kPa 范围内，－200 kPa 的物理意义是:正应力大小为200 kPa，方向是外法线方向(很可能使肺泡产生内爆裂)(图7(f))。由此可见，胸部在受球形冲击波作用，当t=1.1 ms 时，肺组织表面达到最大正应力值1 302 kPa，且肺组织中比较大的正应力集中出现在正对爆炸源的前表面，这部分的肺泡很可能会压垮、破裂，成为肺损伤较为严重的区域之一。

图8 是最大切应力单元的τxy、τyz、τzx3 个方向的切应力变化曲线，其位置处于肋骨与肺组织接触的边界上。在t=1.9 ms 之前3 个方向的切应力都为零。当t=2.33 ms 时，切应力达到第1 个峰值，分别是τxy=－227 kPa、τyz=－78.4 kPa、τzx=206 kPa。当t=4.35 ms 时，切应力达到第2 个峰值，分别是τxy=－64.8 kPa、τyz=90.8 kPa、τzx=321 kPa(上述切应力中的负号表示切应力的方向指向坐标轴的负方向)。由此可知，当t=4.35 ms 时，此单元最大切应力为τzx=321 kPa。由于缺乏高应变率下生物组织的实验拉伸数据，特别是对微观组织，所以只能参考Yamada 的静态拉伸张应力数据(见图9)。由图9可见，肌肉最大静态拉伸破坏应力为140 kPa，因此通过计算得到的胸内组织已经处于破坏状态。实际上，由于肺中的毛细管和肺泡壁只有一个细胞厚度且强度较低，因此对强应力波效应几乎无防卫能力［1］，研究表明，应力波是肺创伤的主要原因。肺中的复杂波是由胸腔和其他生物材料的反射和聚焦产生的，而且剪切波可造成较硬的细支气管剪切肺组织。这些现象已得到实验证实，如肋膜破口和在肺泡及小叶内静脉壁之间的撕裂等［3］。由此可见，通过数值分析，可得到肺泡及胸内组织最可能撕裂的位置和程度。此外，通过与Bowen 损伤曲线比较，可知肺损伤处于Dc与D1损伤状态之间，而应力分析的结果也表明存在损伤。

图8 切应力变化曲线Fig.8 Shear stress curves

图9 肌肉组织静态拉伸张应力应变曲线Fig.9 Static tensile stress-strain curve of muscle

总之，通过模拟人体胸部在自由爆炸空间中受冲击波作用的力学响应过程，并采用不同性质的冲击波，可以获得对人体胸部的作用过程(包括最大正应力位置和最大剪应力位置)，从而判断损伤的位置，并可根据解剖学结论得到验证。我们已经知道，爆炸肺损伤主要是由于肺泡的破裂引起的，而肺泡尺寸的量级是10－4m，因此微观软组织的破坏机理是一个需要进一步深入研究的问题。

5 结 论

通过数值模拟自由空间爆炸场中冲击波与三维人体胸部作用的力学响应过程，得出以下结论:

(1)基于CT 图像，利用Mimics 软件可以重建人体胸部三维模型。三维模型经过一定的简化，并给各个生物组织赋予合理的材料模型和参数，可以利用LS-DYNA 有限元程序进行计算;

(2)在球形冲击波作用下，计算结果中比较大的正应力集中出现在肺部的前表面，最大应力值为1 302 kPa，稍大于入射冲击波的超压峰值1 160 kPa，这部分区域是肺损伤较为严重的区域之一;

(3)肺部最大切应力位于肋骨与肺组织接触边界上，肺部最大切应力峰值321 kPa，大于肌肉组织最大静态破坏应力140 kPa，可以导致组织的撕裂，引起肺损伤，这部分区域也是肺损伤较为严重的区域之一;

(4)参照Bowen 曲线的损伤标准，模型的损伤状态处于Dc与D1之间，而计算结果表明切应力引起肺部破裂，引起肺损伤，这与对应的Bowen 曲线损伤标准结论一致。

［1］ D’yachenkoa A I，Manyuhina O V.Modeling of weak blast wave propagation in the lung［J］.Journal of Biomechanics，2006，39(11):2113-2122.

［2］ Argyros G J.Management of primary blast injury［J］.Toxicology，1997，121(1):105-115.

［3］ Geer A D.Numerical modeling for the prediction of primary blast injury to the lung［D］.Ontario，Canada:University of Waterloo，2006.

［4］ Axe J R，Abbrecht P H.Analysis of the pressure-volume relationship of excised lungs［J］.Annals of Biomedical Engineering，1985，13(2):101-117.

［5］ Bowen I G，Fletcher E R，Richmond D R.Estimate of man’s tolerance to the direct effects of air blast［R］.DASA-2113，1968.

［6］ Mukherjee S，Chawla A，Karthikeyan B，et al.Finite element crash simulations of the human body:Passive and active muscle modeling［J］.Sadhana，2007，32(4):409-426.

［7］ Chang H T.The development，validation and comparison of a finite element human thorax model for automotive impact injury studies［D］.Lowa，America:The University of Lowa，2001.