CN112464335A - 高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法 - Google Patents

Abstract

本发明公开了一种高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法，首先建立高大空间复杂建筑结构及周边环境的三维模型；然后分析危险源房间及周边房间的破损、超压情况；以在不同时间节点下以不同颜色包围球表示冲击波等效压力，以不同色块表示高大空间复杂建筑结构破坏区域。本发明能够增强损伤评估的可操作性和可视性。

Description

高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法

技术领域

发明涉及工程规划的技术领域，具体涉及一种高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法。

背景技术

四氧化二氮、偏二甲肼具有较高的爆炸危险性，是场区各类设施的重大危险源。为了保障人员、产品及地面设施设备的安全，需要针对场区可能存在的爆炸事故进行数值仿真与安全评估。为了增强厂房内外意外爆炸事故对建筑物及附属房间人员的损伤评估的可操作性和可视性，需要一种高大空间复杂结构内危险品爆炸可视化仿真分析方法。

发明内容

有鉴于此，本发明提供了一种高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法，能够增强损伤评估的可操作性和可视性。

本发明采用的技术方案如下：

高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法，包括以下步骤：

步骤1，建立高大空间复杂建筑结构及周边环境的三维模型；

步骤2，将高大空间复杂建筑结构及周边环境划分为若干个子空间；所有子空间的波阵面都将经历由未扰动区、波前区、波阵区最后到波后区的过程，将每个波前区子空间再划分为子网格，并将所述建筑结构的墙面与子网格的面重合；

步骤3，判断波前区子空间的边界是空气还是墙面，若边界是空气，则采用透射边界条件；若边界是墙面，则采用反射边界条件，且对于墙面边界，若在冲击波和墙面相互作用中墙面没有破坏，则波阵区子空间的冲击波不传递到相邻子空间，若墙面被破坏，计算墙面与冲击波相互作用时墙面所受到的载荷，冲击波强度根据所述载荷进行衰减，然后传递给相邻子空间；所述波阵区子空间变为波后区子空间，所述相邻子空间变为新的波阵区子空间，衰减后的冲击波由新的波阵区子空间传递到下一相邻子空间，以此类推，直到冲击波传递到周边环境或不传递；

步骤4，所有子空间由未扰动区、波前区、波阵区变为波后区后，在不同时间节点下以不同颜色包围球表示冲击波等效压力，以不同色块表示高大空间复杂建筑结构破坏区域。

进一步地，所述高大空间复杂建筑结构进行分级层次建模，所述周边环境进行三维示意建模。

进一步地，所述墙面是否破环根据能量原则进行初步估算，同时提取出采用反射边界条件的墙面上的压强，将压强转化为节点力输出含载荷信息的文件，并使用仿真软件进行验证计算，若墙面破环结果有误则修正。

进一步地，求解危险品爆炸后空气作用力时，对欧拉运动微分方程在时间方向上采用一阶的迎风格式进行离散：

|λ_max|_i+1/2＝|u|+a

|λ_max|_j+1/2＝|v|+a

|λ_max|_k+1/2＝|w|+a

其中，ε＝0.5

U为爆炸初始能量，是守恒变量，F、G、H分别为守恒变量在三个方向的通量，n为迭代次数，x、y、z坐标方向分别为东西方向、高度方向、南北方向，下标i、j、k为子网格节点上的通量大小，u、v、w分别为东西方向、高度方向、南北方向三个方向的速度，a为节点上的流体声速，ε为对数之基数，λ_max为中间变量，无实际含义。

进一步地，墙面与冲击波相互作用时墙面所受到的载荷为

其中，c_wall为墙面材料的波速，ρ_air为空气密度，c_air为空气的声速，e_break为墙面单位体积的破坏能，

σ_y为材料在弹性阶段的屈服应力，σ_s为材料发生破坏时的应力，E为材料在弹性阶段的杨氏模量，E_t为材料进入塑性后的模量。

进一步地，冲击波强度衰减后传递的能量的计算方法为：

其中，l_react＝tc_air＝h_wallc_air/c_wall为相互作用过程中扰动在空气中传递的距离，l为子网格的边长，e₀为冲击波的初始能量，e_final为冲击波衰减后传递的能量；

在得到e_final后，若e_final<0，则物理量信息不用传递给相邻的子空间；若e_final>0，则墙面发生破坏，边界上的物理量需要经过衰减后传递给相邻的子空间，能量由e₀衰减为e_final。

进一步地，冲击波等效压力的显示根据空间某点所处压力范围区间选取对应颜色包围球来表示，渲染时颜色包围球由爆心向四周逐渐传播开去，球体出现地方即表示被超压所覆盖；

高大空间复杂建筑结构破坏区域的显示根据墙面某块区域所受破坏程度范围区间选取对应颜色色块来表示，渲染时机与超压传播到来时刻相对应，超压指e_final>0。

有益效果：

