一种大规模并行计算的自适应扩大计算域的数值模拟方法与流程

本发明涉及一种大规模并行计算的自适应扩大计算域的数值模拟方法，属于数值模拟领域。

背景技术：

爆炸毁伤实验测试能够快速并准确得到冲击波压力值，但是冲击波速度在几千米每秒，近场压力一般高于十个大气压，对传感器、示波器等测量技术要求较高。国外学者利用实验手段对不同药量、不同位置处的压力值进行测量及标定，拟合出了广泛应用的henrych经验公式(j.henrych,thedynamicsofexplosionanditsuse,elsevier,newyork,1979)、mills经验公式(c.mills,1stint.conf.hazardprot.edinburgh,uk,1987)、brode经验公式(h.brode,phys.fluids,1959,2:217-229)，但是这些经验公式仅在一定假设下成立，研究结果存在很大差距；经验公式中带有实验确定的系数，不同的研究者所得出的系数各不相同，不适用于实际工程问题；针对大药量的炸药或不同边界条件，在选用公式需考虑实际情况。此外，实验测试的不足在于：(1)测量仪器难以在高速、高压下得到正确的响应。(2)近场压力极高，并产生高温火球，对传感器耐冲击、耐高温的性能要求较高，且易损坏。(3)多次进行实验误差较大，具有人为、环境等多种不确定因素。

随着数值模拟技术的发展与高性能计算机的出现，数值仿真已经成为爆炸毁伤研究的主要技术手段。技术的优势在于：(1)相比于实验研究，数值模拟非常安全且成本较低。(2)研究爆炸毁伤机理的效率较高，操作简单。而爆炸毁伤问题涉及高速、高压、高温、相变，气体、液体和固体等多种介质间相互耦合甚至混合，材料不但会发生严重变形和破碎，还会融化甚至汽化。国外学者进行了一些管道内爆炸的数值模拟研究，如yanez等人(j.yanez,internationaljournalofhydrogenenergy,2011,36:2613-2619)采用三种方法进行了rut管道内爆炸研究，由于管道尺寸较大，因此采用了比爆炸胞格大的网格宽度，无法得到准确的模拟结果，同时这三种普通方法计算耗时较长；heidari(a.heidari,internationaljournalofhydrogenenergy,2011,36:2538-2544)采用较粗的网格(网格宽度5cm)与20cpu+2200000并行求解rut管道爆炸问题，由于计算方法与计算资源的限制，同样无法获得准确的计算结果。因此，在进行大规模高精度并行计算时，无法高效、准确的得到数值模拟结果。本发明提出一种大规模并行计算的自适应扩大计算域的数值模拟方法，能够在计算资源一定的条件下，对大范围爆炸流场快速求解，三维求解速度远远快于普通方法。综上，本发明提出的方法具有重要的技术背景。

技术实现要素：

本发明的目的是提供一种大规模并行计算的自适应扩大计算域的数值模拟方法；该方法在计算资源一定的条件下，既满足高精度，又能够高效的对大规模爆炸流场问题进行模拟，显著提高效率，进而解决大规模高精度数值模拟技术领域的工程问题。

所述大规模高精度数值模拟技术领域包括爆炸毁伤、航空航天、机械工程领域。

本发明的目的是通过下述技术方案实现的：

本发明提出一种大规模高精度并行计算的自适应扩大计算域方法。具体而言，该方法是一种基于本质无震荡(weno)有限差分方法、水平集方法(level-set)和真实虚拟流体方法(rgfm)的算法的自适应网格放大方法，能够大幅提高大规模问题的计算效率。weno有限差分格式利用间断处权重为零，将欧拉方程中的通量求解转换为离散求解，在保证精度的情况下，同时降低了数值震荡；level-set方法能够隐式追踪复杂多介质界面的位置演化；rgfm方法利用虚拟流场处理界面状态，保持了界面的压力、速度连续特性；自适应网格放大方法能够在一定的计算资源情况下，能够对米、公里尺度的大范围流场进行高精度计算，显著提高计算效率。通过该方法，能高效高精度的模拟爆炸流场传播过程。

