基于Level-Set的多介质流体动力学数值分析

　　目前大多数商业软件如Ls-dyna等都是一阶或者二阶精度,程序库里也不能满足实际问题中需要的状态方程及本构方程,而且程序源代码不公开, 因此,针对各自学科需要开发适用的程序是当务之急.关于多介质流体动力学的数值计算受到了国内外研究者的重视,已有许多很成熟的计算方法不断出现 [1-2].水下爆炸的计算涉及到多介质高压比高密度比,特别是对界面的处理,一直都是高精度数值计算的难点.作者在国内外一些研究分析方法基础上,提出 了处理这类问题的方便有效的方法,并且对多介质流体动力学进行了编程.

　　作者阐述了程序内部采用的计算方法及计算步骤,通过一些典型的算例对程序进行了数值模拟,计算结果很好地反映了其物理本质.对1 kg的TNT炸药水下爆炸问题的参数进行了求解,通过与经验公式的计算结果的对比表明,该程序能够满足工程的实际需要,可以用来解决水下爆炸的能量输出特 性问题.

　　1　数值计算程序

　　1·1　直角坐标下的守恒型欧拉控制方程

　　令ρ,p和e分别表示密度、压力和单位质量的热力学能. u和v分别表示x和y方向的速度分量,E=e+(u2+v2)/2和H=E+p/ρ分别表示单位质量的总能和总焓,二维直角坐标下的Euler方程可简写成Ut+F(U)x+G(U)y=0,即

　　1·2　状态方程

　　①水的状态方程

　　采用刚性状态方程

　　式中:γ=7·15;p0为大气压力;β=3·31×108Pa;ρ0为水的初始密度.

　　定义热力学能

　　②爆轰产物状态方程

　　采用JWL状态方程

　　式中:A,B,R1,R2,ω为与炸药相关的常数;ρe为装药密度。

　　1·3　计算方法

　　①方程的差分离散

　　采用Steger-Warming流通量分裂方法对流通量进行分裂,给定p=p(ρ,e),ξ= F/ U的特征值矩阵为Λ=diag(λ1,λ2,λ3,λ4)=diag(u-c,u,u,u+c),相对应的左右特征向量分别为Lk和Rk.因此,

　　对Fx和Gy在空间上采用张涵信提出的三阶精度的ENN有限差分格式进行差分离散[3],而在时间上采用三阶精度的Runge-Kutta法进行差分离散。

　　②内边界处理

　　采用Level-Set对内边界进行追踪,然后利用改进的GFM对边界进行处理[4-7].

　　Level-Set方程为

　　其中为距离函数,满足>0,流体1=0,边界上<0,流体2.Level-Set离散采用了五阶的weno格式.

　　重新初始化方程为

　　为了确定界面接触间断两侧流体的状态,需要沿着界面的法向求解一维Riemann问题,界面的法线方向N由Level-Set函数给出,用N= /| |表示单位法向量.沿着法向定义每种流体的界面边界条件,在界面附近的带状区域内求解Riemann问题.通过定常外插方程It±N I=0把流体1的状态从界面的一侧外推到另一侧,这样在界面附近带状区域内的网格点上存在真实值和外推值.通过这些真实值和外推值在这些网格点的法线方向 上构造Riemann问题.如图1所示,求解UIL=Unext,UIR=Unb的Riemann问题得出ρIL,ρIR,uI,pI.这样就可以得到流 体1在b点的ghost流体状态为ρb=ρIR,pb=pI,ub=(uIN)N+ua-(ubN)N.