1、本发明利用虚拟仿真与数值计算技术实现高大空间复杂建筑结构内危险品爆炸分析，可视化分析建筑物内冲击波超压分布范围以及结构破坏效果，验证防护墙结构对近距离内爆炸事故的防护效果，增强了损伤评估的可操作性，最后给出建筑物破坏和人员安全分析以及防护建议。

2、本发明简化爆炸冲击载荷及其与结构相互作用的演化过程，分成两步完成流体区域和结构的模拟，突破了科学计算可视化中大数据量的存储管理加载显示等一系列技术难点。

3、本发明求解危险品爆炸后空气作用力时的求解方程具有较强的鲁棒性，能够处理超强间断的问题，对于爆炸波刚开始传播时的超高压具有较好的计算稳定性。

附图说明

图1为波阵面区域演化过程示意图；

图2为波阵面跟踪算法信息记录传递示意图；

图3为高大空间复杂建筑结构CAD模型示意图；

图4为建筑模型和爆炸点及监测点示意图；

图5为爆炸分析数据可视化流程示意图；

图6为超压分布与结构破坏情况三维展示示意图。

具体实施方式

下面结合附图并举实施例，对本发明进行详细描述。

本发明提供了一种高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法，步骤如下：

步骤1，建立高大空间复杂建筑结构及周边环境的三维模型。

考虑到实际厂房结构非常复杂，包含柱、梁、板等基本结构，以及门、窗、沟、洞等各种开口，还有各种装饰元素，连接关系也多种多样，建立完整的几何模型具有很大难度，即便能够完成建模，庞大的网格数量也无法用于数值计算。考虑到厂房地面以上的主要承载部件包括柱、梁、楼板、剪力墙、砌块墙等等，门、窗等构件及小尺寸开口不会对厂房的结构和承载能力造成很大影响，所以几何模型主要包括柱、梁、楼板、剪力墙、砌块墙、屋面、大门等要素。建模处理方法如下：

步骤101：建立各段(层)的图层群，包括建筑构件、平台、火箭、爆炸源、监测点图层，其中建筑构件指柱、梁、楼板、剪力墙、砌块墙、屋面、大门。

步骤102：根据建筑构件的材料，建立“1_普通钢筋砼”、“2_高强钢筋砼”、“3_砌块”、“4_Q235钢”和“5_Q345钢”等五种典型材料名称，并将每个图层赋予对应的材料信息。

步骤103：建立存储柱、梁截面尺寸的文件，包含截面说明、截面编号、截面宽、截面高共4个元素。

步骤104：分段(层)建立几何模型，如图3所示，在每段(层)中一般按照柱->梁->剪力墙->砌块墙->大门->楼板的顺序进行。每建立一个柱、梁时，需要赋予截面和材料信息，并将截面编号加入到步骤103的文件中。

步骤105：使用二次开发的Autolisp程序对几何模型进行处理，获取几何信息、检查重叠、剖分几何体，保证相邻几何体共用节点。处理完毕后输出剖分后的dwg文件以及txt格式的数据文件。

步骤106：3D Max建模。按照输出剖分后的dwg文件，等比例开展点线面建模，场区模型精度控制在分米级，建筑外观模型精度控制在厘米级，建筑房间、构件模型精度控制在毫米级。

步骤107：纹理贴图。不满足要求的纹理不能正常显示，因此为了场景驱动的要求，需对用到的纹理图像作以下处理：

步骤107-1：转换纹理格式。将纹理文件转换成3D Max中支持的文件格式：TGA、PNG、BMP、JPG；

步骤107-2：标准化纹理大小。将纹理文件长、宽全部限制为2的幂次单位的象素大小；

步骤107-3：透明纹理。如果在RGB文件中包含了a值，那么文件的后缀名必须是RGBA才能正常显示。

步骤108：模型简化。模型的复杂度是影响实时渲染的关键因素，因此在不影响视觉效果的前提下，应尽量减少模型复杂度：

步骤108-1：去除冗余面。这里的冗余多边形主要是指在实体外部观察模型时不可见的部分。由于场景浏览时它们处于不可见的位置，去除它们并不影响实体的视觉效果，但消除它们可以很大程度上降低场景的复杂度。

步骤108-2：“合并”面片。在影响不大的情况下，造型中相近的面片应予以合并。

步骤2：开展爆炸分析，重点分析危险源房间及周边房间的破损，超压情况。

(1)基本思路为：

首先，由于爆炸产生的冲击波本质上是一个波传播的问题，冲击波第一次到达结构并和结构的相互作用是最强的，建设第一次相互作用后，由于耗散，反射波的强度会显著降低，其破坏作用不会超过冲击波首次与结构的作用，所以忽略反射波等作用。

其次，由于只考虑冲击波的首次作用，所以只需考虑模拟冲击波波前的演化过程。可以将整个空间分为若干个子空间，每个子空间尺寸为1m×1m×1m。冲击波从爆源开始，每一步向外扩展一层子空间，将每个波前子空间再细分为10×10×10至20×20×20的计算子网格，使得网格精度满足冲击波计算的基本要求，通过依次计算每层子空间集从而实现对波前演化的模拟。

