基于伪弧长移动网格算法的爆炸与冲击多介质问(3)
式中:G1和G2为给定的正定对称矩阵;则相应的Euler-Lagrange方程为
将G1和G2都取为MI,得到:
式中:M为弧长自适应控制函数,
参数α1、α2是非负常数,w为W的子集。在计算中,使用Gauss-Seidel迭代方法对(24)式进行求解:
式中:1≤i≤Nξ,1≤j≤Nη,Nξ、Nη分别为伪弧长空间ξ轴、η轴方向上的网格数;上标[z]和[z+1]表示旧网格与新网格;z=0,1,2,….Gauss-Seidel迭代过程如下:
在计算过程中,为保持计算网格的平滑性,通常会对函数M做光滑性处理:
2.2.2 物理量更新
进一步考虑新网格x[z+1]和旧网格x[z]之间物理量的更新方案。新旧网格之间物理量的重构要满足物理量的守恒性:
令Dk为经过一次迭代后边所扫过的面积大小(见图2),进而得到如下守恒插值重构公式:
图2 新旧网格变化示意图Fig.2 Transform of old grid to new mesh
式中:为物理量W在区域Dk的积分,k=1,2,3,4.
对(31)式采用(32)式近似求解:
关于面积Dk,以D1为例进行计算:
当为逆时针排列时,D1为正;当为顺时针排列时,D1为负。
2.3 人工界面压缩技术
2.3.1 体积分数重构
以一维情形为例,介绍人工界面压缩技术。对于多维情形,可以采用维分裂的思想逐维计算。
定义介质1在网格单元上的体积分数为
式中:z1,i(x)用如下双曲正切函数[20]表示:
下面定义界面单元。当网格单元的体积分数满足以下两个条件时该网格单元即为界面单元:
将zi(x)代入(34)式,可得体积分数的左右重构值为
式中:β为常数(本文取β=1.6);
2.3.2 其他物理量重构
下面对界面单元处的物理量进行重构。假设在重构时间步内网格单元的各介质密度ρ1、ρ2保持不变,从而有如下重构公式:
1)质量
式中:ωl为权重系数,本文计算中取为常数。
2)动量
3)能量
式中:
综合以上论述,将MM-PALM算法计算步骤总结如下:
1)初始化。在计算空间与物理空间中给出初始网格分布计算求得网格物理量
2)根据Gauss-Seidel迭代公式(28)式,由旧网格坐标获取新的网格坐标
3)通过守恒差值公式(31)式,计算新网格的物理量
4)利用(20)式进行2阶MUSCL重构。
5)根据界面单元判断条件(37)式,确定界面单元网格。
6)对各维度界面单元网格分别进行重构。
7)利用(18)式计算数值通量。
8)利用Runge-Kutta公式(22)式计算tn+1时刻的物理量
9)如果tn+1≤te(te为终止时间),则令重复步骤2.
3 数值算例
3.1 一维激波管问题
本文算例是经典的Air-helium激波管问题[21],介质1为空气(l=1),介质2为氦气(l=2)。计算区域为x∈[0,1]并划分为200个计算网格。状态方程采用理想气体状态方程,其中参数γ1=1.4,γ2=1.667.本算例中所有参数均为无量纲,计算初始条件为
对PALM算法和MM-PALM算法采用不同的伪弧长自适应函数。PALM算法采用如下控制函数:
NM-PALM算法采用如下控制函数:
MM-PALM参照解为2 000个计算网格的计算结果。图3为t=0.2时刻的计算结果图,图4为相应的冲击波区域局部放大图。从图3和图4中可以看出,在冲击波区域,PALM算法与多介质伪弧长算法无明显差别,计算结果都好于固定网格计算结果。图5为体积分数随时间的变化曲线。从图3(d)和图5中可以看出,PALM算法较固定网格算法有更好的界面分辨率,但与MM-PALM算法相比,其界面分辨率仍有非常大的差距。从图5中还可以看出,随着计算时间的增加,固定网格算法与PALM算法的界面区域都变得越来越模糊,而MM-PALM算法则将界面单元保持在5个网格之内。图6所示为网格自适应变化图。从图6中可以看出:PALM算法与MM-PALM算法的网格都在冲击波区域加密,计算结果都好于固定网格算法;只有PALM算法在多物质界面区域计算网格加密,因此PALM算法的界面分辨率要好于固定网格。尽管PALM算法在界面单元区域进行了网格加密,但与加入了人工界面压缩技术的MM-PALM算法相比,仍旧难以保持精确的界面分辨率。由此算例可以看出,MM-PALM算法不仅能够很好地处理冲击波区域的强间断问题,更重要的是可以保持多介质之间的尖锐界面。
文章来源:《爆炸与冲击》 网址: http://www.bzycjzz.cn/qikandaodu/2021/0709/1295.html
上一篇:冲击载荷下螺栓预紧力对应力波影响分析
下一篇:三维多物质欧拉界面处理的并行算法研究