发明内容

针对现有技术中爆炸反应流预测方法存在计算效率较低、守恒性缺失等问题，本发明公开的一种基于广义黎曼解法器的爆炸反应流高精度预测方法要解决的技术问题是：实现基于广义黎曼解法器的高精度数值通量构造方法，提高高精度预测方法的计算效率；同时有效解决非守恒问题，实现对爆轰和爆燃过程的精确预测，进而解决爆炸力学领域相关工程技术问题。所述爆炸力学领域相关工程技术问题包括武器弹药设计与毁伤评估以及爆炸事故的防灾减灾。

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

本发明公开的一种基于广义黎曼解法器的爆炸反应流高精度预测方法，是通过两步四阶Lax-Wendroff方法以及守恒界面方法的结合实现的。利用水平集方法求解爆轰界面和爆燃界面的拓扑变化，快捷且精确地追踪界面演化过程。利用虚拟流体方法在爆轰界面和爆燃界面附近的计算单元内构造虚拟流体，将多介质问题转化为单介质问题。利用守恒界面方法构造爆炸反应流守恒型离散格式，通过求解含爆轰波和爆燃波的黎曼问题，构造沿着界面的守恒量交换项。结合两步四阶Lax-Wendroff方法以及守恒界面方法，构造基于广义黎曼解法器的高精度数值通量求解方法。该方法能够提高爆炸反应流高精度预测方法的计算效率；同时有效解决多介质界面非守恒问题，实现对可压缩反应流体爆轰和爆燃过程的精确预测。

本发明公开的一种基于广义黎曼解法器的爆炸反应流高精度预测方法，包括如下步骤：

步骤1：建立直角坐标系，定义x方向以及y方向的空间步长，在计算区域内设置(M+1)*(N+1)个计算单元，记为单元(i,j)，其中i＝0,1,…,M，j＝0,1,…,N。

步骤2：根据初始条件定义初始时刻计算区域内所有计算单元的水平集函数以及未燃气体和已燃气体的流体物理量，包括密度、x方向的运动速度、y方向的运动速度以及压力。

在水平集方法中，水平集函数为一种符号距离函数，用φ表示。令φ＞0和φ＜0的区域分别表示未燃气体以及已燃气体两种流体的分布区域，φ＝0的位置表示爆轰或爆燃界面。在初始时刻，采用以下公式定义水平集函数φ的数值：

其中Γ⁰表示初始时刻界面的分布曲线，d(x,y)表示计算单元中心(x,y)到初始界面Γ⁰的距离，Ω₁、Ω₂分别表示未燃气体以及已燃气体分布的区域。

根据初始条件，定义初始界面两侧未燃气体以及已燃气体的密度、x方向的运动速度、y方向的运动速度以及压力。

步骤3：建立可压缩反应流体的控制方程组，并且设置边界条件。

步骤3.1：建立可压缩反应流体的控制方程组。

应用CJ模型建立可压缩流体的控制方程组，即不考虑化学反应率，爆轰以及爆燃界面视为跳跃间断，化学反应物沿着界面立即转化为产物，并且释放热量。未燃气体以及已燃气体满足欧拉方程：

其中ρ、u、v、p、E分别表示密度、x方向的运动速度、y方向的运动速度、压力以及单位质量气体的总能。总能表示内能和动能的总和，即其中e表示单位质量气体的内能。未燃气体以及已燃气体的状态方程满足：

p＝(γ-1)ρ(e-q), (2)

其中γ表示比热比，q表示由于化学反应释放的热量参数。

步骤3.2：设置可压缩反应流体的控制方程组边界条件。

在高阶精度数值通量的构造过程中，需要使用到周围计算单元的物理量。因此，为了更新计算区域边界附近的守恒量，需要将计算区域往外延拓三层。

左右边界：

对于来流边界条件，延拓单元的守恒量设置为：

其中i＝1,2,3，j＝0,1,2,…,N，U＝(ρ,ρu,ρv,ρE)^T。U_-1,j、U_-2,j、U_-3,j分别表示计算区域左边界往左延拓的第一、第二以及第三层的守恒量，U_M+1,j、U_M+2,j、U_M+3,j分别表示计算区域右边界往右延拓的第一、第二以及第三层的守恒量，的具体数值需要根据实际工况进行设置。

对于出流边界条件，延拓单元的守恒量设置为：

