爆炸冲击波作用下颅脑损伤机理的数值模拟研究*

2020-04-09栗志杰由小川柳占立杜智博

爆炸与冲击 2020年1期

栗志杰，由小川，柳占立，杜智博，张仡，杨策，庄茁

（1. 清华大学航天航空学院，北京 100084；2. 陆军军医大学大坪医院，重庆 400038）

现代局部战争战伤分析显示，由爆炸导致的单兵伤亡约占总伤亡人数的70%，其中冲击波直接作用引起的伤害约占60%[1]。同时，国际上爆炸恐怖袭击日益猖獗，严重危害民众的生命安全。而作为人体中枢神经系统的头部，是爆炸冲击波的重要靶器官。据相关报告显示，美军中受到创伤性脑损伤（traumatic brain injury, TBI）影响的军人为总服役人数的20%，这其中由于爆炸冲击波所引起的脑损伤(blast-induced traumatic brain injury, b-TBI)占比高达40%～60%[1]。爆炸冲击波引起颅脑损伤的途径主要包括冲击波对头部的直接作用、飞起破片对头部的冲击作用、头部与其他物体的撞击作用、冲击波对头部的电磁热等作用[2-3]。随着爆炸伤的日渐频发，爆炸冲击波作用下的颅脑致伤机理与对应的防护策略成为研究热点。关于颅脑致伤机理，当前的研究主要集中在冲击波对头部直接作用所造成的颅脑损伤。

虽然，经过大量的动物实验、病理生理学分析以及数值模拟研究，建立了颅脑碰撞损伤与头部加速度之间的联系，并确定了相应的损伤准则[4-6]，而相关理论并不能应用于爆炸冲击波直接作用下的颅脑损伤分析。Cernak 等[7]在生物激波管实验中，将老鼠进行固定，保证头部在冲击波条件下为稳定状态，但是实验结果表明，老鼠的脑组织仍然出现了器质性病变。费周等利用不同尺寸的微型球形炸药，对大鼠进行不同损伤程度的b-TBI 实验，并从神经功能、宏观病理学、组织病理学及各种生物标志物表达水平方面，建立了可靠且可重复的b-TBI 模型及损伤特征指标[8]。该实验将大鼠放置在泡沫垫块上，并将球形炸药悬空挂于大鼠头部上面，大鼠头部在爆炸冲击波作用下是相对稳定的，但是仍然出现了脑组织损伤。这些动物实验均表明，爆炸冲击波对颅脑的损伤可以通过直接作用产生。而爆炸冲击波对头部直接作用的时间通常在1～3 ms，是应力波主导的物理现象。因此，基于头部加速度建立的颅脑碰撞损伤理论无法解释上述实验现象。

针对爆炸冲击波直接作用下的颅脑损伤，学者们做了大量的研究工作，并提出了多种致伤机理。Cloots 等[9-11]将微观体积单元（critical volume element, CVE）与常规有限元单元进行结合，建立了脑组织的多尺度模型，从而得到脑组织细胞层次的力学响应，并结合生物细胞的力学损伤实验，建立了相应的损伤准则。Radovitzky 等[12]、Jean 等[13]发展了脑组织的多尺度生物力学模型，研究了爆炸冲击波作用下颅脑在组织、细胞等不同尺度的生物响应特征，并建立了评估b-TBI 风险的动物与人之间的尺度律。Goeller 等[14]、Salzar 等[15]提出脑脊液空穴致伤机制，通过激波管实验与数值模拟，对脑脊液的空穴效应进行了研究，研究结果表明，脑脊液在冲击波作用下会发生空化与湮灭，形成局部点高压，从而造成脑组织损伤。Franck[16]基于类脑软材料，研究了微型空穴效应在b-TBI 当中的作用，并利用实时成像技术对不同应变率下神经元细胞的损伤模式进行了探究。Bhattacharjee[17]提出主动脉冲压致伤理论，冲击波通过压缩主动脉血液，从而使其形成高压液体，并经毛细血管对脑组织形成冲击，造成脑损伤。同时，Courtney 等[18]设计了动物实验，验证了冲击波可以通过作用躯干部位间接对脑组织造成损伤。Moss 等[19]通过数值模拟研究发现，颅骨在爆炸冲击波作用下会发生局部弯曲，从而使脑组织压力出现较大的压力梯度，颅骨的弯曲变形与颅内压力的分布密切相关。但是，颅骨局部弯曲变形与脑组织压力之间的内在联系有待于进一步确立，相应的颅脑损伤机理仍需进一步研究。