一种大规模并行计算的自适应扩大计算域的数值模拟方法，包括如下步骤：

步骤1：根据需要模拟的爆炸流场尺度以及计算资源(cpu核数和内存大小)，确定最终计算域尺寸，x、y、z三个方向并行分区数量，给定自适应扩大计算域方法中的扩大次数、扩大比率，以及每个分区中的网格数；给定计算域各cpu的序号；确定初始计算域l＝0时x、y、z方向的网格宽度、相应计算域的尺寸及各计算网格点的坐标。

步骤1.1：确定最终计算域尺寸[πx1,πx2]*[πy1,πy2]*[πz1,πz2]，x、y、z三个方向并行分区数量cutx、cuty、cutz，以及每个分区中的网格数nx×ny×nz，其中nx表示x方向网格数量，ny表示y方向网格数量，nz表示z方向网格数量。并行使用的cpu总量为cutx×cuty×cutz，爆炸流场以计算域中心点o(xo,yo,zo)向x-、x+、y-、y+、z-、z+方向扩大，x、y、z三个方向扩大的次数largex、largey、largez，扩大比率r，可以得到网格数总量为cutx*cuty*cutz*nx*ny*nz。

步骤1.2：给定计算域各cpu的序号。

定义myid为cpu序号(从0开始)；posx、posy、posz表示该cpu处于x方向的第posx个、y方向的第posy个、z方向的第posz个cpu。posx的取值范围是1至cutx，posy的取值范围是1至cuty，posz的取值范围是1至cutz。

则cpu序号myid与posx、posy、posz的关系为：

posx＝mod(myid,cutx)+1

posy＝myid/cutx+1(1)

posz＝myid/(cutx×cuty)+1

其中mod为取余符号。

步骤1.3：初始计算域l＝0时计算区域的x、y、z方向的网格宽度分别为

初始计算域l＝0时计算域尺寸ω0为：

其中r^largex为r的largex次方，r^largey为r的largey次方，r^largez为r的largez次方，其中πx1为最终计算域x方向最左侧坐标,πx2为最终计算域x方向最右侧坐标,πy1为最终计算域y方向最左侧坐标,πy2为最终计算域y方向最右侧坐标,πz1为最终计算域z方向最左侧坐标,πz2为最终计算域z方向最右侧坐标。

因此，初始计算域内每个网格的x、y、z方向坐标(xi,yjzk)为：

其中i,j,k为该网格在此cpu中的x方向第i个，y方向第j个，z方向第k个网格。

步骤2：定义初始计算区域内初始状态t＝0时刻的用于区分爆轰产物与空气交界面位置的level-set函数；对爆轰产物与空气两区域物理量赋初始值；建立描述爆轰产物与空气的守恒型欧拉控制方程组；建立爆轰产物与空气两种材料的状态方程函数。

步骤2.1：定义最小(初始)计算区域内初始状态t＝0时刻的用于区分爆轰产物与空气交界面位置的level-set函数，利用符号距离函数给出任意时刻多物质界面的位置并实时追踪level-set函数的演化过程，t时刻的界面γ(t)位于level-set函数φ(xi,yj,zk,t)为零的位置：

γ(t)＝{(xi,yj,zk)|φ(xi,yj,zk,t)＝0}(5)

t＝0时刻，爆轰产物与空气交界面位置可由计算域中每个计算网格点的位置与界面位置的距离求得：

其中，x、y、z分别为计算网格点的物理空间坐标，d(γ(0),(xi,yj,zk))表示计算网格点(xi,yj,zk)到初始界面γ(0)的距离，s爆轰产物与s空气分别为物质交界面两侧的区域。

步骤2.2：给定步骤2.1得到的爆轰产物与空气两个区域的初始状态。

给定爆轰产物的初始密度ρ1、三个方向速度u1、v1、w1、压力p1；或者爆轰产物的初始密度ρ1、三个方向速度u1、v1、w1、单位质量的内能e1；