其中i＝1,2,3，j＝0,1,2,…,N。即将计算区域内最靠近左边界或者右边界的那一层守恒量直接赋值给延拓单元的守恒量。

对于固壁边界条件，延拓单元的守恒量设置为：

其中i＝1,2,3，j＝0,1,2,…,N。即将计算区域内最靠近左边界或者右边界的第一、第二和第三层的质量、y方向动量以及总能分别赋值给往外延拓的第一、第二和第三层的对应物理量；并且将计算区域内最靠近左边界或者右边界的第一、第二和第三层的x方向动量的相反数分别赋值给往外延拓的第一、第二和第三层的对应物理量。

上下边界：

对于来流边界条件，延拓单元的守恒量设置为：

其中i＝0,1,2,…,M，j＝1,2,3。U_i,N+1、U_i,N+2、U_i,N+3分别表示计算区域上边界往上延拓的第一、第二以及第三层的守恒量，U_i,-1、U_i,-2、U_i,-3分别表示计算区域下边界往下延拓的第一、第二以及第三层的守恒量，的具体数值需要根据实际工况进行设置。

对于出流边界条件，延拓单元的守恒量设置为：

其中i＝0,1,2,…,M，j＝1,2,3。即将计算区域内最靠近上边界或者下边界的那一层守恒量直接赋值给延拓单元的守恒量。

对于固壁边界，延拓单元的守恒量设置为：

其中i＝0,1,2,…,M，j＝1,2,3。即将计算区域内最靠近上边界或者下边界的第一、第二和第三层的质量、x方向动量以及总能分别赋值给往外延拓的第一、第二和第三层的对应物理量；并且将计算区域内最靠近上边界或者下边界的第一、第二和第三层的y方向动量的相反数分别赋值给往外延拓的第一、第二和第三层的对应物理量。

步骤4：计算更新未燃气体和已燃气体守恒量的时间步长。

计算满足CFL条件的时间步长：

其中CFL表示CFL系数，取值范围为(0,1)；Δx和Δy分别表示x方向和y方向的空间步长；u_i,j和v_i,j分别表示计算单元(i,j)内流体沿x方向和y方向的运动速度；c_i,j表示计算单元(i,j)内流体的声速。

步骤5：根据平流方程更新水平集函数，并进行重新初始化，求解爆轰或爆燃界面的演化过程。

根据如下式所示的平流方程更新水平集函数：

φ_t+μφ_x+νφ_y＝0, (4)

其中μ和ν分别表示水平集函数沿着x方向和y方向的运动速度，需要通过以下方程组进行定义：

其中D表示爆轰或者爆燃界面的运动速度，N_x和N_y分别表示单位法方向向量N沿着x方向和y方向的分量。

应用三阶TVD Runge-Kutta方法对方程(4)进行时间离散，水平集函数在x方向和y方向的偏导数φ_x和φ_y需要根据修正的Godunov方法求得。通过下式能够求得空间导数φ_x和φ_y：

其中S(φ)表示水平集函数φ的符号函数。

其中

为了提高计算效率，只需要在界面附近的一至四层计算单元内，通过方程(4)更新水平集函数。

为了使得水平集函数维持距离函数的性质，需要通过以下方程进行重新初始化：

步骤6：根据水平集函数的数值定义爆轰或爆燃界面的几何参数。

令Γ表示爆轰或者爆燃界面，分别表示计算单元(i,j)的左边界、右边界、下边界、上边界被界面Γ切割的比例，具体数值能够通过水平集函数的算符运算计算获得。对于φ＞0的流体，A通过以下方程进行定义：

其中

对于φ＜0的流体，界面切割比例A通过方程A^-＝1-A⁺求得，其中上标-和+分别表示φ＜0和φ＞0的流体对应的界面切割比例几何参数。

各流体的体积分数α也能通过水平集函数的算符运算计算获得，对于φ＞0的流体，α通过以下方程进行定义：

其中

对于φ＜0的流体，体积分数α则通过方程α^-＝1-α⁺求得，其中上标-和+分别表示φ＜0和φ＞0的流体对应的体积分数几何参数。

步骤7：在界面附近的计算单元内根据延拓方程构造虚拟流体。

根据延拓方程，构造虚拟流体的物理量状态：

V_τ±N·▽V＝0, (12)