然后，为了反应出墙面破坏后冲击波继续传播的特性，需要将结构的每个墙面都与子网格的面重合，即结构节点坐标需要按米m取整。

计算时，结构面按照固壁边界条件进行处理。在对波阵面进行模拟后，根据能量法则定性判断该墙面是否被破坏，若破坏，则对冲击波的强度进行衰减，在冲击波下一步扩展过程中传递到相邻的子空间；若没有破坏，则该波阵面不会向下传递。计算中对墙面的破坏仅使用能量原则进行初步估算，可能存在较大误差，为此在对冲击波模拟的同时，提取出模拟为固壁边界条件(反射边界条件)的墙面上的压强，将其转化为节点力并进行记录，计算结束后输出含载荷信息的K文件，使用LS-DYNA进行验证计算，若墙面破环结果有误则修正。，使得流体区域和结构的模拟分成两步完成，分别进行计算。在施加爆炸载荷后结构的模拟由LS-DYNA完成，技术可行，计算量可接受。

计算过程涉及：波阵面演化的计算流体力学格式，墙面的破坏模型，以及流体的数据信息在不同子空间的传递。

(2)波阵面演化的计算流体力学格式

求解危险品爆炸后空气作用力时，对空气的模拟求解如下所示的三维欧拉Euler运动微分方程：

其中，

U为爆炸初始能量，是守恒变量，u、v、w分别为东西方向、高度方向、南北方向三个方向的速度，ρρuρvρwρE分别为质量、三个方向的动量和能量，F,G,H为守恒变量在三个方向的通量。

在空间上，将Euler方程采用中心型的鲁沙诺夫Rusanov格式进行离散：

其中，

|λ_max|_i+1/2＝|u|+a

|λ_max|_j+1/2＝|v|+a

|λ_max|_k+1/2＝|w|+a

ε＝0.5

下标i、j、k为子网格节点上的通量大小，a为节点上的流体声速即空气传播速度，ε为对数之基数，λ_max为中间变量，无实际含义。

在时间方向上采用一阶的迎风格式进行离散，所以最终的求解格式为：

该格式具有较强的鲁棒性，能够处理超强间断的问题，对于爆炸波刚开始传播时的超高压具有较好的计算稳定性。

(3)边界条件处理

主要有两种边界条件：透射边界条件和光滑固壁边界条件(反射边界条件)。在求解区域外部布置一层虚拟的求解节点，根据不同的边界条件给虚拟节点赋值不同的物理量。若边界是空气，则采用透射边界条件；若边界是墙面，则采用反射边界条件，且对于墙面边界。

对于透射边界条件，直接将边界节点的物理量映射到虚拟节点即可，即：

ρ_ghost＝ρ_bd；ρu_ghost＝ρu_bd；ρv_ghost＝ρv_bd；ρw_ghost＝ρw_bd；ρE_ghost＝ρE_bd

其中，下标ghost表示虚拟节点上的物理量，下标bd表示边界节点的物理量。

对于光滑固壁边界条件，将沿边界法向的速度映射为零，其余物理量保持不变即可。具体而言，对x_min和x_max处的边界有：

ρ_ghost＝ρ_bd；ρu_ghost＝0；ρv_ghost＝ρv_bd；ρw_ghost＝ρw_bd；ρE_ghost＝ρE_bd

对y_min和y_max处的边界有：

ρ_ghost＝ρ_bd；ρu_ghost＝ρu_bd；ρv_ghost＝0；ρw_ghost＝ρw_bd；ρE_ghost＝ρE_bd

对z_min和z_max处的边界有：

ρ_ghost＝ρ_bd；ρu_ghost＝ρu_bd；ρv_ghost＝ρv_bd；ρw_ghost＝0；ρE_ghost＝ρE_bd

这种边界条件的提法可以准确描述冲击波在求解边界处反射和透射的现象。

(4)波阵面传播

(a)子空间分类

按照1m×1m×1m的划分方法把整个空间划分为若干个子空间，由于简化之后的问题只考虑波阵面所在的子空间内的演化，所以首先要提取需要计算的子空间。按照冲击波传播的规律，波阵面只会向外扩大，所以每一个子空间只会被波阵面通过一次，为此可以将波阵面分为以下四种状态：未扰动区(波阵面距离该子空间尚且遥远，暂时完全不用考虑)、波前区(波阵面进入该子空间，将要被进行计算)、波阵区(已经进行过计算的子空间，但是其相邻的子空间尚未被计算)、波后区(已经进行过计算的子空间，且其相邻的子空间也已经计算完毕)。

所有区域都将经历由“未扰动区”，“波前区”，“波阵区”最后到“波后区”的过程。以平面问题为例，如图1所示：初始爆炸物所在的子空间处于“波前区”状态，其余子空间均为“未扰动区”状态；对“波前区”子空间进行计算，计算完毕后，该子空间变为“波阵区”状态，与它相邻的子空间变为“波前区”状态；对所有“准备”子空间进行计算并变成“波阵”状态，与“准备”子空间相邻的“未扰动区”子空间变为“波前区”状态，原来的“波阵区”状态变为“波后区”状态；重复上述过程，各个需要计算的子空间的状态变更为“波后区”。

