0 引言

临近空间飞行器、跨大气层飞行器以及返回舱等高超声速飞行器，在其再入过程中，流动状态涵盖从稀薄流到连续流，涉及到多尺度、跨流域问题，在近连续、连续流区还会出现黏性干扰、湍流等复杂流动现象，严重影响飞行器的气动力和气动热特性。对其进行准确模拟具有非常高的挑战性。

在高空稀薄流区，受限于连续介质假设，基于Euler和Navier-Stokes方程的传统CFD方法不再适用。直接模拟蒙特卡罗方法(DSMC)从分子平均自由程或者分子平均碰撞时间尺度出发，刻画宏观气体流动的稀薄效应^[1]。传统的离散速度坐标法(DOM或DVM)基于算子分裂^[2]，直接在高维空间离散Boltzmann模型方程或其模型方程。对于三维复杂飞行器再入各流域绕流问题，需要在三维位置空间和三维速度空间离散求解分布函数，计算量太大。为求解高超声速飞行器面临的复杂跨流域绕流问题，亟需发展高效的多尺度计算方法。

香港科技大学徐昆等发展的气体动理学格式(GKS)^[3-4]基于BGK类方程的解析解来计算网格界面上的数值通量，耦合了分子间相互碰撞与分子的自由运动，离散尺度仅取决于求解精度的要求，不再受限于平均自由程和平均碰撞时间^[3]。对连续及近连续流，可以将分布函数简化，得到逼近N-S方程的GKS-NS方法，避免了离散速度坐标，极大地减少了存储量和计算量。GKS-NS在多种流动，特别是高速黏性流动中表现优异^[5]。

为了将GKS应用于实际再入问题，本文主要开展了以下几个方面的工作。首先，耦合常用湍流模型对GKS-NS进行拓展，使之能有效模拟高超声速飞行器飞行高度不太高时面临的复杂高雷诺数高速湍流问题。其次，发展非结构网格上的高精度多维GKS-NS格式，以对高超声速复杂流动问题进行精细模拟。最后，对结合速度空间的离散发展起来的统一气体动理学格式(UGKS)进行拓展及优化，使之能够模拟跨流域多尺度流动问题。

1 GKS简介

气体动理学格式基于Boltzmann方程的模型方程，BGK方程。它具有如下形式：

$ \begin{array}{l} \frac{{\partial f}}{{\partial t}} + \mathit{\boldsymbol{u}}\frac{{\partial f}}{{\partial x}} = \frac{{g - f}}{\tau }\\ g = \rho {\left( {2\pi RT} \right)^{{{\rm - }}\left( {K + 3} \right)/2}}{{{\rm e}}^{{{\rm - }}\left[ {{{\left( {\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{U}}} \right)}^2} + {\mathit{\boldsymbol{\xi }}^2}} \right]/(2RT)}} \end{array} $

(1)

其中，f=f(x, u, ξ, t)为气体分子速度分布函数，u为分子平动速度，ξ为内自由度，内自由度数为K。g为平衡态分布函数，ρ、U、T、E分别为气体的密度、速度、温度和总能，τ为分子平均碰撞时间，R为气体常数。对f 和BGK方程取矩，可以得到宏观守恒量Q及守恒方程：

$ \begin{array}{l} \frac{{\partial \mathit{\boldsymbol{Q}}}}{{\partial t}} + \nabla \cdot F = 0, \\ \mathit{\boldsymbol{Q}} = {\left( {\rho , \rho \mathit{\boldsymbol{U}}, \rho E} \right)^{{\rm T}}} = \int {\mathit{\boldsymbol{f\psi }}} {{\rm d}}\mathit{\Xi }, \\ {\mathit{\boldsymbol{F}}_\alpha } = \int {{u_\alpha }} \mathit{\boldsymbol{f\psi }}{{\rm d}}\mathit{\Xi }, \alpha = 1, 2, 3, \\ \mathit{\boldsymbol{ \boldsymbol{\varPsi} }} = {\left[ {1, \mathit{\boldsymbol{u}}, \left( {{\mathit{\boldsymbol{u}}^2} + {\mathit{\boldsymbol{\xi }}^2}} \right)/2} \right]^{{\rm T}}}{{\rm 。}} \end{array} $

(2)

相比于宏观方程，BGK方程通过分布函数的变化来刻画分子的自由运动和相互碰撞的过程，具有更丰富的物理内涵^[6-7]。

BGK方程(1)有如下形式的演化解：

