基于伪弧长移动网格算法的爆炸与冲击多介质问(2)
本文首先简单介绍了五方程模型及Mie-Grüneisen形式状态方程。然后介绍了控制方程演化算法、PALM算法和人工界面压缩技术,并总结了多介质伪弧长移动网格(MM-PALM)算法的计算步骤。通过一维激波管问题和二维水下爆炸问题的数值算例表明(MM-PALM)算法是可行的,能够精确捕捉爆炸冲击波阵面和精确描述多物质界面。
1 五方程模型
1.1 五方程模型
对于高速、非黏、可压缩两相流,可以采用如下五方程模型进行描述:
式中:t为时间变量;W为物理量,W=(z1ρ1,z2ρ2,ρu,ρE,z1)T,z1、z2为各介质体积分数,对于两种介质的体积分数,z1∈[0,1],z2=1-z1,ρ1、ρ2为各介质的密度,ρ为流体的混合密度,u为速度矢量,E为单位质量流体所具有的总能量;H(W)为通量,H(W)=(z1ρ1u,z2ρ2u,ρu?u+pI,(ρE+p)u,z1u)T,p为流体的压力,I为单位张量矩阵;S(W)为源项,
对模型(1)式,可通过增加额外的输送方程和质量守恒方程来扩展到多种介质。流体的总密度和总能量由各介质的密度ρl和内能el(l=1,2)计算得到:
1.2 流体状态方程
将可压缩多介质流体中各介质的状态方程写成统一的Mie-Grüneisen形式:
式中:Γ(ρ)、eref(ρ)、pref(ρ)为关于密度ρ的函数,Γ(ρ)为Grüneisen系数,eref(ρ)和pref(ρ)为通过实验得到的材料内能和压力值。参数Γ(ρ)、eref(ρ)、pref(ρ)通过不同的取值,(4)式可代表多种不同的材料,可以是气体、液体或者固体等。根据不同材料的Mie-Grüneisen状态方程,五方程模型(1)式可以对不同物理问题进行建模。
为更方便地表示Mie-Grüneisen状态方程,将状态方程(4)式改写为
式中:
进一步根据等压假设,即各介质在界面两边的压力相等(p=pl),可得混合流体压力为
式中:ξ混合参数,
混合声速对于柯朗-弗里德里希斯-列维(CFL)条件十分重要,经典力学中对声速的定义如下:
对Mie-Grüneisen形式的状态方程(4)式,各介质的声速(10)式可写成如下形式:
式中:Γ′l、(pref)′l、(eref)′l是状态方程(4)式关于ρ的导数。因此定义流体的平均混合声速为
式中:yl为质量分数,
常用的两种Mie-Grüneisen形式的状态方程如下:
1)理想气体的状态方程为
式中:γ为比热比。对应于(4)式中的参数Γ(ρ)=γ-1,pref(ρ)=0,eref(ρ)=0.
2)JWL状态方程为
式中:ρ0、A1、A2、R1、R2、e0均为材料参数。
2 数值算法
2.1 控制方程演化
2.1.1 空间离散格式
本文采用有限体积的方法,空间离散采用2阶MUSCL方案,时间离散采用3阶TVD Runge-Kutta方法。将模型(1)式改写为如下形式:
式中:F(W)和G(W)分别为x轴、y轴方向的通量。设网格Ki,j(1≤i≤Nx,1≤j≤Ny)上的物理量均值为为x轴、y轴方向的网格总数),将(15)式在Ki,j上积分,得到下式:
式中:dσ为面积微元。利用Green-Gauss定理,可以将(16)式转化为
式中:|Ki,j|表示网格Ki,j的面积;?表示网格Ki,j的第k条边;ni,j为边?Ki,j单位外法向量;ds为边界长度微元;δ(W)=(F(W),G(W))。定义Lax-Friedrichs通量格式为
式中:Wl、Wr分别为单元格左右两边的物理量值;n为单元格的单位外法向量;c0=max(δ′(Wl)·n,δ′(Wr)·n),δ′(Wl)、δ′(Wr)分别为δ(W)对Wl、Wr求导数。则(17)式可以写成如下半离散格式:
式中:和为网格Ki,j边界的内外逼近值。为了获得2阶空间精度,在网格单元的边界处进行2阶重构,如图1所示。
图1 二维网格单元计算示意图Fig.1 Calculation diagram of 2D mesh unit
2阶重构形式如下:
式中:ψ(*,*)为非线性限制函数,当ψ(*,*)≡0时,半离散系统(19)式在空间上为1阶精度。为实现2阶精度,在计算中利用van Leer’s限制器进行2阶重构:
式中:a、b为任意值;ε为常数,ε=10-16.
2.1.2 时间离散格式
时间离散采用3阶TVD Runge-Kutta格式:
式中:为每一子步插值;w(n)、w(n+1)分别为tn、tn+1时刻的物理量值;L(*)由(19)式求得。
2.2 PALM算法
2.2.1 网格点移动
令x=(x,y)和ζ=(ξ,η)分别为物理空间Ωp和计算空间Ωc的坐标,〈xi-1,j-1,xi-1,j,xi,j,xi,j-1〉表示网格单元Ki,j的4个顶点。则物理空间坐标与计算空间坐标的一一对应关系为(x,y)=(x(ξ,η),y(ξ,η))。通过变分方法可知,其最小能量泛函为
文章来源:《爆炸与冲击》 网址: http://www.bzycjzz.cn/qikandaodu/2021/0709/1295.html
上一篇:冲击载荷下螺栓预紧力对应力波影响分析
下一篇:三维多物质欧拉界面处理的并行算法研究