以及空气的初始密度ρ2、三个方向速度u2、v2、w2、压力p2；或者空气的初始密度ρ2、三个方向速度u2、v2、w2、单位质量的内能e2。

步骤2.3：构建描述爆轰产物与空气的守恒型欧拉控制方程组。

守恒型欧拉控制方程组为：

其中，

u＝[ρρuρvρwρe]^t，

f＝[ρuρu²+pρuvρuwu(ρe+p)]^t，

g＝[ρvρuvρv²+pρvwv(ρe+p)]^t，

h＝[ρwρuwρvwρw²+pw(ρe+p)]^t，

且e为单位质量的总能量，若计算爆轰产物，将式(7)中的ρ、u、v、w、p、e用ρ1、u1、v1、w1、p1、e1代入；若计算空气，将式(7)中的ρ、u、v、w、p、e用ρ2、u2、v2、w2、p2、e2代入。

步骤2.4：选取材料的状态方程函数，进而确定材料在变形、运动下的物理状态。

由于式(6)中仅有五个方程，但共有六个未知数，因此需要选取状态方程对整个控制方程组进行封闭。爆轰产物的单位质量的内能e1与压力p1的jwl状态方程关系表示为：

其中ρ0是炸药的初始密度，a、b、r1、r2、ω是状态方程参数。空气的理想气体状态方程关系表示为：

p2＝(γ-1)ρ2e2(9)

其中γ是多方气体指数。

步骤3：设置各级区域的边界条件。

由于计算域在不同时刻的尺寸不同，因此每一级的计算域边界条件不相同。针对空间离散格式的不同，每一个计算网格点需要周围n个点的物理量计算通量，故在计算域边界处需要n层虚拟点。

(1)对于流出边界，在每一级计算域边界处，有

①x方向且当posx＝1时，有

i＝1,2,...,n,j＝0,1,...,ny,k＝0,1,...,nz；

②x方向且当posx＝nx时，有

i＝1,2,...,n,j＝0,1,...,ny,k＝0,1,...,nz；

③y方向且当posy＝1时，有

i＝0,1,...,nx,j＝1,2,...,n,k＝0,1,...,nz；④y方向且当posy＝ny时，有

i＝0,1,...,nx,j＝1,2,...,n,k＝0,1,...,nz；

⑤z方向且当posz＝1时，有

i＝0,1,...,nx,j＝0,1,...,ny,k＝1,2,...,n；

⑥z方向且当posz＝nz时，有

i＝0,1,...,nx,j＝0,1,...,ny,k＝1,2,...,n

(2)对于固壁有滑移边界和对称边界，

①x方向且当posx＝1时，有

φ(x-i,yj,zk,t)＝φ(xi,yj,zk,t)，

i＝1,2,...,n,j＝0,1,...,ny,k＝0,1,...,nz；

②x方向且当posx＝nx时，有

φ(xnx+i,yj,zk,t)＝φ(xnx-i,yj,zk,t)，

i＝1,2,...,n,j＝0,1,...,ny,k＝0,1,...,nz；

③y方向且当posy＝1时，有

φ(xi,y-j,zk,t)＝φ(xi,yj,zk,t)，

i＝0,1,...,nx,j＝1,2,...,n,k＝0,1,...,nz；

④y方向且当posy＝ny时，有

φ(xi,yny+j,zk,t)＝φ(xi,yny-j,zk,t)，

i＝0,1,...,nx,j＝1,2,...,n,k＝0,1,...,nz；

⑤z方向且当posz＝1时，有

φ(xi,yj,z-k,t)＝φ(xi,yj,zk,t)，

i＝0,1,...,nx,j＝0,1,...,ny,k＝1,2,...,n；

⑥z方向且当posz＝nz时，有

φ(xi,yj,znz+k,t)＝φ(xi,yj,znz-k,t)，

i＝0,1,...,nx,j＝0,1,...,ny,k＝1,2,...,n

所述流出边界包括如无限域空气；

所述固壁有滑移边界包括刚性地面；

所述对称边界包括1/2模型、1/4模型或1/8模型；

