一种冲击波作用下结构毁伤算法研究*

周 朗，徐春光

（中山大学航空航天学院，广东 深圳 518107）

对于常规弹药而言，爆炸冲击波是造成结构破坏的主要因素，研究爆炸冲击波传播与演化规律，对准确评估弹药的毁伤效能具有重要意义。针对自由大气内爆炸冲击波的传播规律，研究人员构建了多个经验公式[1-4]，并应用于工程实际中，是研究爆炸冲击波传播与演化的重要手段，但对于建筑物等复杂空间，经验公式不能反映冲击波的演化过程，与实际物理过程存在比较明显的差异，且无法给出冲击波演化的细节，也不能反映冲击波与结构的相互作用过程及影响。

随着数值模拟技术和计算机技术的发展，研究人员利用数值模拟的方法对冲击波与障碍物存在相互作用的问题进行研究。刘伟等[5]在考虑地面因素的影响条件下，研究了不同当量TNT 装药空爆下冲击波超压规律，获得了有意义的结果。宁建国等[6]基于JWL 状态方程和瞬时爆轰模型，研究了存在防护挡墙情况下，爆炸流场的传播及演化规律，所得超压计算结果比已有经验公式计算结果偏高。王飞等[7]采用ALE 方法，对冲击波绕过障碍物后的超压分布情况进行了数值模拟，数值模拟结果能较好的反映超压变化的规律。张晓伟等[8]对爆炸冲击波作用于隔爆墙的绕射进行实验和数值仿真，得到墙后冲击波超压分布规律，并给出了墙后近地面超压峰值的工程计算公式。文献[9-12]对冲击波在坑道、隧道内的传播规律进行了研究。上述工作对冲击波与结构的相互作用进行了较为详细的研究，但研究中均假设结构或障碍物为刚体，没有考虑结构的变形和运动，难以描述结构在冲击波作用下的变形及破坏。

对于考虑冲击波作用下结构变形、运动及其对冲击波传播影响的问题，多采用国外软件进行数值模拟，譬如LS-DYNA、AUTODYN、ABAQUS 等软件。邓荣兵等[13]利用LS-DYNA 对爆炸冲击波与幕墙结构的作用过程进行了研究，分析了幕墙在冲击波作用下的动态响应特性和内外层玻璃的破坏情况，与实验结果符合较好。董湘乾[14]利用ABAQUS 对爆炸冲击波与混凝土框架结构的作用过程进行了模拟，分析了多层框架结构在冲击波作用下的动态响应。张秀华等[15]利用LS-DYNA 对爆炸冲击波与框架结构的作用过程进行了模拟，分析了爆炸冲击波的传播特性,结果表明合理的设置泄压口能一定程度的减小冲击波对结构的破坏。游俊[16]利用ABAQUS 对冲击波与破片联合作用下砌墙体进行了数值模拟，对结构动力响应进行分析，得到了冲击波传播以及与结构相互作用导致结构毁坏的损伤机理。LS-DYNA 等商业软件是目前此类问题数值模拟的主要手段，国内未见具有类似功能的算法软件，同时软件详细处理模型、方法也缺乏完整的介绍。

基于上述问题，本文发展一种冲击波作用下结构毁伤算法，并对爆炸冲击波与建筑物结构的相互作用过程进行模拟，研究冲击波作用下建筑物的毁伤范围、建筑物破坏以后冲击波的传播特性等，验证该算法的适用性。

1 冲击波与结构相互作用机理分析

构建如图1 所示的建筑物墙体简化模型，采用冲击波冲量[17]作为结构是否破坏的指标，分析冲击波与结构相互作用机理。图1 中将包含建筑物墙体等结构的整个计算空间全部划分为网格单元，共有1～4 个网格单元，计算每个单元受到的比冲量，判断该单元是否被破坏：

⑴ 如果载荷不足以破坏结构，建筑物单元按照刚性壁面假设处理，冲击波发生反射，形成比入射波高很多的超压；

⑵ 如果冲量超过破坏载荷则该单元毁坏，那么它就变成可以在冲击波作用下向前运动的碎块。