基于上述研究现状，本文中针对颅骨局部弯曲变形这个致伤机制，开展进一步的深入研究。首先，建立充分反映颅脑生理结构的三维头部数值模型，并利用头部碰撞试验数据对头部模型的有效性进行验证。然后，基于欧拉-拉格朗日耦合法（Euler-Lagrangian coupling method，CEL）建立爆炸冲击波-头部流固耦合模型，对头部在正面冲击工况下的动态响应过程进行数值模拟，确立颅骨局部弯曲变形与脑组织压力之间的内在联系，初步揭示脑组织的损伤机理。

1 三维头部数值模型与验证

1.1 三维头部有限元模型

在颅脑结构中，脑组织包括大脑、小脑与脑干3 部分，其表面有3 层被膜包裹，由内向外依次为软脑膜、蛛网膜和硬脑膜。软脑膜薄而富含血管，紧贴脑表面并深入脑的沟裂中。蛛网膜是半透明结缔组织薄膜，它与软脑膜之间是蛛网膜下腔，内含脑脊液（cerebral spinal fluid, CSF）、蛛网膜小梁和较大血管。硬脑膜呈套状包被脑，在颅腔内分化为大脑镰与小脑幕，两者将颅腔分为左、右与下3 个腔体，限制脑组织相对于颅骨的过大滑移。颅骨与硬脑膜外侧紧密结合，其内、外为刚度与密度较大的骨密质，中间为骨松质，为典型的类三明治结构。图1 显示了人体头部的解剖结构。

图 1 头部解剖结构示意图[20]Fig. 1 Schematic diagram for anatomical structures of head [20]

基于颅脑的核磁共振切片，建立了三维头部有限元模型，该模型包括脑组织（大脑、小脑与脑干），膜结构（软脑膜、蛛网膜、硬脑膜、大脑镰与小脑幕），包围脑组织的脑脊液以及颅骨、皮肤与颈部。该头部模型共有303 588 个单元与530 954 个节点，单元尺寸为2～3 mm，充分考虑了人体头部的生理结构与相关细节，如图2 所示。

1.2 本构模型与材料参数

对于大脑、小脑与脑干等脑组织，采用黏超弹性本构模型进行表征[21]。利用Mooney-Rivlin 模型表征其超弹性行为，对应的应变能函数定义为：

图 2 三维头部有限元模型Fig. 2 3D finite element model of human head

表 1 黏超弹性材料参数Table 1 Material properties for hyper-viscoelastic model

脑脊液与水的性质较为接近，在蛛网膜下腔内流动；该腔是由蛛网膜小梁支撑的腔体，能够传递剪切变形，如图1 所示。为了便于建模，头部模型中的脑脊液是物理脑脊液、蛛网膜小梁与较大血管的综合体。对于脑脊液CSF 的本构模型，文献[22-23]中直接采用剪切模量相对很小的线弹性模型来进行表征。为了真实反映脑脊液与蛛网膜小梁的力学特性，模型中采用Mie-Grüneisen 状态方程表征脑脊液的不可压缩特性，同时采用较小的剪切模量（22 kPa）表征蛛网膜小梁的剪切传递特性[24-26]。Mie-Grüneisen 状态方程表示为：

表 2 线弹性材料参数Table 2 Material properties for liner elastic model

1.3 接触关系与边界条件

