爆炸沖擊波對肺損傷的數(shù)值模擬*

周 杰，陶 鋼，王 健

(南京理工大學(xué)能源與動力工程學(xué)院，江蘇 南京210094)

爆炸沖擊波對人體的損傷主要集中在含有空氣的生物器官:耳膜、肺及腸胃［1］。G.J.Argyros［2］指出，沖擊波對肺的損傷是使人致命的關(guān)鍵因素。沖擊波與人體胸部作用時，誘發(fā)生物組織產(chǎn)生應(yīng)力波，應(yīng)力波在肺組織中傳播時，會導(dǎo)致肺泡破裂、撕裂、內(nèi)爆裂［3-4］。肺在受到應(yīng)力作用后的損傷位置、損傷程度和損傷機理一直受到人們的關(guān)注，有限元方法已成為研究損傷問題最主要的方法。

A.I.D＇yachenko 等［1］建立了一維有限元肺模型，并選用粘彈性材料模型研究了沖擊波傳播時造成的損傷。A.D.Geer［3］建立了人體胸部二維有限元模型，根據(jù)軟組織和硬組織的特性，選擇合理的材料模型，研究了肺在沖擊波作用下的損傷過程。

本文中，基于CT 圖片，建立人體胸部三維有限元模型。模型簡化成由肌肉、骨骼及肺3 部分組成。采用LS-DYNA 有限元程序中流固耦合方法，對人體胸部在自由空間爆炸場中受沖擊波作用的力學(xué)過程進行數(shù)值模擬。利用計算獲得的沖擊波入射超壓峰值和正壓持續(xù)時間，參照Bowen 曲線［5］，驗證模型的損傷狀態(tài)。通過觀察肺部應(yīng)力變化過程，分析肺部前表面應(yīng)力的變化規(guī)律，得到損傷最嚴重的區(qū)域，從而了解損傷的力學(xué)機制。

1 幾何模型的建立

本文中CT 圖像源自網(wǎng)絡(luò)。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 軟件，根據(jù)CT 圖像灰度值的不同，分別提取人體胸部的肌肉、骨骼和肺，獲得ASII 點云數(shù)據(jù)。采用CATA 軟件對點云數(shù)據(jù)分別進行去噪、濾波和光滑等處理，去除了不存在的奇點和多余的拐點，使點云的曲率變化變小和二階幾何連續(xù)。通過Geomagic 軟件，將各區(qū)域的點云數(shù)據(jù)按照給定的精度擬合成曲面，再根據(jù)胸部肌肉、骨骼和肺各部分的曲面曲率變化情況劃分成多個區(qū)域，最后將各區(qū)域曲面擬合生成CAD 幾何模型。由于對原始數(shù)據(jù)進行了去噪、濾波和光滑處理，因此生成的人體胸部幾何模型中肌肉、骨骼和肺的幾何關(guān)系也發(fā)生了一定的變化，三者之間會出現(xiàn)間隙或重疊等不匹配的情況，再利用Solidworks 軟件進行幾何模型修復(fù)，得到封閉的CAD 幾何模型。

圖1 CT 圖像Fig.1 CT image

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

2 有限元模型的建立

基于人體CT 圖像建立人體三維有限元模型。為了提高工作的可行性和計算效率，建立的有限元模型只保留了人體胸部的基本特征。因為心臟和肌肉材料的靜態(tài)力學(xué)特征大致相同，因此用肌肉代替了心臟。采用映射方法，對人體胸部模型進行六面體結(jié)構(gòu)化網(wǎng)格劃分，由于骨骼的外形十分復(fù)雜，在CAE 軟件處理過程中困難較大，因此進行了大量簡化。人體的胸部網(wǎng)格模型是由肌肉、骨骼和肺3 個模型組成的，如圖3 所示。由于肌肉體積較大，選用的最大單元尺寸為5 mm，而骨骼體積較小且形狀不規(guī)則。這里采用的最大單元尺寸為2 mm，肺的結(jié)構(gòu)比較規(guī)則，選用的最大單元尺寸為5 mm，整個模型網(wǎng)格數(shù)為43 477。胸部與頸部的接觸平面、胸部與腹部的接觸平面設(shè)定為透射邊界，即應(yīng)力波可透過接觸平面而不發(fā)生反射。

