原始变量守恒形式控制方程的时间准确性分析*

张怀宝，王光学，王靖宇

(中山大学 物理学院， 广东 广州 510006)

传统的可压缩流体求解器通常采用以守恒变量为自变量的守恒形式的控制方程，其中为代表的国内外比较著名的求解器有HOSTA[1]，TRIP[2]，OpenCFD[3]，FUN3D[4]，OpenFOAM[5]等。采用有限体积方法，积分形式下的控制方程一般可以表示为：

(1)

式中，Q为守恒变量，F为无黏通量，Fd为黏性通量，S为源项，t为时间，n为单元控制面的法向矢量。通过链式法则，在时间项上引入雅克比矩阵J=∂Q/∂P，可以将上述守恒变量守恒形式的控制方程转换为以原始变量P为自变量的守恒形式的控制方程[6-7]，即

(2)

对于一维Euler方程，上述各个变量的定义为：

(3)

在某些具体的数值运算中，都需要用到原始变量，例如半点处的变量重构多采用原始变量[8]；另外，对黏性流体的求解，动力黏度和传热系数都是温度的函数[9]，而黏性应力的计算需要采用速度梯度；边界条件，例如无滑移边界，一般假设壁面压强的法向梯度为零[10]。而采用守恒变量形式的控制方程时，需要在每一个时间步推进后，将更新的守恒变量转换成原始变量，进行相关计算。特别地，当涉及热化学非平衡问题时，温度是总能的隐函数，直接求取温度比较困难，实际一般采用迭代方法，时间耗费较大，而采用原始变量形式，无须该转换过程[11]。另外，采用原始变量形式可以大大简化时间隐格式中雅克比矩阵的构造，推导解析形式的雅克比矩阵相对容易，而实现数值形式的雅克比矩阵更为直接[11]。

本文采用的原始变量类型的控制方程基于压力、速度和温度，一方面该组变量的选取与Weiss等[6]在构造低速预处理矩阵时的选取保持一致，因此现有的计算平台可以直接采用其发展的预处理方法，并将计算能力扩展到低马赫数问题；另一方面，温度的计算相对于其他变量，尤其针对热化学问题，较为重要也较为困难[12-14]；最后采用速度作为原始变量，对应动量方程，较为直观。在对低速预处理系统的相关研究中，Turkel[15]指出，预处理后的控制方程尽管对于定常问题满足守恒性质，但是对于非定常问题则不成立，原因是预处理矩阵的引入破坏了时间准确性。然而从数值分析发现，即使不对引入的雅克比矩阵进行低速预处理，基于原始变量的控制方程在离散之后，对于非定常问题也不能得到准确的数值解。数值试验也发现，如果流场光滑，非定常的数值结果误差或许并不明显，但是在梯度较大或者存在强间断时，流场解与准确解相比会出现较大差异。为了消除该数值方法带来的误差，获得准确的非定常解，可以构造相应的双时间步方法，虚拟时间项采用原始变量，而物理时间项采用守恒变量，在相邻的两个物理时间步内作为定常问题求解。尽管虚拟时间推进过程中存在误差，但是方程最终会收敛到物理时间上的守恒形式。计算结果表明，该双时间步方法与守恒变量守恒形式下的双时间步方法计算结果相同，两种方法等价。

1 数值方法

1.1 守恒变量守恒形式控制方程的离散

积分公式(1)可以表示为：

(4)

式中：

R=-A(F-Fd)·ndA+∭VSdV

(5)

R为右端项(或者残差项)。 对某一个控制单元i，右端项在空间上离散，可以表示为：

(6)

式中，Aj为单元i的第j个控制面的面积，nj为该面的法向矢量，N为单元i的控制面总个数，Vcell为单元i的体积。 不考虑网格的运动，时间项的离散采用一阶向后差分，即

(7)

(8)

第n+1时间步的流场解为：

(9)

1.2 原始变量守恒形式控制方程的离散

类似地，积分公式(2)可以表示为：

(10)

式中：

R=-A(F-Fd)·ndA+∭VSdV

(11)

即右端项与式(5)相同，该项的空间离散可以采用式(6)表示。 不考虑网格的运动，时间项的离散采用一阶向后差分，即

(12)

(13)

第n+1时间步的流场解为：

(14)

1.3 非定常解的时间准确性

不难看出，尽管J=∂Q/∂P，但是Q≠JP，而是Q=JP+M。其中：

(15)

因此原始变量方法得到的守恒变量值为：

(16)

(17)

(18)

1.4 基于原始变量的双时间步方程

(19)

对式(19)进行一阶显式离散，得：

(20)

(21)

(22)

应用近似关系式

(23)

双时间步方程的线性方程组最终表示为：

(24)

2 激波管问题与结果分析

非定常流动算例采用典型的一维激波管问题，其初值条件见表1。驱动段与非驱动段的分隔位置位于x=0，气体常数R=287.06，比热比γ=1.4，各个物理量均进行无量纲化处理。

表1 激波管问题初值条件

