CN113158404A - 一种用于计算网状坑道内爆炸冲击波传播的方法 - Google Patents

Abstract

本发明涉及一种计算网状坑道内爆炸冲击波传播的方法，其包括：建立坑道网络图和初始节点的初始化、根据预设的三个指标建立最不利传播路径、计算最不利路径的衰减系数、计算节点冲击波超压以及计算工程离散点冲击波超压，本发明提供了基于最小累计距离、最小累计节点度、最小累计排除空气体积三个指标计算最不利传播路径的方法，在计算衰减系数时考虑爆炸初始段、口部传入、直线传播、拐点、分叉、断面典型工况条件下爆炸冲击波衰减影响，提高了最终计算的精度和可靠度，同时保证了较低的算法复杂度。

Description

一种用于计算网状坑道内爆炸冲击波传播的方法

技术领域

本发明领域为防护工程领域，具体为一种用于计算网状坑道内爆炸冲击波传播的方法。

背景技术

坑道结构是地铁、矿井、综合管廊、地下隧道、人防工程、军事防护工程等地下工程设施中较为常见的一种结构形式，一直以来，爆炸冲击波在坑道中的传播规律是地铁反恐防爆、矿井安全管理、人防建设、军事防护等关注的核心问题之一，与无约束爆炸冲击波比，坑道内爆炸冲击波由于受到坑道四周壁面阻挠以及约束作用而不断沿着坑道传播，爆炸冲击波的超压峰值、作用时间等效应参数的衰减规律相较自由场中传播要更为复杂，研究表明，不同的坑道构型对爆炸冲击波传播衰减规律影响存在较大差异，其中复杂坑道构型，尤其是交叉成环、成网坑道中的爆炸冲击波传播由于涉及多组冲击波相互间的扰流、对流或汇聚，其传播衰减规律比在简单坑道构型中也更为复杂；

目前，针对坑道结构爆炸冲击波传播规律影响的研究主要集中在简单坑道构型上，例如直通式出入口、单项式出入口、穿廊式出入口、直角拐弯构型、直角分叉构型等，反映爆炸冲击波传播规律的核心参数主要是冲击波超压峰值；采用的计算方法主要是小比例模型坑道试验、原型坑道试验、数值模拟计算等，并已在实际工程设计和论证中被广泛使用。但是，针对简单坑道构型的爆炸冲击波传播计算方法无法很好支持对复杂网状坑道的计算，原因如下：

1、小比例模型坑道试验或原型坑道试验方法的建设周期长、投入成本高、通用性弱，一般仅用于部分坑道段或关键技术验证，不可能对完整的网状坑道开展试验；数值模拟计算方法的建模工作量大、仿真时间长，不便于开展快速计算；

2、传统爆炸冲击波入射超压峰值计算方法一般仅考虑爆炸点在坑道口部外，很少涉及在坑道内部爆炸的情况，现有技术中也缺乏针对坑道内部不同位置处爆炸情况的系统研究；

3、某些已有算法在计算坑道内某点的爆炸冲击波入射超压峰值时，将爆炸冲击波传播的最不利路径等价为爆炸点至坑道内该点的最短距离路径进行求解，实际是仅考虑了爆炸冲击波沿直线传播的距离衰减因素，忽略了爆炸冲击波在坑道拐点、分叉、界面突变等构型处的衰减作用，这些位置的衰减作用一般比距离衰减因素影响更大，因此得出的结果与理论值会存在较大偏差。

发明内容

针对现有技术中存在的不足，本发明的目的在于提出一种用于计算网状坑道内爆炸冲击波传播的方法，基于最小累计距离、最小累计节点度、最小累计排除空气体积三个指标，提出多指标综合最不利传播路径算法，同时通过对爆炸初始段、口部传入、直线传播、拐点、分叉、断面等典型工况条件下爆炸冲击波衰减影响分情况考虑计算综合衰减系数考虑了坑道连接处衰减、坑道截面变化等因素对冲击波传播的衰减作用，同时，考虑了不同衰减因素的影响，进一步提高了坑道内冲击波大小的计算精度，也能够保证较低的算法复杂度，实现较高的计算效率。

一种用于计算网状坑道内爆炸冲击波传播的方法，其具体包括如下步骤：

步骤一：构建网状坑道的坑道网络图；

步骤二：根据爆炸点位置对所述坑道网络图更新，进行初始节点参数初始化，确定爆炸冲击波入射超压峰值初始计算节点；

步骤三：确定所述坑道网络图中的任意一个未计算爆炸冲击波入射超压峰值的坑道网络内部节点作为下一个计算节点，并转入步骤四；

步骤四:预设表征衰减因素的全局静态指标，根据所述全局静态指标，均按照网络最短路方法计算由所述初始计算节点到达所选定内部节点的冲击波最不利传播路径；

所述全局静态指标包括指标一、指标二以及指标三，所述指标一为传播累计距离，所述指标二为传播累计节点度，所述指标三为传播累计体积；

步骤五：根据所述步骤三中得出的三条最不利传播路径以及步骤二中更新的坑道网络图，计算三条最不利传播路径各自对应的综合衰减系数，并确定最大综合衰减系数和全局最不利传播路径，根据所述步骤二中的初始计算节点爆炸冲击波入射超压峰值和所述全局最不利传播路径，计算坑道网络节点的爆炸冲击波全局最不利入射超压峰值；

步骤六:判定是否还有内部节点未计算爆炸冲击波入射超压峰值，若有节点未计算爆炸冲击波入射超压峰值，则返回步骤三，重新确定计算节点；若没有内部节点未计算爆炸冲击波入射超压峰值，则进行步骤七；

步骤七：计算坑道段(v_i,v_j)中任意位置处离散点v_x的冲击波超压峰值P_x：

