AAM  >> Vol. 8 No. 2 (February 2019)

    Mie-Grüneisen多介质混合计算模型的Godunov型格式
    The Godunov Type Scheme for Multi-Component Mie-Grüneisen Mixture Model

  • 全文下载: PDF(1138KB)    PP.350-356   DOI: 10.12677/AAM.2019.82039  
  • 下载量: 157  浏览量: 279   国家自然科学基金支持

作者:  

吴宗铎,严 谨,董智惠:广东海洋大学,海洋工程学院,广东 湛江;
张 建:广东省航道事务中心,广东 广州;
刘柄宏:上海海洋大学,水产与生命学院,上海

关键词:
Godunov类型格式多介质Mie-Grüneisen方程Godunov Type Scheme Multi-Component Mie-Grüneisen Equation

摘要:

本文的主要目的是为了推导Mie-Grüneisen多介质混合模型下的Godunov类型格式。Godunov格式在求解多介质问题的间断面时,通常能得到使数值解保持无振荡。为了适应Mie-Grüneisen多介质混合模型,将二维矩形的Godunov类型格式在有限体积算法下进行了重新构造。数值实验表明,Godunov格式在Mie-Grüneisen混合模型下能很好地相结合。

The aim of this article is the deduction of Godunov type scheme for Mie-Grüneisen mixture model. For multi-component problem with discontinuous sections, the Godunov schemes always keep the solution non-oscillation. To adapt to Mie-Grüneisen mixture model, the Godunov fluxes are recon-structed under finite volume method in 2D rectangular grids. According to numerical test, the Go-dunov scheme can couple well with Mie-Grüneisen mixture model.

1. 引言

可压缩流场中的多介质Riemann问题,多年以来一直是研究者们所关注的一个焦点问题 [1] 。流体的运动特性可由欧拉方程描述,而流体的物理性质则由状态方程来描述。相比于其它形式的状态方程,Mie-Grüneisen方式无论是对流体 [2] 还是固体 [3] ,均能很好地描述其力学和热力学特性,因而得到广泛的应用。而薛克明提出的Mie-Grüneisen多介质混合模型 [4] ,则为Mie-Grüneisen方程下的多介质问题提供了很大的便利 [5] 。但是由于Mie-Grüneisen多介质混合模型目前并未得到广泛的重视,这使得它可供参考的数值离散格式仍然比较少 [6] 。而Godunov型格式则是计算流体力学有限体积方法中比较常见的一种格式,因为它能给带有强间断面的问题提供稳定的数值解 [7] 。相比于GFM等数值算法,它省去了其它的数值处理界面的过程 [8] ,仅自身的数值格式便具备较好的分辨率。对于Mie-Grüneisen多介质混合模型来说,其方程系统同时包括Euler守恒方程以及其它热力学参数的非守恒方程,不仅涉及的方程式较为复杂,而且具有很强的非守恒性,Godunov型格式的构造形式将显得非常复杂。

本文在Mie-Grüneisen多介质混合模型下,按照二维有限体积差分方法进行了数值离散化,离散后确定了守恒变量与通量的计算形式,并按特定的时间步长和均匀的空间步长进行数值迭代。该过程中,对Godunov格式下的通量做了细致的推导。为了验证该格式的计算稳定性,利用该模型模拟了一个复杂的二维流场的相互作用问题。

2. 计算方法

2.1. 控制方程组

这里考虑由多种流场介质共存的二维可压缩流场,那么整个流场的运动可以由欧拉方程来描述:

{ ρ t + x ( ρ u ) + y ( ρ v ) = 0 ( ρ u ) t + x ( ρ u 2 + p ) + y ( ρ u v ) = 0 ( ρ v ) t + x ( ρ u v ) + y ( ρ v 2 + p ) = 0 ( ρ E ) t + x ( ρ E u + p u ) + y ( ρ E v + p v ) = 0 (1)

式中,r为密度,u、v分别为横向和纵向速度,E为能量,p为压力。欧拉方程(1)需要结合状态方程才能完成求解。这里使用Mie-Grüneisen状态方程:

p p r e f ( ρ ) = ρ Γ ( ρ ) ( e e r e f ( ρ ) ) (2)

这里,G、pref、eref分别是由流体物理特性决定的热力学参数 [4] 。由于G、pref、eref通常都具有复杂的函数表达式,可以考虑利用非守恒型方程来对它们进行分别求解 [5] :

{ t ( 1 Γ ) + u x ( 1 Γ ) + v y ( 1 Γ ) + ρ ϕ u x + ρ ϕ v y = 0 t ( p r e f Γ ) + u x ( p r e f Γ ) + v y ( p r e f Γ ) + ρ φ u x + ρ φ v y = 0 t ( ρ e r e f ) + u x ( ρ e r e f ) + v y ( ρ e r e f ) + ρ ψ u x + ρ ψ v y = 0 (3)

其中f、j、y代表1/G、pref/G、reref对密度r的导数项。

对于多介质流场,需要对界面进行捕捉。假设质量分数为YI,界面的运动可以通过以下方式来求解:

( ρ Y I ) t + ( ρ Y I u ) x + ( ρ Y I v ) y = 0 , i = 1 , 2 , , m 1 (4)

将守恒方程(1)、非守恒方程(3)和界面运动方程(4)联立起来,即可得到整个Mie-Grüneisen多介质混合模型的求解方程组。

2.2. 有限体积算法

Mie-Grüneisen多介质混合模型的有限体积格式可以表示成如下所示的形式:

W j n + 1 W j n Δ t + F i + 1 / 2 , j F i 1 / 2 , j Δ x + G i , j + 1 / 2 G i , j 1 / 2 Δ y = 0 (5)

这里Δt为时间步长,Δx、Δy为空间步长。W代表变量,分别用F、G表示x、y方向的通量。角标i和j分别代表x、y方向的节点序号,n为时间步长的序号。变量和通量各自的表达式分别为:

W = { ρ ρ u ρ v ρ E 1 / Γ p r e f / Γ ρ p r e f ρ Y I } F = { ρ u ρ u 2 + p ρ u v ( ρ E + p ) u ρ ϕ u ρ φ u ρ ψ u ρ Y I u } + u { 0 0 0 0 1 / Γ p r e f / Γ ρ p r e f 0 } G = { ρ u ρ u v ρ v 2 + p ( ρ E + p ) v ρ ϕ v ρ φ v ρ ψ v ρ Y I v } + v { 0 0 0 0 1 / Γ p r e f / Γ ρ p r e f 0 }

Godunov型的通量F和G,其表达形式如下 [7] :

{ F i + 1 / 2 , j = 1 2 [ F L i + 1 / 2 , j + F R i + 1 / 2 , j R ^ | Λ ^ | L ^ ( W L i + 1 / 2 , j W R i + 1 / 2 , j ) ] G i , j + 1 / 2 = 1 2 [ G L i , j + 1 / 2 + G R i , j + 1 / 2 R ^ | Λ ^ | L ^ ( W L i , j + 1 / 2 W R i , j + 1 / 2 ) ] (6)

这里,WL、WR代表通量所在的网格左右两侧的变量,FL、FR和GL、GR代表两侧变量形成的x、y方向通量。 L ^ R ^ Λ ^ 则代表了Jacobi矩阵的左右特征向量和特征值对角矩阵的平均值。Jacobi矩阵为通量对变量的偏导数,即 F / W G / W

2.3. 通量相关矩阵

对于左右两侧的变量及通量,表达形式如下:

W K = { ρ K ρ K u K ρ K v K ρ K E K ( 1 / Γ ) K ( p r e f / Γ ) K ( ρ p r e f ) K ( ρ Y I ) K } F K G K = { ( ρ θ ) K ( ρ θ ) K u K + p K ζ x ( ρ θ ) K v K + p K ζ y ( ρ K E K + p K ) θ K θ ( 1 / Γ ) K + ρ ϕ θ K θ ( p r e f / Γ ) K + ρ φ θ K θ ( ρ p r e f ) K + ρ ψ θ K ( ρ Y I ) K θ K } ( K = L , R ) (7)

这里,ζx和ζy的数值分别定义为:

ζ x = { 1 F 0 G , ζ y = { 0 F 1 G

其中 θ = u ζ x + v ζ y 。(7)中,未注明角标的参数θ、r、f、j、y则代表节点i,j处的数值。这样Jacobi矩阵A可以通过(7)式求偏导数得到:

A = [ ζ x ζ y K ζ x u θ ( 1 Γ ) u ζ x + θ u ζ y Γ v ζ x Γ ζ x p Γ ζ x Γ ζ x Γ ζ x K ζ y v θ v ζ x Γ u ζ y ( 1 Γ ) u ζ y + θ Γ ζ y p Γ ζ y Γ ζ y Γ ζ y ( K H ) θ H ζ x Γ u θ H ζ y Γ v θ ( Γ + 1 ) θ p Γ θ Γ θ Γ θ ϕ θ ϕ θ φ θ φ θ ψ θ ψ θ Y I θ Y I θ ] (8)

这里,H为熵: H = E + p / ρ 。且有 K = Γ ( u 2 + v 2 ) / 2 。由矩阵A可以用偏微分关系求得特征值对角矩和左右特征向量矩阵:

Λ = [ θ c θ θ θ + c θ θ θ θ ]

这里c代表声速,计算式为:

c = Γ ( H 1 2 ( u 2 + v 2 ) p ϕ + φ ψ ) (9)

L = Γ 2 c 2 [ θ ξ + u 2 + v 2 2 u ζ x ξ v ζ y ξ 1 p 1 1 u 2 v 2 + κ 2 u 2 v 2 2 p 2 2 ( ζ y u ζ x v ) κ ζ y κ ζ x κ θ ξ + u 2 + v 2 2 u + ζ x ξ v + ζ y ξ 1 p 1 1 ϕ ( u 2 + v 2 ) 2 ϕ u 2 ϕ v 2 ϕ 2 ϕ p + κ 2 ϕ 2 ϕ φ ( u 2 + v 2 ) 2 φ u 2 φ v 2 φ 2 φ 2 φ + κ φ 2 φ ψ ( u 2 + v 2 ) 2 ψ u 2 ψ v 2 ψ 2 ψ 2 ψ 2 ψ + κ ψ Y ( u 2 + v 2 ) 2 Y u 2 Y v 2 Y 2 Y 2 Y 2 Y κ ]

R = [ 1 1 0 1 u ζ x c u ζ y u + ζ x c v ζ y c v ζ x v + ζ y c H θ c 1 2 ( u 2 + v 2 ) ζ y u + ζ x v H + θ c p 1 1 ϕ ϕ 1 φ φ 1 ψ ψ 1 Y Y 1 ]

这里 k = 2 c 2 / Γ ξ = c / Γ

2.4. Roe平均值计算

在方程组(6)中, L ^ R ^ Λ ^ 中的参数,需要根据左右两侧的状态取平均值。这里采用最为常用的一种Roe平均值,计算方式如下:

u ^ = ρ L u L + ρ R u R ρ L + ρ R , v ^ = ρ L v L + ρ R v R ρ L + ρ R , H ^ = ρ L H L + ρ R H R ρ L + ρ R

原始的Roe格式未提供p和G的处理方法,这里采用如下计算方式:

1 Γ ^ = ρ L ( 1 / Γ L ) + ρ R ( 1 / Γ R ) ρ L + ρ R , p ^ Γ ^ = ρ L ( p L / Γ L ) + ρ R ( p R / Γ R ) ρ L + ρ R

参数f、j、y为节点处的数值,因此不需要对网格左右侧取平均值。声速c可根据(9)式完成求解。

3. 数值算例

这里再考虑一个钼和液态MORB(Mid-Ocean Ridge Basalt)相互作用的问题。初始时刻,整个流场计算域长宽都为1 m。MORB分布在[0.4, 0.6] × [0, 0.5]的一块方形区域,其余区域均为钼。左侧[0, 0.3 m]区域内的钼为高压态,除此以外都为正常态的钼。钼和液态MORB的物理特性均由Mie-Grüneisen状态方程(2)来表示,其中的参数为:

Γ = Γ 0 ( V V 0 ) q , p r e f = p 0 + c 0 2 ( V 0 V ) [ V 0 s ( V 0 V ) ] 2 , e r e f = e 0 + 1 2 ( p r e f + p 0 ) ( V 0 V )

这里V为相对体积,V = 1/r。流场中各介质的物理状态为:

{ MORB : ρ = 2.26 × 10 3 kg / m 3 , p = 1 × 10 5 pa , u = 0 , : ρ = 9.961 × 10 3 kg / m 3 , p = 1 × 10 5 pa , u = 0 : ρ = 11.042 × 10 3 kg / m 3 , p = 3 × 10 10 pa , u = 543 m / s MORB : ρ 0 = 2.26 × 10 3 kg / m 3 , c 0 = 2100 m / s , s = 1.68 , Γ 0 = 1.56 , q = 1.0 : ρ 0 = 9.96 × 10 3 kg / m 3 , c 0 = 4770 m / s , s = 1.43 , Γ 0 = 0.18 , q = 1.0

Figure 1. Density contours of Godunov scheme (left) and HLLC scheme (right)

图1. Godunov格式(左)和HLLC格式(右)的密度云图

Figure 2. Pressure contours of Godunov scheme (left) and HLLC scheme (right)

图2. Godunov格式(左)和HLLC格式(右)的压力云图

当高压态的钼以冲击波的形式向右推进时,上半部分没有受到其它介质的影响,一直沿直线传播;而下半部分冲击波冲击到MORB块上,与MORB块产生相互作用。图1图2为50 μs时刻的密度与压力的分布状态。此时,一部分冲击波进入MORB块中继续传播,一部分返回钼中形成稀疏波。冲击波为强间断面,在图中可以观察到明显的波面。而稀疏波属于弱间断面,在图中则显示为大范围的密度变化。冲击波的波面与稀疏波的边界组合到一起,形成一个近似的圆形。

为了验证计算结果的准确性,本文将Godunov型格式的计算结果,与文献 [6] 中HLLC的计算结果进行了比较。比较的结果可以发现,无论是密度和压力,两者的计算结果均吻合的比较好。但是由于HLLC格式中,为了保证守恒方程与非守恒方程的波面运动速度一致,对差分格式上的进行了改动,这使得两种格式的计算结果中,冲击波波面的位置与稀疏波区域边界的位置出现了小幅度的差异。

4. 结论

本文针对Mie-Grüneisen多介质混合模型结构复杂、非线性强的特点,提供了Godunov型的计算格式。经过推导后的Godunov型格式,可以同时考虑到模型中的守恒与非守恒关系式,且求解过程不需要借助其它的数值手段。通过和其它格式的数值结果对比,证明Godunov型格式可以得到稳定无振荡的数值解。

基金项目

国家自然科学基金资助项目(11702066);广东省自然科学基金(2017A030313275);广东海洋大学“创新强校工程”省财政资金支持项目(GDOU2016050258, GDOU2017052623)。

文章引用:
吴宗铎, 张建, 严谨, 刘柄宏, 董智惠. Mie-Grüneisen多介质混合计算模型的Godunov型格式[J]. 应用数学进展, 2019, 8(2): 350-356. https://doi.org/10.12677/AAM.2019.82039

参考文献

[1] Toro, E.F. (2009) Riemann Solvers and Numerical Methods for Fluid Dynamics. Springer-Verlag, Berlin Heidel-berg.
[2] 王东红, 赵宁, 徐爽. 适用于Mie-Grüneisen状态方程的一种通用激波限制器[J]. 南京航空航天大学学报, 2014, 46(2): 232-238.
[3] 李晓杰, 闫鸿浩, 王小红, 等. Mie-Grüneisen状态方程的物理力学释义[J]. 高压物理学报, 2014, 28(2): 227-231.
[4] Shyue, K.-M. (2001) A Fluid-Mixture Type Algorithm for Compressible Mul-ticomponent Flow with Mie-Grüneisen Equation of States. Journal of Computational Physics, 171, 678-707.
https://doi.org/10.1006/jcph.2001.6801
[5] 吴宗铎, 宗智. Mie-Grüneisen状态方程下多介质守恒型欧拉方程组的数值计算[J]. 计算物理, 2011, 28(6): 803-809.
[6] 吴宗铎, 严谨, 宗智, 等. 一种基于HLLC算法的Mie-Grüneisen多介质混合模型[EB/OL]. 计算物理. http://kns.cnki.net/kcms/detail/11.2011.o4.20181219.1145.004.html, 2018-12-24.
[7] 任玉新, 陈海昕. 计算流体力学基础[M]. 北京: 清华大学出版社, 2006.
[8] 刘晓慧, 李肖. ALE框架下虚拟流体方法模拟移动边界的一维可压缩多介质流动问题[J]. 应用数学进展, 2018, 7(12): 1638-1644.