在结构单元已经破坏的情况下，实际物理过程为结构碎块在冲击波作用下飞散。对于建筑物类脆性材料，由于结构碎块四周没有固支，碎块的变形很小，可以将结构碎块作为刚体处理，则冲击波作用下结构碎块的飞散问题是动边界绕流问题。对于此类问题，可以采用非结构动网格技术进行模拟，即将流体力学方程与6 自由度弹道方程耦合求解，求解流体力学方程获得碎块表面压力，进而得到碎块受到的气动力，将其作为弹道方程的输入，求解弹道方程获得碎块的新位置，并据此对计算网格进行变形移动，获得新的计算网格，重复这一过程即可获得气动力作用下的碎块轨迹。Löhner 等[18]和Baum 等[-19]对大规模破片飞散过程进行数值模拟，模拟破片数量达到上百枚；爆炸波作用下的钢筋混凝土墙发生变形、断裂、碎裂过程[20]也得到了较好模拟。王巍等[21-23]对非结构动网格技术开展了大量研究，并将其应用于多种场景下多体运动问题的模拟，如飞行器头罩分离[21]，助推器分离[22]，子母弹抛壳[23]等，计算精度得到了确认。本文采用上述发展的非结构动网格技术及程序，对如图2 所示单块碎块在冲击波作用下的运动过程进行模拟，分析冲击波与结构相互作用机理。

计算模型如图2 所示，封闭空间一侧墙体上设置一洞口，洞口内存在一个堵块，堵块与墙体之间留有缝隙。其中堵块的尺寸为1.14 m×0.18 m×0.6 m，质量为50 kg；洞口尺寸为1.18 m×0.18 m×0.67 m；封闭空间尺寸为5.8 m×6.0 m×3.3 m。

图3 为计算得到的堵块运动位置和压力等值线。从图中可以清晰看出，在碎块完全脱离墙体之前，碎块运动的稀疏波引起超压降低；在碎块刚刚脱离墙体时，气体透过碎块与墙体之间的缝隙传播，由于缝隙面积小于洞口入口处面积，缝隙处气流速度达到声速，形成壅塞，此时冲击波既有绕射，也有反射，流动非常复杂；在碎块脱离墙体较远后，碎块和墙体之间缝隙形成的气流壅塞消失，爆炸点所在空间内部高压气体在洞口处形成壅塞，且由于碎块运动速度小于气体速度，大量气体绕过碎块，在绕过碎块后规整成球面波。

根据以上分析，图4 给出了冲击波与结构碎块作用机理，冲击波与结构碎块作用过程中主要存在如下物理效应：(1) 碎块运动会产生稀疏波效应，如图4(a)示意；(2) 碎块刚脱离墙体后的壅塞效应，如图4(b)示意；(3) 绕射波重新规整为平面波的效应，如图4(c)示意。从建立数学模型的角度看，能够刻画以上机理的模型就可以模拟冲击波与结构碎块作用过程。

上述结果表明，在建筑物结构破坏以后，利用动网格技术能够模拟冲击波与结构碎块的相互作用过程，能够获得较为真实的数值模拟结果，并且可以直接利用现有的动网格计算模块进行计算。但在实际应用中，对于整栋楼的数值模拟，结构单元的破坏会产生成百上千，甚至上万的结构碎块。利用动网格技术处理，需要进行网格变形计算，并且由于每个碎块独立运动，使得每时每刻都可能会出现局部网格变差的情况，就需要将计算程序停止，进行局部重构和插值，计算量巨大，大大降低了计算效率。因此，尽管非结构动网格技术能够很好地处理冲击波推动结构碎块飞散的问题，技术上也不存在不可克服的难点，但对计算机硬件要求太高，不利于工程应用，需要更高效率的计算方法。

2 模型建立

从上一节分析可知，冲击波作用下的结构毁伤过程，存在结构未破坏形成的冲击波刚性壁面反射、结构完全破坏后形成的冲击波透射以及结构破坏过程中因碎块运动形成的稀疏波和壅塞效应。对于未被破坏结构，按照刚性壁面假设处理，不是本文关注的重点。对于已破坏结构单元，冲击波在碎块下游的绕射波重新规整为平面波效应有比较好描述，不需要特别的处理。因此，碎块引起的稀疏波和壅塞两种效应是建模的关键。

2.1 运动稀疏波效应的模型

采用任意拉格朗日-欧拉方法（ALE）坐标描述三维无量纲可压缩非定常流动的Euler 方程：

式中：S是控制体V的表面，n为其外法向单位矢量；Q为流体守恒变量；Fc(Q,xc) 表示相对流通矢量，xc是网格运动速度。