準確模擬爆炸沖擊波對人體胸部的作用過程，關(guān)鍵在于選取合理的人體胸部生物組織的材料模型和參數(shù)。對于模型中的胸部軟組織部分:肌肉組織采用粘彈性模型［6］，肺組織采用連續(xù)介質(zhì)的Fung-Vawter 模型［7］。硬組織骨骼則采用線彈性模型［7］。

3 計算模型和方法

模擬自由空間爆炸場中沖擊波對人體胸部損傷效應(yīng)的計算模型如圖4 所示:人體胸部正面朝向爆炸源，為了獲得一定的沖擊波壓力，模型距離爆炸源0.7 m;邊界的幾何中心處加入一小塊TNT 炸藥，運用點起爆方式獲得球形爆炸沖擊波，炸藥的網(wǎng)格數(shù)為344，藥量為70 g;設(shè)定長方體空氣域為1.2 m×1.4 m×0.45 m，網(wǎng)格數(shù)量為515 656，為模擬無限大空間效應(yīng)，6 個邊界均設(shè)為透射邊界，這些邊界允許流體介質(zhì)流出空氣域。

圖4 胸部損傷效應(yīng)的計算模型Fig.4 Calculation model of chest damage effect

采用流固耦合計算方法，對計算區(qū)域進行單元離散后，利用LS-DYNA 顯式積分法求解。計算域內(nèi)所有網(wǎng)格采用八結(jié)點映射網(wǎng)格單元，空氣和炸藥為多物質(zhì)流體ALE 網(wǎng)格，人體胸部生物組織為固體Lagrangian 網(wǎng)格，流固耦合采用罰函數(shù)約束方式追蹤結(jié)構(gòu)和流體位置間的相對位移，計算界面力并分布到流體結(jié)點上實現(xiàn)耦合。

圖5 Bowen 損傷曲線Fig.5 Bowen damage curves

4 結(jié)果與分析

4.1 肺的損傷評估

Bowen 損傷曲線是基于爆炸沖擊波對動物的損傷實驗數(shù)據(jù)、根據(jù)相應(yīng)的準則把動物的損傷狀態(tài)折合成人在自由空間爆炸場中的損傷狀態(tài)［5］，如圖5(人體胸部正面朝向爆炸源)所示。通過沖擊波入射的超壓峰值和正壓持續(xù)時間，參照Bowen 損傷曲線就可以評估肺臨界損傷(Dc)到99%致命損傷(D99)范圍內(nèi)的損傷狀態(tài)，D1、D99表示在受爆炸沖擊波作用后的24 h 內(nèi)的死亡率分別為1%和99%。

根據(jù)計算獲得的沖擊波入射超壓峰值pmax=1 160 kPa、正壓持續(xù)時間t=0.3 ms，如圖6 所示。參照圖5 的Bowen 損傷曲線獲得Dc、D1、D10和D50損傷狀態(tài)下的超壓值。正壓時間0.3 ms 時，肺的臨界損傷狀態(tài)的超壓值為900 kPa。1%致命的損傷狀態(tài)超壓值為2 500 kPa，10%致命的損傷狀態(tài)超壓值為2 900 kPa，50%致命的損傷狀態(tài)的超壓值為3 600 kPa。由此可以得出本文計算模型中肺處于Dc與D1損傷狀態(tài)之間。

4.2 力學(xué)響應(yīng)與損傷效應(yīng)分析

