高阶气体动理学格式在湍流数值模拟中的应用

曹贵瑜, 潘 亮, 徐 昆

(1. 南方科技大学 前沿与交叉科学研究院, 深圳 518055; 2. 北京师范大学 数学科学学院, 数学与复杂系统教育部重点实验室, 北京 100875;3. 香港科技大学 数学系, 香港)

0 引 言

湍流在自然现象和工程应用中无处不在[1]，湍流问题的研究对于航空、航天、天气预报等领域的国家战略需求和学科发展具有重要意义。目前，难以采用理论解析方法直接研究湍流问题，因此通常需要借助实验手段和数值模拟进行研究。

随着数值计算方法和超级计算机的发展，直接数值模拟[2-4](Direct Numerical Simulation，DNS)逐渐成为研究湍流机理的重要手段。对于不可压缩湍流的直接数值模拟，谱方法[3](Spectral Method)、伪谱法[5](Pesudo-Spectral Method)以及格子玻尔兹曼方法[6](Lattice Boltzmann Method，LBM)得到了充分的发展。但是，以上方法都不能应用于带间断的可压缩湍流模拟。目前，高阶有限差分法[7-8]被广泛应用于可压缩湍流的直接数值模拟，但由于其数值不稳定性，仍难以适用于超声速湍流马赫数流动。为了在间断区域避免非物理振荡、光滑流动区域保持高精度，Wang等结合WENO重构和紧致有限差分格式发展了高精度混合格式[9]，并成功应用于超声速湍流马赫数的各向同性湍流直接数值模拟中。鉴于有限体积格式的物理守恒性，且能适用于复杂网格和复杂几何外形，二阶有限体积格式被广泛应用于可压缩工程湍流的计算中[10]。但是，对可压缩湍流的直接数值模拟却鲜有报道。在可压缩湍流场中通常包含激波、边界层、涡以及它们之间的相互作用等复杂的流动结构。二阶精度的数值方法具有较大的数值耗散和数值色散，难以捕捉精细的流场结构。近年来，高阶有限体积格式逐渐兴起，克服了二阶格式的精度不足问题，逐渐成为可压缩湍流直接数值模拟的重要研究工具。

在过去的三十年里，气体动理学格式[11-12](Gas-Kinetic Scheme, GKS)得到了系统性的发展，目前已成功应用于从低速到高超声速的全速域流动问题的数值模拟。在有限体积格式的框架下，基于Bhatnagar-Gross-Krook(BGK)方程[13]的积分解直接构造网格界面上的分布函数，求矩得到相应的数值通量进而更新下一时刻网格单元内的宏观守恒量。相比于传统的Riemann求解器，可以同时考虑从自由分子流到连续流的多尺度演化过程，不仅能准确给出流场光滑区域的解，而且能在间断区域很好地捕捉激波。基于Chanman-Enskog(C-E)展开，可以同时处理无黏和黏性问题，避免额外的黏性项离散。此外，气体动理学格式是真正的多维格式，可以将切向信息包含在数值通量中，这对于三维流动和间断流动十分重要。近年来，高阶精度气体动理学格式(High-Order Gas-Kinetic Scheme, HGKS)发展也极为迅速。基于WENO空间重构[14-15]和速度分布函数高阶展开，发展了单步三阶的气体动理学格式[16-17]。该格式可避免使用传统高精度格式中界面上通量的Gauss积分和Runge-Kutta时间离散。但是，高阶分布函数的形式十分复杂，给发展更高阶数值格式和三维实战编程带来巨大的困难。相比于时间空间解耦的Riemann求解器，时空耦合的气体动理学求解器可以提供通量的时间导数。基于两步四阶时间离散和高阶精度空间重构[18]，我们发展了具有四阶时间精度的气体动理学格式[19]。与原有的单步高阶气体动理学方法相比，两步四阶时间离散方法仅需要二阶GKS通量，通量形式得到大大简化，大大降低计算量。目前，结构网格和非结构网格上的两步四阶气体动理学格式都得到了充分发展，并对无黏流和层流的标准算例都进行了验证[20-21]。该格式不仅具有更高的数值精度和更好的稳定性，而且有处理复杂流动问题的能力。