基于颅脑的生理结构并参考相关文献[23-28]，三维头部数值模型中的各部分组织之间的接触关系定义为绑定连结。同时，采用状态方程的脑脊液实现了脑组织与颅骨之间的有限滑动，符合颅脑的生理结构。而对于头部模型中颈部的边界条件，至今没有明确统一的结论。Ganpule 等[23]直接采用固定边界条件，对头部碰撞实验[29]进行数值模拟，两者的结果具有相对较好的吻合性。Chen等[27]对比分析了自由边界条件与固定边界条件两种工况下的模拟结果，虽然均与实验结果存在一定的差异，但是自由边界对应的模拟结果更能真实反映头部的动态响应过程。躯体对颅脑动态响应的影响是与颈部的肌肉状态密切相关的。当肌肉处于紧绷状态，颈部较为“刚硬”，躯体对颅脑动态响应的影响较大；而当肌肉处于放松状态时，颈部则较为柔软，这种影响作用相对较小。通常情况下，颈部肌肉是处于放松状态的。本文中基于显式有限元程序，对包含躯体的完整模型与自由边界条件下的头部模型在撞击载荷作用下的动态响应过程进行了数值模拟。其中，颈部采用较小的弹性模量（3.54 MPa），表征颈部肌肉的放松状态。图3 对比显示了两种数值模型对应的前额处与枕部处脑组织的压力曲线。由图3 可知，采用自由边界条件的头部模型，其模拟结果能够较好地吻合完整模型的结果。这表明：当颈部处于放松状态时，躯体对颅脑动态响应过程的影响在碰撞等冲击载荷条件下（通常为几个毫秒）基本上是可以忽略的，此时的颅脑动态响应过程呈现出明显的波动效应。同时，两种数值模型对应的前额处脑组织压力在t=7 ms 时开始产生差异，并且呈现出不断扩大的趋势。在这个过程中，颅脑的动态响应会由应力波主导转化为整体结构响应，受到躯体对颈部约束作用的影响。

图 3 前额处与枕部处脑组织的压力曲线Fig. 3 Pressure curues of brain tissue at the forehead and occiput

1.4 头部模型验证与数值模拟

Nahum 等[29]进行的头部碰撞实验，已成为头部数值模型有效性验证的经典实验。该实验利用运动块体撞击前额部位，利用传感器记录特征部位（前额、颅顶、枕部与颅后窝处）的颅内压力变化过程。由于该头部碰撞实验没有提供与头部发生碰撞块体的材料力学参数，参考Ganpule 等[23]进行头部数值模型验证时所使用的方法：将实验中测量所得到的碰撞接触力转化为平均压力作用于前额部位，作用面积为1 470 mm2。对应的压力曲线如图4 所示。基于上述头部碰撞模型，可获取颅内压力的动态响应过程，与头部碰撞实验数据进行对比分析。为减小采样位置对模拟结果与实验数据对比分析的影响，在前额、颅顶、枕部与颅后窝处各取3 个特征点（P1、P2 与P3），特征点位置如图5 所示。

图 4 碰撞压力曲线Fig. 4 Impact pressure curve

本文中提取特征点处颅内压力的模拟结果，与头部碰撞实验数据进行对比，两者具有较好的吻合性，如图6 所示。由碰撞实验结果可知，前额与颅顶处的颅内压力在整个碰撞过程中基本上均为正压，表明这两处脑组织处于受压缩状态。同时，前额处的颅内压力要明显高于颅顶处的相应值，这与头部简化模型的理论解吻合[30]。相比而言，枕部与颅后窝处的颅内压力在碰撞前期阶段为负压，表明脑组织处于拉伸状态，相应的压力幅值较为接近。但是，这两处的颅内压力在后期阶段会发生“反转”，由负压转变为正压，正压幅值相比于负压幅值较小。基于三维头部数值模型所得到的模拟结果，能够较好地预测颅内压力在碰撞条件下的演化过程，与实验数据吻合较好，从而验证了该头部数值模型的有效性。