選擇緊靠胸廓前的空氣單元中的沖擊波的波形，由于沖擊波是從低阻抗材料空氣向高阻抗的胸廓傳播的，沖擊波在胸廓表面會產(chǎn)生反射，所以波形會呈現(xiàn)曲線上升的趨勢。參考理想沖擊波的波形［3］，可以獲得入射沖擊波的最大超壓峰值為1 160 kPa，正壓持續(xù)時間為0.3 ms。如圖6 所示，在0.46 ms 之前肺表面的壓力為環(huán)境壓力，應(yīng)力波還沒有作用到肺組織;當t=0.46 ms 時應(yīng)力波開始作用到肺表面，隨后壓力迅速增大;當t=1.0 ms時肺表面的超壓峰值達到最大值1 302 kPa，隨后壓力迅速下降到環(huán)境壓力以下;在隨后的時間內(nèi)，壓力振動衰減到環(huán)境壓力。由此可知，在球形沖擊波作用下，肺組織中的最大正應(yīng)力峰值大于入射沖擊波的超壓峰值。

圖6 沖擊波和肺組織中最大受壓單元壓力變化曲線Fig.6 Incident shock wave and stress wave of the lung tissue maximum pressure element

圖7 典型時刻的肺的應(yīng)力云圖Fig.7 Stress contours of the lung in certain time

圖7 是典型時刻肺的應(yīng)力云圖。t=0.45 ms時，應(yīng)力波還沒有作用到肺組織上，肺的應(yīng)力大小與環(huán)境壓力相同(圖7(a));當t=0.61 ms 時，應(yīng)力波最初作用到肺組織表面區(qū)域的正應(yīng)力峰值達到1 230 kPa(圖7(b));當t=1.0 ms 時，應(yīng)力波的正應(yīng)力峰值集中分布在肺的前表面，此時正應(yīng)力峰值為1 302 kPa 達到最大，正應(yīng)力峰值集中的部分也可能是肺組織損傷最為嚴重的區(qū)域(圖7(c))。為了更好地觀察肺組織內(nèi)部應(yīng)力變化的情況，圖7(d)～(e)中對肺部做了切片處理。當t=2.59 ms 時，肺的正應(yīng)力峰值下降到575 kPa，是因為應(yīng)力波在人體胸部傳播受到人體生物組織的阻尼作用產(chǎn)生了衰減(圖7(d));當t=4.14 ms 時，肺的正應(yīng)力峰值為534 kPa，雖然肺組織對應(yīng)力波傳播有阻尼作用，但是由于胸后肋骨對應(yīng)力波有反射作用，所以在肺后半段與肋骨接觸的部分也會出現(xiàn)比較高的正應(yīng)力值(圖7(e));當t=5.0 ms 時，應(yīng)力峰值下降到184 kPa，且大部分區(qū)域的值都比較低，在－200 ～120 kPa 范圍內(nèi)，－200 kPa 的物理意義是:正應(yīng)力大小為200 kPa，方向是外法線方向(很可能使肺泡產(chǎn)生內(nèi)爆裂)(圖7(f))。由此可見，胸部在受球形沖擊波作用，當t=1.1 ms 時，肺組織表面達到最大正應(yīng)力值1 302 kPa，且肺組織中比較大的正應(yīng)力集中出現(xiàn)在正對爆炸源的前表面，這部分的肺泡很可能會壓垮、破裂，成為肺損傷較為嚴重的區(qū)域之一。