其中，

为坑道段(v_i,v_j)在v_i点的全局最大入射冲击波超压峰值；

为坑道段(v_i,v_j)在v_j点的全局最大入射冲击波超压峰值；|v_iv_j|表示坑道段(v_i,v_j) 长度。

进一步地，所述步骤二中，根据爆炸点位置对所述坑道网络图进行更新，确定初始计算节点时，其过程包括：

Step1：确定爆炸点是否为口部外点；

若爆炸点为口部外点，则原坑道网络图不变，初始计算节点确定为距离爆炸点最近的坑道口部节点，计算该初始计算节点的冲击波入射超压峰值ΔP₀，

其中：ΔP₀为入口处的冲击波入射超压峰值；W为爆炸物TNT当量；X₀为爆炸点到坑道口部无遮挡直线距离；

若爆炸点不是口部节点，则转入Step2；

Step2：确定爆炸点是否位于坑道内网络节点上；

若爆炸点在坑道内网络节点上，则原坑道网络图不变，将初始计算节点确定为该爆炸点所在网络节点，并依次计算到达与所述网络节点相邻的其他各节点冲击波入射超压峰值；

(a)当爆炸点正好位于坑道段拐点处时，则按照以下公式分别计算与爆炸点邻接的两个坑道网络节点的入射冲击波超压峰值ΔP₀，

其中：δ_X为爆炸点的初始爆距，所述初始爆距为冲击波传播方向上最接近爆炸点的坑道网络节点距离；S为坑道截面积，W为爆炸物TNT 当量；

(b)当爆炸点正好位于坑道段三岔口或十字交叉口时，则按照以下公式计算与爆炸点邻接的三个或四个坑道网络节点的入射冲击波超压峰值ΔP₀，

其中，当所述爆炸点位于三岔口时，S′＝3·S/2，位于十字叉口时， S′＝4·S/2，S为坑道截面积，δ_X为爆炸点与相邻各节点的初始爆距；

若爆炸点不在坑道内网络节点上，则进行Step3；

Step3、确定爆炸点是否在坑道内部某坑道段上；

若爆炸点在坑道内部某坑道段上，则需要在坑道网络图中新增节点为初始计算点，其中；

设新增节点v_n+1所在坑道段为(v_a,v_b)，更新原坑道网络图G(V,E)为 G′(V′,E′)，将所述坑道段中心轴线交叉点及出入口端点构成的点集 V＝{v_i|i＝1,2,…,n}更新为V′＝V∪{v_n+1}，坑道段边集E＝{(v_i,v_j)|v_i,v_j∈V,i≠j}更新为E′＝E∪{(v_a，v_n+1)}∪{(v_n+1，v_b)}-{(v_a，v_b)}；初始计算节点确定为新增节点v_n+1且其与爆炸点重合，计算与爆炸点邻接的两个坑道网络节点的入射冲击波超压峰值ΔP₀，

其中：δ_X为爆炸点的初始爆距，所述初始爆距为冲击波传播方向上最接近爆炸点的坑道网络节点距离；S为坑道截面积，W为爆炸物TNT当量。

进一步地，所述指标一用下式表示：

其中：L(·)为坑道段长度函数，L(R_i，R_i+1)表示冲击波传播路径R中相邻的第 i个和第i+1个节点所在坑道段的坑道长度；R为由两两相邻的坑道网络节点所构成的连通的冲击波传播路径，为一组坑道网络节点序列；R_i，R_i+1为路径R中沿传播方向的第i和第i+1个节点；|R|为路径中节点个数；R1为I₁(R)最小值对应的冲击波最不利传播路径，以指标一为准求解冲击波最不利传播路径，将坑道段长度表示坑道网络的边权，然后采用网络最短路方法求解最不利传播路径；

所述指标二用以下公式表示：

其中：Deg(·)为坑道网络节点度函数，Deg(R_i)表示与路径R中第i个节点相连的坑道段的条数，可用网络邻接矩阵中该点所在行的非0项总个数来计算；R_i为路径R中沿传播方向的第i个节点；R2为I₂(R)最小值对应的冲击波最不利传播路径；

以指标二为准求解冲击波最不利传播路径，先计算坑道网络图中所有节点的度数，再用坑道段两端节点度的平均值表示坑道网络的边权，从而转化为网络最短路方法求解最不利传播路径；

所述指标三用以下公式表示：

其中，L(·)为坑道段长度函数，L(R_i，R_i+1)表示路径R中相邻的第i个和第i+1 个节点所在坑道段的坑道长度；S_i，i+1为传播路径中相邻的第i个和第i+1个节点所在坑道段的通道截面积，R3为经优化得到的累计体积最小的路径；

先计算坑道网络所有坑道段长度和通道截面积，再用求得各坑道段的体积表示坑道网络的边权，从而转化为网络最短路方法求解最不利传播路径。

进一步地，所述综合衰减系数是冲击波从初始计算节点开始，按照某一冲击波传播路径R，经过坑道网络内部衰减到达计算节点后的超压峰值相对初始节点超压峰值的比例。

进一步地，计算所述综合衰减系数时，根据爆炸点所处的坑道位置，分别采用相应的综合衰减系数计算方法，其过程包括：

若爆炸点为坑道口部外点时，所述综合衰减系数为，

其中：n_i表示第i种传播衰减在某次冲击波传播路径R中重复出现的次数， i＝1，2，3，4，5，n_i为整数；ξ₁为口部传入衰减系数；ξ₂为直线传播衰减系数；ξ₃为拐点传播衰减系数；ξ₄为分叉传播衰减系数；ξ₅为断面变化衰减系数；

若爆炸点位于坑道内部时，则不存在口部传入衰减情况，所述综合衰减系数为，