步骤4：选取第l级计算时间步长。

选取第l级计算时间步长δtl：

其中，参数cfl为介于0到1之间的常数，δxl为第l级计算区域中x方向的网格宽度；δyl为第l级计算区域中y方向的网格宽度；δzl为第l级计算区域中z方向的网格宽度；u、v、w为当前计算网格点的三个方向速度，若处于爆轰产物区域，则取为u1、v1、w1，若处于空气区域，则取为u2、v2、w2；c为当前计算网格点的声速。

步骤5：利用level-set函数求解每一级区域中物质界面的法向量，求解爆轰产物与空气界面处的法向黎曼问题，并采用真实虚拟流体方法确定界面附近n层网格的物理状态，更新全流场物理状态。

采用level-set方程对流场中界面位置进行推进，对于计算域内每一个计算网格点，均采用下式推进，为了简化公式，将φ(xi,yj,zk,t)简写为φ；

若计算网格点处于爆轰产物区域，将式(11)中的u、v、w用u1、v1、w1代入；若计算网格点处于空气区域，将式(11)中的u、v、w用u2、v2、w2代入。

每一个计算网格点的法向量n求解公式为：

先将每个计算网格点的速度按该点法向量进行分解，分解为法向速度与切向速度，再利用靠近界面附近的爆轰产物、空气的密度、法向速度、压力物理量状态，根据界面两侧法向速度连续、压力连续以及守恒定律，利用黎曼求解器给出界面处的法向速度、压力以及界面两侧爆轰产物与空气的密度，然而法向速度与切向速度为局部坐标系下的物理状态，需要保持界面附近的计算网格点的法向速度不变，对速度分量逆向变换得到全局坐标系下的速度。在界面两侧爆轰产物与空气区域内分别建立一层20个网格宽度的带状区域，带状区域内的计算网格点为虚拟点；虚拟点处的密度、三个方向速度及压力通过延拓得到。延拓方程为：

其中i表示虚拟点处物理量，所述虚拟点处物理量包括密度、三个方向速度和压力。最后通过利用高精度weno格式对爆轰产物与空气分别求解。

步骤6：判断当前计算域是否需要扩大。

在计算一段时间步后，初始计算域中的爆轰产物向外膨胀，即将到达初始计算域边界，采用冲击波(物理量突变)及稀疏波(物理量缓慢上升或下降)分别进行分析。判断是否需要扩大；如果判定为不需要扩大，则继续运行步骤3，若判定即将扩大，则执行步骤7-9；

当冲击波即将传播至达初始计算域边界时，所处的位置xm为临界点；若临界点满足如下阈值条件时，即扩大计算域，

|im-i0|＞icr(14)

其中，i可代表密度、压力；i0为空气处的常值密度、常值压力，im为临界点的密度、压力；icr为临界阈值，一般取为较小的数，使得当冲击波到达临界点时，计算域自动扩大，防止冲击波流出计算域边界导致重要信息丢失。

当稀疏波即将传播至达初始计算域边界时，所处的位置xn为临界点；当该点满足如下阈值条件时，计算域即将扩大，

其中i可代表密度、压力；in为临界点的密度、压力；ic'r为临界梯度阈值，一般取为较小的数、且与网格宽度有关。如果仍然采用冲击波阈值(如式14)计算，在有限计算步长内，在边界点xn处会产生“阶跃”，从而引入误差。因此，采用式(15)来判断是否扩大，能够防止稀疏波流出计算域边界导致重要信息丢失。

步骤7：由于计算过程中，爆轰产物向外传播，压力与密度向外扩散，第l级计算域边界处满足步骤6中的扩大阈值条件，此时l’＝l+1，即为扩大后计算域，需要给定扩大后计算域的计算域、相应的网格宽度以及每个计算域内每个计算网格点的坐标。

将计算域分级为l＝0,1,...,maxlarge，其中maxlarge为largex、largey、largez中的最大值，其中l＝0代表初始计算域，l＝maxlarge代表最终计算域。因此，扩大后第l’级计算区域的x、y、z方向的网格宽度分别为

