爆炸冲击波致颅脑冲击伤数值模拟研究

张文超，王舒，梁增友，覃彬，卢海涛，陈新元，卢文杰

（1. 中北大学 机电工程学院，山西，太原 030051；2. 瞬态冲击技术重点实验室，北京 102202；3. 中国人民解放军93320 部队，黑龙江，齐齐哈尔 161000；4. 河南北方红阳机电有限公司，河南，南阳 473000）

爆炸冲击脑损伤（blast-induced traumatic brain injury，bTBI）是近年来国外战争中战斗人员遭受的一种较为严重又多发的创伤，它被认为是伊拉克和阿富汗战争中最常见的军事伤害，超过50%的创伤性脑损伤是由爆炸引起的[1-2]. 美国五角大楼和退伍军人事务部估计，在过去的20 年里，全世界有40.8 万名军人遭受了不同程度的创伤性脑损伤[3]. 爆炸冲击脑损伤可导致记忆缺失、注意力不集中、情绪紧张等问题，还会导致睡眠障碍、疲劳、头晕、易怒、头痛和癫痫发作，严重影响士兵作战能力.

在研究早期，普遍认为爆炸冲击脑损伤是由于冲击波作用头部使头部获得速度撞击其他物体而产生[4]，后来通过动物实验发现，即使固定头部，冲击波仍能造成脑损伤[5-6]. 针对颅脑冲击伤问题，很多学者展开了研究. GOELLER 等[7]通过实验及数值模拟得出当冲击波正面作用头部时，冲击波通过颅骨、脑脊液(CSF)和脑组织传播，会在挫伤处产生负压导致空化，颅骨变形是引起空化的重要因素. MOSS 等[8]通过三维流体动力学模拟表明，爆炸冲击波直接作用于头部会导致颅骨弯曲，在脑组织中产生的机械载荷与在伤害性撞击中产生的载荷相当. GRUJICIC 等[9]通过数值模拟，将两次爆炸冲击波作用头部与头部钝击进行对比，结果表明即使在爆炸冲击波作用下头部米塞斯应力没有超过头部钝击时的脑损伤阈值，爆炸致使的颅内压也可能导致轻症外伤性 脑 损 伤(mild traumatic brain injury,mTBI). SÄLJÖ等[10]通过动物实验，发现爆炸冲击波作用下猪大脑中的颅内压峰值与头部外空气中的爆炸压力峰值相当，结果表明头骨几乎不能保护大脑免受爆炸的伤害. 韩建保等[11]建立多相非线性脑组织形变仿真计算模型，通过转动惯性载荷下多相非线性脑组织的变形状况计算结果表明，脑组织各相交界处应变较大，最大剪应力出现在大脑额叶侧. 栗志杰等[12]针对颅骨局部弯曲变形这一致伤机制对爆炸冲击脑损进行研究，确立了颅骨局部弯曲变形与脑组织压力之间的内在联系.

基于以上研究现状可知国内外对爆炸冲击波作用头部的作用过程、致伤机理及颅脑在冲击波作用下的动态响应特性尚未有一个统一的认识，对冲击波与头部相互作用过程中的反射叠加增强效应、颅脑易发生创伤区域等研究成果较少，这些研究内容对bTBI 的防护及治疗有很大意义，需要进行深入研究. 本文在此背景下做了爆炸冲击波作用头部仿真研究，建立了具有典型颅脑结构的三维头部有限元模型，对所建模型进行了撞击实验验证；利用三维多物质ALE 流固耦合方法模拟冲击波与头部相互作用过程，观察冲击波正面冲击头部时空气压力场变化规律；同时对脑组织压力进行分析，找出脑组织压力变化特性及易发生创伤区域，可为以后进一步的研究提供指导意义.

1 模型介绍及验证

1.1 有限元模型及接触关系