计算采用中山大学计算流体力学研究中心开发的非结构网格计算平台MulPhy。针对激波管问题，该求解器只求解Euler方程，采用Median-Dual格点格式离散控制体，通量格式采用AUSM+-up[20]，控制方程可以采用原始变量或者守恒变量作为自变量，时间推进采用一阶显式格式，线性求解器使用GMRES[21]方法。计算网格为二维，网格数量为600×20，分别对应x与y方向，且数值实验表明，该网格尺度下空间上的数值误差与时间上的数值误差相比，可以忽略不计。时间推进的步长为Δt=1×10-7，总共推进1000个时间步，对应物理时间t=1×10-4。双时间步计算中，每两物理时间步之间的收敛残差降到初始残差的10个量级以下。图1(a)～(c)分别是压力、速度、温度的计算结果。在每张图中，四组计算结果分别对应解析解，守恒变量方法(式(9))数值解，原始变量方法(式(14))数值解，原始变量+双时间步方法(式(24))数值解。从结果对比中可以看出，采用原始变量方法与守恒变量方法得到的压力和速度曲线略有差别，但差别不大，而温度曲线差别较为明显。特别地，在激波强间断附近，原始变量方法会产生较大误差。从图1(c)可以看出，激波后的温度出现明显的上冲，而采用守恒变量方法得到的数值解与解析解保持一致。图2为温度对比放大图，由图2可以看出，激波的预测位置，也出现了偏差。进一步研究表明，随着时间的继续推进，激波位置的误差会继续增大。

(a) 压力曲线 (a) Pressure curve

(b) 速度曲线 (b) Velocity curve

(c) 温度曲线 (c) Temperature curve图1 三种数值解与解析解的原始变量对比Fig.1 Comparison of the primitive variables using three numerical approaches against analytical solution

(a) 激波与接触间断附近 (a) In the vicinity of the shock and contact discontinuity

(b) 激波附近 (b) In the vicinity of the shock图2 三种数值解与解析解的温度对比放大图Fig.2 Contrast magnification of the temperatures using three numerical approaches against analytical solution

将守恒变量方法的数值解作为准确解，计算原始变量方法数值解的误差范数，结果见表2。其中相对误差范数的参考量采用驱动段的初值条件，由于初速度为零，涉及速度的参考量均采用驱动段的初始声速a=347.3。可以看出，采用原始变量方法，各个原始变量(p,u,T)与各个守恒变量(ρ,ρu,ρE)均有不同程度的误差。

表2 激波管问题数值结果的误差范数统计

(25)

的解。

因此将守恒变量方程(1)引入双时间步方法并增加一组数值解，对两组采用双时间步方法的数值解进行比较，温度曲线如图3所示，可以看出两曲线完全重合；再次对各个流场变量的误差范数进行统计，L2范数都在1.0×10-10以下，表明两种方法等价。

图3 基于原始变量的双时间步方法与基于守恒变量 的双时间步方法的温度计算结果对比Fig.3 Comparison of the temperatures predicted by primitive-variable based and conservative-variable based approaches both together with dual-time stepping technique

3 双马赫反射问题与结果分析

为进一步验证基于原始变量的双时间步方法能够恢复准确的数值解，对Woodward等[22]首先提出的双马赫反射问题进行计算。双马赫反射问题的初始条件和边界条件如图4所示，初始时刻，Ma=10的激波与平板呈60°，并向平板方向运动，激波与平板壁面相交产生激波反射过程。运动激波前后的气体状态分别为：

其中：u为x方向速度，v为y方向速度。

二维计算域为[0,4]×[0,1]，采用均匀网格，网格量为1921×480，分别对应x与y方向，网格收敛性试验表明，该网格尺度满足空间精度要求。时间推进的物理时间步长为Δt=4×10-5，总共推进5000个时间步，对应物理时间t=0.2。双时间步计算中采用定虚拟时间步Δτ=4×10-5，并发现每两物理时间步之间迭代20个虚拟时间步后结果即已经收敛，实际计算采用30个虚拟时间步。通量格式采用AUSMPW+[23]。

图4 双马赫反射问题的初始条件和边界条件Fig.4 Initial and boundary conditions of double Mach reflection

数值计算首先采用基于原始变量的双时间步方法(式(24))，得到t=0.2时流场的密度云图，如图5所示(图中30条等势线表示的密度范围为2～22)。可以看出，随着时间推进，流场演化出丰富的流动现象，数值方法能够清晰地捕捉到马赫反射后形成的三叉点、入射激波、马赫杆、反射激波和接触间断以及近壁面结构。

图5 双马赫反射问题t=0.2时密度云图Fig.5 Density contours at t=0.2 for the double Mach reflection

作为对比，采用基于守恒变量方程(1)的双时间步方法重新计算，得到另一组数值解，然后对两组数值解的密度变量进行误差范数统计，L2范数值为8.14×10-12，表明尽管该问题物理现象非常复杂，两种方法的计算结果仍然相等，证明两种双时间步方法等价。

4 结论

1)将守恒变量守恒形式控制方程的时间项引入原始变量作为自变量，非定常数值解无法保证时间准确性。

2)针对性地构造基于原始变量的双时间步方法，能够恢复准确的数值解，该方法等价于基于守恒变量的双时间步方法。

3)通过一维激波管问题与双马赫反射的数值实验，上述两个结论均得到验证。基于原始变量的双时间步方法可以推广到一般的非定常流动问题中。