近几年，GKS也不断被国内外研究者应用于湍流的数值模拟。对于低雷诺数湍流问题，二阶GKS和基于WENO重构的空间高分辨率GKS已被用作直接数值模拟工具[22-23]。对于高雷诺数工程湍流问题，通过适当地给定BGK方程湍流弛豫时间[24]能引入湍流黏性，可以将气体动理学格式与涡黏湍流模式耦合，实现不可解网格上RANS模拟。目前，耦合各类RANS模型、LES模型和混合湍流模型[25-30]的二阶GKS和高分辨率GKS已被国内外多个课题组成功应用于高雷诺数工程湍流问题的计算中。考虑到湍流数值模拟对空间和时间都要求高分辨率，例如直接数值模拟需要解析Kolmogorov时间尺度和空间尺度上的所有湍流结构，将高精度高数值稳定性的两步四阶HGKS推广到湍流直接数值模拟和RANS计算中十分必要。鉴于此，本文回顾HGKS在湍流直接数值模拟和RANS模拟中的应用。

1 高阶气体动理学格式

对于有限体积方法，关键过程是通过数值通量更新每个控制体内的守恒变量。GKS的数值通量是基于BGK方程的积分解。无外力场的三维BGK方程为:

式中碰撞不变量矢量

相空间微元为dΞ=dudvdwdξ1…dξN。在连续流区域，对BGK方程做C-E展开，

f=g-τDug+τDu(τDug)+…

(4)

式中Du=∂/∂t+ui∂/∂xi。根据气体动理论，零阶展开f=g对应推导Euler方程；一阶展开f=g-τDug对应推导NS方程。这就是BGK方程求解NS方程的基础。GKS就是通过求解BGK方程式(1)达到求解NS方程的目的。

BGK方程对碰撞不变量ψ求矩，并在有限控制体上离散：

式中，|Ωijk|为有限控制体体积。以x方向上数值通量为例，

式中，xi+1/2,jm,kn=(xi+1/2,yjm,zkn)T，(yjm,zkn)是在界面yj×zk上的高斯点坐标，ωmn是高斯点权重。对于HGKS，核心任务就是获得界面高斯点上的分布函数f(xi+1/2,jm,kn,t,u,ξ)，进而通过取矩求得界面通量以更新下一时刻单元内的宏观守恒量。BGK方程局部形式解提供了界面上的分布函数：

f(xi+1/2,jm,kn,t,u,ξ)=

(8)

式中：x′=xi+1/2,jm,kn-u(t-t′)是分子运动轨迹，f0是界面初始速度分布函数，g是界面上沿分子运动轨迹对应的平衡态分布函数。可以看到，这个积分解包含了时间和空间的所有信息。经过简单计算，界面分布函数显式表达为：

f(xi+1/2,jm,kn,t,u,ξ)=(1-e-t/τ)g0+

e-t/τgl[1-(τ+t)(alu+blv+clw)-τAl)]H(u)+

e-t/τgr[1-(τ+t)(aru+brv+crw)-

τAr)](1-H(u))

(9)

其中H(…)为Heaviside函数，gl和gr对应界面两侧的初始平衡态分布函数，g0是界面上初始平衡态分布。g0通过如下求得。

〈aku+bkv+ckw+Ak〉=0

为了获得时空上的高阶精度，高阶有限体积格式需要高阶空间重构和多步时间离散。基于高精度空间重构，可以完全确定界面上的分布函数如式(11)，对其求矩可以得到界面上时空耦合的数值通量。为了获得高阶空间精度，在结构网格的计算中，采用经典的五阶WENO重构，具体细节见参见文献[14-15]。对于二维和三维问题，采用一维一维重构方式，来获得界面高斯点值和多维空间梯度值，进而获得时空耦合的数值通量。为获得时间高阶精度，基于守恒律方程的Lax-Wendroff形求解器，Li等发展了两步四阶的时间离散方法[18-19]。对于时空耦合的演化过程，可以利用数值通量和数值通量的时间导数值实现高阶时间精度。对于半离散有限体积格式式(5)，可以采用以下的两步格式进行离散：