提出如下模型：如图5 所示，假设结构单元2、3 已经被破坏，在计算网格位置不变情况下，根据流固界面耦合算法、流场物理量和结构单元质量等参数，可以计算得到结构单元2、3 的网格运动速度 →xc，将其反映在流体控制方程中，计算流场参数，即可模拟结构单元2、3 运动引起的泄压过程。在该计算过程中，尽管结构单元2、3 的网格位置没有变化，但其单元运动速度真实的反映在流体控制方程中，从而在一定程度上模拟了结构单元运动引起的泄压作用。

2.2 碎块壅塞效应的计算模型

利用文献[24]提出的虚拟网格透气技术解决壅塞效应这一问题，图6 为虚拟网格透气技术原理图，图中单元1～3 为描述建筑物墙体结构的网格单元，单元4 为冲击波一侧流场中的气体网格单元。完全气体假设下的控制方程：

式中：Φ 的形式可简写为：

对于计算区域中单元4 物理量的计算可以表述为如下半离散形式：

式中：下标k用以表示不同表面。

对于图6 所示情况，令单元4 与2 的边界面为Sb，计算中采用如下形式：

式中：Φb为边界面上的对流项，Sb表示单元4 和2 边界投影面积。由于单元2 为结构单元，通常计算中式(5)中的Sb只能表示为固壁边界，完全不通气，不能模拟结构单元2 破坏后，冲击波透过破坏单元传播的过程，因而也就不能模拟碎块破裂后的壅塞效应。

为了实现对结构单元破坏后边界面由不透气逐渐变为透气，并最终达到完全透气问题的模拟，文献[24]引入如下计算模型：

式中：Sf、Sw分别是边界面Sb通气部分和不通气部分的面积，即将边界面分为通气和不通气两部分，并计算对流项Φb：

式中：Φf是作为内部单元面处理时的对流项，可按一般式处理；Φw是作为固壁处理时的对流项。其对表面的通量积分为：

式中：U为流体相对于网格的速度；at为网格运动的法向速度；nx、ny、nz是n的三个分量;Sw为面积矢量。

对于具有虚拟挡板表面的单元，应该针对该表面通气和壁面两种状态分别求解变量梯度，计算Φf和Φw时分别采用相应的梯度进行基本变量的重建。

将式(7)代入式(5)，则有如下形式：

分析式(9)，可以看出当Sf=0，即图5 中单元4 与2 之间边界面完全为固壁条件时，式(9)退化为式(5)，边界面不透气；当Sw=0 ，即单元4 与2 之间边界面完全为流场内部面时，式(9)退化为式(4)，边界面完全通气；当处于边界面状态处于两者之间时，单元4 与2 之间的边界面可以实现按比例通气。

根据第1 节中图4 的分析结果，定义系数к为如下形式：

式中：S1为碎块与墙体之间缝隙环带的面积，可由碎块运动速度、碎块周长及运动时间等计算获得；S2为碎块的在来流方向的法向投影面积，固定不变，如图7 所示。将к代入式(9)，S2代替其中Sb，则有如下形式：

当碎块处于图4(a)状态时к=0，单元4 与2 之间边界面完全为固壁条件，完全不通气；当碎块处于图4(b)状态时，к表示缝隙环带面积与结构在来流方向的法向投影面积比值，模拟了缝隙对通气量的限制；当碎块处于图4(c)状态时к=1，单元4 与2 之间边界面变为完全通气状体，此时结构碎块对绕射后形成的第一道冲击波强度影响很小。式(11)即为网格位置不变条件下，结构对冲击波壅塞效应计算模型。

3 结果与分析

3.1 与非结构动网格模拟结果对比确认

对本文所提方法的计算过程，采用如图6 所示的计算模型，模型的形状、参数与图2 中的模型一致，但洞口处没有堵块，不存在动网格，而是设置一个“虚拟挡板[24]”。计算过程中，当冲击波到达洞口后，采用本文提出的计算模型，对冲击波的传播过程进行模拟，监测与图2 中模型相同位置处超压变化，并进行对比，分析验证计算模型的准确性。其中，具体监测点位置A、B如图8 所示，在墙壁中心位置；给定高压气团球心位置时，坐标定义如图9 所示，坐标原点在堵块中心。