其中：n_i表示第i种传播衰减在某次冲击波传播路径中重复出现的次数， i＝2,3,4,5,n_i为整数；ξ₂为直线传播衰减系数；ξ₃为拐点传播衰减系数；ξ₄为分叉传播衰减系数；ξ₅为断面变化衰减系数；

进一步地，所述步骤五中，最大综合衰减系数指以初始计算节点为传播路径起点且到达坑道内部节点为传播路径终点的冲击波传播路径R1，R2， R3所对应的综合衰减系数的最大值，其中R1为I₁(R)最小值对应的冲击波最不利传播路径，R2为I₂(R)最小值对应的冲击波最不利传播路径；R3为经优化得到的累计体积最小的路径。

进一步地，所述综合衰减系数的最大值为

其中，R1为I₁(R)最小值对应的冲击波最不利传播路径，R2为I₂(R)最小值对应的冲击波最不利传播路径；R3为经优化得到的累计体积最小的路径；

坑道网络内部节点v_q的爆炸冲击波全局最不利入射超压峰值为

ΔP_p为起始点v_p的入射超压峰值。

与现有技术相比，本发明的技术效果在于：

本发明的基于多指标综合最不利传播路径的网状坑道内冲击波传播算法，本发明基于最小累计距离、最小累计节点度、最小累计排除空气体积三个指标，提出多指标综合最不利传播路径算法，同时通过对爆炸初始段、口部传入、直线传播、拐点、分叉、断面等典型工况条件下爆炸冲击波衰减影响分情况考虑计算综合衰减系数考虑了坑道连接处衰减、坑道截面变化等因素对冲击波传播的衰减作用，同时，考虑了不同衰减因素的影响，进一步提高了坑道内冲击波大小的计算精度，也能够保证较低的算法复杂度，实现较高的计算效率。

尤其，本发明在计算综合衰减系数时，分情况考虑，考虑了爆炸初始段、口部传入、直线传播、拐点、分叉、断面等典型工况条件下爆炸冲击波衰减影响因素，对每种情况提出了对应的衰减系数计算方式，进一步提高了最终计算结果的准确性与可靠性。

尤其，本发明设置最小累计距离、最小累计节点度、最小累计排除空气体积三个指标，根据所述指标生成最不利路径，并且本发明还提出了用于计算网络最短路、任意坑道拐点夹角、任意坑道分叉夹角和坑道中间离散点插值的关键方法，计算更为精确、可靠，进一步提高了最终计算结果的精度。

附图说明

图1为本发明实施例的用于计算网状坑道内爆炸冲击波传播的工程化算法的流程图；

图2为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的网络图更新步骤图；

图3为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的网络最短路方法流程图；

图4为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的冲击波口部传入情况及对应衰减系数示意图；

图5为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的冲击波在坑道段拐弯处的情况及对应衰减系数示意图；

图6为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的冲击波在坑道段三分叉工况一传播的情况及对应衰减系数示意图；

图7为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的冲击波在坑道段三分叉工况二传播的情况及对应衰减系数示意图；

图8为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的冲击波在坑道段四分叉工况传播的情况及对应衰减系数示意图；

图9为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的冲击波从小截面坑道传入大截面坑道的情况及对应衰减系数示意图；

图10为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的应用实例坑道路线图；

图11为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的应用实例坑道截面图；

图12为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的应用实例坑道距离权网络图；

图13为根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法所提供的应用实例坑道内冲击波超压分布图。

具体实施方式

以下结合附图，对本发明上述的和另外的技术特征和优点作更详细的说明。

下面参照附图来描述本发明的优选实施方式。本领域技术人员应当理解的是，这些实施方式仅仅用于解释本发明的技术原理，并非旨在限制本发明的保护范围。

需要说明的是，在本发明的描述中，术语“上”、“下”、“左”、“右”、“内”、“外”等指示的方向或位置关系的术语是基于附图所示的方向或位置关系，这仅仅是为了便于描述，而不是指示或暗示所述装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

此外，还需要说明的是，在本发明的描述中，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域技术人员而言，可根据具体情况理解上述术语在本发明中的具体含义。

参阅图1所示，其为本发明实施例所提供的一种用于计算网状坑道内爆炸冲击波传播的工程化算法流程图。具体地，如图1至图13所示，根据本发明实施例的用于计算网状坑道内爆炸冲击波传播的方法，其具体包括如下步骤：

步骤一：构建网状坑道的坑道网络图；也就是说，对所述坑道网络图进行初始化；

具体地，根据网状坑道的空间地理位置和结构形式构建所述网状坑道的坑道网络图，基于网络图理论将所述网状坑道的坑道段抽象为所述坑道网络图的边，将所述网状坑道的出入口端点以及相邻坑道段所在中心轴线的交叉点抽象为所述坑道网络图的节点，根据所述网状坑道的坑道网络节点连接关系建立网络节点邻接矩阵，其中：

设G(V,E)是一个完全无向图，其用以代表坑道网络图。

V＝{v_i|i＝1,2,…,n}表示坑道段中心轴线交叉点及出入口端点构成的点集，

其中，v_i表示第i个坑道网络节点，n代表坑道网络节点数量。

E＝{(v_i,v_j)|v_i,v_j∈V,i≠j}代表坑道段边集，v_i表示第i个坑道网络节点，v_j表示第i个坑道网络节点。

(v_i,v_j)表示坑道网络图内以v_i及v_j为两端节点的坑道段，根据坑道网络的连接情况和坑道网络节点的邻接关系得到该网状坑道的邻接矩阵，用 A＝[a_ij]_n×n表示。

更进一步地，若v_i与v_j为某一坑道段的两端节点，则a_ij＝1；否则，a_ij＝0。