所建头部有限元模型皮肤、颅骨几何数据来源于中国人民解放军第三军医大学可视化人体切片数据集（CVH）. 基于皮肤、颅骨几何模型，利用Geomagic 对几何模型进行优化，利用hypermesh 对几何进行头部有限元网格划分. 建立的头部有限元模型包括皮肤、颅骨、面颅、眼、脑脊液、脑组织等组织器官及硬脑膜、软脑膜等膜结构；膜结构采用四边形壳单元划分，其余部分均采用6 面体实体单元，头部模型共有168 229 个节点，152 069 个单元，模型单元尺寸4～5 mm，头部模型总质量为4.6 kg；皮肤颅骨及脑组织间采用AUTOMATIC_SURFACE_TO_SURFACE_TIEBREAK 接触，膜结构与其他结构间采用AUTOMATIC_SURFACE_TO_SURFACE 接 触. 有 限元模型如图1.

图1 头部有限元模型Fig. 1 Finite element model of head

爆炸空气环境用60 cm×60 cm 立方体区域表征，头部放置空气域中心，冲击波与头部模型发生作用部分进行网格加密，以保证耦合界面不出现压力泄露；其他部分采用渐变网格，网格逐渐变大从而减少计算成本. 空气域中心区域单元尺寸4 mm，此时相应的计算结果收敛. 头部与空气域相对位置、冲击波作用方向、空气域网格示意图如图2.

图2 头部与空气域相对位置Fig. 2 Relative position of head and air domain

冲击波作用头部过程采用三维多物质ALE 流固耦合方法[13]来实现，流固耦合算法用于定义空气和头部之间的相互作用，此算法使空气给头部施加压力载荷，而头部则相当于空气的边界条件，用于约束空气的运动. 有限元模型中空气单元采用Euler 网格描述，头部单元采用Lagrangian 网格描述，利用*CONSTRAINED_LAGRANGE_IN_SOLID 关键字进行罚耦合算法的定义，定义空气网格边界参考压力为标准大气压101.325 kPa，这样冲击波负压段压力才会回升，与实际情况相符，如图3.

图3 空气域网格示意图Fig. 3 Schematic diagram of air domain grid

1.2 本构模型与材料参数

头部模型中脑组织材料接近不可压缩，是一种质地很软的生物软组织，其表现为不可压缩性、非线性、各向异性及黏弹性. 大量脑组织实验证明脑组织的变形只取决于其剪切模量，本文脑组织采用黏弹性本构模型[14](*MAT_VISCOELASTIC)对材料性能进行表征，其剪切 弹性模量的计算公式为

式中：G0为 短效剪切模量；G∞为 长效剪切模量； β为衰减系数.

脑脊液采用弹性流体本构模型[15]（*MAT_ELASTIC_FLUID）表征.

头部模型皮肤、颅骨、面颅、硬脑膜、软脑膜均采用线弹性模型进行表征，各部分材料参数如表1～表3 所示.

表1 弹性模型材料参数Tab. 1 Material parameters of elastic model

表2 黏弹性模型材料参数Tab. 2 Material parameters of viscoelastic model

表3 弹性流体模型材料参数Tab. 3 Material parameters of elastofluid model

1.3 头部有限元模型验证

基于头部三维有限元模型进行冲击波作用头部仿真研究，首先需要对头部模型有效性进行验证.NAHUM 等[16]做了正面冲击作用下颅内压变化的人体死尸头部碰撞实验，它已成为头部数值模型有效性验证的经典实验. 该实验中一个以恒定速度移动的刚性物体撞击前额部位，利用传感器记录头部特定部位颅内压变化过程；由于实验没有提供撞击物体的材料力学性能，参考RUAN 等[17]验证头部模型所使用的方法：将实验中测量所得到的碰撞接触力载荷转化为平均压力作用于前额部位，作用面积为1 480 mm2，通过数值仿真与实验结果进行对比分析，从而验证头部模型的有效性.

Nahum 实验中测量了额骨、顶骨、枕骨及小脑窝处颅内压变化过程，数值模拟中分别在相应位置设置3 个观测点，观察观测点颅内压变化情况，与实验数据进行对比. 观测点位置示意图如图4 所示，头部加载区域如图5 阴影所示.

图4 观测点位置示意图Fig. 4 Schematic diagram of observation point position

图5 撞击区域示意图Fig. 5 Schematic diagram of impact area

Nahum 实验中撞击物撞击头部输出撞击力时程曲线如图6，转化为平均压力曲线如图7. 数值模拟结果与实验数据对比如图8.

图6 实验撞击力时程曲线Fig. 6 Time history curve of experimental impact force