则第l级计算域尺寸ωl’为：

其中πx1为计算域x方向最左侧坐标,πx2为计算域x方向最右侧坐标,πy1为计算域y方向最左侧坐标,πy2为计算域y方向最右侧坐标,πz1为计算域z方向最左侧坐标,πz2为计算域z方向最右侧坐标。

第l’级计算域内每个计算网格点的x、y、z方向坐标(xi,yjzk)为：

其中i,j,k为该网格在此cpu中的x方向第i个，y方向第j个，z方向第k个网格。

步骤8：对扩大后计算域(第l’级)中各cpu所有网格点的密度、三个方向速度、压力、level-set函数赋值。将扩大后范围内的的网格分为三类：新计算域点，即扩大后得到的心的网格点；原计算域网格点，即扩大前和扩大后均存在于原计算域且能够重合的网格点；原计算域新点，即扩大前和扩大后均存在于原计算域且不重合的网格点；

针对原计算域网格点，可以采用第l级计算域直接赋值的方法。

针对原计算域新点，由于空间离散格式为2n-1阶精度，因此在这些计算网格点处也需要给出2n-1阶精度的插值点，保证整体计算域一致的精度。因此，采用2n-1阶weno插值对这些点进行赋值。

针对新计算域点，可直接赋值为常温常压空气的密度、速度、压力。

步骤9：采用自适应扩大计算域mpi并行计算模块进行多cpu联合求解。

采用以上步骤，能够对单cpu的串行程序进行自适应扩大计算域快速高效求解。但是为了能够计算大规模问题，还需要引入并行计算。由于mpi并行中，每个cpu占据一块独有的内存空间，需要将扩大前区域(第l级计算域)所有cpu中的所有网格点的密度、三个方向速度、压力、level-set函数值共享至扩大后区域(第l’级计算域)的每一个cpu中，并且需要寻找第l’级计算域与第l级计算域中距离多个相临近的点。

利用并行分块搜索方法，能够大幅减少并行自适应扩大计算域模块的计算时间，提高并行效率。该并行分块搜索方法为：由于边界条件的不同，如步骤3中的流出边界与固壁有滑移边界，采用不同的并行分块搜索方法。

针对流出边界，可按对称轴将其分为4块(二维)或8块(三维)，在每一块区域内重复步骤8，得到第l’级各cpu所有网格点的密度、三个方向速度、压力、level-set函数值；

针对固壁有滑移边界和对称边界，对第l’级的部分区域赋值，得到扩大后的第l’级各cpu所有网格点的密度、三个方向速度、压力、level-set函数值，这样也可以减少3/4(三维7/8)的计算量。所述部分区域为：第l’级区域与第l级区域相重合的部分；此时计算域已扩大，生成了新的x、y、z方向的尺寸，重新设置各级区域的边界条件，即不断重复步骤3-9；直至达到规定时间，循环结束；

还包括步骤10：利用步骤1至步骤9进行基于自适应扩大计算域的高精度大规模数值模拟，降低冲击波耗散，能够高精度、高效率的解决大规模及超大规模的实际工程问题；

所述大规模高精度数值模拟技术领域包括高速杀伤武器、航空航天、机械工程领域。

预测所述大规模高精度数值模拟技术领域中的大规模计算过程，包括激波打气泡，炸药空中爆炸、近地面爆炸、水下爆炸、聚能射流，爆炸冲击波与建筑物相互作用，多组分混合气体化学爆炸以及超新星爆炸。

有益效果：

1、本发明公开的一种大规模并行计算的自适应扩大计算域方法，能够在计算资源一定的条件下，既满足高精度，又能够高效的对大规模爆炸流场问题进行求解，显著提高求解效率，进而准确预测激波打气泡，炸药空中爆炸、近地面爆炸、水下爆炸、聚能射流，爆炸冲击波与建筑物相互作用，多组分混合气体化学爆炸以及超新星爆炸等大规模高精度数值模拟技术领域的工程问题。