根据邻接矩阵得到坑道网络的边权W,W＝{w_ij|(v_i,v_j)∈E,w_ij＞0,i≠j}。

其中，w_ij是坑道段属性的参数，其包括坑道段(v_i,v_j)的长度、坑道段两端节点度的平均值和坑道段的体积.

当冲击波在坑道网络传播时，所形成的冲击波的传播路径R，R＝[v_p,…,v_q]_1×k。

其中，R为冲击波传播路径，v_p为爆炸起始节点且p为该节点在坑道的网络图中的序号，v_q为传播末节点且q为该节点在坑道的网络图中的序号， k为冲击波传播路径R中节点数量。

优选地，所述冲击波传播路径R必须满足如下条件：

所有节点属于所述坑道段中心轴线交叉点及出入口端点构成的点集V，以表示冲击波只在坑道网络中传播；

任意两个连续节点构成的边属于所述坑道段边集E且对应网状坑道网络的邻接矩阵A中的元素值为1，表示传播路径中的节点在坑道网络图中必须相互连接；

任意一个节点只出现1次，表示传播路径中没有内环流。

步骤二：根据爆炸点位置对所述坑道网络图更新，进行初始节点参数初始化，确定初始计算节点并计算该初始计算节点的入射冲击波超压峰值；

具体而言，请参阅图2所示，其为本发明实施例所提供的网络图更新步骤图，所述步骤二中根据爆炸点位置对所述坑道网络图更新，确定初始计算节点并计算坑道网络图基础参数的过程包括：

Step1、确定爆炸点是否为口部外点，若爆炸点为口部外点，则原坑道网络图不变，初始计算节点确定为距离爆炸点最近的坑道口部节点，并且按照以下公式计算该初始计算节点的冲击波入射超压峰值ΔP₀，若爆炸点不是口部节点，则进行Step2；

其中：ΔP₀为入口处的冲击波入射超压峰值；W为爆炸物TNT当量；X₀为爆炸点到坑道口部无遮挡直线距离，

Step2、确定爆炸点是否位于坑道内网络节点上，若爆炸点在坑道内网络节点上，则原坑道网络图不变，将初始计算节点确定为该爆炸点所在网络节点，并依次计算到达与所述网络节点相邻的其他各节点冲击波入射超压峰值，若爆炸点不在坑道内，则进行Step3；特殊地，

当爆炸点正好位于坑道段拐点处时，则按照以下公式分别计算与爆炸点邻接的两个坑道网络节点的入射冲击波超压峰值ΔP₀，

其中：δ_X为爆炸点的初始爆距，所述初始爆距为冲击波传播方向上最接近爆炸点的坑道网络节点距离；S为坑道截面积，W为爆炸物TNT 当量；

当爆炸点正好位于坑道段三岔口或十字交叉口时，则按照以下公式计算与爆炸点邻接的三个或四个节点的入射冲击波超压峰值ΔP₀，

其中，当所述爆炸点位于三岔口时，S′＝3·S/2，位于十字叉口时， S′＝4·S/2，S为坑道截面积，δ_X为爆炸点与相邻各节点的初始爆距。

Step3、确定爆炸点是否在坑道内部某坑道段上，若爆炸点在坑道内部某坑道段上，则需要在坑道网络图中新增节点为初始计算点。假设新增节点v_n+1所在坑道段为(v_a,v_b)，则更新原坑道网络图G(V,E)为G′(V′,E′)，其中，将所述坑道段中心轴线交叉点及出入口端点构成的点集V＝{v_i|i＝1,2,…,n}更新为V′＝V∪{v_n+1}，坑道段边集E＝{(v_i,v_j)|v_i,v_j∈V,i≠j}更新为 E′＝E∪{(v_a,v_n+1)}∪{(v_n+1,v_b)}-{(v_a,v_b)}。初始计算节点确定为新增的节点v_n+1(与爆炸点重合)，分别按照以下公式计算与爆炸点邻接的两个坑道网络节点的入射冲击波超压峰值ΔP₀，

其中：δ_X为爆炸点的初始爆距，所述初始爆距为冲击波传播方向上最接近爆炸点的坑道网络节点距离；S为坑道截面积，W为爆炸物TNT当量。

步骤三：确定所述坑道网络图中的任意一个未计算超压冲击波峰值的坑道网络节点为下一个计算节点，并进行步骤四；

步骤四：设计三组用以表征衰减因素的全局静态指标，根据所述三组全局静态指标，均按照网络最短路方法计算由所述初始计算节点到达所选定内部节点的冲击波最不利传播路径；所述三组全局静态指标包括：指标一、指标二以及指标三，所述指标一为传播累计距离，所述指标二为传播累计节点度，所述指标三为传播累计体积，

所述指标一用下式表示：

其中：L(·)为坑道段长度函数，L(R_i，R_i+1)表示冲击波传播路径R中相邻的第 i个和第i+1个节点所在坑道段的坑道长度；R为由两两相邻的坑道网络节点所构成的连通的冲击波传播路径，为一组坑道网络节点序列；R_i，R_i+1为路径R中沿传播方向的第i和第i+1个节点；|R|为路径中节点个数；R1为I₁(R)最小值对应的冲击波最不利传播路径，以指标一为准求解冲击波最不利传播路径，将坑道段长度表示坑道网络的边权，然后采用网络最短路方法求解最不利传播路径；

所述指标二用以下公式表示：

其中：Deg(·)为坑道网络节点度函数，Deg(R_i)表示与路径R中第i个节点相连的坑道段的条数，可用网络邻接矩阵中该点所在行的非0项总个数来计算；R_i为路径R中沿传播方向的第i个节点；R2为I₂(R)最小值对应的冲击波最不利传播路径：