图7 仿真加载压力时程曲线Fig. 7 Time history curve of simulated loading pressure

由图8 可以看出数值模拟中观测点颅内压变化趋势、颅内压峰值与实验结果吻合较好. 由不同部位观测点颅内压-时间历程曲线知，头部前额受到撞击时，皮肤、颅骨受力获得较大加速度向后运动，同时颅骨将压力传递给脑组织，此时前额脑脊液、脑组织受压，颅内压升高；由于惯性力差异脑组织与颅骨发生相对位移，与前额相对的枕部、脑后窝处脑组织与颅骨远离，颅内压降为负值. 随着压力载荷的增加前额处脑组织继续受压，同时压力增大，不同特征部位颅内压值也相应增大. 随着压力载荷减小，颅骨、脑组织相对位移变小，颅内压逐渐减小. 数值仿真结果与撞击实验数据吻合性较好，验证结果表明头部模型可以正确反映头部受到物理冲击时动态响应特性，能够用于进一步的研究.

图8 仿真与实验结果对比Fig. 8 Comparison of simulation and experimental results

2 爆炸冲击波载荷和边界条件

模拟了60 g TNT 距离头部1 m 处爆炸时冲击波与头部相互作用过程. 空气参考压力为标准大气压101.325 kPa，在欧拉域迎爆面设置压力加载边界，冲击波载荷从此边界输入. 冲击波载荷输入曲线如图9，为了避免压力波在欧拉域边界反射对头部产生不真实的二次冲击，在欧拉域四周及出口设置非反射边界[18].

图9 迎爆边界冲击波输入载荷Fig. 9 Shock wave input load curve at detonation boundary

空气冲击波超压在正压段大致按指数规律衰减，一些经验公式可以描述此衰减过程，其中FRIEDLANDER[19]方程较接近实际且又简单易于计算

冲击波即将与头部接触时,取冲击波波阵面中心单元观测冲击波超压时间历程曲线，其与典型冲击波曲线衰减规律相符，此时冲击波压力为227 kPa，正压持续时间为730 us. 观测单元冲击波超压时间历程曲线如图10.

图10 作用头部冲击波载荷Fig. 10 Acting shock wave loads on the head

3 爆炸冲击波作用头部数值模拟与分析

基于上述头部有限元模型及爆炸冲击波加载，利用LS-DYNA 显示动力分析软件对爆炸冲击波正面作用头部进行数值模拟，分析爆炸冲击波与头部相互作用规律，找出脑部易受损区域.

3.1 空气压力场演化规律及冲击波反射叠加增强效应

爆炸冲击波作用头部时空气压力场演化过程如图11 所示.t=350 μs 时，冲击波到达头部前方，此时头部前方冲击波压力峰值为227 kPa；t=390 μs 时，冲击波与面部相互作用，由于面部皮肤波阻抗较大，冲击波与面部皮肤相互作用时很大部分冲击波发生反射，与后部压力波叠加而使超压增强，此时与面部接触冲击波峰值压力为326 kPa；t=470 μs 时，冲击波作用于前额，同样因为冲击波反射叠加增强使得峰值压力增大为387 kPa；t=700 μs 时，颅顶曲面结构使冲击波出现边界层分离现象，此时颅顶冲击波压力为194 kPa，冲击波在颅顶部位反射叠效增强效应较弱；t=1 110 μs 时，冲击波绕过头部在枕部发生二次汇聚，对头颅枕部发生冲击，汇聚点压力峰值最大，较上下两侧未汇聚点压力峰值增大约20～40 kPa，冲击波汇聚点峰值压力为183 kPa.

图11 爆炸冲击波压力云图Fig. 11 Pressure nephogram of explosion shock wave

取1～12 不同位置特征点，观测空气压力场时间历程曲线，更好地说明冲击波反射叠加增强现象，观测点位置如图12.

图12 空气压力场观测点Fig. 12 Observation point of air pressure field

在前额部位空气域中选3 点观测压力变化，2 号位置和3 号位置距头部较远，由图13 压力曲线可看到2、3 观测点处空气压力反射叠加使压力曲线产生二次峰值；1 点距头部较近，无法看到冲击波压力衰减段，压力直接叠加增强为322 kPa，较无头部阻挡时1 号点相同位置227 kPa 增强1.4 倍.