针对上述两个模型，模型1（堵块模型），模型2（虚拟透气计算模型），对不同参数设置的算例进行计算，算例情况见表1。

表1 算例参数设置Table 1 Example parameter setting

图10～图13 为不同算例的监测点A、B压力变化，可以看出两种模型情况下得到的监测点压力峰值及规律均符合较好，表明本文提出的计算模型具有较高的准确度，能够用于此类问题的模拟。

3.2 典型建筑物作用过程模拟结果分析

利用本文构建的爆炸冲击波与建筑物结构作用计算模型，对冲击波作用下建筑物的毁伤及冲击波传播过程进行了模拟，获得了建筑物毁伤效果和冲击波超压分布。

3.2.1 计算模型

如图14 所示，构建了一个大楼模型。为了提高计算效率，横向宽度只包含3 个房间。其中，楼房左侧的3、4、5 层建立结构计算单元，即墙体、楼板等结构部件内部建立网格单元，计算中根据结构单元毁伤载荷判据计算结构单元是否破坏，若破坏则变为透气单元，楼房的其他部分不参与结构计算。计算中装药当量、结构毁伤载荷根据相关设计参数确定。

3.2.2 计算结果与分析

图15 为建筑物结构毁伤计算结果，给出了结构单元破坏的过程。从计算结果来看，406 房间左右隔墙、上下楼板和走廊隔墙完全崩塌；405、407 房间外侧墙体大部分崩塌；4 层走廊楼板大部分崩塌；306 房间下楼板和506 房间上楼板局部崩塌，较为真实地反映了试验中建筑物的毁伤效果（其中，中间的房间是406，左边405、右边407，上边506、下边306）。计算中对每个单元独立判断其是否破坏，只与超压分布有关，由于房间冲击波反射、相交影响，使得超压分布非常不均匀，导致出现墙体破坏也非常不均匀，且计算模型只考虑结构是否破坏透气，未考虑结构间的支撑及运动情况，出现了四周均破坏而中间单元未破坏的情况，但计算结果仍具有较高的参考价值。

图16 为爆炸流场演化云图，与图12 对照，可以看出冲击波传播与结构破坏特征一致，并且由于“虚拟挡板”对通气量的限制，刚开始透过破坏单元传播的冲击波强度较弱，表明本文计算模型较好的模拟了结构碎块对冲击波的阻碍效应。

图17 为406 房间监测点冲击波超压计算结果，由于没有毁伤效应试验超压测量数据，将计算结果与无限空气中爆炸冲击波超压经验公式计算结果进行对比，其中TNT 爆炸当量为120 kg。表2 为经验公式计算结果，经验公式采用Baker 公式[25]，具体如下：

表2 监测点超压经验公式计算结果Table 2 The numerical calculation results of shock wave overpressure under different monitoring point

正反射超压公式如下：

式中：w为TNT 药量，kg；r为距离，m；Δp为超压，kg/cm2。

比较图14 与表2 计算结果，可以看出在爆炸房间内，冲击波反射超压数值计算结果与经验公式结算结果符合较好，表明本文采用的爆炸流场计算方法具有较好的计算精度。

4 结 论

本文针对爆炸冲击波与建筑物结构作用过程，在分析冲击波与结构碎块作用机理的基础上，构建了一种计算模型，能够在计算网格不动的条件下模拟冲击波与结构的作用过程，主要结论如下：

(1)根据非结构动网格技术计算冲击波推动结构碎块飞散模拟中，非定常流场的变化规律，分析得出冲击波与碎块之间相互作用机理主要分为碎块运动会产生稀疏波效应、碎块刚脱离墙体的壅塞效应、绕射波重新规整为平面波效应。

(2)在机理清晰认识的基础上，本文采用“虚拟运动”模拟稀疏波效应、“虚拟网格透气技术[1]”处理碎块和墙体之间的壅塞效应，发展了一种新的模拟爆炸冲击波与建筑物作用过程的模型。通过算例对比，该模型计算结果与非结构动网格技术计算结果符合较好，表明该模型能够用于冲击波作用下结构破坏以及冲击波传播过程的模拟，具有良好的工程应用价值。

(3)建立建筑物模型，采用本文构建的爆炸冲击波与建筑物结构作用计算模型 ，模拟了对建筑物的毁伤及冲击波传播过程，获得了建筑物毁伤效果和冲击波超压分布。