以指标2为准求解冲击波最不利传播路径，先计算坑道网络图中所有节点的度数，再用坑道段两端节点度的平均值表示坑道网络的边权，从而转化为网络最短路方法求解最不利传播路径；

所述指标三用以下公式表示：

其中，L(·)为坑道段长度函数，L(R_i，R_i+1)表示路径R中相邻的第i个和第i+1 个节点所在坑道段的坑道长度；S_i，i+1为传播路径中相邻的第i个和第i+1个节点所在坑道段的通道截面积，R3为经优化得到的累计体积最小的路径；

先计算坑道网络所有坑道段长度和通道截面积，再用求得各坑道段的体积表示坑道网络的边权，从而转化为网络最短路方法求解最不利传播路径。

具体而言，请参阅图3所示，其为本发明实施例所提供的网络最短路方法流程图，本领域技术人员应当明白的是，网络最短路方法已经被提出，本实施例采用网络最短路方法中的一种，其具体过程参阅图3所示，本实施例提出的网络最短路方法中的“距离”除表示实际长度距离外，还表示坑道网络的边权值，其包括，在所述指标一中表示坑道段距离所表示的边权、在所述指标二中坑道网络节点度平均值所表示的边权以及所述指标三中坑道段体积所表示的边权。

步骤五：根据所述步骤三中得出的三条最不利传播路径以及步骤二中更新的坑道网络图，计算三条路径对应的综合衰减系数，并确定最大综合衰减系数和全局最不利传播路径，根据所述步骤二中的初始计算节点爆炸冲击波入射超压峰值和所述全局最不利传播路径，计算坑道网络节点的爆炸冲击波全局最不利入射超压峰值。

所述综合衰减系数是指冲击波从初始计算节点开始，按照某一冲击波传播路径R，经过坑道网络内部多种因素衰减到达计算节点后的超压峰值相对初始节点超压峰值的比例。特别地，对于爆炸点相对坑道所处位置的不同，采用不同的综合衰减系数计算方式，其包括，

若爆炸点为坑道口部外点时，综合衰减系数按照以下公式计算，

其中：n_i表示第i种传播衰减在某次冲击波传播路径R中重复出现的次数， i＝1,2,3,4,5,n_i为整数；ξ₁为口部传入衰减系数；ξ₂为直线传播衰减系数；ξ₃为拐点传播衰减系数；ξ₄为分叉传播衰减系数；ξ₅为断面变化衰减系数。

若爆炸点位于坑道内部时，则不存在口部传入衰减情况，综合衰减系数按照以下公式进行计算，

其中：n_i表示第i种传播衰减在某次冲击波传播路径中重复出现的次数， i＝2,3,4,5,n_i为整数；ξ₂为直线传播衰减系数；ξ₃为拐点传播衰减系数；ξ₄为分叉传播衰减系数；ξ₅为断面变化衰减系数。

所述步骤五中的最大综合衰减系数指以初始计算节点为传播路径起点到达坑道内部节点为传播路径终点的所有可能的冲击波传播路径R1，R2，R3，所对应的综合衰减系数的最大值，R1为I₁(R)最小值对应的冲击波最不利传播路径，R2为I₂(R)最小值对应的冲击波最不利传播路径；R3为经优化得到的累计体积最小的路径；

所述综合衰减系数的最大值为

其中，R1为 I₁(R)最小值对应的冲击波最不利传播路径，R2为I₂(R)最小值对应的冲击波最不利传播路径；R3为经优化得到的累计体积最小的路径；

坑道网络内部节点v_q的爆炸冲击波全局最不利入射超压峰值为

ΔP_p为起始点v_p的入射超压峰值。

具体而言，本发明实施例所提供的冲击波口部传入情况及对应衰减系数示意图，口部传入衰减是指爆炸冲击波到达坑道口部时，受口部壁面阻挡、坡面朝向以及坑道入口方向影响发生突然性衰减。

口部传入衰减系数ξ₁按照以下公式进行计算：

其中，α₁为坑道轴线方向(指向内侧)与冲击波阵面传播方向形成的口部夹角。

确定所述口部夹角α₁时，根据爆炸点v₀坐标、最接近爆炸点的口部点v_k坐标，以及与口部点v_k直接相邻的节点v_k+1坐标，按照下式计算：

直线传播衰减是指爆炸冲击波在坑道中沿直线传播时，受空气阻力、重力、壁面反射与壁面凸起扰流等因素综合作用影响发生连续性衰减，直线传播衰减系数ξ₂按照以下公式计算：

其中，C_f为坑道壁面摩擦系数；D为坑道等效直径，

S为坑道截面积；X为传播路径中某一计算节点与坑道段另一端节点间的距离，且满足 X≥6D；δ₀为坑道壁面凸起高度，可取平均值。

参阅图5所示，其为本发明实施例所提供的冲击波在坑道段拐弯处的情况及对应衰减系数示意图，拐点传播衰减是指爆炸冲击波到达坑道段两两连接的拐弯处时，受壁面反射扰动而发生突然性衰减，传播方向角度变化越大，衰减幅度越大。拐点传播衰减系数ξ₃按照以下公式计算：

其中，α₃为相邻两条坑道轴线相交所形成的夹角。

确定所述相邻两条坑道轴线相交所形成的夹角α₃时，首先基于邻接矩阵判断坑道网络图节点是否为拐点，即按照节点度等于2(邻接矩阵中该点所在行的元素之和，等于2表明只有两个坑道网络图节点与其相邻)来筛选；然后根据传播路径R中拐点v_k坐标、与其相邻的2个节点v_k-1和v_k+1坐标，按照下式计算：