圖8 是最大切應(yīng)力單元的τxy、τyz、τzx3 個方向的切應(yīng)力變化曲線，其位置處于肋骨與肺組織接觸的邊界上。在t=1.9 ms 之前3 個方向的切應(yīng)力都為零。當t=2.33 ms 時，切應(yīng)力達到第1 個峰值，分別是τxy=－227 kPa、τyz=－78.4 kPa、τzx=206 kPa。當t=4.35 ms 時，切應(yīng)力達到第2 個峰值，分別是τxy=－64.8 kPa、τyz=90.8 kPa、τzx=321 kPa(上述切應(yīng)力中的負號表示切應(yīng)力的方向指向坐標軸的負方向)。由此可知，當t=4.35 ms 時，此單元最大切應(yīng)力為τzx=321 kPa。由于缺乏高應(yīng)變率下生物組織的實驗拉伸數(shù)據(jù)，特別是對微觀組織，所以只能參考Yamada 的靜態(tài)拉伸張應(yīng)力數(shù)據(jù)(見圖9)。由圖9可見，肌肉最大靜態(tài)拉伸破壞應(yīng)力為140 kPa，因此通過計算得到的胸內(nèi)組織已經(jīng)處于破壞狀態(tài)。實際上，由于肺中的毛細管和肺泡壁只有一個細胞厚度且強度較低，因此對強應(yīng)力波效應(yīng)幾乎無防衛(wèi)能力［1］，研究表明，應(yīng)力波是肺創(chuàng)傷的主要原因。肺中的復(fù)雜波是由胸腔和其他生物材料的反射和聚焦產(chǎn)生的，而且剪切波可造成較硬的細支氣管剪切肺組織。這些現(xiàn)象已得到實驗證實，如肋膜破口和在肺泡及小葉內(nèi)靜脈壁之間的撕裂等［3］。由此可見，通過數(shù)值分析，可得到肺泡及胸內(nèi)組織最可能撕裂的位置和程度。此外，通過與Bowen 損傷曲線比較，可知肺損傷處于Dc與D1損傷狀態(tài)之間，而應(yīng)力分析的結(jié)果也表明存在損傷。

圖8 切應(yīng)力變化曲線Fig.8 Shear stress curves

圖9 肌肉組織靜態(tài)拉伸張應(yīng)力應(yīng)變曲線Fig.9 Static tensile stress-strain curve of muscle

總之，通過模擬人體胸部在自由爆炸空間中受沖擊波作用的力學(xué)響應(yīng)過程，并采用不同性質(zhì)的沖擊波，可以獲得對人體胸部的作用過程(包括最大正應(yīng)力位置和最大剪應(yīng)力位置)，從而判斷損傷的位置，并可根據(jù)解剖學(xué)結(jié)論得到驗證。我們已經(jīng)知道，爆炸肺損傷主要是由于肺泡的破裂引起的，而肺泡尺寸的量級是10－4m，因此微觀軟組織的破壞機理是一個需要進一步深入研究的問題。

5 結(jié) 論

通過數(shù)值模擬自由空間爆炸場中沖擊波與三維人體胸部作用的力學(xué)響應(yīng)過程，得出以下結(jié)論:

(1)基于CT 圖像，利用Mimics 軟件可以重建人體胸部三維模型。三維模型經(jīng)過一定的簡化，并給各個生物組織賦予合理的材料模型和參數(shù)，可以利用LS-DYNA 有限元程序進行計算;

(2)在球形沖擊波作用下，計算結(jié)果中比較大的正應(yīng)力集中出現(xiàn)在肺部的前表面，最大應(yīng)力值為1 302 kPa，稍大于入射沖擊波的超壓峰值1 160 kPa，這部分區(qū)域是肺損傷較為嚴重的區(qū)域之一;

(3)肺部最大切應(yīng)力位于肋骨與肺組織接觸邊界上，肺部最大切應(yīng)力峰值321 kPa，大于肌肉組織最大靜態(tài)破壞應(yīng)力140 kPa，可以導(dǎo)致組織的撕裂，引起肺損傷，這部分區(qū)域也是肺損傷較為嚴重的區(qū)域之一;

(4)參照Bowen 曲線的損傷標準，模型的損傷狀態(tài)處于Dc與D1之間，而計算結(jié)果表明切應(yīng)力引起肺部破裂，引起肺損傷，這與對應(yīng)的Bowen 曲線損傷標準結(jié)論一致。

［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，F(xiàn)letcher 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.