KT格式数值求解几类典型流体力学问题*
2018-04-20蔡振宇李明军
蔡振宇, 李明军
(湘潭大学 数学与计算科学学院,湖南 湘潭 411105)
KT格式[1-3]是一种新的高精度中心格式,于2000年首先由Kurganov和Tadmor共同提出.KT格式利用波在当地传播的最大速度来估计黎曼扇的宽度,从而将计算区域分为光滑区域和非光滑区域,并对两部分分别进行积分求解.KT格式除了具有中心格式所具有的无需求解黎曼求解器,算法简单的优点以外,该格式还具有数值耗散小,稳定性高,捕捉激波能力强等特点.在已有中心格式的基础上,KT格式精确地利用波在当地传播的最大速度来计算非光滑区域.在非光滑区域内使用与当地速度相关的更小网格,精确地估计黎曼扇的宽度,对黎曼扇进行积分近似求解.KT格式最大的特点是无需求解黎曼求解器,所以算法简单.该格式的数值粘性与Ο(Δx2r-1)同阶,与时间步长没有关系.所以当Δt→0时,该格式可以写成半离散格式,它的数值耗散更小,分辨率高.
本文利用KT格式分别求解一维、二维激波管问题,并首次利用该格式实现双马赫反射问题计算[4-5].通过观察、分析计算结果,并将一维激波管问题计算结果与其他几种经典格式(Steger-Warming,van Leer,AUSM)进行比较,发现KT格式能较好地模拟这几类流体问题.通过比较发现,KT格式在间断处展现出更高的精度,具有一定的优越性.本文结果很好地验证了KT格式数值耗散小、分辨率高、捕捉激波能力强等相关理论.
1 控制方程
激波管问题和双马赫问题所用的控制方程都是欧拉方程,它包含了质量守恒,动量守恒和能量守恒定律[9].二维欧拉方程形式为:
(1)
式中,U=(ρ,ρu,ρv,ρe0),F=(ρu,ρu2+p,ρuv,(ρe0+p)u),G=(ρv,ρuv,ρv2+p,(ρe0+p)v).在这里,ρ是密度,u是x方向的速度,v是y方向的速度,p是压力,e0是单位质量总能量e0=e+(u2+v2)/2,e是单位质量内能.
2 数值方法
2.1 空间离散[10-12]
为了阐述更加简单,便于理解,在这部分我们以标量守恒方程为研究对象:
(2)
(3)
在这里使用积分中值定理可得
(4)
若Δt项趋近于零,则
可得:
(5)
2.2 时间离散[5,8]
本文采用三阶Runge-Kutta显式方法对时间进行离散、推进求解.其具体表达形式为:
u(1)=un+ΔtnC[un],u(l+1)=ηlun+(1-ηl)(u(l)+ΔtnC[u(l)]),un+1=u(3),l=1, 2.
3 结果与讨论
3.1 一维激波管问题计算结果及其分析
对于一维激波管问题,在[0,1]区间内划分了500个均匀的网格.在边界处,采用了固壁边界条件.初始条件如下:
p=105,ρ=1,u=0,x≤0.4;p=104,ρ=0.125,u=0,0.4
图1显示了用KT格式计算一维激波管问题时,所得流体的密度、速度和压力的分布情况.在每个图中对比了数值解和精确解.结果显示,KT格式对流体中经典物理量的真实分布情况可以很好的近似.
3.2 二维激波管问题计算结果及其分析
对于二维激波管问题,在[0,1] ×[0,0.1]矩形区域内,划分了200 ×10个均匀的网格.在边界处采用了固壁边界条件.初始条件如下:
p=1,ρ=1,u=0,v=0,x≤0.5;p=0.1,ρ=0.125,u=0,v=0,0.5
如图2(a)、(b)、(c)所示,分别展现了在二维激波管问题中,流体的密度、速度和压力的分布情况.容易看到,KT格式的计算结果可以很好地反应真实的物理流动,无论是膨胀波、接触间断还是冲激波都被比较准确地捕捉到.
4 结 论
本文研究了KT格式在模拟几种典型流体力学问题时的性能.KT格式利用波在当地传播的最大速度来估计黎曼扇的宽度,将光滑区域与非光滑区域分开分别进行求解,能更准确地捕捉激波位置.通过对比KT、Steger-Warming格式、van Leer格式和AUSM格式在计算一维激波管问题的数值模拟结果,可以得知,KT格式具有二阶精度.在间断(膨胀波、接触间断、冲激波)出现的位置,相较于其他几种格式,具有较小的数值耗散.通过观察二维激波管问题的数值计算结果,可以看到,KT格式可以很好地模拟此类问题.由于KT格式无需求解黎曼求解器,而且能写成半离散格式,所以算法简单,编写程序也相对简单.
[1]KURGANOV A,TADMOR E.New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[J]. J Comput Phys, 2000, 160:241-282.
[2]GOTTLIEB S, SHU C W. Total variation diminishing runge-kutta schemes[J]. Mathematics of Computation, 1998, 67(221): 73-85.
[3]JAMESON A, SCHMIDT W, TURKEL E. Numeical solution of the euler equations by finite volume methods using runge kutta time stepping schemes[C]//AIAA, Fluid & Plasma Dynamics Conference, 1981, 1259(11): 2004-4325.
[4]LIEPMANN H W,ROSHKO A.Elements of gasdynamics[M]. New York:Dover Publications, 1957.
[5]张德良. 计算流体力学教程[M]. 北京:高等教育出版社, 2010.
[6]WARMING R F,BEAM R M. On the construction and application of implicit factored schemes for conservation laws[C].SIAM-AMS Proceeding, New York, 1977.
[7]STEGER J L. Coefficient matrices for implicit finite difference solution of the inviscid fluid conservation law equations[J]. Computer Methods on Applied Mechanics and Engineering, 1978, 13:175-188.
[8]LIOU M S,STEFFEN C J. A new flux splitting scheme[J]. Journal of Computational Physics, 1993, 87: 23-29.
[9]LI H,BEN-DOR G. Analysis of double-Mach-reflection wave configurations with convexly curved Mach stems[J]. Shock Waves, 1999, 9(5): 319-326.
[10]CHERTOCK A, CUI S, KURGANOV A,et al. Well-balanced positivity central-upwind scheme for the shallow water system with friction terms[J]. International Journal for Numerical Methods in Fluids, 2015, 78(6): 355-383.
[11]BELJADID A, MOHAMMADIAN A,KURGANOV A. Well-balanced positivity preserving cell-vertex central-upwind scheme for the shallow water flows[J]. Computers & Fluids, 2016, 136: 193-206.
[12]DO S, HA Y, KANG M,et al. Application of a multi-dimensional limiting process to central-upwind schemes for solving hyperbolic systems of conservation laws[J]. Journal of Scientific Computing, 2016, 69(1): 274-291.