复杂结构内爆炸冲击波传播的大规模计算方法研究.pdf

复杂结构内爆炸冲击波传播的大规模计算方法研究王成1，王雪桃1，李涛1，张志刚2,3，孙伟福1（1.北京理工大学 爆炸科学与技术国家重点实验室，北京100081；2.中煤科工集团重庆研究院有限公司，重庆400039；3.煤矿灾害防控全国重点实验室，重庆400039）摘 要：为研究爆炸冲击波在矿井巷道等复杂地下密闭空间内的传播过程及衰减规律，提出了一种基于数据流的分块续算方法.该方法解决了计算资源的大量浪费问题，实现了大尺寸复杂地下密闭空间内冲击波传播的大规模、高效率、高精度数值模拟.方形坑道内爆过程的数值结果与实验结果基本一致，证明了该计算方法的准确性.标准模型算例全尺寸计算与分块续算数值模拟结果的误差在 0.33%以内，且续算方法在给定算例下相比普通全尺寸计算节约了 88.68%的计算核时，证明该方法不仅能保证计算过程的准确性，而且还极大地提高了计算效率.同时，对一个包含长方体型工房和圆柱形分叉泄压洞的大尺寸复杂结构内的爆炸问题进行了数值模拟，结果表明：冲击波超压在泄压洞内是以波动衰减的趋势逐渐减小的，且其衰减速度比在自由场的小；T 型分流时管道间采用直角交汇或缩小管道直径在一定程度上有利于冲击波超压的衰减.数值模拟结果表明，基于数据流的分块续算方法可用于地下空间等大尺度复杂结构防爆抗爆设计的理论指导.关键词：爆炸冲击波；复杂结构；数值模拟；分块续算中图分类号：O382.1 文献标志码：A 文章编号：1001-0645(2024)01-0039-12DOI：10.15918/j.tbit1001-0645.2023.078Study on Large Scale Calculation Method of Explosive Shock WavePropagation in Complex StructureWANG Cheng1，WANG Xuetao1，LI Tao1，ZHANG Zhigang2,3，SUN Weifu1（1.State Key Laboratory of Explosion Science and Technology,Beijing Institute of Technology,Beijing 100081,China;2.Chongqing Research Institute,China Coal Research Institute,Chongqing 400039,China;3.State Key Laboratoryof Coal Mine Disaster Dynamics and Control,Chongqing 400039,China）Abstract：In order to study the propagation process and attenuation law of explosion shock wave in complex un-derground confined space such as mine tunnel,a block continuation calculation method based on data flow wasproposed.This method solved the problem of huge waste of computational resources and realized large-scale,high-efficiency and high-precision numerical simulation of shock wave propagation in large-size complex under-ground confined spaces.The numerical results of the implosion process in the square trench were basically con-sistent with the experimental results,which proved the accuracy of the calculation method.The error between thefull-scale calculation of the standard model case and the block continuation numerical simulation was within0.33%,and the continuation method saved 88.68%of the computational core time compared with the ordinaryfull-scale calculation for a given case,which proved that the method not only ensured the accuracy of the calcula-tion process,but also greatly improved the computational efficiency.At the same time,the explosion of a largecomplex structure containing a rectangular room and a cylindrical bifurcated pressure relief cavity was numeric-ally simulated.Simulation results show that:shock wave overpressure in the pressure relief cavity decreased 收稿日期：2023 03 31作者简介：王成（1972）男，教授，博士生导师，E-mail:.第 44 卷第 1 期北 京 理 工 大 学 学 报Vol.44No.12024 年 1 月Transactions of Beijing Institute of TechnologyJan.2024gradually with a trend of fluctuating decline,and its decay rate was smaller than that in the free field;T-shapedmanifold with right-angle intersection between pipes or reduced diameter of the pipe was to a certain extent con-ducive to the attenuation of shock wave overpressure.Numerical simulation results show that the data flow-basedblock continuation calculation method can be used for theoretical guidance of explosion-proof and blast-resistantdesign of large scale complex structures such as underground spaces.Key words：blast wave；complex structure；numerical simulation；block continuation 爆炸冲击波是矿井等复杂密闭建筑结构发生爆炸事故产生的主要灾害效应之一，冲击波不仅会造成人员伤亡，还会对周围结构及建筑物造成大规模破坏.空中自由场等开放空间发生爆炸事故时，由于爆炸冲击波未受任何约束，通常以指数形式快速衰减.然而，矿井巷道等密闭空间内爆炸产生的冲击波由于受到壁面等结构的约束作用，在其扩散过程中会发生多次的反射、折射等效应，并随之发生多次叠加作用，使得冲击波超压峰值较多且难以衰减至常压，在很大程度上提高了爆炸压力和能量.因此，爆炸冲击波在复杂地下密闭空间内的传播规律是矿井、巷道等制定防爆及抗爆安全的主要科学依据，在爆炸事故的快速应急响应及防护救援等过程中具有极其重要的理论意义和实际应用价值.关于类似巷道、矿井等复杂密闭结构内的爆炸问题，陈海天等1、庞伟宾等2 3及其他研究者4 6通过实验对其进行了理论研究，并基于自由场爆炸冲击波压强预测公式，得到了直边墙拱形隧道、方形坑道、T 型坑道以及不同分叉角度的坑道内爆炸冲击波的衰减规律.但由于巷道在弯曲处和交叉处的冲击波流场分布极其复杂，理论研究无法对其传播过程进行直观的反映，所以仅通过实验对巷道内爆问题进行深入细致的研究存在许多实际困难.因此，近些年来数值模拟技术在结构内爆问题中得到了广泛应用.张奇7、屈康康等8及其他研究者们9 18利用Autodyn、LSDYNA 等有限元软件对长直坑道、L 形坑道、S 形坑道、分叉管道、变径管道等不同形状巷道内的冲击波传播规律进行了研究.樊保龙等19、王成等20对密闭巷道内瓦斯的爆炸过程进行了数值分析.胡洋等21、宋诗谦等22、张晓伟等23及吴文溢等24对密闭空间内 TNT 爆炸效应进行了数值模拟分析.张学博等25提出了冲击波传播的分段接力数值模拟法，并实现了矿井大尺度爆炸问题的二维模型模拟冲击波传播.杨科之等26通过自主计算程序对简单巷道内爆流场进行了数值计算，并得出冲击波在巷道内衰减的计算公式，但是 Autodyn 等商业软件无法对炸药放置于地面或近壁面的爆炸过程进行高效计算，也无法对近、中场的爆炸问题进行准确计算27.另外，目前现有的自主编制爆炸流场的程序更多地是对简单、小尺寸的爆炸问题进行计算，若对大尺寸的复杂密闭空间结构进行数值模拟，则存在一定的难度和局限性.对于类似于地下矿井等大尺度的复杂密闭空间结构，若在计算时整体计算域采用较粗的网格，则会导致计算结果不准确；若全流场采用较为精细的网格尺寸，则会导致由于网格数目过多而无法正常计算，即使利用高性能的计算设备和大规模的并行处理器，也会浪费大量的内存和计算资源.针对大尺寸复杂结构内爆数值计算浪费大量计算资源的问题，在此提出了一种基于数据流的分块续算方法，并通过文献实验结果以及标准模型算例验证了该计算方法的准确性和有效性，实现了复杂地下密闭空间内爆问题的大规模、高效率、高精度数值模拟.同时，以复杂地下密闭空间结构中的爆炸问题为研究对象，采用基于数据流的分块续算方法对一个包含长方体形工房、圆柱形泄压洞的复杂结构内的爆炸冲击波传播过程进行了数值模拟，并给出了冲击波在巷道内的传播规律.1 数学模型及算法 1.1 数学模型由于炸药在爆轰传播过程中，冲击波的传播速度很快，因此分子扩散、热传导及黏性效应在模型计算中可以忽略.控制方程组为由质量守恒、动量守恒和能量守恒组成的非定常可压缩三维 Euler 方程组：Ut+Fx+Gy+Hz=0,U=uvwET,F=uu2+puvuwu(E+p)T,G=vuvv2+pvwv(E+p)T,H=wuwvww2+pw(E+p)T（1）pu、v、w式中：为密度；为压强；为速度矢量分量；40北 京 理 工 大 学 学 报第 44 卷Ee为单位质量的总能量，是单位质量的比内能 与动能的总和，满足E=e+12(u2+v2+w2)（2）为了使上述控制方程组封闭，需引入描述压强、密度和比内能之间关系的状态方程.爆轰产物采用JWL 状态方程，空气介质采用理想气体状态方程.JWL 状态方程为：p=A(1R10)exp(R10)+B(1R20)exp(R20)+e（3）peABR1R20式中：为压强；为密度；为单位质量的比内能；、和为爆轰产物的材料常数.理想气体状态方程为：p=(1)e（4）pe式中：为压强；为密度；为比内能；为多方气体指数，对于空气 一般取为 1.4.1.2 数值计算方法在复杂结构内爆炸冲击波流场的数值模拟过程中，为确保冲击波的追踪精度，采用 5 阶 WENO（wei-ghted essentially non-oscillatory）高精度有限差分格式28进行空间离散，采用 3 阶精度的 TVD Runge-Kutta 方法进行时间离散.针对多介质爆炸流场问题，采用RGFM（real ghost fluid method）29和 Level-set30 31 相结合的物质界面处理方法，通过 Level-set 方法对界面进行隐式追踪，通过求解法向黎曼问题确定物质界面附近流场相互接触作用之后的流动状态，并利用 RGFM 方法对界面处流场进行处理，实现了复杂地下密闭空间结构内爆轰产物及爆炸冲击波传播过程的高精度数值模拟.另外，通过采用一种基于 WENO格式的自适应扩大计算域高精度算法32，提高了多介质界面追踪与处理的精度和分辨率，加快了计算效率.在此基础上，引入 MPI（message passing interface）并行计算方法，对自适应扩大计算域进行并行计算，即可实现高精度大规模爆炸问题的数值模拟.obstacleTNTair在该大规模计算方法中，复杂结构内爆炸问题的数值模型一般包括两部分区域：计算域和障碍物区域，如图 1 所示.其中，在计算域中设置炸药区域、空气域两种介质以及刚性壁面区域，爆轰产物采用 JWL 状态方程计算，空气采用理想气体状态方程.障碍物区域不参与计算，但由于矿井、巷道等地下空间的壁面一般为钢筋混凝土墙，所以在障碍物与空气/爆轰产物的交界处给定固壁无滑移边界条件，用于模拟刚性壁面的反射过程.障碍物区域obstacle空气域air炸药区域TNT刚性壁面图 1 复杂结构数值模型示意图Fig.1 Diagram of numerical model of complex structure 2 分块续算方法由于地下空间或矿井等大规模复杂结构爆炸问题的数值研究都是百米、千米量级，而炸药爆炸形成的爆轰波阵面厚度通常在微米量级；另外，三维的爆炸问题存在时间和空间上的多尺度特性，若计算区域采用粗网格划分会带来较大的累积误差，导致计算结果不准确甚至不收敛；若采用精细网格则只能对小尺寸的简单结构进行模拟，大规模复杂密闭空间结构由于网格数目巨大，即使采用并行技术也很难满足对其全尺寸快速计算的需求，对其进行全尺寸网格划分更会浪费大量的 CPU 计算资源.以图 2所示的二维结构的爆炸问题为例，若对其进行全尺寸计算时需要 72 个 CPU 并行计算，且其中有 36 个处于障碍物区域不参与计算（图中阴影区域所示），极大地浪费了计算资源，如图 2（a）所示.为了减少大量计算资源的浪费，提高计算效率，实现对大规模复杂结构内爆炸问题的全流场数值模拟，提出了一种基于数据流的分块续算方法，通过合理划分模型并对每个部分进行单独计算可对 CPU 进行有效利用，如图 2（b）所示.对该二维结构模型利用分块续算方法计算时，共需要 44 个 CPU 但只有 6 个处于障碍物区域，与全尺寸计算方法相比，节约了 39%的 CPU核数，有效利用了计算域，减少了计算资源的浪费.另外，若全尺寸计算需 72 n 核时，则分块续算共需22 n 核时，计算核时相比于全尺寸计算节约了 69%.若以三维结构来看，全尺寸计算浪费的 CPU 计算资源更多，因此与全尺寸计算相比，数据流分块续算方法极大地减少了计算资源的浪费，对复杂结构进行合理划分即可有效实现其内爆流场的数值模拟.以下简要叙述数据流分块续算方法的流程：第 1 期王成等：复杂结构内爆炸冲击波传播的大规模计算方法研究41 划分模型.i=1,2,3,在对矿井、巷道等复杂地下空间大尺度爆炸问题进行模拟时，可根据实际情况将计算模型进行简化，并根据模型几何形状和爆炸流场的流体动力学特征将模型各部分进行合理划分，即将整个计算模型划分为若干部分（PART-i，），使得在对每个部分进行数值计算时尽可能少地浪费计算资源.数据储存.通常情况下，将放置炸药的模型部分划分为PART-1 对其先进行计算.首先根据具体的爆炸流场尺度以及计算资源，确定该部分的计算域尺寸、每个方向的 CPU 数以及每个分区中的网格数，并给定计算域各 CPU 的序号.另外，在 PART-1 模型的出口处设置一个数据储存区域，记为，如图 2（b）所示.在计算若干时间步后，该部分模型计算域中的爆炸冲击波向外传播，当冲击波传播至 PART-1 中 区域的 80%85%时，程序停止计算并将这一时刻 区域内各 CPU 所有计算网格点的物理量（包括密度、3个方向的速度、总能量、压力等）保存至文件中，即数据储存.数据赋值.在进行下一部分即 PART-2 的计算时，同样先确定该部分的计算域尺寸，并按 PART-1 的网格大小对其进行网格划分，使得相邻 PART 间能够进行“映射”数据流的传递.在开始计算前，先更新时间步，让两部分程序的时间同步.另外，在 PART-2 中设置一块与 PART-1 完全重叠的区域，在计算该部分模型时，将上一部分模型 PART-1 保存的区域 的数据作为该部分的初始条件，在开始计算时读入此文件并重新启动程序，将区域 内各 CPU 所有网格点的物理量直接映射在该部分区域 内相对应的计算网格点上，然后继续求解欧拉方程进行计算，即数据赋值.若进行计算的整个大规模模型被划分为很多部分，则根据流场中爆炸冲击波的传播过程对所有 PART依次进行数据储存、数据赋值的数据流传递，即可实现对整个模型的爆炸流场进行数值模拟.分块续算方法由于在相邻 PART 间分割处都会设置重叠区域，且各部分模型之间是通过数据储存、数据赋值的方式来进行数据流传递的，因此，相对于全尺寸计算，这种计算方法增加了一定的计算量.但是，相邻 PART 间的这种直接映射的数据流传递方式也保证了数据的完整性，区域 在一定程度上可以代表上一个 PART 计算得到的结果，使得接续计算时产生的误差可以忽略（下文对此进行了验证）.更为重要的是，分块续算方法解决了大尺度复杂结构由于网格数目过大而无法计算的难题，更解决了计算资源的大量浪费问题，实现了复杂地下密闭空间内爆炸冲击波传播的大规模数值模拟，在矿井、巷道等复杂地下空间的防爆抗爆设计中具有广阔的应用前景.3 验证性算例 3.1 计算方法验证为了验证该计算方法与其数值计算结果的准确性，结合文献 6 对方形坑道内爆炸冲击波传播的实验结果，设计了两端开口的 TNT 方形坑道内爆的计算模型，并利用该计算方法对其爆炸流场中的冲击波衰减过程进行数值计算.该计算模型示意图如图 3(a)全尺寸计算示意图数据储存数据赋值PART-1PART-2(b)分步续算示意图t时刻区域内所有网格点的物理量图 2 二维结构的数值计算示意图Fig.2 Diagram of numerical calculation of 2-D structure42北 京 理 工 大 学 学 报第 44 卷所示，方形坑道截面长宽为 1.2 m1.8 m，坑道总长为 41 m.圆柱形 TNT 装药的密度为 1.57 g/cm3，质量为 3.25 kg.柱状装药距离地面 0.9 m，且悬挂于距离坑道端面 5 m 的位置.在坑道一侧壁面的爆心所在截面上分别距离爆心 6、10、16、20、28、36 m 处设置超压测点，其布置如图 3 所示.6 m10 m16 m20 m5 m28 m36 m1.2 m1.8 m1.2 mTNT图 3 方形坑道数值计算测点布置示意图Fig.3 Diagram of measuring point layout for numerical calculation of square tunnel 表 1 列出了文献 5 所得实验值以及数值模拟所得典型测点处的超压峰值.从表中可以看出，20 m处冲击波超压的数值计算结果与实验结果基本一致，而越靠近爆炸中心的测点，超压峰值的偏差相对越大.因为在实验中距离炸药越近的地方冲击波压力越高，传感器受震动会产生较多杂信号，而且灵敏度（分辨率）不同，会影响读数.但总体而言，利用该数值计算方法预测的峰值压力与实验所测得的压力值基本一致，相对误差在 9.55%以内，证明该计算方法能够较为准确地预测结构内爆问题.表 1 典型测点处的文献实验及数值模拟超压峰值Tab.1 Literature experiments and numerical simulations ofoverpressure peaks at typical measurement points测点位置/m超压峰值/MPa误差/%数值计算文献5实验值60.6880.6289.55100.5110.4767.35160.2870.3004.33200.2340.2350.43280.1810.1904.74360.1610.1572.55 3.2 分块续算方法验证由于数据流分块续算方法在相邻 PART 间用于“数据储存、数据赋值”的 区域只保存了 PART-i中爆炸冲击波的主要流场数据，实际上 区域以外的爆炸流场对后续冲击波的传播也有一定的影响.为了研究 区域保存的流场数据能否作为下一部分 PART-（i+1）的初始条件，实现 PART 之间的冲击波续算过程，设计了如图 4 所示的标准模型算例对分块续算方法的准确性进行了验证.该计算模型由一个 4 m4 m3 m 的长方体和半径为 1 m、长 10 m 的管道组成，长方体内放置一块药量为 50 kg 的 TNT 正方体状装药，位于距离地面高0.5 m 处.整个计算域有两种介质：空气介质和 TNT炸药，空气介质采用理想气体状态方程，TNT 炸药采用 JWL 状态方程.另外，在计算域长管道中心轴线方向每隔 1 m 设置一个冲击波超压峰值测点，用于记录冲击波超压在管道内随距离的衰减过程.为验证数据流分块续算方法的准确性，对该计算模型内的爆炸问题进行两次数值计算即全尺寸计算和分块续算.全尺寸计算时，计算模型不进行拆分直接对整个计算域划分网格进行计算，如图 4 所示，网格尺寸 dx=dy=dz=0.01 m，网格总数约为 2.94 亿，采用 600 个 CPU 并行计算，总的计算时长为 19 h；分块续算时，将整个计算模型划分为 PART-1 和 PART-2两部分，并在 PART-1 和 PART-2 中分别设置一个长为 2 m 的区域，用于接续计算时数据流的传递，如图 4（b）所示.续算时，网格尺寸不变，网格总数约为1.17 亿，共采用 240 个 CPU 并行计算，总的计算时长为 8.5 h 左右.与全尺寸计算方法相比，分步续算方法节约了 60.00%的 CPU 核数，极大地减少了计算资源的浪费；另外，计算时长相比于全尺寸计算降低了55.26%左右，即分步续算方法在给定算例下相比于普通全尺寸计算，节约了 88.68%的计算核时，极大地提高了计算效率.TNT12345678910(a)全尺寸计算模型PART-212345678910PART-1TNT12(b)分块续算模型图 4 标准模型算例示意图Fig.4 Schematic diagram of standard model calculation example 表 2 列出了全尺寸计算和分块续算得到的管道内典型测点处冲击波的超压峰值.可见分块续算得到的超压峰值与全尺寸计算时的基本一致，相对误第 1 期王成等：复杂结构内爆炸冲击波传播的大规模计算方法研究43差在 0.33%以内.表明分块续算方法可实现基于数据流的数据传递，并能够保证数据在传递过程中的准确性.另外，对于复杂地下密闭空间等大规模结构，分块续算方法可以有效利用计算资源，节约计算时长，提高计算效率，实现爆炸流场中冲击波传播的大规模、高效率、高精度数值模拟.表 2 全尺寸计算和分块续算时典型测点处的超压峰值Tab.2 Peak overpressure at typical measuring points duringfull-scale calculation and block continuation测点编号超压峰值/MPa误差/%全尺寸计算分块续算117.56917.5690.00216.19316.1930.0039.2419.2390.0246.8616.8610.0057.1997.2020.0464.9924.9930.0275.6135.6130.0084.1004.1010.0294.7874.7990.25103.5993.5870.33 4 爆炸冲击波在复杂结构中的数值模拟 4.1 计算模型设置针对复杂地下密闭空间结构中的爆炸问题，采用上文提出的基于数据流的分块续算方法，设计了一个包含长方体型工房及圆柱形泄压洞的大规模复杂结构对爆炸流场中的冲击波传播及衰减过程进行数值模拟，其计算模型示意图如图 5 所示.整个计算模型由工房、泄压洞 A、B、C 及泄压洞B 和 C 出口的自由场区域组成，工房内放置一块质量为 392 kg 的方形炸药，TNT 装药密度为 1.58 g/cm3.其中，工房长宽高为 14 m5 m5 m，泄压洞 A 和 B的直径均为 2 m，泄压洞 C 直径为 1 m.泄压洞 A 长为 10.5 m，泄压洞 B 包括 7 m 的水平直段部分和15 m 的向上倾斜 45部分，泄压洞 C 长为 25 m，整个模型的壁面都设置为固壁无滑移边界条件.另外，在泄压洞 B 和 C 的出口处分别设置一块 4 m4 m2 m及 16 m8 m2 m 的自由场区域，该部分边界条件设置为无反射出流.由于该计算模型较为复杂，为了减少计算资源的浪费，根据模型具体几何形状特征将其划分为 6 个部分（PART-1、PART-2、PART-3.1、PART-3.2、PART-4.1、PART-4.2），并利用基于数据流的分块续算方法对其爆炸流场进行数值计算，如图 6 所示.在该计算模型中，用于数据储存、数据赋 PART-1 PART-2PART-4.1PART-4.2PART-3.1PART-3.2侧视图工房自由场自由场泄压洞A泄压洞B泄压洞C图 5 大规模复杂结构计算模型示意图Fig.5 Diagram of large-scale complex structure calculation model PART-1PART-2PART-3.1PART-3.2PART-4.1PART-4.2计算工房部分及1.6 m的工房出口部分,储存区域所有网格点的物理量计算泄压洞A和T型管道部分,储存和区域所有网格点的物理量计算泄压洞B的向上倾斜45部分,储存区域所有网格点的物理量计算泄压洞B尾部及其出口处的自由流场区域计算泄压洞C部分,储存区域所有网格点的物理量计算泄压洞C尾部及其出口处的自由流场区域图 6 数据流传递过程示意图Fig.6 Diagram of data flow transmission process44北 京 理 工 大 学 学 报第 44 卷值的 区域（即区域）设置为 1.6 m 长，整个模型各 PART 间数据流的传递过程如图 6 所示.4.2 计算结果及分析 4.2.1 工房内部及工房出口处的冲击波传播规律在工房内前后左右每个侧壁面的三等分点处分别设置 2 个测点，其高度与炸药中心平齐，监测每个测点处的超压时间曲线.通过数值计算可以发现，距离炸药最近的测点 2、4 处的反射超压高达281.1 MPa，主要是由于高压冲击波直接撞击至壁面处形成反射波（如图 7（b）所示），而在撞击壁面后冲击波在向外传播的同时压力快速衰减，因此测点 1、3 处的超压峰值仅有 29 MPa.(a)0.417 8 ms(c)1.030 ms(e)2.084 ms1211223344345566567778(b)0.721 8 ms(d)1.650 ms(f)2.677 ms11122233344455566677788888图 7 PART-1 典型时刻压力云图Fig.7 Pressure distributions at typical time of PART-1 在 1.030 ms 时，冲击波传播至泄压洞 A 内，由于工房与泄压洞 A 的交汇处存在 T 字型拐角，致使此处产生极强的稀疏效应，密度与压强均接近于 0，如图 7（c）所示；冲击波在 1.333 6 ms 时传播至泄压洞 A右侧壁面，并在 1.650 ms 时刻从右侧反射至管道中心轴线区域，因此使得泄压洞 A 始段约 0.8、1.0 m 处的压力并不高，反而在泄压洞 A 约 1.21.6 m 处的压力升高.在 2.084、2.677 ms 时冲击波传播至工房右、左两侧壁面，并产生较高压力的反射波，分别高达 54、45 MPa.4.2.2 泄压洞 A 内部冲击波衰减规律图 8 给出了泄压洞 A 内部冲击波超压峰值随距离的衰减规律.在泄压洞 A 入口 0 m 处，管道内中心轴线的超压峰值为 9.78 MPa，随后爆炸冲击波以波动衰减的趋势逐渐减小，经过约 10 m 长的泄压洞A 后，冲击波波动衰减至约 1.2 MPa 左右，与入口处相比降低了约 88%.因为冲击波从工房传播出来时并不是平面波，从而导致冲击波在泄压洞内的壁面上来回反射，因此在泄压洞 A 内部的冲击波并不像空中自由场爆炸一样随着爆距的增加而快速衰减，而是以波动衰减的趋势逐渐降低，其波动衰减过程如图 9（a）9（d）所示.01234567891002468101214p/MPaR/m图 8 泄压洞 A 内超压峰值衰减过程Fig.8 Attenuation process of overpressure peak in pressurerelief tunnel A第 1 期王成等：复杂结构内爆炸冲击波传播的大规模计算方法研究45(a)2.189 4 ms(b)3.069 4 ms(c)4.263 0 ms(d)5.717 4 ms(e)7.424 0 ms(f)9.268 4 ms(g)10.270 8 ms(h)11.292 6 ms图 9 PART-2 典型时刻冲击波压力云图Fig.9 Pressure distributions at typical time of PART-2 4.2.3 泄压洞 A、B、C 交汇处 T 型通道冲击波传播规律由图 9 可以看出，冲击波在约 7.424 0 ms 时刻传播至泄压洞 A 侧 R15 圆弧处，在 9.268 4 ms 时撞击至泄压洞 B 侧 R15 圆弧处，并产生较高反射压力，此后，该反射冲击波向泄压洞 A、B 和泄压洞 C 内传播，使得 T 型区域的超压峰值有一定程度的升高.由于冲击波在 T 型区域处进行无规则反射，致使该区域超压较不稳定，为研究冲击波在 T 型通道处的分流情况，此处取距离三管道中心轴线交点2 m 处的测点进行冲击波分流分析.通过读取相关测点处的超压峰值可知，泄压洞 A 内距离交点 2 m处的超压峰值为 1.39 MPa，泄压洞 B 的超压峰值为0.93 MPa，泄压洞 C 的超压峰值为 0.89 MPa.即在 T型通道分流处，与泄压洞 A 同等直径的泄压洞 B 的冲击波超压峰值约占分流前的 67%，而呈直角且直径只有泄压洞 A1/2 的泄压洞 C 的超压峰值约占分流前的 64%.由此可见，分流时管道间采用直角交汇和缩小管道直径在一定程度上有利于冲击波超压的衰减.4.2.4 泄压洞 B 和 C 内部冲击波衰减规律从 T 型区域开始，随着冲击波向外传播，泄压洞C 内的冲击波超压值出现了明显减小，其衰减过程如图 10（b）所示.冲击波超压峰值从泄压洞 C 内 0 m的 1.22 MPa 衰减至 24 m 处的 0.19 MPa，降低了 84%，且其衰减速度逐渐减慢.由于泄压洞 C 前端部分位于 T 型区域内，因此冲击波在 03 m 范围内有一定的反射，超压峰值有所升高，在 2.6 m 处达到最高值1.8 MPa.距离 16 m 处开始波动衰减消失，超压峰值以缓慢的速度逐渐降低，说明在泄压洞 C 内 16 m 左右,冲击波经过多次壁面反射后已基本形成平面波.冲击波于 T 型通道分流后，在泄压洞 B 内经过 7 m的水平直段及 15 m 的向上倾斜 45部分的传播，其在斜段部分的衰减过程如图 10（a）所示.爆炸冲击波46北 京 理 工 大 学 学 报第 44 卷的超压峰值从斜坡起始处的 0.82 MPa 衰减至 15 m处的 0.20 MPa，降低了 76%，而且其衰减速度逐渐减慢.同样，冲击波在斜坡管道内传播时，其超压峰值也是一个波动衰减的过程.024681012141600.250.500.751.001.25(a)泄压洞B斜段部分p/MPaR/m0481216202400.51.01.52.0(b)泄压洞Cp/MPaR/m图 10 泄压洞 B 和 C 内部超压峰值衰减规律Fig.10 Attenuation law of overpressure peak inside relief tunnels B and C 4.2.5 泄压洞 B 和 C 出口自由场区域冲击波衰减规律泄压洞 B 和 C 外部自由场区域冲击波超压峰值测点的分布示意图如图 11（a）和 11（b）所示，其所对应的超压峰值分别如图 11（c）和 11（d）所示，可以看出，与管道内冲击波衰减规律不同，冲击波在自由场区域是持续衰减的，且其衰减速度比较快.爆炸冲击波从泄压内传出后，分别于 1 m 和 5 m 处“整合”为平面波.2 m测点19测点1020(左斜)测点2131(右斜)(b)泄压洞C出口处测点布置示意图测点120测点2141(左斜45)测点4262(右斜45)2 m(c)泄压洞B出口超压峰值衰减过程p/MPap/MPaR/m01234500.050.100.150.20轴线方向右斜方向左斜方向(d)泄压洞C出口自由场超压峰值衰减过程R/m02468101200.050.100.150.20轴线方向右斜45方向左斜45方向(a)泄压洞B出口处测点布置示意图图 11 泄压洞出口测点分布及超压衰减过程Fig.11 Distribution of measuring points and the process of overpressure attenuation 第 1 期王成等：复杂结构内爆炸冲击波传播的大规模计算方法研究474.3 区域取值对数值计算结果的影响上述采用基于数据流的分块续算方法对爆炸冲击波在复杂结构内的传播过程进行了数值计算，为了在一定程度上能够减少计算量，用于数据储存、数据赋予的 区域不宜取得过大.上述计算模型中用于数据传递的 区域取值为 1.6 m，在此通过对 区域进行适当的扩大，研究 区域的取值对分块续算结果的影响大小，进而验证数据流分块续算方法的准确性.选用 PART-3.1 和 PART-3.2 两部分模型，对 区域的取值对计算结果的影响进行研究.用于数据传递的 区域扩大为原来的 3 倍，即在计算 PART-3.1 时，保存所有物理量的区域由原来的 1.6 m 扩大为 4.8 m，在 PART-3.2 中也设置一块与 PART-3.1中区域重叠的区域.图 12 给出了 区域为1.6 m 和 4.8 m 时，泄压洞 C 出口区域超压峰值衰减曲线比较图.由图中对比结果可以看出，冲击波在轴线方向、左斜及右斜方向的超压峰值衰减过程基本一致，因此，对上述大规模复杂结构进行分块续算时，进行数据传递的 区域取为 1.6 m 是合理的，对整体爆炸流场的数值计算结果影响极小.02468101200.050.100.150.2000.050.100.150.20(a)轴线方向024681012 区域为1.6 m 区域为4.8 m 区域为1.6 m 区域为4.8 m(b)左斜方向R/mR/mp/MPa 区域为1.6 m 区域为4.8 m02468101200.050.100.150.20R/m(c)右斜方向p/MPap/MPa图 12 区域为 1.6 m 和 4.8 m 时泄压洞 C 出口区域超压峰值衰减过程比较Fig.12 Comparison of peak overpress