$ \begin{array}{l} f\left( {\mathit{\boldsymbol{x}}, t, \mathit{\boldsymbol{u}}, \mathit{\boldsymbol{\xi }}} \right) = \frac{1}{\tau }\int_0^t g \left( {\mathit{\boldsymbol{x}}', t', \mathit{\boldsymbol{u}}, \mathit{\boldsymbol{\xi }}} \right){{{\rm e}}^{{{\rm - }}\left( {t{{\rm - }}t'} \right)/\tau }}{{\rm d}}t' + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{{\rm e}}^{ - t/\tau }}{f_0}\left( {\mathit{\boldsymbol{x}}{{\rm - }}\mathit{\boldsymbol{u}}t, \mathit{\boldsymbol{u}}, \mathit{\boldsymbol{\xi }}} \right), \\ \mathit{\boldsymbol{x'}} = \mathit{\boldsymbol{x}}{{\rm - }}\mathit{\boldsymbol{u}}\left( {t - t'} \right) \end{array} $

(3)

其中，f₀是初始分布函数，g是随时间演化的平衡态分布。上述演化解自动耦合了分子的自由运动和相互碰撞，这是确保基于该演化解的GKS成为真正多尺度方法的关键^[3]。

GKS采用有限体积法，通过演化解(3)来计算单元界面上的数值通量，从而对宏观量和分布函数进行更新：

$ \mathit{\boldsymbol{Q}}_{ijk}^{n+1}=\mathit{\boldsymbol{Q}}_{ijk}^{n}+\frac{1}{{{\mathit{\Omega }}_{ijk}}}\ \mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\nolimits_{{S}_{ijk}} \int_{{{t}^{n}}}^{{{t}^{n+1}}}{\mathit{\boldsymbol{F}}{\rm d}t\cdot {\rm d}\mathit{\boldsymbol{S}}} $

(4)

$ \begin{align} & f_{ijk}^{n+1}=f_{ijk}^{n}+\frac{1}{{{\mathit{\Omega }}_{ijk}}}\mathop{{\int\!\!\!\!\!\int}\mkern-21mu \bigcirc}\nolimits_{{S}_{ijk}} \int_{{{t}^{n}}}^{{{t}^{n+1}}}{\mathit{\boldsymbol{u}}\mathit{f}{\rm d}t\cdot {\rm d}\mathit{\boldsymbol{S}}{\rm +}} \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \frac{1}{{{\mathit{\Omega }}_{ijk}}}\iiint_{{{\mathit{\Omega }}_{{\rm i}jk}}}{\int_{{{t}^{n}}}^{{{t}^{n+1}}}{\left( g{\rm -}f \right)/\tau {\rm d}t{\rm d}\mathit{\Omega }}} \\ \end{align} $

(5)

其中，S_ijk为离散物理空间上网格单元(i, j, k)的界面面积，Ω_ijk为单元体积。演化解中的平衡态分布g可根据格式精度需要由单元界面中心关于时、空的Taylor展开逼近。

UGKS同时更新宏观量(4)和分布函数(5)，适用于从连续流到稀薄流区的全流域。并且，由于通过演化解来计算通量，可以刻画出不同网格尺度上、不同时间间隔内分子的自由运动和相互碰撞的物理过程，网格尺度、时间步长可以不受分子平均自由程和平均碰撞时间的限制，因而是真正的多尺度方法，模拟再入过程中的跨尺度、跨流域流动问题时具有很高的求解效率。

在近连续流区，初始分布函数可以直接由一阶Chapmann-Enskog展开逼近，

$ {{f}_{0}}={{\left[ g\text{-}\tau \left( \frac{\partial g}{\partial t}+\mathit{\boldsymbol{u}}\cdot \frac{\partial g}{\partial \mathit{\boldsymbol{x}}} \right) \right]}_{t=0}} $

(6)

此时演化解可完全由平衡态分布及其时、空导数确定，因此在宏观量的更新过程中，式(4)，通量F仅与Q及其空间导数有关^[4]。更为重要的是，这种情况下无需直接更新分布函数，即式(5)，因而避免了速度和内自由度空间(u, ξ)的离散，极大地减少了计算量。这样建立的GKS-NS格式仍然耦合处理了无黏与黏性项，保证了其在高超声速黏性模拟中的优异性能。

2 GKS的拓展及应用 2.1 高超声速湍流模拟

GKS-NS由于内含黏性并且耦合处理无黏与黏性项，在高超声速黏性流动模拟中表现优异。在工程高雷诺数问题中，为了高效地模拟湍流场，可进一步基于扩展BGK方程^[9]发展耦合湍流模式的拓展GKS^[10-11]。扩展BGK方程为：

$ \frac{\partial f}{\partial t}+\mathit{\boldsymbol{u}}\cdot \frac{\partial f}{\partial \mathit{\boldsymbol{x}}}=\frac{g\text{-}f}{{{\tau }_{\text{eff}}}} $

(7)

τ_eff=τ+τ_t为流动的有效松弛时间，其中湍流松弛时间通过湍流黏性系数定义τ_t=μ_t/(p+2k/3)。通过μ_t，GKS可直接与常用的湍流模式耦合，在积分尺度上描述湍流流动。τ_t定义中k为流场湍动能，通常在高超声速流动中不能忽略。在该框架下建立的拓展GKS-NS中，扩展BGK方程的多尺度演化过程退化为格式的自适应耗散机制，使格式可在间断处保持强鲁棒性，同时在光滑流场区保持较高的分辨率。另外，大部分湍流模式中需要求解湍流量的输运方程。本文中通过Li^[8]等发展的带有标量输运的GKS-NS格式，以耦合方式求解湍流量输运方程，湍流量的数值精度和耗散机制与守恒量保持一致。

拓展GKS中采用了多种湍流模型，包括零方程模型B-L、一方程模型S-A、两方程模型k-ε、k-ω、SST以及Menter转捩模式等。这里给出结合含有可压缩修正的k-ω-LS模式对高超声速压缩拐角流动的计算结果。来流马赫数为9.22，斜坡角度为θ=24°^[12]。在斜坡拐角较大时，分离区的大小以及壁面热流的准确预测对数值计算非常具有挑战性。图 1中为GKS的计算结果及其与已有实验测量及数值结果的比较。当计算采用相同湍流模式时，拓展GKS的计算结果比已有数值结果^[13]与实验符合得更好。其中，热流峰值位置基本与实验一致，热流峰值大小与实验的误差约为41%，精度高于参考计算结果。拓展GKS在高超声速模拟中的优势来自于其内涵的自适应耗散机制，而湍流量与守恒量耦合的求解方法也有助于提高计算精度。

2.2 基于CPR框架的高精度GKS

高精度格式能够以相对较少的计算量得到精度更高的数值结果，从而提高计算效率。近年来，Huynh和Wang等人所发展的通量重构方法(FR/CPR)^[14-15]为DG、SV/SD等经典的高精度方法提供了一个统一框架。相比于传统的DG方法，FR/CPR方法摆脱了数值积分的计算量，并且通过特定求解点和通量点的选取，能够有效减少重构的计算量，因而FR/CPR方法具有很高的效率。而考虑到通量演化，基于气体动理学格式所构造的高精度GKS-NS通量求解器耦合了无黏通量和黏性通量，并且包含了通量随时间的演化，因而只需要单步推进即可同时达到时空高阶精度。此外，GKS-NS具有真正的多维特性，可以直接构造出通量在单元内的空间分布。因此，结合高效的FR/CPR框架和高精度GKS-NS通量求解器，发展适用于飞行器复杂外形的非结构网格上的新型高精度格式具有很好的应用前景。

以二维问题为例，FR/CPR框架的表达式如下

$ \frac{\partial {{\mathit{\boldsymbol{Q}}}_{ij}}}{\partial t}+\nabla \cdot \mathit{\boldsymbol{F}}({{\mathit{\boldsymbol{Q}}}_{ij}})+{{\delta }_{ij}}=0 $

(8)

其利用单元内设置的若干求解点可以插值得到守恒量和通量在单元内的分布。考虑到单元界面两侧可能存在间断，还需在每个单元界面上设置若干通量点，计算耦合相邻单元影响的通量值，并由此对单元内的通量分布进行修正，得到修正项δ_ij，进而更新每个求解点上的守恒变量Q_ij。

为了计算每个求解点和通量点上的通量值，需要分别在三角形单元内和单元界面上构造出通量分布。首先对气体分子速度分布函数进行高阶泰勒展开，通过重构得到的守恒量分布及其空间导数，以及相容条件可以计算出展开式中的各项系数。高精度GKS通量求解器所具有的多维特性使得所有求解点以及通量点上的通量值均可以同时获得，而传统的通量求解器则需要逐点计算通量。

为了计算含间断的可压缩流动问题，格式中引入了紧致型WENO限制器^[16]。该限制器能够较好地抑制间断附近的振荡，同时保证光滑区域的高精度。多种典型算例表明，新发展的时-空三阶精度CPR-GKS在可压缩流动中兼具高精度和较好的激波捕捉性能，并且计算效率能和原CPR格式相当。本文采用CPR-GKS计算了Re=200的黏性激波管问题。它包含了强激波和边界层的相互作用，不仅要求格式具有良好的稳健性，同时还对格式的精度和分辨率有较高的要求，因而是验证格式在非定常黏性流动中捕捉激波和复杂流动结构的经典算例。图 2为计算得到的密度等值线图，其中网格尺度为1/300，计算时间为t=1。由图中可以看出，格式较好地捕捉到了激波和边界层干扰产生的复杂流场，尤其是对涡结构的捕捉展现了格式的高精度。计算出的中间涡结构的高度约为0.172，与有限体积GKS得到的0.171吻合较好^[18]。该有限体积GKS采用了五阶重构，时间精度为四阶，且网格尺度为1/500。由此也进一步表明了基于内点重构的CPR框架的高精度。

2.3 适合大规模并行化计算的UGKS

UGKS能够有效模拟跨流域稀薄流动，但需要同时求解式(4)和式(5)。除了需要离散物理空间，还需要对速度空间进行离散。一般地，模拟三维稀薄流动问题，相比于直接求解宏观方程，如Navier-Stokes方程所适用的连续流问题，计算规模将有量级上的增加。对于高超声速流动，由于分布函数更宽，速度空间的离散规模将进一步增加。因此，为了模拟高超声速跨流域稀薄流动，迫切需要发展适合UGKS的大规模高效并行算法。

一方面，UGKS算法需要同时对物理空间和速度空间进行离散，核心部分为气体分子速度分布函数在物理空间上的重构、演化和投影，在演化和投影过程中涉及到分布函数在速度空间上的积分求矩。考虑到UGKS算法本身在离散物理空间和速度空间上计算过程的差异，在并行算法的设计中可以对二者进行区分。另一方面，考虑到在实际三维问题中，需要离散三维速度空间，特别是在非定常高超声速稀薄流模拟中，其存储量尤其巨大，因此有必要对速度空间进行并行分块。

本文发展的UGKS并行程序适合物理空间分块搭接的结构网格，且物理空间和速度空间同时并行分块。为了加速定常问题收敛，还引入了当地时间步长、隐式格式^[24]等措施。在数据交换中，采用了减带宽技术^[17]优化了分块界面两侧不同物理分块之间的宏观物理量和分布函数的数据交换，有效减少了数据交换的搜索和等待时间。针对速度空间并行分块策略，采用逻辑网格MPI分组算法^[19-20]，划分不同子通信域，将分块编号排列如图 3，物理空间分N_x块，数据传递在子域(Col_common)中进行，速度空间分N_v块，在子通信域(Row_common)中归约，分块编号矩阵保持了仅对物理空间进行分块时的大小，速度空间同时进行分块并没有对数据传递中分块编号的搜索和等待造成更大的负担。该并行策略能充分考虑到UGKS算法中物理空间和速度空间上的算法差异，并能显著提高大规模并行计算的效率。对其进行并行效率测试，其加速比呈线性增长^[21]，显示出了UGKS并行程序的高效并行特性。

基于本文发展的大规模高效并行UGKS程序，对多个典型稀薄流动，如三维方腔流动、圆球绕流、卫星体飞行器绕流进行了验证^[21-22]。在此基础上，对典型返回舱定常绕流问题进行了模拟。模拟真实飞行器复杂外形下的高超声速稀薄流动，尤其是捕捉多尺度精细流动结构，对格式的稳健性和计算效率提出了很高的要求。在此算例中，来流迎角20°，马赫数23.5，努森数0.27。总计算规模为140亿，其中速度空间离散点数为2万。物理空间和速度空间均分块并行，使用了5040个进程进行计算，总计算时间为25万CPU小时。图 4和图 5所示为返回舱周围的温度和压力分布云图。

3 结论

本文介绍了气体动理学格式GKS在飞行器再入过程中面临的高超声速跨流域问题中的拓展和应用。

研究工作主要包括三个方面：1)基于拓展BGK方程发展了耦合湍流模式的拓展GKS，格式在典型高超声速流动模拟中具有良好的表现；2)基于CPR框架，结合紧致型WENO限制器发展了非结构网格上的高精度动理学格式CPR-GKS，通过黏性激波管等典型算例验证了其在可压缩黏性流动中的高精度和较好的激波捕捉性能；3)发展了适合复杂分块结构网格上跨流域稀薄流动大规模高效并行计算的UGKS方法，初步对典型飞行器的再入进行了精细模拟。

本文研究揭示了GKS在再入飞行器面临的高超声速湍流及跨流域稀薄流动模拟中的优异性能，具有广阔的应用前景。

致谢: 感谢江南计算技术研究所的徐金秀高级工程师和浙江省水利河口研究院的于普兵博士的帮助。清华探索100高性能计算平台、并行科技高性能计算支持服务为本文工作提供部分计算机时。