可以看出，“波后区”子空间的范围在不断扩大，最后当所有需要计算的子空间全部变成“波后区”后则认为计算结束，完成对空气部分的模拟。

由于需要进行计算的“波前区”子空间是一层一层向外扩展的，所以把每一次扩展称为一个循环，每一个循环内完成对所有当前“波前区”子空间的计算并更新各个子空间所处的状态使得波阵面能传播出去。

(b)“波前区”子空间内的计算

每一个循环内的主要工作量在于对所有“波前区”子空间进行计算。由于所有子空间的尺寸均为1m×1m×1m，所以首先将子空间在空间上进行网格划分为子网格，并且所有子空间需要按照同样的尺寸进行网格划分。之后再计算网格上按照第(2)节波阵面演化的计算流体力学格式的最终求解格式进行求解，得到波阵面在该“波前区”子空间内的演化情况。

各个“波前区”子空间的初始状态为常温常压，即p＝0.1Mpρ＝1.225kg/m³，波阵面的信息是通过子空间边界传送过来，在上一个循环中，“波阵区”子空间已经计算完成，得到该子空间内的波阵面的演化过程。为了将这个演化过程传递给与其相邻的“波前区”子空间，需要把图2中“波阵区”子空间边界上所有节点的物理量演化信息记录下来。

为了避免所有时间步信息记录储存量过大的问题，将“波阵区”子空间的求解时间等分为100份，记录该边界所有节点上这100个时刻的物理量的值。在计算“波前区”子空间时，这个边界上记录的信息会传递给该子空间作为一个和时间相关的边界条件进行处理，这样便能够让波阵面的信息在“波前区”子空间内继续传播。由于各个子空间网格划分的方式是一致的，所以“波阵区”子空间和“波前区”子空间的公共边界处的节点是重合的，故“波阵区”子空间记录的物理量信息可以直接施加到“波前区”子空间的边界上，

“波前区”子空间的边界条件会在下一个循环中变为时间相关的边界条件，它将会为与其相邻的“未扰动区”子空间(下一个循环成为“波前区”子空间)提供波阵面的信息，如图2的绿色边界部分。

随着冲击波的扩散，波速降低，格子的计算时间变长，所以当前计算的“波阵区”格子的时间相关边界条件结束点会早于当前格子的计算时间，后续边界条件将按最后一个时间点的数值进行指数衰减。

(c)初始子空间的计算

初始子空间由于包含有炸药，属于对多介质的问题，需采用LS-DYNA进行计算。首先确定房间内的爆炸物当量，之后建立用于计算的k文件，将其提交给LS-DYNA进行计算后提取出单元格六个边界面上各个节点的密度和压强信息(提取从爆炸波传递到边界开始后的100个时刻的物理量的值)作为下一个循环中准备子空间的时间相关边界条件。在采用LS-DYNA进行爆炸计算时，网格可能划分得较密，所以此时其边界上的节点位置和相邻的子空间的节点位置不一定一致。为此，在决定初始时间相关边界条件时，首先寻找出与各个节点最接近的LS-DYNA上的计算节点，然后将该LS-DYNA节点上的物理量信息赋值给时间相关边界上的对应节点。

(d)波阵面强度衰减

在对每个“波前区”子空间进行计算时，若边界是空气，则采用透射边界条件，若边界是墙面，则采用反射边界条件。对于空气边界，物理量直接穿过该边界传递给相邻的子空间，但是对于墙面边界，若在冲击波和墙面相互作用中墙面没有破坏，则不会有信息传递到相邻的子空间，若墙面被破坏，则冲击波的强度会有一定的衰减然后再传递给相邻的子空间。由于在对空气进行模拟时不会考虑到其和墙面结构的整体作用，所以这里只能根据能量守恒的原理定性地判断墙面是否被破坏以及波阵面强度的衰减量。

在冲击波和墙面的相互作用中，设参与相互作用的空气质量为m_air，密度为ρ_air，相互作用之前的压强为p₀，相互作用结束之后压强衰减为p₁。由于相互作用的时间非常短，所以参与相互作用的空气的体积来不及发生变化，故作用前后空气的密度保持不变，都为ρ_air。

在相互作用中，能量转化的形式主要是空气的内能转化为墙面的变形能，若墙面最终破坏则空气的部分内能会被墙面破坏所吸收，从而导致压强的减小。墙面单位体积的破坏能e_break由材料的本构方程决定：

其中σ_y为材料在弹性阶段的屈服应力，σ_s为材料发生破坏时的应力，E为材料在弹性阶段的杨氏模量，E_t为材料进入塑性后的模量。因此，面积为S的墙面发生破坏时吸收的能量为：

E_break＝S_wallh_walle_break

其中h_wall为墙面的厚度。

另一方面，对空气而言，由能量守恒有：

m_aire₀-m_aire₁＝S_wallh_walle_break

e₀为冲击波的初始能量，e₁为墙面破坏后冲击波的能量，即：