2、本发明公开的一种大规模并行计算的自适应扩大计算域方法，采用weno格式对欧拉方程空间导数项进行离散，采用tvdrunge-kutta格式对欧拉方程时间导数项进行离散，相比传统方法具有更高的精度，能够降低冲击波及稀疏波传播过程中的耗散。

3、本发明公开的一种大规模并行计算的自适应扩大计算域方法，采用levelset与rgfm相结合的黎曼求解方法处理多物质界面的相互作用的强间断问题，具有物质界面高分辨率的间断特性；相比于经典的gfm方法，能有效地处理复杂流动过程，能有效抑制界面附近的非物理震荡现象，保证数值模拟过程的稳定运行。

4、本发明公开的一种大规模并行计算的自适应扩大计算域方法，利用冲击波/稀疏波扩大阈值条件、扩大后网格高精度赋值方法以及并行自适应扩大模块，高精度、高效的求解大规模爆炸问题，与传统并行方法相比，极大地提高了计算效率，在处理大规模及超大规模工程问题时，具有明显优势。

附图说明

图1为算例的数值模型及观测点示意图；

图2为算例的每一级计算域三维模型示意图；

图3为冲击波阈值条件示意图；

图4为稀疏波阈值条件示意图；

图5为自适应扩大计算域赋值方法示意图；

图6为流出边界条件下并行搜索方法示意图；

图7为对称/固壁有滑移边界条件下并行搜索方法示意图；

图8为算例t＝0.025ms时的扩大前、扩大后压力云图；

图9为算例在0.05ms、0.1ms、0.2ms、0.5ms、1.0ms、5.0ms、15ms、25ms时刻的三维压力云图；

图10为算例在距离爆炸中心3m、5m、7m、9m、11m处的压力-时间曲线数值模拟与实验对比图。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

实施例1：

本实施例公开一种大规模并行计算的自适应扩大计算域方法，为了更好的说明本发明的目的和优点，下面结合附图与实施例对本发明作进一步说明。

本实例是在2kgtnt炸药放置于离地1.6米处的近地面爆炸问题中的应用。

本实施例公开的算例中的初始几何模型在附图1与附图2中呈现，不同时刻，计算得到的爆炸冲击波、地面反射波及二者叠加形成的马赫波的传播效果在附图8和附图9中呈现，典型位置压力-时间曲线与实验的对比结果在附图10中呈现。

首先给出步骤1、2中所需的信息：如附图1和附图2所示，本实施例采用1/4模型计算，此时x-与y-方向为对称边界，x+、y+、z+为流出边界，z-在计算域未接触至地面时(第三级计算域之前即l＝0,1,2)为流出边界，计算域在贴近地面z＝-1.6之后(l＝3,4,5)为固壁有滑移边界。取largex＝largey＝largez＝5，即5级扩大，扩大比率取为r＝2。炸药放置于距离地面1.6米处，因此在设置计算域及其边界条件时，需要考虑地面的影响。最终计算域尺寸设定为[0,12.0]×[0,12.0]×[-1.6,11.2]，x、y、z方向并行分区数量cutx＝cuty＝6、cutz＝4，每个cpu中的网格数为nx×ny×nz＝50×50×80，即共有144个cpu参与运算，网格总量高达2880万。设置计算停止时刻为40ms。爆炸中心选取为o(0,0,0)，在计算域未接触至地面时(第三级计算域之前即l＝0,1,2)，向x+、y+、z-、z+方向扩大，计算域在贴近地面z＝-1.6之后(l＝3,4,5)，仅向x+、y+、z+扩大。

计算域的cpu序号为myid，取为0-143，cpu序号myid与posx、posy、posz的关系为：

posx＝mod(myid,6)+1

posy＝myid/6+1

posz＝myid/36+1

初始计算域网格宽度为δx0＝δy0＝δz0＝0.00125m。

初始计算域的尺寸为ω0＝[0,0.375]×[0,0.375]×[-0.2,0.2]

初始计算域内每个网格的x、y、z方向坐标(xi,yjzk)为

xi＝i·δxl+(posx-1)·50·δxl