为了对比分析头部动态响应在碰撞载荷与爆炸冲击波载荷作用下的差异，提取了头颈部在特征时刻的压力云图，如图7 所示。由于颅骨与脑组织的惯性力差异，在碰撞载荷作用下，颅骨加速度明显大于脑组织加速度。因此，前额处（撞击处）颅骨与脑组织之间产生相对径向压缩位移，使得脑组织处于压缩状态，对应的压力为正；而枕部处（对撞处）两者之间产生相对径向拉伸位移，脑组织处于拉伸状态，对应的压力为负。随着撞击载荷的持续作用，颅内正、负压的作用范围持续扩大，发展成为靠近前额的前半部分脑组织受压而靠近枕部的后半部分脑组织受拉的受力状态。随着外部载荷的减弱和头颈部之间相对变形的增大，头部运动会在某一特定时刻发生转向，颅骨与脑组织的惯性差异使得两者之间的相对位移发生转变，致使脑组织压力的状态发生“反转”。但是，“反转”后的脑组织压力，其幅值则相对较小。因此，颅骨与脑组织间的相对径向位移决定脑组织压力，撞击处前半部分脑组织主要处于压缩状态，而对撞处后半部分脑组织则主要处于拉伸状态。

图 6 颅内压力模拟结果与实验结果的对比Fig. 6 Comparison of intracranial pressure between experiments and calculations

图 7 前额撞击时脑组织压力云图Fig. 7 Nephogram of brain pressure for the forehead collision

2 爆炸冲击波-头部流固耦合模型

为了获得头部在爆炸冲击波下的动态响应规律，需要建立爆炸冲击波-头部流固耦合模型。头部在常规爆炸冲击波作用下的变形相对较小，可采用拉格朗日单元进行网格划分。而空气需要模拟爆炸冲击波的产生、传递以及与头部相互作用的过程，对应区域应当采用欧拉网格，避免单元过度变形。爆炸冲击波与头部在接触边界上的相互作用，需要通过求解流固耦合方程进行确定。通用有限元软件ABAQUS 通过非线性瞬态程序与欧拉-拉格朗日耦合法（Euler-Lagrangian coupling method，CEL），对整体系统的三大守恒偏微分方程（动量、质量与能量）同时进行求解。欧拉区域能够对导致严重网格扭曲的高度动态事件进行模拟（如爆炸冲击波），并为拉格朗日区域提供了压力边界条件；拉格朗日区域可完全或部分位于欧拉区域内，拉格朗日表面为欧拉区域提供边界条件，该边界条件不允许在其表面的法线方向上有流动产生。同时，增强浸入边界方法提供欧拉区域和拉格朗日区域之间的耦合作用。因此，固定欧拉网格与基于增强浸入边界法的流固界面模型，能够模拟爆炸冲击波与头部之间耦合作用，获得头部的动态响应规律。

2.1 有限元离散

模拟爆炸冲击波的欧拉区域采用正方体形状，边长为600 mm，这个尺寸的选择充分考虑了冲击波与头部的作用时间，防止冲击波边界反射对头部产生非真实的二次冲击作用。经过验证的三维头部数值模型，在爆炸冲击波-头部流固耦合模型中可以直接使用，并放置在欧拉区域的中心位置。欧拉区域使用缩减积分六面体欧拉单元（EC3D8R）进行网格划分，并采用分块划分单元网格的方法：与头部发生接触的区域采用较小的单元尺寸，利于模拟爆炸冲击波与头部之间的流固耦合作用；没有与头部发生作用的区域则采用相对较大的单元，利于节约计算成本。同时，针对单元尺寸对模拟结果的影响进行参数化研究。研究表明：当中心区域的单元尺寸为3 mm 时，相应的模拟结果能够收敛，与文献[23]结论基本一致。因此，将单元尺寸定义为3 mm，欧拉区域共有8 524 800 个单元。欧拉区域采用理想气体状态方程来表征空气，空气密度为1.18 kg/m3，气体常数为287.04 J/(kg·K)，比热容为1 000 J/(kg·K)。最终，建立的爆炸冲击波-头部流固耦合模型如图8 所示。

2.2 荷载与边界条件

根据我国工程设计规范中的设计准则，爆炸冲击波的超压值（以大气压力为参考零值）会随时间呈现指数型衰减。其理论表达式为：

图 8 爆炸冲击波-头部的流固耦合模型Fig. 8 Fluid-solid coupling model of explosive blast wave-head