其中Q*为t*=tn+Δt/2时刻的守恒量。可以证明，由式(12)和式(13)给出的Qn+1具有四阶时间精度。和四阶Runge-Kutta方法相比，仅需要一个中间步就可实现四阶精度，提高了数值格式的计算效率[20]。更重要的是，两步四阶格式有着和二阶格式相当的数值稳定性。至此，基于二阶GKS通量，五阶WENO空间重构和两步时间离散的两步四阶HGKS简述完毕。

2 数值算例

本节中，选取经典的无黏算例、层流算例、低雷诺数湍流直接数值模拟算例和高雷诺数工程湍流RANS模拟算例，来验证HGKS计算复杂流动问题的能力。

2.1 Titarev-Toro问题

第一个算例是Titarev-Toro 问题[31]，初值条件如下：

这个算例描述激波和高频振荡波的相互作用。左边界采用无反射边界条件，右边界给定固定的初值波形。采用1000个均匀网格进行计算，并基于两步四阶、两步五阶、三步五阶气体动理学格式和WENO重构测试该算例[20]。在图1中给出Titarev-Toro问题在t=1.8时刻的数值解和精确解。相比与传统的基于Runge-Kutta时间离散和Riemann求解器的数值格式相比，HGKS能很好地捕捉流场中的高频波。这也证明了时间空间耦合的数值通量对于高阶格式的重要性。

图1 Titarev-Toro问题，t=5的密度数值解和精确解[20]

2.2 双马赫反射问题

双马赫反射问题是可压缩无黏流动的经典算例[32]。该算例的计算域为[0,4]×[0,1]。初始时刻x=1/6处有一道与x轴呈60°的马赫数为10的激波，激波前后的初值条件如下：

(ρ,U,V,p)=(1.4,0,0,1)

从x=1/6开始壁面上采用反射边条件，其余部分的下边界和左侧边界一样采用来流条件，右侧边界采用出口条件，并根据激波的移动位置给出上边界的边界条件。图2中给出WENO-Z格式采用1440×480均匀网格时的局部密度分布。数值结果表明HGKS有很好的健壮性，并且能够很好地分辨剪切层的不稳定性。

图2 双马赫反射问题，t=0.2的数值解[19]

2.3 平板边界层问题

二维平板边界层问题是低速黏性流的经典算例。算例中，来流马赫数Ma=0.1，雷诺数Re=1×105。壁面上采用无滑移条件，平板之前采用对称条件，其它边界根据Riemann不变量给定无反射条件。网格数为80×40的非均匀网格。图3中给出无量纲的U和V，数值结果和Blasius 参考解吻合得很好。HGKS仅用五个网格就能解析边界层速度型。

图3 平板边界层问题，无量纲的U和V，数值结果和Blasius 参考解[19]

2.4 低速槽道湍流

槽道湍流是研究壁湍流的典型算例。低速槽道湍流计算条件：马赫数Ma=0.1，壁面摩擦雷诺数Reτ=180。计算域为[0,2π]×[-1,1]×[0,π]，G1和G2两套计算网格分别为963和1283。初始速度场采用泊肃叶理论解叠加幅值为当地速度10%的白噪声随机场。通过外加质量力驱动流动, 保证流量为常数。流向和展向给定为周期边界条件，法向壁面为等温壁。

图4(a)中给出了归一化的流向平均速度分布。其中，作为参考解的谱方法计算网格为129×192×160。LBM作为模拟不可压缩流动流行的介观方法，其计算网格为200×400×200，DUGKS则采用1283计算网格[33]。可以看出HGKS、LBM、DUGKS和参考解谱方法[2]计算结果都符合得很好。图4(b)显示了平均雷诺剪切应力曲线。与谱方法结果相比，G2网格上HGKS和DUGKS结果明显优于LBM结果，尤其是在近壁区域。