yj＝j·δyl+(posy-1)·50·δyl

zk＝k·δzl+(posz-1)·80·δzl-0.2

tnt炸药中心位于(0，0，0)，质量为2kg，其半径为0.0664m，周围其他物质为常温常压空气。初始时刻(l＝0计算域)每个点的levelset符号距离函数(炸药定义为负，空气定义为正)为

炸药初始状态：密度ρ1＝1630kg/m³、速度为u1＝v1＝w1＝0、压力p1＝8600mpa；

空气初始状态：密度ρ2＝1.225kg/m³、速度为u2＝v2＝w2＝0、压力p2＝0.101mpa；

建立如式(7)的守恒型欧拉控制方程组，并依照炸药与爆轰产物所在的区域将两种物质的密度、速度、压力代入式(7)中。

炸药采用jwl状态方程，其参数为：a＝373800mpa、b＝3747mpa、r1＝4.15、r2＝0.9、ω＝0.35，空气采用理想气体状态方程，其参数为γ＝1.4。

由此，完成步骤1、2。

接着，循环进行步骤3至步骤5。由于本实施例采用1/4模型计算，并考虑到地面的存在，因此设置x-与y-方向为对称边界，x+、y+、z+为流出边界，z-在计算域未接触至地面时(第三级计算域之前即l＝0,1,2)为流出边界，计算域在贴近地面z＝-1.6之后(l＝3,4,5)为固壁有滑移边界。第l级计算时间步长δtl由式(10)给出。利用式(11)与式(12)可以对流场中界面位置进行推进，并利用高精度weno格式与真实虚拟流体方法更新全流场物理状态。由此，完成步骤3-5。

同时在步骤5后进行步骤6中的判断，如果计算域边界处压力达到以下条件时，采用步骤7-9自动扩大。

当冲击波即将传播至达初始计算域边界时，当边界右侧且仅距离一个网格点(附图3中的点xm-1)时，若临界点满足如下阈值条件时，即扩大计算域，

其中超压δpn(＝pn-0.101)为x+,y+,z-和z+方向的超压，px+,py+,pz-和pz+为计算域右侧、后侧、底面、上面的压力。

当稀疏波即将传播至边界右侧且仅距离一个网格点(附图4中的点xm-1)时，当该点满足如下阈值条件时，计算域即将扩大，

其中，由于地面的存在，底面压力的判断p^*z-仅在第0，1，2级计算域时有效，底面压力可写做

步骤7中，扩大后的计算域网格宽度为

δx1＝δy1＝δz1＝0.0025m

δx2＝δy2＝δz2＝0.005m

δx3＝δy3＝δz3＝0.01m

δx4＝δy4＝δz4＝0.02m

δx5＝δy5＝δz5＝0.04m

扩大后计算域的尺寸为

ω1＝[0,0.75]×[0,0.75]×[-0.4,0.4]

ω2＝[0,1.5]×[0,1.5]×[-0.8,0.8]

ω3＝[0,3.0]×[0,3.0]×[-1.6,1.6](贴近地面，此时及以后不向z-方向扩大)

ω4＝[0,6.0]×[0,6.0]×[-1.6,4.8]

ω5＝[0,12.0]×[0,12.0]×[-1.6,11.2](最终计算域)

第l＝1,2,3级计算域中每个网格的x、y、z方向坐标(xi,yjzk)为：

xi＝i·δxl+(posx-1)·50·δxl

yj＝j·δyl+(posy-1)·50·δyl

第l＝4,5级计算域中每个网格的x、y、z方向坐标(xi,yjzk)为：

xi＝i·δxl+(posx-1)·50·δxl

yj＝j·δyl+(posy-1)·50·δyl

zk＝k·δzl+(posz-1)·80·δzl-1.6

针对步骤8的赋值方法，下面分为三类点进行分析，如附图5中所示的第l+1级计算域中原计算域网格点(蓝色圆圈)、原计算域新点(蓝色叉)以及新计算域点(红色圆圈)。