基于上述理论公式，可获得炸药当量为7 kg且距离爆心3.8 m 处的压力曲线。将理论结果与实爆实验结果进行对比，两者具有较好的吻合性，如图9 所示。

本文中，采用这个典型爆炸冲击波压力曲线作为欧拉区域的输入载荷，直接作用在入口处。欧拉场的出口处采用无反射边界条件，冲击波可以直接从这里逸出，不再对欧拉场产生影响。而欧拉场的环向4 个边界则采用法向速度约束边界条件，防止冲击波从侧面逸出，模拟平面型爆炸冲击波与头部的相互作用，可真实反映现场环境。

图 9 爆炸冲击波压力曲线的理论结果与实验结果对比Fig. 9 Comparison of pressure curves of explosion blast between theoretical and experimental results

2.3 接触与初始条件

由爆炸冲击波-头部流固耦合模型的理论基础可知，两者之间的相互作用需要借助基于增强浸入边界法的流固界面模型来实现。在程序中，通用接触（general contact）采用增强浸入边界法能够较好地处理拉格朗日单元与欧拉单元之间的接触问题。因此，通用接触是进行流固耦合分析的必备条件。

对于初始条件，系统中有两个方面需要进行合理设置。（1）欧拉区域中空气的初始温度与参考压力。这里设定初始时刻的空气温度为303 K，以标准大气压力作为参考值，即欧拉区域的压力初始值为零。（2）空气材料的初始分布。这里可以使用体积分数工具（volume fraction tool）来设定空气在离散欧拉场中的分布，保证每个位置点处只有一种材料。

3 结果与分析

基于上述爆炸冲击波-头部流固耦合模型，可以获得头部在爆炸冲击波作用下的动态响应规律，进而揭示颅脑的损伤机理。

3.1 流场分析

爆炸冲击波从正面作用于头部时，流场压力的整个演化过程如图10 所示。在t=0.375 ms 时，平面冲击波与头部即将发生接触，此时的冲击波压力峰值为170 kPa。随后，冲击波与头部发生相互作用，反射叠加效应使其压力峰值急剧上升；当t=0.450 ms 时，前额处的压力峰值达到600 kPa，是入射冲击波峰值的3.5 倍。在t=0.720 ms 时，冲击波经过前额进入颅顶，曲面型的颅顶结构使其在该处发生边界层分离现象。因此，冲击波对颅顶处的压力峰值相对较小。当t=0.960 ms 时，由于头部阻挡作用而分散的冲击波会在后脑部位重新汇集，该处的压力值迅速上升，从而对头部造成二次冲击。相比于面部首次冲击的压力峰值，二次冲击的对应峰值出现了显著下降。

图 10 正面冲击时流场压力分布Fig. 10 Pressure distribution of flow field in blast frontal impact

图 11 正面冲击时眼部的流场压力分布Fig. 11 Pressure distribution of flow field around the eyes in blast frontal impact

3.2 脑组织压力

对于正面冲击，先提取脑组织特征点处的压力时程曲线，如图12 所示。由图12 可知，各特征点处的脑组织压力均呈现出明显的波动特性，并且随着时间的推移，这种波动特性会逐步弱化。同时，从前额至颅顶处，特征点1～4 的脑组织压力在前期阶段呈现出明显的周期特性，其中前额处的频率高达8 kHz；而特征点5～8 则没有明显的周期特性。在爆炸冲击波下，脑组织压力的高频波动特性与碰撞载荷下的脑组织响应是完全不同的，这就决定了针对两种不同工况需要采用不同的颅脑防护策略。结合脑组织压力的动态演化过程与损伤判断准则可知，多处脑组织的压力将超过重度损伤阈值（235 kPa）[31]，造成脑组织出现多处重度损伤。

图 12 正面冲击时特征点处脑组织的压力曲线Fig. 12 Pressure curves of brain tissue at feature points in blast frontal impact