为了确定参与相互作用的空气的质量，首先需要估算出相互作用的时间。由于墙面的破坏是由于一个一个强的应力波在墙中传播导致的，所以相互作用时间可以近似认为是波穿过墙面所用的时间，即：

其中

是墙面材料的波速。在冲击波穿过墙面的同时，扰动也会以c_air在空气中传播，c_air为空气的声速，被扰动传过的空气则认为是参与相互作用的空气，故参与相互作用的空气质量为：

将其带入到压强衰减的表达式中得到：

一般情况下，固体的声速会远远大于气体中的声速，而且墙面也相对较薄，所以在相互作用过程中扰动不容易穿过一个计算网格。而在子空间传递边界物理量时传递的是边界网格内的物理量，所以需要把边界网格内的参与相互作用的空气和没参与相互作用的空气的内能按体积份额进行平均后再传递给相邻的子空间，即：

其中，l_react＝tc_air＝h_wallc_air/c_wall为相互作用过程中扰动在空气中传递的距离，l为计算网格即子网格的边长尺寸，e_final为冲击波衰减后传递的能量。在得到e_final后，若e_final<0，则说明冲击波提供的能量不足以使墙面发生破坏，所以墙面保持不变，物理量信息不用传递给相邻的子空间；若e_final>0，则说明冲击波提供的能量能够使墙面发生破坏，边界上的物理量需要经过衰减后传递给相邻的子空间。对密度而言，相互作用前后密度保持不变，所以密度不用衰减；对三个方向的速度而言，在相互作用中，先假设了墙面是固壁，所以三个方向的速度都重置为零；对能量而言，由于经过了破坏的吸能过程，内能由e₀衰减为e_fina^l。

(e)结构载荷施加

当“波前区”子空间的边界是墙面时，需要确定当其和冲击波相互作用时墙面上所受到的载荷。由于子空间的节点与结构墙面上网格节点同样不一定完全匹配，所以与之前一样，首先确定与结构网格点最接近的一个子空间网格节点，然后便认为该结构网格节点在此时刻受到的压强即为最接近的子空间节点的压强。最后将结构受到的压强减去标准大气压后乘以结构网格的面积便得到结构网格点在该时刻的节点力。

在对冲击波传播的整个过程中，所有的结构墙面都会被遍历一次，当冲击波模拟计算结束后，便得到了冲击波与结构相互作用的载荷，最后将其提交给LS-DYNA计算即可得到最终的结构响应。

在利用软件进行仿真操作时，利用Fortran90语言编写实现计算程序，具备数据读入，计算、输出监测点压力，整个空间的最大压强分布，墙体破坏程度估算等信息，同时可生成含载荷信息的K文件，用于后续使用LS-DYNA程序计算结构的响应。

步骤201：输入每个子空间的信息，包括其子空间的对角坐标，子空间的每个面的属性(包括墙面材料，厚度等信息)。另外输入各个材料的参数，得到材料声速c_wall和单位体积破坏能e_break。

步骤202：输入需要输出压强的观测点，如图4所示。

步骤203：分配每个子空间计算冲击波演化时所需要的内存。

步骤204：输入结构网格，主要是壳单元和梁单元的节点信息。

步骤205：输入采用LS-DYNA模拟初始格子爆炸时所采用的网格。

步骤206：输入用LS-DYNA模拟得到的边界处密度和压强的演化信息。

步骤207：分配初始的6个时间相关边界条件的内存并根据LS-DYNA的结果得到这些时间相关边界条件的具体信息。

步骤208：根据当前循环结束后的各个子空间状态确定下一个循环各个子空间的状态并返回下一个循环中有多少个“波前区”子空间等待被计算。

步骤209：将上一个循环中的透射及经过衰减后的反射边界条件上的物理量演化信息赋值给这一个循环中的时间相关边界条件，之后释放掉上一个循环中时间相关边界条件所用的内存。

步骤210：根据这个循环中的“波前区”子空间确定出需要储存演化过程的边界面并分配其所需要的内存。这些边界面上的信息之后会作为时间相关边界条件传递给下一次循环。

步骤211：将上一个循环中的“波阵区”子空间转化为“波后区”子空间。

步骤212：确定子空间的各个面上的边界条件类型(透射边界条件或反射边界条件)并将其储存。

步骤213：根据子空间的各个面上的边界条件类型对子空间进行计算。

步骤213-1：确定完成冲击波演化模拟所需要的物理时间。

步骤213-2：确定冲击波在演化的起始时间以及终止时间。

步骤213-3：根据“波前区”子空间的对角坐标生成其计算网格。

步骤213-4：给定计算网格上节点初始常温常压状态。

步骤213-5：根据子空间各个边界条件的类型确定出当前时刻各个边界上的物理量，包括透射边界条件和反射边界条件。

步骤213-6：计算时间步长。

步骤213-7：储存透射和反射边界面上的物理量信息。

步骤213-8：计算网格节点上x方向上的通量F。

步骤213-9：计算网格节点上y方向上的通量G。