参阅图6-图7所示，其为本发明实施例所提供的冲击波在坑道段三分叉工况一传播的情况及对应衰减系数示意图以及冲击波在坑道段三分叉工况二传播的情况及对应衰减系数示意图，分叉传播衰减是指爆炸冲击波到达由三段坑道或四段坑道同时相交形成的分叉点时，受坑道分流、反射扰动而发生突然性衰减，分叉越多、分叉后传播方向变化越大，衰减幅度越大，所述分叉传播衰减系数ξ₄需要根据坑道出现三分叉和四分叉情况时的不同工况展开计算，且假设分叉处坑道轴线所在直线有且只有2条，即为“T”字型或“十”字型，不存在“丫”字型或“K”字型分叉，

当坑道出现三分叉时根据工况一以及工况二采用不同的公式计算所述综合衰减系数，所述工况一为传出坑道方向全部改变，所述工况二为传出坑道方向部分改变；

计算所述综合衰减系数时，

继续参阅图6所示，当冲击波传播路径某节点处出现三分叉工况一情况时，则按照以下公式计算分叉传播衰减系数：

其中，α₄为坑道轴线所在2条直线所夹钝角；

为冲击波传出方向与冲击波传入方向夹角呈钝角的冲击波分叉传播衰减系数；

为冲击波传出方向与冲击波传入方向夹角呈锐角的冲击波分叉传播衰减系数。

继续参阅图7所示当冲击波传播路径某节点处出现所述三分叉工况二情况时，则按照以下公式计算分叉传播衰减系数：

若α′₄∈(0°，90°]，则

若α′₄∈(90°，180°)，则

其中，α′₄为与冲击波传入方向不一致的传出方向与冲击波传入方向的夹角；

为冲击波传出方向与冲击波传入方向一致的分叉传播衰减系数；

为冲击波传出方向与冲击波传入方向不一致的分叉传播衰减系数。

参阅图8所示，其为本发明实施例所提供的冲击波在坑道段四分叉工况传播的情况及对应衰减系数示意图，当冲击波传播路径某节点处出现所述四分叉情况时，则按照以下公式计算分叉传播衰减系数：

其中，α″₄为坑道轴线所在2条直线所夹的锐角或直角；

为冲击波传出方向与冲击波传入方向一致的分叉传播衰减系数；

为冲击波传出方向与冲击波传入方向呈钝角的分叉传播衰减系数；

为冲击波传出方向与冲击波传入方向呈锐角的分叉传播衰减系数。

确定坑道轴线所在2条直线所夹钝角α₄、与冲击波传入方向不一致的传出方向与冲击波传入方向的夹角α′₄以及坑道轴线所在2条直线所夹的锐角或直角α″₄时，首先基于邻接矩阵判断坑道网络图中该节点是否为拐点，即按照节点度等于3或4来筛选，所述节点度为邻接矩阵中某坑道网络图节点所在行的元素之和，节点度等于3表明有3个坑道网络图节点与该节点相邻，等于4表明有4 个坑道网络图节点与该节点相邻。

对于路径R中节点度等于3的分叉点，假设某一分叉点为v_k，与其相邻的节点分别为v_k-1和v_k+1，满足v_k-1,v_k,v_k+1∈R，且还有1个不属于路径R但也与v_k相邻的节点v_x，则该分叉点夹角α₄或α′₄，以及分叉衰减系数ξ₄按以下不同情况确定：

如果

与

不共线且

与

不共线，则属于所述三分叉工况一，首先按照以下公式计算

其中，若

若

则

如果

与

或

共线，则属于三分叉工况二，直接按照以下公式计算α′₄：

其中，如果

与

共线，则

如果

与

共线，则

对于路径R中节点度等于4的分叉点，假设某一分叉点为v_k，与其相邻的传播节点分别为v_k-1和v_k+1，满足v_k-1,v_k,v_k+1∈R，还有2个不属于路径R但也与v_k相邻的节点v_x1、v_x2，则该分叉点夹角α″₄按不同情况计算：

如果

与

共线，则首先按照以下公式计算

其中，若

则

若

则

且

如果

与

不共线，则首先按照以下公式计算

其中，若

则

若

则

参阅图9所示，其本发明实施例所提供的冲击波从小截面坑道传入大截面坑道的情况及对应衰减系数示意图，断面变化衰减是指当爆炸冲击波从小截面坑道传入大截面坑道时，由于扩散作用强度会减小；而当冲击波从大截面坑道传入小截面坑道时由于空气压缩作用强度会变大(此时衰减系数大于1)。强度变化大小取决于爆炸冲击波传出坑道截面积与传入坑道截面积的比值，因此，根据下式所述确定衰减系数ξ₅；

其中，D₁为冲击波传入坑道等效直径，

S₁为冲击波传入坑道截面积；D₂为冲击波传出坑道等效直径，

S₂为冲击波传出坑道截面积。

步骤六:判定是否还有节点未计算爆炸冲击波入射超压峰值，若有节点未计算爆炸冲击波入射超压峰值，则返回步骤三，重新确定计算节点；若没有节点未计算爆炸冲击波入射超压峰值，则进行步骤六。

步骤七：按照最不利原则，用以下线性插值公式计算坑道段(v_i,v_j)中任意位置处离散点v_x的冲击波超压峰值P_x：

其中，

为坑道段(v_i,v_j)在v_i点的全局最大入射冲击波超压峰值；

为坑道段(v_i,v_j)在v_j点的全局最大入射冲击波超压峰值；|v_iv_j|表示坑道段(v_i,v_j) 长度。

本领域技术人员应当明白，实际情况中坑道段中间任意位置处的实际冲击波超压峰值应小于该坑道段节点两端冲击波超压峰值在该处的线性插值，因此，才可以按照最不利原则，用所述线性插值公式计算坑道段(v_i,v_j) 中任意位置处离散点v_x的冲击波超压峰值P_x。