其中V表示需要延拓的物理量，包括密度、x方向的运动速度、y方向的运动速度和压力；N表示沿着界面的单位法方向向量，通过方程求得；τ表示人工时间步长，定义为：方程(12)中，+号用于计算φ＞0的流体在虚拟单元的物理量；-号用于计算φ＜0的流体在虚拟单元的物理量。

应用三阶TVD Runge-Kutta方法对方程(12)进行时间离散，应用迎风格式计算空间导数以及

步骤8：在爆轰或爆燃界面附近的计算单元内，构造并求解含爆轰波或爆燃波的黎曼问题，计算得到爆轰波或者爆燃波两侧的物理量状态，包括密度、运动速度、压力以及比内能。

在界面附近的计算单元内，根据真实流体和虚拟流体的状态，沿着法方向构造黎曼问题的初始条件。根据质量、动量以及能量守恒方程，爆轰波或爆燃波的波前和波后的流体满足以下关系式：

其中下标0和1分别表示爆轰波或爆燃波波前未燃气体和波后已燃气体的对应物理量，w表示流体相对于爆轰波或者爆燃波的运动速度。激波波前和波后的流体状态需要满足激波关系式；稀疏波波前和波后的流体状态需要满足黎曼不变量为常数。

当发生CJ爆轰时，黎曼解由非反应波(激波或稀疏波)、接触间断、稀疏波和CJ爆轰波组成；当发生强爆轰时，黎曼解由非反应波、接触间断和强爆轰波组成。当发生爆燃时，黎曼解由非反应波、接触间断、爆燃波和非反应波组成。根据激波、稀疏波、接触间断以及爆轰波或爆燃波关系式，求解各波系两侧的物理量状态，包括密度、运动速度、压力以及比内能。

通过迭代法求解含爆轰波或爆燃波的黎曼解中各波系两侧的物理量状态，作为优选，选用二分法，二分法求解过程虽然比较繁杂，但具有很强的鲁棒性。

步骤9：根据含爆轰波以及爆燃波的黎曼解，构造未燃气体和已燃气体沿着界面的守恒量交换项。

对于爆轰情形，根据含爆轰波的黎曼解，计算爆轰波波速：

其中ρ_r和u_r分别表示爆轰波波前未燃气体的密度和速度；和分别表示爆轰波波后已燃气体的密度和速度。

对于未燃气体，界面项通过下式计算得到：

其中ρ₁、u₁、v₁、p₁、E₁分别表示爆轰波波前未燃气体的密度、x方向的运动速度、y方向的运动速度、压力和单位质量气体的总能；表示计算单元(i,j)内界面截断的长度，通过下式求得：

已燃气体的界面项直接取值为未燃气体界面项的相反数，由此保证沿着爆轰界面是满足质量、动量和能量守恒的。

爆燃的情形与爆轰略有不同，由于含爆燃波的黎曼解由四个波组成，爆燃波的侧边多了一个非反应波，因此在爆燃波波速已知的条件下，才能求解爆燃问题。

对于未燃气体，界面项可以通过下式计算得到：

其中ρ₀、u₀、p₀、E₀分别表示黎曼解中爆燃波波前未燃气体的密度、流体运动速度、压力和单位质量气体的总能；u₁、v₁、u_r分别表示黎曼问题的初始条件中未燃气体的x方向运动速度、y方向运动速度和法方向运动速度；D表示爆燃波波速。

已燃气体的界面项直接取值为未燃气体界面项的相反数，由此保证沿着爆燃界面是满足质量、动量和能量守恒的。

步骤10：根据修正的有限体积方法，构造控制方程组的守恒型离散格式。

在切割单元(i,j)内，对未燃气体以及已燃气体分别进行修正的有限体积离散，得到控制方程组的离散形式：

其中界面的几何参数α_i,j、根据步骤6求得；流体沿着界面的守恒量交换项X(Γ_i,j)根据步骤7-9求得；F、G分别表示计算单元边界沿着x方向和y方向的数值通量。