步骤213-10：计算网格节点上z方向上的通量H。

步骤213-11：在完成对子空间的计算后输出其边界前面上的载荷信息。

步骤214：对墙面边界，根据衰减模型

对冲击波的强度进行衰减。

步骤215：确定下一个循环中各个时间相关边界条件的起始时间。

步骤216：把这个循环的“波前区”子空间变更为“波阵区”子空间。

步骤217：判断是否所有子空间均已经“波后区”。

步骤3：利用流体作用于结构的载荷，将带载荷的结构模型提交给LS-DYNA进行计算，得到厂房在冲击作用下的响应。分析建筑物内冲击波超压的分布范围，重点关心人员伤亡和建筑物破坏的阈值，根据《爆炸基本原理》，冲击波超压低于0.01MPa时，人员基本无损伤，超过0.075MPa时人员死亡，中间分为轻伤、中伤和重伤三等；根据《爆破安全规程》(GB6722-2003)，冲击波超压低于0.009MPa时，建筑物次轻度破坏，墙体基本无损，超过0.076MPa时建筑完全破坏，中间分为轻度破坏、中等破坏、次严重破坏、严重破坏等，因此应对0.009～0.076MPa的范围进行重点研究，给出建筑物破坏和人员安全分析以及防护建议。

步骤301：打开LS-DYNA软件，并点开Solver求解器菜单。

步骤302：在下拉列表中选择start LS-DYNA Analysis命令按钮。

步骤303：弹出Start Input and Output对话框，在Input File 1＝中点击browse，选择K文件。

步骤304：点击Run。

步骤305：求解器会对K文件进行求解计算。

步骤4：基于Quest3D软件平台编写爆炸分析数据可视化程序，涉及三维场景、仿真数据的动态加载、空间坐标变换和三维渲染等，流程如图5所示。主要是采用模型数据与程序分离方法，将所有模型数据进行外部存储并且支持动态加载和卸除，突破了高大空间复杂结构模型和超压、结构破坏数据的大数据量存储管理加载技术难题，有效减少了系统复杂性，采用不同时间节点的基于不同颜色包围球表示的冲击波等效压力和基于不同色块表示的建筑破坏区域的方法，实现了数据到可视化的直观映射，在确保冲击波毁伤情况展示效果的前提下，避免了动态替换结构破坏构件模型所带来的高成本、低效率的问题。

步骤401：构建用于Quest3D程序动态加载的外部模型、数据文件。

步骤401-1：构建外部模型文件。

(1)3dsMax模型导出

从3dsMax导出时记录场区位置信息，坐标归零后导出，导出后按原坐标在Quest3D中进行坐标位移转换。导出时记录的位置信息包括：场区编码、场区名称、位置坐标(X，Y，Z)、旋转角度，其中，X坐标指向东西方向，Y坐标指向高度方向，Z坐标指向南北方向，旋转方向为自顶向下看，顺时针为正。

建筑外观模型坐标原点选取在建筑外墙的零平面西南角点。在3dsMax中记录坐标原点位置信息，将坐标归零后进行导出，导出后按原坐标在Quest3D中进行坐标位移转换。导出时记录的位置信息包括：建筑编码、建筑名称、位置坐标(X，Y，Z)、旋转角度、建筑尺寸(长、宽、高)，其中，X坐标指向东西方向，Y坐标指向高度方向，Z坐标指向南北方向，旋转方向为自顶向下看、顺时针为正，建筑长为东西方向、按最大规则形状计算，建筑宽为高度方向，按最大规则形状计算，建筑高为南北方向，按最大规则形状计算。

建筑房间模型坐标原点选取在顶视图中房间内部内沿左下角点。记录房间坐标原点与建筑坐标原点相对坐标值，不进行坐标归零，直接导出，导入至Quest3D后房间位于建筑内部，不需进行坐标转换，但此时需要使用检验程序进行坐标检查。导出时记录的位置信息包括：房间编码、房间名称、位置坐标(X，Y，Z)、房间尺寸(长、宽、高)，其中，X坐标指向东西方向，Y坐标指向高度方向，Z坐标指向南北方向，房间长为东西方向、按最小规则形状计算，房间宽为高度方向，按最小规则形状计算，房间高为南北方向，按最小规则形状计算。

建筑构件模型与建筑房间模型导出方法一样。

(2)Quest3D程序中创建模型位置信息表

根据模型导出时记录的位置信息，在Quest3D程序中构建模型位置信息表，表字段主要包括：模型名称(Name)、位置坐标(Position<Vector>)、旋转角度(Rotation<Vector>)和尺寸大小(Size<Vector>)等信息，Quest3D程序运行时将对模型位置进行自动赋值，减少人为操作难度。

(3)Quest3D程序中构建工况-模型关联关系

将3dsMax导出的模型文件(*模型.cgr)导入至Quest3D程序中，在程序中将出现文件中包含的模型数据，创建模型实例并按照各工况仿真需要用到的场区、建筑外观、房间、构件类模型，建立工况与模型实例间的关联关系，然后在程序中删除导入的模型数据，并按照模型文件物理存储地址重新建立外部模型文件中的模型数据与程序中对应的模型实例间的映射关系，至此，可供动态加载的模型文件构建完成。