(a) 平均流向速度

图5给出了用壁面摩擦速度无量纲化的三个方向的脉动速度均方根值。由图可见，流向速度脉动峰值明显大于另外两个方向的脉动，且在槽道中心区域湍流脉动趋于各向同性。从图5还可以看出，LBM在近中心线区域表现较好而在近壁区欠佳，DUGKS在近壁区域表现较好而在近中心线区域欠佳。HGKS在近壁区域和近中心线区域都与谱方法结果更接近，优于二阶LBM和DUGKS。此外，与LBM相比，HGKS结果是在较大计算域的粗网格分辨率获得的。由于LBM需要等距网格，因此网格被限制为一个极小的值以求解黏性底层，即在相同尺寸的计算域中，LBM所需的网格数将远超HGKS。

(a) 流向脉动速度均方根

2.5 超声速均匀各向同性衰减湍流

可压缩均匀各向同性衰减湍流常用来检验湍流直接数值模拟中数值格式的鲁棒性。算例R1和R2的初始泰勒微尺度雷诺数分别为Reλ=72和Reλ=120，初始湍流马赫数均为Mat=2.0，动力黏性系数按幂律指数为0.76给出。计算域为[0,2π]×[0,2π]×[0,2π]，R1计算网格为3843，R2计算网格为5123，HGKS网格的分辨率和时间步长以参考准则[35]设定。初始速度脉动场按无散条件给定，能谱为:

式中，A0=0.00013确定初始湍动能，κ是波数，κ0=8确定能谱峰值。密度和压力场为均匀分布。三个方向均采用周期边界条件。

图6(a)给出了无量纲时间t/τto=0.5和t/τto=1.0(τto为大涡翻转时间)速度场胀量θ概率密度分布函数。所有胀量PDF都显示出负尾巴，这是由激波结构导致的可压缩各向同性湍流中最显著的流动结构。图6(b)给出了无量纲时间t/τto=0.5和t/τto=1.0的胀量和x方向速度分量在x=0和z=0截线分布。截线预示，较强的可压缩区域和高膨胀区域频繁且随机出现在湍流场中。为了更直观地显示流场结构，图7给出了无量纲时刻t/τto=0.5时θ/〈θ〉*云图。其中〈θ〉*是流场胀量均方根。通常，θ/〈θ〉*≤-3的强压缩区域被识别为小激波[36]。这些随机分布的小激波和高膨胀区域导致流场中的剧烈变化的空间梯度和时间梯度，给高阶数值格式稳定性带来巨大的挑战。很少有高阶数值格式能模拟超声速湍流马赫数的均匀各向同性衰减湍流。HGKS在高达Mat=2.0的高湍流马赫数的DNS成功应用证实了HGKS优秀的鲁棒性，为高可压缩湍流提供了有力的计算工具。

(a) 速度场胀量θ概率密度分布函数

(b) 胀量θ和x方向速度分量在x=0和z=0截线分布

图7 超声速各向同性衰减湍流算例R1在无量纲t/τto=0.5时间θ/〈θ〉*云图[34]

2.6 亚声速NACA0012翼型湍流

亚声速NACA0012翼型湍流是NASA湍流模拟网站[37]提供的标准湍流模型验证算例。自由来流条件如下：马赫数Ma=0.15，雷诺数Re=3.0×106，翼型攻角15°，弦长c=1.0，升力系数和阻力系数计算中参考面积为1.0。 计算域和边界条件与NASA网站保持一致，其中G3和G4两套计算网格设置分别为225×65和897×257。粗网格G3被应用于隐式HGKS(Implicit HGKS, IHGKS)和二阶隐式GKS(Implicit GKS, IGKS)模拟，参考解为CFL3D在密网格G4计算结果。对于RANS计算，式(1)中分子弛豫时间τ调整为:

其中μt为湍流模型[30]提供的湍流黏性系数。本文耦合k-ωSST 湍流模型求得湍流黏性系数，进而调整BGK方程中的湍流弛豫时间，时间离散采用LUSGS隐式算法。

图8中显示了NACA0012翼型周围的压力系数分布和翼型上壁面摩擦系数分布。粗网格G3上，相比于二阶隐式GKS的结果，隐式HGKS的压力系数和摩擦系数与实验数据[38]和CFL3D在密网格G4上参考解更加接近。阻力系数CD对翼型周围的压力分布和壁面摩擦力分布非常敏感。 如表1所示，隐式HGKS的升力系数CL和阻力系数CD在粗网格G3上求解结果非常接近密网格G4上CFL3D的参考解。 特别是粗网格G3上隐式HGKS的CD和参考解的差异在3个阻力数(0.0001)之内，达到了工程湍流模拟的要求。 对于二阶隐式GKS，升力系数是可以接受的，而它过高地预测了阻力系数达到40个阻力数。 该算例证实了隐式HGKS高精度RANS模拟湍流场的能力。

(a)

表1 NACA0012升力系数CL和阻力系数CD

2.7 跨声速ARA M100翼身湍流

三维跨声速ARA M100翼身湍流来流条件：马赫数Ma=0.8027，基于弦长雷诺数Relc=1.31×107(弦长lc=0.245 )，翼身攻角2.873°。采用CFL3D第6版网站[39]提供的计算域和网格，其中C-O型网格设置为321×57×49。 图9显示了ARA M100机翼机身，其中黑色部分为机翼，绿色部分为机身。

图9 ARA M100翼身构型[34]

图10(a)给出了翼尖附近截面Z/b=0.935的马赫数的云图，该云图显示激波-边界层相互作用，证实了当前隐式HGKS在激波捕捉方面的鲁棒性。CFL3D网站的实验数据、隐式HGKS、二阶隐式GKS和基于SA模型的CFL3D的在翼根附近典型截面Z/b=0.123压力系数Cp的结果比较如图9(b)所示。所有方法在该截面压力系数Cp都与实验数据吻合较好。以上结果表明相比于数值离散的误差，湍流模型的误差在跨声速三维复杂RANS模拟中占主导地位。该算例验证了隐式HGKS的鲁棒性和模拟三维工程湍流的能力。

(a)

3 结 论

本文回顾了高阶气体动理学格式及其在湍流数值模拟中的应用研究。时空耦合的气体动理学求解器可以提供通量的时间导数。基于两步四阶时间离散和高精度空间重构，我们发展了具有四阶时间精度的气体动理学格式。该格式不仅具有更高的数值精度和更好的稳定性，而且有处理复杂流动问题的能力。数值结果表明，高阶气体动理学格式可以为可压缩湍流的数值模拟提供有力的计算工具。未来，将使用高阶气体动理学格式研究更具有挑战性的可压缩湍流问题，例如超声速湍流边界层和激波边界层相互作用等。希望为可压缩强间断湍流的流动机理研究探索新方向。

同时值得指出的是，气体动理学格式的主要优势是给出了在网格边界上气体分布函数的演化解。除了给出对应的数值通量外，在网格边界上的其它物理量也能够显示地给出来，比如网格边界上在这一时间步长结束时的守恒量。所以格式除了更新网格里面的守恒量的值，每个网格内的守恒量的梯度也可以根据高斯定律由边界上的值积分直接确定。这一性质对构造紧致高阶气体动理学格式提供了可靠的网格局部的动力学基础，从而避免了用很远处和本网格没有直接动力学关系的量来做高阶重构。这方面的研究在过去几年得到了飞速发展[21]。这一基于局部高阶动力学演化解的指导方向是发展下一代高阶计算流体力学格式的可行之道。

致谢:感谢广州天河超算中心提供计算资源。