KT格式数值求解几类典型流体力学问题*

蔡振宇,　李明军

(湘潭大学 数学与计算科学学院，湖南 湘潭 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.