步骤401-2：构建外部数据文件。构建系统外部存放的仿真数据文件(*数据.cgr)，包括超压数据表、结构破坏数据表、工况-仿真节点表、工况-仿真节点-超压数据索引表和工况-仿真节点-结构破坏数据索引表。其中，超压数据表、结构破坏数据表即为步骤2中计算所得数据，按照仿真节点的数量构建相应数量的超压数据表、结构破坏数据表，分别命名为Blockpressure_i、Breakratio_i，分别用于集中存储在第i个仿真节点时所有工况的超压和结构破坏数据，超压数据表包括序号、超压球中心点坐标(X、Y、Z)和超压值，结构破坏数据表包括序号、毁伤面端点1坐标(X、Y、Z)、毁伤面端点2坐标(X、Y、Z)、毁伤面端点3坐标(X、Y、Z)、毁伤面端点4坐标(X、Y、Z)、毁伤率；工况-仿真节点表命名为FileNumInEachLogic，用于存储各工况的仿真节点数量，包括工况序号、该工况仿真节点数量，如不同的爆炸当量、爆炸源不同的放置位置；工况-仿真节点-超压数据索引表命名为DataNumInEachLogic&File，用于存储各工况各仿真节点中包含的超压数据量，包括工况序号、仿真节点序号和数据量；工况-仿真节点-结构破坏数据索引表命名为DataNumInEachLogic&File_BR，用于存储各工况各仿真节点中包含的结构破坏数据量，包括工况序号、仿真节点序号和数据量。

步骤402：选择工况。设计导航界面，按照“场区—建筑单体—单体房间—工况名称”组织结构对所有工况进行分类索引，帮助快速定位并进入相应的工况模块。

步骤403：动态加载模型、数据。

步骤403-1：动态加载模型。按照所选工况，依据程序中建立的工况-模型关联关系，按照存放地址快速查找外部存放的模型文件，并按需加载其中包含的模型至内存中，释放之前工况模型所占据的内存。

步骤403-2：动态加载数据。(1)动态加载外部数据文件，向Quest3D程序中即时拷贝工况-仿真节点表和工况-仿真节点-超压数据索引表；(2)根据工况序号，从工况-仿真节点表中获取该工况的仿真节点数M，在系统中即时建立M个临时超压数据表、结构破坏数据表，分别命名为BlockpressureL_i、BreakratioL_i，(i＝1..M),分别用于集中存储在第i个仿真节点时该工况的超压和结构破坏数据，临时超压数据表包括序号、超压球中心点坐标(X、Y、Z)和超压值，临时结构破坏数据表包括序号、毁伤面端点1坐标(X、Y、Z)、毁伤面端点2坐标(X、Y、Z)、毁伤面端点3坐标(X、Y、Z)、毁伤面端点4坐标(X、Y、Z)、毁伤率。(3)根据工况序号，从工况-仿真节点-超压数据索引表中获取该工况各仿真节点中包含的数据量，换算成在超压数据表、结构破坏数据表中的存储位置，并依据存储位置抽取该工况所有的超压、结构破坏数据写入临时超压数据表、结构破坏数据表中。

步骤404：按照仿真节点序号，开始仿真循环。

步骤405：按当前仿真节点序号i，从相应的临时超压数据表、结构破坏数据表BlockpressureL_i、BreakratioL_i中抽取所有超压、结构破坏数据，并写入内存。

步骤406：开展室内爆炸超压与结构破坏效果三维可视化展示。

室内爆炸超压显示会根据空间某点所处压力范围区间选取对应颜色球体来表示，渲染时颜色球体会由爆心向四周逐渐传播开去，出现球体的地方即表示被超压所覆盖；结构破坏显示根据结构面某块区域所受破坏程度范围区间选取对应颜色色块来表示，其渲染时机与超压传播到来时刻相对应，可形象直观地展示超压对结构产生了一定的破坏作用，倘若超压穿过墙体继续传播，即表示该墙体已被彻底破坏，此时墙面色块颜色代表的破坏率应为1。超压指e_final>0。

步骤406-1：根据压力值所处的范围区间，换算得出超压球体的颜色。

步骤406-2：根据四个端点坐标，计算毁伤色块的中心点坐标、大小及朝向，根据破坏率所处范围区间换算色块颜色。

给定毁伤色块四个端点P1、P2、P3、P4的坐标(X1,Y1,Z1),(X2,Y2,Z2),(X3,Y3,Z3),(X4,Y4,Z4)，中心点P坐标(X,Y,Z)计算方法：(X,Y,Z)＝((X1,Y1,Z1)+(X2,Y2,Z2)+(X3,Y3,Z3)+(X4,Y4,Z4))/4；

毁伤色块的大小(XSize，YSize，ZSize)计算方法：

XSize＝if(((X1-X2+X1-X3+X1-X4)＝0)？0:1),1；else 0.01；