实施例二，提供一种用于计算网状坑道内爆炸冲击波传播的工程化算法的应用案例，参阅图10-12所示，其为本发明实施例所提供的应用实例坑道路线图、应用实例坑道截面图以及应用实例坑道距离权网络图；在本实施例所提供的应用案例中，假设一个实验坑道，所述实验坑道由11个坑道段构成，其中的数字标号1-10为坑道网络节点，特殊的，0号标号为爆炸点，坑道段1-2长度为90m，坑道段2-3长度为53m，坑道段3-4长度为127m，坑道段4-5长度为4-3m，坑道段5-6的长度为55m，坑道段4-7为41m，坑道段5-8为68m，坑道段3-8为127m，坑道段7-8为48m，坑道段7-10为50m，坑道段8-9为50m；β1为爆炸冲击波入射角度，β2～β7为坑道段夹角，其中β1＝60°，β2＝135°，β3＝135°，β4＝105°，β5＝90°，β6＝60°，β7＝90°；爆炸物当量为1000kgTNT，继续参阅图 11所示，其为本实施例实验坑道其中，坑道穹顶的高度H2为1m，坑道本体高度的高度H1为2m坑道本体的宽度D为2m。

通过本发明实施例一所提出的算法，计算得到的各节点最不利冲击波超压值如下表1所示；

表1不同方法得到各节点最不利超压值(MPa)

节点	指标一	指标二	指标三	综合最不利超压峰值
					1	0.2581	0.2581	0.2581	0.2581
2	0.0716	0.0716	0.0716	0.0716
					3	0.0169	0.0169	0.0169	0.0169
4	0.0018	0.0018	0.0018	0.0018
					5	0.0005	0.0005	0.0005	0.0005
6	0.0001	0.0002	0.0001	0.0002
					7	0.0003	0.0003	0.0003	0.0003
8	0.0015	0.0015	0.0015	0.0015
					9	0.0002	0.0002	0.0002	0.0002
10	0.0001	0.0001	0.0001	0.0001

表中，各节点最不利值为取三个指标中的最大值，例如，在6号节点中，综合最不利超压峰值在指标2中取得最大值0.0002MPa，因此其最不利超压值取指标二中的最大值0.0002MPa。

参阅图13所示，其为本发明实施例所提供的应用实例坑道内冲击波超压分布图，按照2米精度对坑道段进行离散化，计算得到坑道网络内部各离散点的冲击波超压值。

至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征作出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。

Claims (7)

1.一种用于计算网状坑道内爆炸冲击波传播的方法，其特征在于，包括：

步骤一：构建网状坑道的坑道网络图；

步骤二：根据爆炸点位置对所述坑道网络图更新，进行初始节点参数初始化，确定爆炸冲击波入射超压峰值初始计算节点；

步骤三：确定所述坑道网络图中的任意一个未计算爆炸冲击波入射超压峰值的坑道网络内部节点作为下一个计算节点，并转入步骤四；

步骤四：预设表征衰减因素的全局静态指标，根据所述全局静态指标，均按照网络最短路方法计算由所述初始计算节点到达所选定内部节点的冲击波最不利传播路径；

所述全局静态指标包括指标一、指标二以及指标三，所述指标一为传播累计距离，所述指标二为传播累计节点度，所述指标三为传播累计体积；

步骤五：根据所述步骤三中得出的三条最不利传播路径以及步骤二中更新的坑道网络图，计算三条最不利传播路径各自对应的综合衰减系数，并确定最大综合衰减系数和全局最不利传播路径，根据所述步骤二中的爆炸冲击波入射超压峰值初始计算节点和所述全局最不利传播路径，计算坑道网络节点的爆炸冲击波全局最不利入射超压峰值；

步骤六：判定是否还有内部节点未计算爆炸冲击波入射超压峰值，若有节点未计算爆炸冲击波入射超压峰值，则返回步骤三，重新确定计算节点；若没有内部节点未计算爆炸冲击波入射超压峰值，则进行步骤七；

步骤七：计算坑道段(v_i，v_j)中任意位置处离散点v_x的冲击波超压峰值P_x：

其中，

为坑道段(v_i，v_j)在v_i点的全局最大入射冲击波超压峰值；

为坑道段(v_i，v_j)在v_j点的全局最大入射冲击波超压峰值；|v_iv_j|表示坑道段(v_i，v_j)长度，|v_iv_x|表示坑道段(v_i，v_x)的长度。

2.如权利要求1所述的用于计算网状坑道内爆炸冲击波传播的方法，其特征在于，所述步骤二中，根据爆炸点位置对所述坑道网络图进行更新，确定初始计算节点时，其过程包括：

Step1：确定爆炸点是否为口部外点；

若爆炸点为口部外点，则原坑道网络图不变，初始计算节点确定为距离爆炸点最近的坑道口部节点，计算该初始计算节点的冲击波入射超压峰值ΔP₀，

其中：ΔP₀为入口处的冲击波入射超压峰值；W为爆炸物TNT当量；X₀为爆炸点到坑道口部无遮挡直线距离；

若爆炸点不是口部节点，则转入Step2；

Step2：确定爆炸点是否位于坑道内网络节点上；

若爆炸点在坑道内网络节点上，则原坑道网络图不变，将初始计算节点确定为该爆炸点所在网络节点，并依次计算到达与所述网络节点相邻的其他各节点冲击波入射超压峰值；

(a)当爆炸点正好位于坑道段拐点处时，则按照以下公式分别计算与爆炸点邻接的两个坑道网络节点的入射冲击波超压峰值ΔP₀，

其中：δ_X为爆炸点的初始爆距，所述初始爆距为冲击波传播方向上最接近爆炸点的坑道网络节点距离；S为坑道截面积，W为爆炸物TNT当量；