针对原计算域网格点，可以采用第l级计算域直接赋值的方法。

针对原计算域新点，由于空间离散格式为5阶精度，因此在这些蓝色叉点处也需要给出5阶精度的插值点，保证整体计算域一致的精度。因此，采用5阶weno插值对这些点进行赋值。

针对新计算域点，可直接赋值为常温常压空气的密度、速度、压力。

针对步骤9的方法，本实施例在计算域未接触至地面时(第三级计算域之前即l＝0,1,2)采用流出边界并行分块搜索方法，如附图6所示，可按对称轴将其分为8，在每一块区域内进行保存数据、搜索距离和赋值，这样能够提高三维问题7/8的效率。当计算域在贴近地面z＝-1.6之后(l＝3,4,5)，采用固壁有滑移边界并行分块搜索方法，如附图7所示，对第l’级计算域中的第l’级区域与第l级区域相重合的部分区域赋值，得到扩大后的第l’级各cpu所有网格点的密度、三个方向速度、压力、level-set函数值，这样也可以减少3/4(三维7/8)的计算量。

这样，近地面爆炸的全物理过程均能够通过自适应扩大计算方法与weno格式、levelset界面追踪与真实虚拟流体界面处理方法联合求解。

结果分析：

附图8为t＝0.025ms时刻的结果图(仅展示x与z方向)，左侧为l＝0的计算域，右侧为l＝1扩大后的计算域，此时压力扰动刚到l＝0计算域的顶端与底端，导致了第一次自适应扩大，从左右两图对比可知，压力结果一致。附图9为t＝0.05ms至t＝25ms的压力云图，右下角的ω1～ω5代表现在该计算域处于哪一级。从图中可以看出计算域在慢慢扩大，并且在0.05ms时刻，压力峰值衰减且冲击波向外传播，第二级计算域的最大压力明显低于第一级。在0.074ms、0.217ms、0.730、6.531ms达到了扩大阈值条件。冲击波在约1.0ms时刻到达地面，并形成地面反射波，而地面反射波的速度较快，慢慢会追上前面的冲击波，叠加形成马赫波。

数值模拟与实验在距离爆炸中心3m、5m、7m、9m、11m处放置观测点与传感器，计算得到的压力-时间曲线与相同药量的实验曲线如附图10所示，数值模拟结果能够体现出冲击波、反射波及马赫波的形成，各点处压力-时间曲线与实验基本一致，证明了该自适应扩大计算域数值模拟技术的有效性。

文献(j.yanez,internationaljournalofhydrogenenergy,2011,36:2613-2619)进行了一系列爆炸流场的数值模拟，采用了三款软件(德国researchcenterkarlsruhe的com3d-v.3.5、俄罗斯russianresearchcentrekurchatovinstitute的b02、英国universityofulster的fluentv6.3.26)模拟rut管道爆炸问题，其精度最高为二阶，均采用结构化正方体网格，网格宽度在0.015-0.045m(大于爆炸问题的胞格尺寸)，分别利用3cpu+3832192网格、单cpu+3832192网格、单cpu+237005网格求解该问题，计算终止时刻约45ms，计算总共耗时分别约为6小时、20小时、60小时。而本实施例采用结构化正方体网格，网格宽度最小为0.00125m，本实施例利用144cpu+41472000网格进行计算，耗时约138小时。本实施例采用了精细的网格进行计算，普通方法需要943718400000网格，而采用本发明的方法仅需要41472000网格，因此若采用上面三款软件求解本实施例的问题，在理想条件下(不考虑并行数据传递等因素)则需要约30782小时、102608小时、1659104小时，远远长于本实施例的138小时。因此本发明与现有模拟方法对比，可以看出本实施例能够快速求解这类爆炸问题，并在一定计算资源条件下计算大规模爆炸问题。

本实施例的结果的冲击波、反射波与马赫波能够对人员、设备进行毁伤，设备受到高压冲击波会发生变形直至破坏，人员受到高压、低压冲击波会产生耳膜破裂、死亡等现象。因此本发明用于本实施例是很有必要的。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。