图13 前额空气域压力时程曲线Fig. 13 Pressure time-history curve of forehead air field

在颅顶上方空气域中取4、5、6、7 这4 个观测点观察压力变化，由图14 压力曲线可知观测点峰值压力无明显增大现象，随着冲击波向后传播，压力峰值逐渐降低，7 号点位置距头部较远，其峰值压力与5号位置基本相同，冲击波反射叠加增强效应较弱.在脑后由上至下取5 个观测点观察压力变化. 由图15 压力曲线可知，冲击波从颅顶、颅底传播至脑后，10 号观测点压力峰值最大，两侧观测点压力峰值相对较小，观测点压力曲线与空气压力云图相呼应，冲击波在脑后汇聚点压力得到一次增强. 冲击波在脑后位置与头部反射叠加增强效应较弱，冲击波在脑后汇聚后继续向后传播.

图14 颅顶空气域压力时程曲线Fig. 14 Pressure time-history curve in the air region of skull top

图15 脑后空气域压力时程曲线Fig. 15 Pressure time history curve of posterior air domain

基于以上可知爆炸冲击波作用头部时，头部迎爆面冲击波反射叠加增强效应较强，冲击波与头部皮肤相互作用可使峰值压力增强1.4 倍，颅顶及头颅后部冲击波反射叠加增强现象不明显. 冲击波绕过头部后会在头部后方一处汇聚，汇聚点冲击波峰值压力得到增强，会对脑后造成二次冲击.

3.2 脑组织压力分析

在过去研究中，人们提出了脑干剪切力、大脑中的压差和撞击对冲部位出现的负压及气穴等各种假说，从不同角度描述了大脑的冲击响应，但他们均是对大脑反应的部分描述，缺乏冲击作用下大脑的总体描述. 基于此，Nahum 等利用脑部有限元模型对一些动物和人体尸体脑部撞击试验以及飞机事故中的头部撞击行为进行数值模拟，得出颅内压峰值压力通常与脑创伤严重程度成正比；这些压力通常会造成脑挫伤或脑表面的小血管出血；通过对颅内压和挫伤的相关综合性分析和实验研究，提出了以脑挫伤和脑出血为基础的脑损伤标准，并将该标准与HIC、GSI 等损伤标准进行对比，认为该标准更合适低瞬态冲击下的颅脑损伤评估. 该颅内压耐受标准指出，当颅内压峰值超过235 kPa 时，会发生严重脑损伤，低于173 kPa 时，会发生轻微或无损伤[20-22]. 本文利用这一标准进一步分析冲击波与头部相互作用规律，找出头部易受损区域. 在前额、颅顶及枕部的脑组织及相邻的空气域中取观测点，观测组织压力变化规律. 观测点位置如图16.取空气中I 号观测点，脑组织中1、2、3 号观测点，观察在冲击波作用下相邻位置脑组织压力变化规律. 前额处空气压力及脑组织压力变化情况如图17，观测点I 处冲击波在490 μs 时作用头部，压力峰值为276 kPa，冲击波作用后相邻位置脑组织压力迅速升高且出现高频波动，波动频率为5.88 kHz；2号观测点压力峰值达229 kPa，压力出现两个较高峰值后，压力波动频率降低，且压力峰值也随时间逐渐降低，1、2 观测点压力峰值均达到颅内压脑损伤阈值175 kPa，可能会出现脑创伤.

图16 压力观测点Fig. 16 Pressure observation point line

图17 前额观测点压力变化Fig. 17 Pressure change of forehead observation point

取空气中II 号观测点，脑组织4、5、6 号观测点，颅顶处空气压力及脑组织压力变化情况如图18. 观测点II 处冲击波在660 μs 时作用头部，压力峰值为201 kPa，此时4、5、6 号观测点处脑组织压力变化规律与前额观测点基本一致，脑组织压力均有一次较大峰值，随着时间推移压力峰值逐渐降低，4 号观测点压力峰值达177 kPa，组织压力初始波动频率为4.76 kHz，随着时间推移压力峰值和压力波动频率逐渐降低.

图18 颅顶观测点压力变化Fig. 18 Pressure changes at cranial observation points