(b)当爆炸点正好位于坑道段三岔口或十字交叉口时，则按照以下公式计算与爆炸点邻接的三个或四个坑道网络节点的入射冲击波超压峰值ΔP₀，

其中，当所述爆炸点位于三岔口时，S′＝3·S/2，位于十字叉口时，S′＝4.S/2，S为坑道截面积，δ_X为爆炸点与相邻各节点的初始爆距；

若爆炸点不在坑道内网络节点上，则进行Step3；

Step3、确定爆炸点是否在坑道内部某坑道段上；

若爆炸点在坑道内部某坑道段上，则需要在坑道网络图中新增节点为初始计算点，其中；

设新增节点v_n+1所在坑道段为(v_a，v_b)，更新原坑道网络图G(V，E)为G′(V′，E′)，将所述坑道段中心轴线交叉点及出入口端点构成的点集V＝{v_i|i＝1，2，…，n}更新为V′＝V∪{v_n+1}，坑道段边集E＝{(v_i，v_j)|v_i，v_j∈V，i≠j}更新为E′＝E∪{(v_a，v_n+1)}∪{(v_n+1，v_b)}-{(v_a，v_b)}；初始计算节点确定为新增节点v_n+1且其与爆炸点重合，计算与爆炸点邻接的两个坑道网络节点的入射冲击波超压峰值ΔP₀，

其中：δ_X为爆炸点的初始爆距，所述初始爆距为冲击波传播方向上最接近爆炸点的坑道网络节点距离；S为坑道截面积，W为爆炸物TNT当量。

3.如权利要求1所述的用于计算网状坑道内爆炸冲击波传播的方法，其特征在于，

所述指标一用下式表示：

其中：L(·)为坑道段长度函数，L(R_i，R_i+1)表示冲击波传播路径R中相邻的第i个和第i+1个节点所在坑道段的坑道长度；R为由两两相邻的坑道网络节点所构成的连通的冲击波传播路径，为一组坑道网络节点序列；R_i，R_i+1为路径R中沿传播方向的第i和第i+1个节点；|R|为路径中节点个数；R1为I₁(R)最小值对应的冲击波最不利传播路径，以指标一为准求解冲击波最不利传播路径，将坑道段长度表示坑道网络的边权，然后采用网络最短路方法求解最不利传播路径；

所述指标二用以下公式表示：

其中：Deg(·)为坑道网络节点度函数，Deg(R_i)表示与路径R中第i个节点相连的坑道段的条数，可用网络邻接矩阵中该点所在行的非0项总个数来计算；R_i为路径R中沿传播方向的第i个节点；R2为I₂(R)最小值对应的冲击波最不利传播路径；

以指标二为准求解冲击波最不利传播路径，先计算坑道网络图中所有节点的度数，再用坑道段两端节点度的平均值表示坑道网络的边权，从而转化为网络最短路方法求解最不利传播路径；

所述指标三用以下公式表示：

其中，L(·)为坑道段长度函数，L(R_i，R_i+1)表示路径R中相邻的第i个和第i+1个节点所在坑道段的坑道长度；S_i，i+1为传播路径中相邻的第i个和第i+1个节点所在坑道段的通道截面积，R3为经优化得到的累计体积最小的路径；

先计算坑道网络所有坑道段长度和通道截面积，再用求得各坑道段的体积表示坑道网络的边权，从而转化为网络最短路方法求解最不利传播路径。

4.如权利要求1所述的用于计算网状坑道内爆炸冲击波传播的方法，其特征在于，所述综合衰减系数是冲击波从初始计算节点开始，按照某一冲击波传播路径R，经过坑道网络内部衰减到达计算节点后的超压峰值相对初始节点超压峰值的比例。

5.如权利要求4所述的用于计算网状坑道内爆炸冲击波传播的方法，其特征在于，所述步骤五中，计算所述综合衰减系数时，根据爆炸点所处的坑道位置，分别采用相应的综合衰减系数计算方法，其过程包括：

若爆炸点为坑道口部外点时，所述综合衰减系数为，

其中：n_i表示第i种传播衰减在某次冲击波传播路径R中重复出现的次数，i＝1，2，3，4，5，n_i为整数；ξ₁为口部传入衰减系数；ξ₂为直线传播衰减系数；ξ₃为拐点传播衰减系数；ξ₄为分叉传播衰减系数；ξ₅为断面变化衰减系数；

若爆炸点位于坑道内部时，则不存在口部传入衰减情况，所述综合衰减系数为，

其中：n_i表示第i种传播衰减在某次冲击波传播路径中重复出现的次数，i＝2，3，4，5，n_i为整数；ξ₂为直线传播衰减系数；ξ₃为拐点传播衰减系数；ξ₄为分叉传播衰减系数；ξ₅为断面变化衰减系数。

6.如权利要求5所述的用于计算网状坑道内爆炸冲击波传播的方法，其特征在于，所述步骤五中，最大综合衰减系数指以初始计算节点为传播路径起点且到达坑道内部节点为传播路径终点的冲击波传播路径R1，R2，R3所对应的综合衰减系数的最大值，其中R1为I₁(R)最小值对应的冲击波最不利传播路径，R2为I₂(R)最小值对应的冲击波最不利传播路径；R3为经优化得到的累计体积最小的路径。

7.如权利要求6所述的用于计算网状坑道内爆炸冲击波传播的方法，其特征在于，所述综合衰减系数的最大值为

其中，R1为I₁(R)最小值对应的冲击波最不利传播路径，R2为I₂(R)最小值对应的冲击波最不利传播路径；R3为经优化得到的累计体积最小的路径；

坑道网络内部节点v_q的爆炸冲击波全局最不利入射超压峰值为

ΔP_p为起始点v_p的入射超压峰值。

Images