YSize＝if(((Y1-Y2+Y1-Y3+Y1-Y4)＝0)？0:1),1；else 0.01；

ZSize＝if(((Z1-Z2+Z1-Z3+Z1-Z4)＝0)？0:1),1；else 0.01。

步骤407：按位置、大小、朝向以及颜色值赋值给球体和色块，主要是Motion-Position(Value Vector)、Motion-Rotation(Value Vector)、Motion-Size(ValueVector)、Surface-Material-Emissive(ValueVector)。

步骤408：渲染该节点的所有球体和色块，如图6所示。

本方法为国内首次实现了高大空间复杂结构内危险品爆炸可视化仿真分析方法，建立了场区高大空间复杂结构的爆炸冲击波超压求解方法，解决了爆炸分析结果空间表达的科学计算可视化难题。经过实际任务检验，通过分析给出的建筑物破坏和人员安全分析以及防护建议，能够科学指导任务实施，确保任务安全。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法，其特征在于，包括以下步骤：

步骤1，建立高大空间复杂建筑结构及周边环境的三维模型；

步骤2，将高大空间复杂建筑结构及周边环境划分为若干个子空间；所有子空间的波阵面都将经历由未扰动区、波前区、波阵区最后到波后区的过程，将每个波前区子空间再划分为子网格，并将所述建筑结构的墙面与子网格的面重合；

步骤3，判断波前区子空间的边界是空气还是墙面，若边界是空气，则采用透射边界条件；若边界是墙面，则采用反射边界条件，且对于墙面边界，若在冲击波和墙面相互作用中墙面没有破坏，则波阵区子空间的冲击波不传递到相邻子空间，若墙面被破坏，计算墙面与冲击波相互作用时墙面所受到的载荷，冲击波强度根据所述载荷进行衰减，然后传递给相邻子空间；所述波阵区子空间变为波后区子空间，所述相邻子空间变为新的波阵区子空间，衰减后的冲击波由新的波阵区子空间传递到下一相邻子空间，以此类推，直到冲击波传递到周边环境或不传递；

步骤4，所有子空间由未扰动区、波前区、波阵区变为波后区后，在不同时间节点下以不同颜色包围球表示冲击波等效压力，以不同色块表示高大空间复杂建筑结构破坏区域。

2.如权利要求1所述的高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法，其特征在于，所述高大空间复杂建筑结构进行分级层次建模，所述周边环境进行三维示意建模。

3.如权利要求1所述的高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法，其特征在于，所述墙面是否破环根据能量原则进行初步估算，同时提取出采用反射边界条件的墙面上的压强，将压强转化为节点力输出含载荷信息的文件，并使用仿真软件进行验证计算，若墙面破环结果有误则修正。

4.如权利要求1所述的高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法，其特征在于，求解危险品爆炸后空气作用力时，对欧拉运动微分方程在时间方向上采用一阶的迎风格式进行离散：

|λ_max|_i+1/2＝|u|+a

|λ_max|_j+1/2＝|v|+a

|λ_max|_k+1/2＝|w|+a

其中，ε＝0.5

U为爆炸初始能量，是守恒变量，F、G、H分别为守恒变量在三个方向的通量，n为迭代次数，x、y、z坐标方向分别为东西方向、高度方向、南北方向，下标i、j、k为子网格节点上的通量大小，u、v、w分别为东西方向、高度方向、南北方向三个方向的速度，a为节点上的流体声速，ε为对数之基数，λ_max为中间变量，无实际含义。

5.如权利要求1所述的高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法，其特征在于，墙面与冲击波相互作用时墙面所受到的载荷为

其中，c_wall为墙面材料的波速，ρ_air为空气密度，c_air为空气的声速，e_break为墙面单位体积的破坏能，

σ_y为材料在弹性阶段的屈服应力，σ_s为材料发生破坏时的应力，E为材料在弹性阶段的杨氏模量，E_t为材料进入塑性后的模量。

6.如权利要求5所述的高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法，其特征在于，冲击波强度衰减后传递的能量的计算方法为：

其中，l_react＝tc_air＝h_wallc_air/c_wall为相互作用过程中扰动在空气中传递的距离，l为子网格的边长，e₀为冲击波的初始能量，e_final为冲击波衰减后传递的能量；

在得到e_final后，若e_final<0，则物理量信息不用传递给相邻的子空间；若e_final>0，则墙面发生破坏，边界上的物理量需要经过衰减后传递给相邻的子空间，能量由e₀衰减为e_final。

7.如权利要求6所述的高大空间复杂建筑结构内危险品爆炸可视化仿真分析方法，其特征在于，冲击波等效压力的显示根据空间某点所处压力范围区间选取对应颜色包围球来表示，渲染时颜色包围球由爆心向四周逐渐传播开去，球体出现地方即表示被超压所覆盖；

高大空间复杂建筑结构破坏区域的显示根据墙面某块区域所受破坏程度范围区间选取对应颜色色块来表示，渲染时机与超压传播到来时刻相对应，超压指e_final>0。

Images