取空气中III 号观测点，脑组织7、8、9 号观测点，脑后空气压力及相邻脑组织压力变化情况如图19.

图19 脑后观测点压力变化Fig. 19 Pressure change of observation point behind brain

III 号观测点冲击波在1 080 μs 时作用头部，压力峰值为186 kPa. 冲击波从正面作用头部，前额先受到冲击波冲击，靠近前额部分脑组织先受压，由于惯性力差异颅骨较脑组织位移较大，脑后组织受拉，所以在冲击波未直接冲击脑后时枕部脑组织压力会先出现负值. Ⅲ号观测点位置冲击波在1 080 μs 绕到脑后对头部进行冲击，但是Ⅲ号观测点附近脑组织在脑后冲击波直接作用前就先出现一次较高压力峰值，本文进一步结合脑组织压力云图（图20）对此现象进行分析.

图20 脑组织压力云图Fig. 20 Pressure cloud of brain tissue

由脑组织压力云图可知，头部受到冲击波作用后，首先靠近前额部位组织受压，同时压缩波在脑组织中向脑后传播，脑组织受压部位不断向脑后扩散，压缩波到达脑后在脑组织-颅骨界面由于波阻抗不匹配发生反射叠加增强，从而在脑组织中出现较大的压力峰值. 冲击波在空气中较压缩波在脑组织中传播较慢，枕部组织中出现一个较大压力峰值后会出现一个第二峰，此由汇聚到脑后的冲击波对脑后进行的二次冲击造成. 第二压力峰值较第一压力峰值较小，说明枕部脑组织所受压力大小受冲击波二次冲击影响较小，脑组织压力需结合压力波在脑组织中传播规律及颅脑受外界压力等因素进行综合分析.

综上可知，头部正面受到冲击波作用时，前额部位脑组织压力峰值较高，颅顶、脑后组织压力峰值逐渐减小，压缩波在脑组织中传播在脑后组织与颅骨界面处可能出现较高压力峰值，在医学救治中需对此位置创伤情况加以重视. 冲击波在脑后汇聚后虽然会对头部再次进行冲击，但是脑组织受二次冲击后未出现较高的压力峰值. 冲击波正面作用头部时，前额、颅顶脑组织压力急剧增大，且出现较高频率的正负压周期性波动，波动频率可高达5.88 kHz，这会使脑组织在局部出现较大压力梯度，增大脑挫伤风险. 随着时间推移，压力峰值逐渐降低. 额叶、顶叶脑组织压力多处超过173 kPa 脑损伤阈值，为易受损区域.

4 结 论

建立了具有典型人体头部结构的三维有限元模型，进行头部撞击数值模拟并与Nahum 头部撞击实验数据进行对比，验证了头部模型的有效性；进一步进行了爆炸冲击波作用头部的数值模拟研究，得出以下结论：

①所建头部有限元模型能正确反映头部在外界冲击下的组织压力动态响应特性，前额受到撞击后，脑组织额叶处于压缩状态颅内压升高，枕叶处于拉伸状态颅内压出现负压，随着压力载荷减小，颅骨、脑组织相对位移变小，颅内压逐渐减小.

②爆炸冲击波正面作用头部后在面部会出现较强的反射叠加增强现象，观测点冲击波压力峰值较自由场爆炸时增大1.4 倍，顶部及枕部冲击波反射叠加增强效应较弱，没有出现较高的二次压力峰值.

③冲击波通过皮肤颅骨等多层介质传入脑组织，由于不同组织介质波阻抗不匹配特性，压力波在不同组织界面会出现反射和透射而出现局部的高压或压力衰减，枕部脑组织与颅骨界面因为压力波反射叠加增强会出现较高的压力峰值且超过脑损伤阈值，在冲击波脑损伤的医学救治方面对此部位需重点观察. 冲击波绕过头部后在头部后方一处进行压力汇聚，汇聚部位会受到冲击波二次冲击，但是二次冲击对组织压力影响较小.

④头部正面受到冲击波作用后，额叶、顶叶脑组织压力会出现高频的正负压周期性波动，且会出现较高的压力峰值，超过脑损伤阈值173 kPa，为易受损区域；为避免头部直接受到冲击，爆炸冲击波的防护装备研发很有必要.