步骤11：结合两步四阶Lax-Wendroff方法以及守恒界面方法，实现对爆轰以及爆燃问题的数值通量高精度构造，并代入修正的有限体积离散格式，更新未燃气体以及已燃气体的守恒量。两步四阶Lax-Wendroff方法能够将四阶精度算法的时间推进步骤减少一半，从而显著提高计算效率，但是主要应用于单介质流体运动问题。因此需要应用水平集方法和虚拟流体方法，精确地捕捉多介质界面位置，并且将多介质问题转化为单介质问题。守恒界面方法是基于水平集方法和虚拟流体方法的联合算法进行发展的。水平集方法能够快捷且准确地处理界面的拓扑变化，精确地捕捉多介质界面的演化过程；虚拟流体方法能够巧妙地将多介质问题转化为单介质问题，从而有效抑制多介质界面问题的非物理振荡；通过求解含爆轰波或者爆燃波的黎曼问题，守恒界面方法能够改善爆炸反应流问题的非守恒性缺陷，提高爆轰和爆燃过程预测结果的可靠性。因此需要将两步四阶Lax-Wendroff方法与守恒界面方法结合，实现对于爆轰以及爆燃过程的高精度预测，同时显著提高高精度预测爆轰以及爆燃过程的计算效率。

利用步骤5，通过水平集方法定义未燃气体以及已燃气体所在区域；利用步骤7，通过虚拟流体方法将多介质爆轰以及爆燃问题转化为关于未燃气体以及已燃气体的单介质问题。

步骤11.1：根据计算单元内的平均守恒量以及计算单元边界处的守恒量和应用HWENO重构技术和广义黎曼解法器得到变量的时间导数和

步骤11.2：计算半个时间步长后的平均守恒量以及计算单元边界处的守恒量和再次应用HWENO重构技术和广义黎曼解法器得到变量的时间导数和

步骤11.3：根据数值通量的时间导数，构造计算单元边界的数值通量，并代入有限体积离散格式，更新未燃气体以及已燃气体的守恒量。

步骤12：根据时间步长Δt判断第n步计算时间tⁿ与终止时刻T的关系，若tⁿ+Δt≤T，则tⁿ⁺¹＝tⁿ+Δt，返回步骤3，更新下一时刻的守恒量Uⁿ⁺¹；否则Δt＝T-tⁿ，返回步骤3，求得终止时刻的守恒量，并且根据状态方程计算终止时刻的质量、x方向的运动速度、y方向的运动速度以及压力，循环结束。

输出终止时刻所有计算单元的水平集函数、质量、x方向的运动速度、y方向的运动速度以及压力，预测爆轰或者爆燃界面的位置及未燃气体和已燃气体所在区域，以及计算区域内流体的质量、x方向的运动速度、y方向的运动速度以及压力。

还包括步骤13：利用步骤1至步骤12进行爆轰和爆燃问题的预测，提高高精度预测方法的计算效率；同时有效解决非守恒问题，实现对爆轰和爆燃过程的可靠预测，进而解决爆炸力学领域相关工程技术问题。

所述爆炸力学领域相关工程技术问题包括武器弹药设计与毁伤评估以及爆炸事故的防灾减灾。

有益效果：

1、本发明公开的一种基于广义黎曼解法器的爆炸反应流高精度预测方法，应用基于广义黎曼解法器的两步四阶Lax-Wendroff方法求解计算单元边界的数值通量，两步四阶Lax-Wendroff方法融合了Runge-Kutta算法的简单性特点以及广义黎曼解法器的时空耦合性质，因此该预测方法能够提高预测爆轰和爆燃过程的计算效率和稳定性。

2、本发明公开的一种基于广义黎曼解法器的爆炸反应流高精度预测方法，应用守恒界面方法得到守恒型计算格式，采用修正的有限体积方法离散控制方程，并且通过求解含爆轰波或者爆燃波的黎曼问题求解守恒量交换项，能够使得沿着界面是满足质量、动量以及能量守恒的。因此该预测方法能够克服传统水平集/虚拟流体方法缺乏守恒性的缺陷，提高预测爆轰和爆燃过程的准确性。

3、本发明公开的一种基于广义黎曼解法器的爆炸反应流高精度预测方法，应用水平集方法追踪爆轰或者爆燃界面，水平集方法能够高效快捷地实现高阶精度，同时能够简单有效地拓展到高维空间。因此，该预测方法能够快速且精确地捕捉爆轰或爆燃界面的拓扑变化过程。

4、本发明公开的一种基于广义黎曼解法器的爆炸反应流高精度预测方法，应用虚拟流体方法构造虚拟单元内未燃气体和已燃气体的状态，能够巧妙地将多介质问题转化为单介质问题。因此，本发明能够有效抑制多介质界面附近的非物理振荡。