基于模拟结果，可提取颅脑在特征时间点处的压力分布云图和颅骨变形云图，如图13 所示。在颅骨变形图中，带有网格的颅骨为初始形态，没有网格的颅骨为当前形态，采用的放大系数为500。由压力云图可以看出，前额处脑组织在冲击波作用下起初处于受压状态（t=0.465 ms），而后转化为受拉状态（t=0.520 ms）；在后续过程中，脑组织压力会在压、拉两种状态之间进行转化。紧邻颅骨内表面处的脑组织，其压力沿着颅骨内表面出现正、负压交替现象，形成较大的压力梯度；而远离颅骨内表面的脑组织，其压力梯度则相对较小。同时，具有典型三明治结构的头部颅骨，利用材料阻抗不匹配特性截断并分散冲击处的能量，从而为脑组织提供防护。

从流场压力的演化过程可知，爆炸冲击波主要作用于前额等面部结构，而颅顶等部位由于存在边界层脱离效应，受到冲击波的作用相对较小，在进行颅骨变形分析时可以忽略这个影响。当t=0.465 ms 时，流固耦合作用使得爆炸冲击波达到最大值，冲击波直接挤压皮肤与颅骨，前额处颅骨的“平面”变形造成该处脑组织处于挤压状态，体积压力为正压。随着前额颅骨变形的进一步增大，当t=0.520 ms 时，该处颅骨的“平面”变形转变为“内凹”弯曲，致使该处与颅骨内表面紧密连结的脑脊液与脑组织处于受拉状态，体积压力为负压。而紧邻前额处的颅骨则相当于发生了“外凸”弯曲，对该处的脑脊液与脑组织进行挤压，体积压力为正压；变形的连续性使得颅顶处颅骨产生较小的“内凹”弯曲，体积压力为负压。当t=0.585 ms 时，冲击波作用的减弱与颅骨的弹性边界使得前额处颅骨的“内凹”弯曲在达到最大值后发生回弹，由“内凹”弯曲转变为“外凸”弯曲，并主导着颅骨其他部位弯曲类型的改变。对比t=0.520 ms 与t=0.585 ms 时的云图可知，各位置处的脑组织的拉压状态发生转变。对比分析后续特征时刻的脑组织压力与颅骨弯曲变形，可以得到相同的结论：理想状态下，颅骨发生“内凹”弯曲，该处脑组织处于受拉状态；颅骨发生“外凸”弯曲，脑组织则为受压状态。而真实条件下，颅骨内表面处脑组织在特定位置处的受力状态同时会受到周边颅骨变形的影响。

图 13 正面冲击时颅脑压力云图与颅骨变形云图Fig. 13 Nephogram of brain pressure and skull displacement in blast frontal impact

本文中进一步提取了前额处颅骨的加速度时程曲线，并将该处特征时刻的加速度标注在该曲线上，如图14 所示。其中实心圆点表示该处脑组织压力为正时的加速度，空心圆点则表示压力为负时的加速度。由图14 可见，爆炸冲击波致使前额颅骨发生局部高频振动，其频率与该处脑组织压力的波动频率相一致。前额颅骨最值加速度与颅骨弯曲变形模式、脑组织状态均具有较好的相关性：颅骨受冲击方向加速度最大时，颅骨为“外凸”弯曲模式，脑组织处于受压状态；而颅骨受冲击方向加速度最小（冲击反方向最大）时，颅骨为“内凹”弯曲模式，脑组织处于受拉状态。

图 14 前额处颅骨加速度曲线Fig. 14 Acceleration curves of forehead skull

对比分析脑组织压力、颅骨局部弯曲变形与加速度的整个演化过程可知，前额处颅骨在爆炸冲击波作用下发生弯曲变形并引起局部振动，该处的局部变形会沿着颅骨进行“传播”，影响着颅骨的变化构型，从而使脑组织在外载作用下由前额处的“单一波源”逐步发展到整个颅骨的“多处波源”，这些波源共同决定脑组织压力的演化过程。颅骨局部弯曲致使脑组织沿着颅骨内表面出现拉、压交替状态，形成较大的压力梯度。同时，颅骨局部弯曲的振动频率较高，使得脑组织压力具有高频波动特性。