[发明专利]一种近场水下爆炸冲击波载荷的高精度数值模拟方法

摘要

本项发明提出了一种近场水下爆炸冲击波载荷的高精度数值模拟方法，属于水下爆炸数值模拟技术领域。本发明采用间断伽辽金方法(RKDG)和三阶TVD Runge‑Kutta方法对计及强可压缩性的非线性欧拉方程进行高精度数值离散；采用虚拟流方法(GFM)对气液界面两侧的物理量进行预处理，有效削减由物理量间断导致的非物理震荡；采用水平集方法(LSM)处理物质界面上发生的复杂拓扑结构变化，实现对运动界面的精准捕捉。本发明提出的方法能够有效模拟近场水下爆炸冲击波载荷的产生及传播过程；数值模拟结果与经典经验公式、实验结果能够较好地吻合，验证了该方法在处理强可压缩、强间断、瞬态强非线性问题上具有一定优势。

主权项

1.一种近场水下爆炸载荷的高精度数值模拟方法，其特征在于：其操作过程为：步骤1：确定水下爆炸的计算区域，建立直角坐标系，将该区域划分为m×n个网格(其中，m表示x方向网格数量，n表示y方向网格数量)；步骤2：定义计算区域内初始状态下的各物理量、相关函数，并建立流体控制方程及物质状态方程，包括：(a)各物质的界面函数，即Level Set函数：其中，Γ(0)表示零等值面，即物质界面，d(x,y,Γ(0))表示计算域内一点(x,y)到物质界面的距离，Ω1、Ω2分别为物质界面两侧的区域；(b)各物质的物理参数：初始状态下炸药的密度、速度、压力，初始状态下水的密度、速度、压力；(c)流体控制方程：根据水下爆炸场的物理特性，建立相适应的无粘、无旋、计及可压缩性的流体控制方程：Ut+▽·F(U)＝S(U)      (2)其中，F(U)＝[f(U),g(U)]，(d)水的状态方程(Tait方程)及相关参数：其中，N、AW、BW是与水的物理性质相关的常数，ρW0表示水在初始状态下的密度，pW、ρW分别表示水的压力和密度；(e)炸药爆轰产物状态方程(JWL方程)及相关参数：其中，AE、BE、R1、R2、ω是与炸药性质相关的常数，ρE0表示炸药在初始状态下的密度，ρE表示爆轰产物的密度，PE表示爆轰产物的密度，E0表示爆轰产物的比内能；步骤3：拓展计算区域，根据不同工况设置边界条件：由于边界处缺少差分格式所必需的网格点，因此要将计算区域向外拓展N层；根据本方法采用的格式，至少需要向外拓展3层网格(即N≥3)；对于自由边界条件，有：UΓ＝U‑1,U‑i＝U‑i+1,i＝2···N，其中，UΓ表示计算区域边界网格上的任意物理量，U‑i表示虚拟网格上的任意物理量；具体而言，计算区域边界外第一层虚拟网格上的物理量等于计算区域边界网格上的物理量；计算区域边界外第二层虚拟网格上的物理量等于第一层虚拟网格上的物理量；计算区域边界外第三层虚拟网格上的物理量等于第二层虚拟网格上的物理量，以此类推；对于固壁边界条件，有：其中，U_Γ表示计算区域边界网格上除速度以外的任意物理量，U_‑i表示虚拟网格上除速度以外的任意物理量，V_Γ表示计算区域边界网格上的速度，V_‑i表示虚拟网格上的速度；具体而言，计算区域边界外第一层虚拟网格上垂直于边界的速度与边界网格上垂直于边界的速度矢量和为零，其他物理量相等；计算区域边界外第二层虚拟网格上垂直于边界的速度与第一层虚拟网格上垂直于边界的速度矢量和为零，其他物理量相等；计算区域边界外第三层虚拟网格上垂直于边界的速度与第二层虚拟网格上垂直于边界的速度矢量和为零，其他物理量相等，以此类推；步骤4：选取CFL参数并计算时间步长：CFL参数为介于0到1之间的常数；当前计算所用的总时间用t＝nΔt表示，n为计算步数，时间步长Δt可由下式得到：Δx为计算区域中x方向的网格步长；Δy为计算区域中y方向的网格步长；u为计算区域中当前网格在x方向的速度分量；v为计算区域中当前网格在y方向的速度分量；c为计算区域中当前网格处的音速；步骤5：为计算区域及虚拟网格上的参量赋值：当t＝0时，将步骤一中给出的计算区域内各物质的物理参数及界面函数值，作为计算区域内各网格的参量初值；当t＞0时，将步骤13中计算得到的计算区域内各物质的物理参量及界面函数值，作为计算区域内各网格参量在下一步计算时的初值；对于虚拟网格，则根据步骤3中确立的方法，通过插值外推对其进行赋值；步骤6：求解物质界面法向：根据Level Set方程，对计算区域内各物质界面的法线方向进行求解，求解公式为：步骤7：求解黎曼问题并对物质界面附近网格赋值：沿物质界面法线方向建立炸药与水之间的黎曼问题并求解出炸药与水之间物质界面处的物理量；用求解黎曼问题得到的物理量，对炸药区域内紧邻物质界面的一层网格进行赋值；同时，通过外推对炸药区域以外的网格进行赋值，具体方法为：在处理物质界面炸药区域一侧的流场时，将物质界面另一侧的网格设置为虚拟流体网格，这些虚拟流体网格上的压力与法向速度与该侧真实单元的一致，而切向速度以及密度则通过炸药区域的真实流体跨物质界面插值外推求得；步骤8：通过RKDG方法离散欧拉方程，求解冲击波超压的数值通量，具体为：对(2)式两边同乘以试探函数Φ(x,y)，并在网格单元K上积分得：对(7)式进行分部积分后，利用高斯求积公式求解式中的线积分和面积分；由于网格单元边界处冲击波超压通量不连续，故采用数值通量代替实际通量用近似解代替精确解U，用数值近似后的试探函数代替试探函数Φ，由此可得：其中，ω表示权重系数，|e|表示积分区域网格单元边界的长度，|K|表示积分区域网格单元的面积；对于(8)式中的数值通量采用(9)式进行求解：其中，α表示Jacobi矩阵的特征值最大的最大绝对值，U⁺(x_el,y_el,t)表示从正方向上逼近网格单元边界的U值，U^‑(x_el,y_el,t)表示从负方向上逼近网格单元边界的U值；对于炸药区域爆轰产物的压力，由步骤一中的(4)式得到；对于(8)式中的近似解采用(10)式进行求解：其中，对于网格单元[x_i‑1/2,x_i+1/2]×[y_j‑1/2,y_j+1/2]，Δx_i＝x_i+1/2‑x_i‑1/2，Δy_j＝y_j+1/2‑y_j‑1/2；将(9)、(10)两式带入(8)式可以得到(2)式的半离散格式，记为：Ut＝R U     (11)其中，R U表示空间导数项的离散算子；步骤9：通过三阶TVD Runge‑Kutta法对(11)式等式左边的时间项进行离散，得到(2)式的完全离散格式，进而求解出下一时刻计算区域内炸药爆轰产物的各物理量；具体公式为：步骤10：与步骤7进行类似操作，再次求解黎曼问题并对物质界面附近网格赋值：沿物质界面法线方向建立炸药与水之间的黎曼问题并求解出炸药与水之间物质界面处的物理量；用求解黎曼问题得到的物理量，对空气区域内紧邻物质界面的一层网格进行赋值；同时，通过外推对空气区域以外的网格进行赋值，具体方法为：在处理物质界面空气区域一侧的流场时，将物质界面另一侧的网格设置为虚拟流体网格，这些虚拟流体网格上的压力与法向速度与该侧真实单元的一致，而切向速度以及密度则通过炸药区域的真实流体跨物质界面插值外推求得；步骤11：在步骤10的基础上，进行与步骤8、步骤9类似的操作；但不同之处在于：步骤8中，对于水区域的压力，由步骤一中的(3)式得到；步骤9中求解得到的是下一时刻计算区域内水的各物理量；步骤12：对Level Set运动方程的空间导数项进行离散，可得：对于炸药与水的界面函数的横向数值通量，可由(14)式求解：其中，表示计算区域内网格单元上炸药与水的界面函数横向上的数值通量，表示该网格单元上处于横向正方向的界面函数数值通量，表示该网格单元上处于横向负方向的界面函数数值通量；表示计算区域内网格单元(i‑r+s,j)上炸药和水的界面函数的横向正方向的通量；表示计算区域内网格单元(i‑r+s,j)上炸药和水的界面函数的横向负方向的通量；具体参数可通过(15)、(16)式求解：其中，表示计算区域内网格单元(i,j)上炸药和水的界面函数值；对于炸药与水的界面函数的纵向数值通量，求解方法相同；经过步骤12，可以得到(13)式的半离散格式，记为：其中，表示空间导数项的离散算子：步骤13：通过三阶TVD Runge‑Kutta方法离散Level Set运动方程的时间导数项，得到(13)式的完全离散格式，进而求解出下一时刻计算区域内各网格单元的界面函数值；步骤14：在之前步骤的基础上，求解并输出计算区域内所有网格单元的压力值：若网格单元处于炸药爆轰产物区域内，则通过步骤一中的(4)式计算求得；若网格单元处于水区域内，则通过步骤一中的(3)式求得；步骤15：判断当前计算所用的总时间t＝nΔt与设定的计算总时间t*的关系：若t＞t*，则结束计算，输出计算区域内网格单元的压力值；若t＞t*，则返回至步骤3，继续运算。经过上述步骤的操作，即可得到水下爆炸计算区域内每一时刻各网格单元处的冲击波压力值。