离散时空直接建模思想及其在模拟多尺度输运中的应用

刘 畅, 徐 昆,*

(1. 香港科技大学 数学系, 香港；2. 香港科技大学 深圳研究院, 深圳 518057)

0 气体动理学理论模型和数值方法的发展简介

历史上人们对于空气动力学的研究是从宏观、近平衡态规律开始，逐步认识到微观、介观及非平衡的气体动理论，进而建立起从微观动理学到宏观动力学的联系。对于宏观连续流域的气体，通常被认为是满足连续介质假设并满足质量、动量和能量守恒定律，在此基础上Bernoulli和Euler在18世纪中叶提出了描述无黏流体的欧拉方程。19世纪中叶，斯托克斯将黏性应力项加入动量方程中，发展了纳维-斯托克斯(N-S)方程，其中应力与应变率的本构关系采用了牛顿的线性本构假设。N-S方程是流体力学中最普遍使用的控制方程，但对于高超声速、非平衡气体流动描述不够准确。在其后两百年的研究中，人们不断尝试将宏观方程的应用范围推广，进而发展出一系列拓展的宏观方程。三大守恒定律方程仍是流体宏观方程中最为可靠的方程，统一气体动理学格式的宏观演化方程也是基于宏观守恒定律而建立起来的，而且这组方程适用于从自由分子流到连续流的整个流域。

微观、介观层面气体动理论的研究可以追溯到1738年，Bernoulli从分子迁移和碰撞的角度解释了恒温气体的压力与密度成正比，19世纪50年代到80年代，Clausius、Maxwell和Boltzmann的工作完成了分子动理论的奠基和全面发展。其中，Clausius提出了分子自由程这一重要概念，即连续两次分子碰撞间分子自由迁移的距离。Maxwell 引入了气体分子速度分布函数的概念，并且得到了平衡态速度分布函数，求得了包含黏性系数、导热系数在内的气体输运系数的表达式，Boltzmann 于1879年提出并证明了H定理，并给出了速度分布函数的积分微分方程，即气体动理学理论中奠基性的Boltzmann方程。由于Boltzmann方程碰撞项较为复杂，研究者们进一步提出了线性化Boltzmann方程[2]，以及一系列的简化动理学模型方程，包括Bhatnagar、Gross和Krook三人提出的BGK模型[3]，以及ES-BGK、Shakhov和Fokker-Planck等模型[4-6]。在动理学中，努森数是衡量气体稀薄程度的参数，也是联系宏观和微观的重要参数，通常努森数Kn定义为气体分子的平均自由程和流场特征长度L的比值

(1)

在20世纪初，Chapman和Enskog研究了Boltzmann方程在小努森数下的渐进性质，从介观层面恢复了欧拉和N-S方程，并得出了宏观输运系数与微观分子碰撞参数之间的关系，建立了宏观和微观间的联系[7]。Chapman和Enskog的理论是研究动理学方程在近连续流域性质的重要方法，也是从Boltzmann方程建立以来最为成功的动理学理论。

钱学森根据努森数将气体的流域分为连续流域、滑移流域、过渡流域和自由分子流域[8]。气体在不同的流域内的流动特性有显著区别，因而通常采用不同的控制方程和计算方法，而基于离散时空直接建模思想的UGKS、DUGKS、UGKWP等多尺度格式能够适用于从无碰撞流域到连续流域的流场流动物理的直接模拟，见图1。在Kn≤10-3的连续流域内，N-S方程能够很好地描述流场的流动特性，连续流域内的计算流体力学的研究主要集中在发展高精度流体数值模拟方法和发展湍流模型两个方向。流体数值模拟方法包括基于数值求解流动控制方程的高阶计算流体力学方法和对流动物理过程进行直接模拟的基于格子气自动机的格子Boltzmann(LB)方法[9-10]。高阶计算流体力学方法的研究包括空间高阶重构和时间高阶演化两个方面，目前较为流行的空间高阶重构方法有加权基本无震荡(WENO)方法[11]，间断伽辽金(DG)方法[12]，spectral volume (SV)[13], spectral difference (SD)[14], 修正重构 (CPR)[15],通量重构(FR)[16]方法等； 在时间精度的研究方面，数值常微分方程时间高精度格式的发展源于20世纪40年代，包括单步的Taylor型方法，和以龙格库塔方法为代表的多步法。而对于数值偏微分方程时间高精度格式的研究，李杰权等首先提出了基于多步多导数的两步四阶格式(two-stage fourth order scheme)[17-19]，综合了单步Taylor型方法和多步法的优点，并强调了在高阶格式构造中时空耦合的重要性[20]。在近期发展中，基于气体动理学格式和两步四阶格式发展的高精度格式在精度和稳定性上都体现出了显著的优势[21-28]，尤其是空间八阶紧致格式能够在捕捉间断的同时，解析声学量级扰动[10, 29-30]，也说明了未来高阶格式的发展要将时间高精度演化和空间高阶重构相结合[20]。 在10-3

图1 根据努森数对气体流域的划分，以及对应流域内的模型方程和数值方法Fig.1 Classification of flow regime based on the Knudsen number, and the corresponding governing equations

在10-210的自由分子流域，非平衡效应显著出现，基于宏观方程的数值方法如N-S方程、矩方法难以准确描述该流域内的非平衡物理，因而在过渡和稀薄流域，通常需要求解动理学方程。对稀薄气体动力学数值模拟主要有两类建模方法，其一为对流动物理过程进行直接模拟的方法，其二为数值离散求解动理学方程的方法。直接模拟的方法中最为流行的是Bird 发展的直接数值模拟蒙特卡洛(Direct Simulation Monte Carlo, DSMC)方法[35]。DSMC 方法通过模拟粒子来代替真实气体粒子，采用算子分裂方法解耦粒子的自由迁移和碰撞。其中粒子迁移过程可以准确跟踪，粒子的碰撞过程采用蒙特卡洛方法计算。DSMC方法对应的物理过程简单，因而容易应用到复杂外形、化学反应和辐射输运等领域[36]。DSMC 对于高超声速问题的模拟具有很高的效率，但由于其解耦了迁移和碰撞，因而时间步长和网格尺寸都需要小于碰撞时间和平均自由程，这大大降低了DSMC方法对于近平衡流动的模拟效率。同时，粒子方法所带来的统计噪音也限制了DSMC方法对于微流动的模拟精度。为了提高DSMC方法在连续流域的计算效率，发展了渐进保持的DSMC方法(AP-DSMC)[37, 38]和矩方程引导DSMC方法[39]等。为了降低DSMC方法的噪音，发展了低噪音DSMC(LVDSMC)方法[40, 41]等。在众多统计学粒子方法中，Jenny发展的粒子Fokker-Planck方法[42]和费飞等人发展的多尺度粒子BGK方法[43]不仅适用于稀薄流域，在连续流域内也有较好的数值表现。

包括离散速度法、谱方法等在内的数值离散求解动理学方程的方法在近几十年也有长足的发展[44-55]。谱方法通过将Boltzmann碰撞项投影到谱空间，将卷积形碰撞项解耦从而提高计算效率，谱方法的计算精度很高，但由于Boltzmann碰撞项积分维度过高，受内存和计算量限制仍需要进一步发展加速算法才能适用于三维工程计算[51-53, 56]。为了降低碰撞项计算量，通常采用离散速度方法求解动理学模型方程，离散速度方法不受统计噪音影响，计算精度高，同时能够通过隐式迭代方法进行加速求解，对于微流动问题的计算有较大优势。近年来，基于离散速度的渐进保持(AP)格式[57-59]、 高阶-低阶(High Order-Low Order, HO-LO)格式的发展拓展了传统离散速度格式的适用流域[60]，并应用到多组分气体、辐射输运等领域。其包括AP格式在内的离散速度方法，其出发点都是从微分动理学方程出发，通过一定的离散方法得到数值格式，而微分动理学方程在粒子迁移和碰撞解耦的建模基础上得到的，其对应于建模尺度小于粒子碰撞时间和平均自由程。完全以离散动理学方程建立起来的数值格式通常求通量时只考虑粒子的自由输运。但对于这类数值格式，解耦了迁移和碰撞过程，为得到精确的物理解，格式需要网格尺度和时间步长小于粒子平均自由程和碰撞时间，对于连续流域的模拟效率和精度都较低[61]。

传统数值求解流动控制方程的方法通常是从微分方程出发，通过离散方法得到对应的数值格式。考虑到偏微分方程是在固定物理尺度上建模得到的，像Boltzmann方程的平均自由程和N-S方程的宏观流体力学尺度。用这种偏微分方程数值解的计算思想只能得到单尺度的计算格式，没有考虑到物理尺度变化对流动物理和方程本身描述的影响。对于复杂工程问题，如高超声速飞行的临近空间飞行器，其周围流场的局部努森数变化可以达到五个数量级(见图2)，传统的计算流体力学方法难以高效准确模拟。同样对于直接模拟方法，如DSMC方法，由于其在建模过程中解耦了迁移和碰撞，因而其适用的时空尺度被限制在平均自由程和碰撞时间量级，对于跨尺度和连续流域的模拟计算效率较低，比如DSMC方法很难适用于在80 km以下对飞行器的流场计算。徐昆提出的离散时空直接数值模拟的建模思想回归建模的出发点，在数值离散控制体上建立数值控制方程，从数值控制体上的守恒律出发，通过微分动理学方程的演化解得到网格界面局部的数值通量[62]。离散时空直接数值模拟的建模思想的核心在于将网格尺度加入建模过程中，耦合粒子迁移、碰撞、加速等过程，以及它们在一个时间步长内的累积效应，实现从稀薄到连续全流域内对流场的流动物理的准确求解。离散时空直接数值模拟的建模思想的关键在于采用了动理学方程的演化解构建数值通量，演化过程在形式上包含了粒子迁移和碰撞效应，在不同尺度上能够准确反映粒子的运动特性，例如在稀薄流域内的自由输运特性和在连续流域内由于碰撞而产生的流体力学波传播的物理特性。离散时空直接数值模拟的建模不仅能够在极限流域给出相应的流动控制方程的解，例如在稀薄流域内的动理学方程解和在连续流域内的N-S/Euler方程的解，更重要的是能够在过渡流域内能够准确地描述流动物理，填补了过渡流域内缺乏合理流动控制方程的空白。由于在网格尺度上直接建模，网格大小变成了影响气体运动模式的动力学量，决定了在不同流域上的气体演化过程。基于直接建模的思想，徐昆与他的合作者提出和发展了统一气体动理学格式(unified gas-kinetic scheme, UGKS)[24, 63-70]。统一气体动理学格式一方面耦合宏观守恒量和微观分布函数的更新，另一方面由网格界面处的演化解构建数值通量，实现了从动理学方程向流体力学方程的自然过渡，能够准确捕捉多尺度复杂流场的流动物理。通过采用隐式加速和多重网格技术[71-73]，统一气体动理学格式能够高效准确应用于三维工程问题的模拟。如图2 、图3所示，统一气体动理学格式能够准确模拟临近空间飞行器周围的多尺度流场[74]。统一气体动理学格式在构造数值格式的过程中，不仅考虑了流场的局部努森数，同时考虑了网格尺度上的流动物理，即网格努森数Knc[75-76]

图2 江定武应用统一气体动理学格式对临近空间飞行器周围流场地模拟得到的流场局部努森数分布[74](Ma∞=4，Re=59 373)
Fig.2 Local Knudsen number of flow around a near space aircraft simulated by Jang[74], with Mach number 4 and Reynolds number 59 373

(2)

其中：τ为局部的碰撞时间，Δt为局部时间步长。统一气体动理学数值通量中的流体力学部分和自由输运部分的所占比例由网格努森数决定，实现了数值格式对网格尺度上的流动物理的自适应(见图3)。基于直接建模思想，研究者们发展了一些统一气体动理学格式的变形格式，包括郭照立等人发展了离散统一气体动理学格式(DUGKS)[77-87]，刘畅和朱亚军等发展了统一气体波粒方法(UGKWP)[88-90]等。在过去十年的研究中，统一气体动理学格式被应用到气体输运、等离子体输运[91]、离散两相流[92]、光子输运[93-96]、中子输运[97]等领域。统一气体动理学格式准确地连接了不同流域的物理过程，比如在光子输运中给出了从光学薄的光子自由传输到光学厚的光扩散过程的连续过渡，在计算精度和计算效率都体现出显著优势。

图3 多尺度格式对于不同网格努森数的自适应Fig.3 Adaptation of the multiscale scheme to local cell Knudsen number

接下来，我们将在第1节中介绍离散时空直接建模的思想，并总结基于离散时空建模构建的多尺度数值格式：统一气体动理学格式、离散统一气体动理学格式、统一气体动理学波粒格式。在第2节中，我们将总结统一气体动理学格式在多尺度气体、等离子体、气固离散两相流、光子输运等领域的应用，以及其隐式加速机制。在第3节中，我们将给出总结和展望。

1 离散时空直接建模方法

离散时空直接建模方法是在数值离散控制体上基于守恒律构建数值控制方程，构建数值控制方程的过程中需要考虑流域变化对控制体界面数值通量的影响。数值通量的构造依据在网格尺度上的演化过程，最为直接的是利用动理学微分方程在网格界面局部的积分解，即考虑了气体演化过程的时间累积效应，这是一种基于偏微分方程演化解的建模方法，而不是传统意义上对偏微分方程的直接离散。在某种意义上说，这种构造数值格式的思想很接近于Gounov方法，即把在网格界面上的演化解直接用于格式构造。而超出Godunov格式的地方包括演化解有时间精度和它的多尺度性质。我们将相空间划分为离散的物理空间∑iΩxi⊂R3，离散的速度空间∑jΩvj⊂Rdv和离散的时间tn∈R+。对于由物理空间控制体Ωxi和速度空间内控制体Ωvj张成的相空间内的有限控制体Ωij=Ωxi⊗Ωvj，在时间区间[tn,tn+1] 内的微观控制方程为

(3)

其中

为速度分布函数f(x,t,v)和碰撞项Q(x,t,v)的控制体平均值。f∂Ωij(t,vj)是控制体界面∂Ωxi上随时间演化的分布函数。对微观控制方程速度空间求矩，可以得到宏观守恒量控制方程

(4)

对应于网格尺度和时间步长，Knc决定了气体通过界面的物理过程。当Δt≤τ时，气体分子以自由输运的方式通过界面，而在Δt≫τ时，每个粒子在一个时间步长内经历了多次碰撞，使其以集体波的形式通过界面。而在中间区域，τ与Δt的比值决定了气体动力学演化。虽然方程(4)演化宏观量，但通过边界的流量是由时间演化的分布函数的积分得到，它的演化从自由分子流到连续流都是精确的，完全超出了直接定义的宏观演化方程，比如Burnett、Grad等矩方程。理论上方程(3)和(4)是两个完全独立的方程，在动力学上又耦合在一起，如何封闭方程(3)和(4)取决于在网格边界上气体演化过程和网格内部在一个时间步长里粒子碰撞的直接建模。

1.1 统一气体动理学格式 (UGKS)

在离散空间直接建模给出的控制方程(3)和方程(4)框架下，基于构建网格界面处的演化解来求得数值通量，从而得到统一气体动理学格式。统一气体动理学格式的有限体积更新框架为

(5)

(6)

其中，ls∈∂Ωxi为控制体界面，界面重心xs处外法向为ns，分布函数的数值通量为

(7)

宏观守恒量的界面通量Fs为

(8)

令tn=0，气体在界面x0处的演化是根据动理学方程的积分解给出

(9)

这个演化解直接连接了自由分子流f0到平衡态g的演化过程，其权重e-t/τ对从动理学尺度到流体力学尺度的建模起了重要的作用。如果想得到在τ时间尺度上更精确的物理过程，两体碰撞的Boltzmann碰撞过程也可以融合在上面的演化解里面[64]。为了得到二阶时间精度的数值通量表达式，我们将初始t=0时刻的速度分布函数f0(x,v)和平衡态分布函数g(x,t,v)在时空展开为

f0(x,v)=f0+xf0·(x-x0),

g(x,t,v)=g0+xg0·(x-x0)+∂tg0t

(10)

其中f0=f0(x0,v)，g0=g(0,x0,v)为初始时刻界面处的初始分布和平衡态分布。平衡态分布以及其空间导数由数值重构的宏观量及宏观量导数求得

(11)

而平衡态分布函数的时间导数由相容性条件得到

(12)

(13)

可以看到，通量包含由平衡态分布贡献的平衡态通量和初始分布函数贡献的自由输运通量两部分构成。 对微观分布函数求矩，我们得到宏观守恒量的数值通量表达式

(14)

其中，时间积分系数C1-5为网格努森数Knc的函数

C1=1-C4

C2=-τC1-C5

C3=Δt/2-τC1

(15)

值得注意的是上述流量通过Knc直接取决于网格大小，即网格大小影响到气体通过界面的演化过程。在连续流域Knc→0，若分布函数初值满足N-S分布

f0(v,x)=[g-τ(∂tg+v·xg)]|t=0

(16)

数值通量将给出与N-S方程一致的通量, 例如分布函数通量收敛到

(17)

对于统一气体动理学格式(UGKS)，即使f0不是N-S分布，在连续流当Knc≪1时，UGKS在网格界面上的分布函数也会回到式(17)所给出的分布，即f0的贡献随时间迅速衰减。在自由分子流域Knc→∞，数值通量(13)将给出与无碰撞动理学方程一致的通量

(18)

在上面两个极限的中间流域，f0和g共同起作用，它们的贡献由网格努森数决定。

1.2 离散统一气体动理学格式(DUGKS)

与UGKS的框架一致， 离散统一气体动理学格式(DUGKS)采用基于有限体积的分布函数演化方程(6)，并通过求矩完成对宏观守恒量的更新。DUGKS与UGKS的主要区别在于数值通量的求解。DUGKS采用梯形积分公式对碰撞项的时间积分进行近似，在时间区间[tn,tn+1/2]内

f(xs,v,tn+1/2)-f(xs-vΔt/2,v,tn)=

(19)

(20)

(21)

与UGKS相同， 在连续流域Knc→0，若分布函数初值满足N-S 分布，DUGKS的数值通量将给出与N-S方程一致的通量(17)，而在自由分子流域Knc→∞，数值通量(20)将给出与无碰撞动理学方程一致的通量(18)。

1.3 统一气体动理学波粒格式(UGKWP)

UGKS和DUGKS都是基于离散坐标方法发展的确定性数值方法，在离散空间直接建模的框架下，我们采用统计学方法，用模拟粒子离散速度空间，发展了统一气体动理学波粒方法[88-90]。 UGKWP方法的关键想法在于，将在流场演化过程中，所有粒子可分为自由输运粒子和碰撞粒子，对于发生碰撞的粒子不需要严格解析其碰撞过程，而可以直接由动理学方程的积分解给出其最终满足的速度分布，从而实现粗粒化多尺度建模。 UGKWP的粒子演化方程为动理学方程积分解，宏观量演化是与UGKS相同的有限体积更新方程。演化过程可以由循环图4表示。

UGKWP的粒子演化方程为

(22)

其中碰撞粒子分布为

[∂tg(x,t,v)+v·xg(x,t,v)]

(23)

在UGKWP的粒子更新过程中，我们只追踪粒子的自由输运过程，并从更新的宏观量中依据碰撞粒子分布g+采样出下一时间步内自由输运的粒子。需要指出的是在UGKWP重采样过程中，我们只需要采样下一时间步自由输运的粒子，其余部分速度分布由解析的速度分布函数表示，因而UGKWP的速度分布由解析部分和粒子部分组成，见图4(c)。 UGKWP的宏观量演化是与UGKS相同的有限体积更新方程

(24)

其中，平衡态通量求解与UGKS平衡态通量求解过程相同，即

(25)

而UGKWP的自由输运通量由解析自由输运通量

(26)

和粒子自由输运通量

(27)

Fig.4 Sampling process of UGKWP in a numerical control volume:(a) Initially, the particles are categorized into collision-less (white) and collisional (black) particles according to their free stream timetf, (b) Stream all particles bytfand keep collision-less particles, (c) Re-sample collisional particles fromg+.

1.4 三种直接建模方法的性质和比较

基于直接建模方法发展的三种数值方法，即UGKS、DUGKS和UGKWP都具有二阶一致渐进保持(unified preserving, UP) 的性质，即格式不仅能够在连续和无碰撞的极限流域收敛到欧拉和无碰撞动理学方程的数值格式，同时在连续流域内，格式能够准确捕捉N-S方程的数值解[64, 98]。 特别的，郭照立在最近的文章中给出了UP的数学定义，即在一定条件下，数值格式对应的修正方程能够保持动理学方程的二阶渐进展开。 以一维为例，我们给出UGKS在Δt/τ≫1条件下，采用空间中心差分重构的修正方程

∂tf+v∂x[g-τ(gt+vgx)]

(28)

对修正方程进行渐进分析，我们将分布函数f展开为

f=f(0)+εf(1)+ε2f(2)+O(ε3)

将时间导数展开为

∂t=∂t0+ε∂t1+ε2∂t2+O(ε3)

其中，∂tk表示εkf(k)的空间梯度对时间导数∂t的贡献[98]。将展开式代入方程(28)，我们可以得到UGKS在连续流域的前两阶渐近展开为

=-τf(2)+O(Δt2,Δx3)

(29)

即UGKS在连续流域给出时间二阶精度的N-S方程数值格式。需要特别指出的是，从UGKS的修正方程(28)可以看出，在小网格努森数(Knc)流域 UGKS的误差主要来源于平衡态通量，这是由于在小网格努森数流域UGKS的平衡态通量占主导。关于DUGKS和UGKWP格式UP性质的论证，参见文献[90, 98]。

UGKS和DUGKS的主要区别在于数值通量系数的不同，尤其是在小网格努森数时差别较大。图5为UGKS和DUGKS通量系数的差别

(30)

可以看到，在小网格努森数流域UGKS的平衡态通量占主导，而DUGKS的连续通量和自由输运通量之比为C1/C4=-2。 因而要提高UGKS在小网格努森数流域的精度，即在连续流的精度，只需要提高平衡态对应的重构精度和时间演化精度，即格式完全有宏观量决定，这个性质使其回到连续流的GKS，虽然这时UGKS还保持着不必要的速度空间网格点。不同的是，要提高DUGKS在连续流域的精度，需要同时提高平衡态和分布函数的重构精度和时间精度， 即它们的结合给出DUGKS在连续流时的数值通量。另一方面由于DUGKS的通量计算相对简洁，因此DUGKS相对于UGKS计算量有一定程度的降低。

基于离散速度的确定性格式UGKS和DUGKS，由于不受统计噪音的影响，因而对于微流动问题的模拟有很高的精度和计算效率。而基于统计学的UGKWP，采用模拟粒子能够高效地实现速度空间的自适应，因而对于高超声速和三维工程问题的计算由较高的计算效率。同时UGKWP只需要采样自由输运粒子，因而UGKWP计算所需粒子数随着网格努森数的降低而以指数减小，计算控制体Ωxi内所需采样的粒子数为

(31)

其中，mp为模拟粒子的参考质量。另一方面，随着网格努森数的降低，UGKWP自由输运通量在总通量中所占比例以指数降低，因此统计噪音也以指数降低。在参考文献[90]中我们证明了UGKWP不仅是多尺度格式同时也是计算复杂度渐进降低和统计噪音渐进降低的格式。 在连续流，由于自由粒子的消失，UGKWP完全回到以前计算N-S方程的气体动理学方法(GKS)[24]。和UGKS比较，UGKWP没有了速度空间，在连续流区内存也大大降低。

总体来讲，基于离散时空直接建模而发展起来的三种格式都具有多尺度性质，不仅对于气体输运，还能够很好地模拟等离子体输运、气固离散两相流、辐射输运、中子输运等等尺度输运问题。

(a)

(b)

图5 统一气体动理学格式和离散统一气体动理学格式平衡态通量和自由输运通量占比(a)、统一气体动理学格式和离散统一气体动理学格式通量系数差别(b)随网格努森数的变化

Fig.5 Proportion of equilibrium flux and free stream flux of UGKS and DUGKS (a)、Flux coefficients of UGKS and DUGKS (b) according to cell Knudsen number

2 统一气体动理学格式的应用

2.1 连续层流及高超声速问题

统一气体动理学格式的首次提出是用于多尺度气体输运问题的计算[50]。在之后的十年中不断发展，在多尺度气体流动、微流动、高超声速气体输运问题的模拟以及三维工程问题的计算中，体现出精度、稳定性以及计算效率等方面的优势[50, 63-64, 71-74, 88-90]。 在这一章节中，我们一方面将通过UGKS在连续层流问题的算例说明UGKS的准确性，另一方面通过UGKWP在高超声速问题的算例说明格式的计算效率。

首先，我们考察二维的Taylor涡，计算域为0

(32)

其中A=B=2π，α=A2+B2，u0=0.01，Kn≪10-4，RT0=0.5网格尺度Δx=2.5×10-2>Kn0.5，CFL数为0.5。由第1节的讨论，UGKS能够在该流域收敛到N-S方程数值格式。图6(a)给出不同格式在时间t=ln2/(μα)的解，我们可以看出二阶UGKS在该流域能够收敛到对N-S求解的二阶GKS格式的数值解。同时，我们只需要将UGKS平衡态通量对应的空间重构提高，就可以提高UGKS 在连续流域的计算精度。另一方面，我们看到由于AP-IMEX 格式[58]的单尺度迎风通量引入过大的数值耗散，因而即使对分布函数采用空间高阶重构，依然无法提高其在N-S流域的精度。

(a) 沿y方向ux的值

(b) uy速度云图

图6 (a) 不同格式对Taylor涡的模拟结果，t=ln2/(μα)，x=0.5;(b) UGKS(实线)和精确解(云图)的uy速度云图,t=ln2/(μα)

Fig.6 (a) Velocity profile of the Taylor vertex att=ln2/(μα),x=0.5; (b) Velocity contour of the Taylor vertex att=ln2/(μα), The solid line is UGKS solution and flood is the analytic solution.

对于高超声速问题，我们给出UGKWP对经典的二维圆柱扰流问题的计算界结果[90]。分子模型采用氩气的VHS模型，来流马赫数为Ma=20、30，以圆柱半径为参考长度的来流努森数为Kn=1.0、1.0×10-4。图7和图8分别给出Ma=20不同努森数下流场温度云图和中轴线上的云图分布，可以看出UGKWP在不同流域都能与参考解吻合。图9给出Ma=30不同努森数下流场温度云图。对于Ma=20,Kn=1.0的算例，UGKWP物理空间网格为径向110网格，周向64网格，粒子参考质量为mp=1.0×10-3，计算时间消耗为36 min，内存消耗为100 MB(i7-8700K CPU)比传统DSMC方法的计算量更低。而对于Ma=20,Kn=10-4的算例，UGKWP物理空间网格为径向150网格，周向100网格，粒子参考质量为mp=4.7×10-3，计算时间消耗为17.2 min，基本与连续流域N-S求解器的计算效率相当。

2.2 多尺度等离子体模拟

统一气体动理学格式还可以应用于多组分气体[99]，以及在外力场作用下的流体[100]。这里我们以等离子体为例来介绍UGKS在多尺度等离子体输运中的应用[91]。 等离子体的流域可以根据等离子体和电磁场的耦合强度参数，即归一化回旋半径rL，划分为双流体流域和磁流体流域，同时可以根据等离子体组分间的耦合强度参数，即努森数，划分为无碰撞等离子体流域和连续流域。我们提出了包含等离子体和电磁场以及等离子体间碰撞的动理学耦合电磁学方程组

(a) Kn=1.0

(b) Kn=1.0×10-4图7 马赫数20的圆柱绕流问题的温度分布云图 (UGKWP 解为黑线，参考解为云图)

Fig.7 Temperature contour of the cylinder flow with Mach 20 (solid line is UGKWP solution and flood is the reference solution)

(a) Kn=1.0

(b) Kn=1.0×10-4图8 马赫数20的圆柱绕流问题的中轴线温度分布 (符号为UGKWP解，实线为参考解)

Fig.8 Temperature profile along the stagnation line with Mach 20 (symbols are UGKWP solution and solid ine is the reference solution)

(a) Kn=1.0

(b) Kn=1.0×10-4图9 UGKWP马赫数30的圆柱绕流问题的温度分布云图Fig.9 Temperature contour of cylinder flow with Mach 30

(33)

其中fα为组分α的分布函数，τα为松弛系数，Fα为洛伦兹力，E、B为电磁场，J为电流。 在数值离散时空控制体上，根据直接建模方法，我们建立等离子体流场微观分布函数的有限体积格式

(34)

和流场宏观量的有限体积格式

(35)

以及电磁场的有限体积格式

(36)

格式的核心在于采用网格界面动理学方程的积分解构建数值通量

(37)

统一气体动理学格式能够捕捉从无碰撞等离子体到磁流体在不同流域的等离子体流动物理，这里我们给出UGKS对朗道阻尼，Orszag-Tang涡问题，和磁重联问题的模拟。 朗道阻尼是一种由电磁场作用引发的无碰撞阻尼现象，由于流场扰动较小，需要格式具有较高的计算精度。 电子的初始扰动为

(38)

其中k=0.5，扰动α=0.01对应线性朗道阻尼，较大扰动α=0.5则会产生非线性效应。我们采用UGKS模拟计算域L=2π/k内的线性和非线性朗道阻尼现象，采用速度空间128个离散点，物理空间128个离散点。如图10所示，UGKS能够准确捕捉无碰撞等离子体的精细结构。

(a) 对线性朗道阻尼问题的模拟

(b) 对非线性朗道阻尼问题的模拟图10 UGKS对线性和非线性朗道阻尼问题的模拟结果Fig.10 The UGKS simulation of linear and nonlinear Landau damping comparing to the theoretical solution

为了验证UGKS在磁流体流域的数值表现，我们计算了Kn=10-4,rLi=10-4的Orszag-Tang问题，计算的初始条件为

ni=ne=γ2,Pi=Pe=γ,By=sin(2x)

ui,x=ue,x=-siny,ui,y=ue,y=sinx

其中γ=5/3，物理空间离散为200×200个网格，速度空间采用28×28的Gauss积分。图11 给出t=3时刻的压力云图以及沿y=0.625π的x方向的压力分布。可以看出UGKS在连续流域能够准确捕捉磁流体方程的解。UGKS不仅能够捕捉极限流域的物理解，同时能够为过渡流域内较为复杂的物理问题提供有效的研究方法，例如对于磁场重联问题，UGKS能够在连续流域恢复双流体模型的结果，在无碰撞流域恢复Vlasov方程的结果，同时能够给出在适中努森数和无量纲回旋半径流域的重联结果，见图12，UGKS为多尺度等离子体物理问题的研究提供有力的数值工具。

图11 (a) Orszag-Tang问题在t=3时刻的压力云图；(b) Orszag-Tan问题在t=3时刻沿y=0.625π的x方向的压力分布(实线)与磁流体方程结果(符号)的比较

Fig.11 (a) Pressure contour of Orszag-Tang vertex problem att=3; (b) Pressure profile alongy=0.625π att=3, solid line is UGKS solution and symbol is reference solution

图12 (a)Kn=10-4,rLi=1的磁重联问题在t=40的磁力线分布；(b) UGKS给出的不同无量纲回旋半径的磁场重连率以及参考解

Fig.12 (a) Magnetic lines of magnetic reconnection problem att=40 withKn=10-4,rLi=1; (b) Reconnection rate of UGKS with different Larmor radius comparing to reference results

2.3 多尺度气固离散两相流模拟

多尺度气固两相流也是一种典型的多尺度输运系统，可以由固体颗粒相的动理学方程和气相的流体力学N-S方程组成的方程组来描述，流域可以由颗粒碰撞频率的强弱划分为无碰撞气固两相流和双流体流域；根据颗粒跟随系数可以划分为颗粒流(granular flow)流域和尘埃流域(dusty flow)。气固离散两相流的控制方程为颗粒相的动理学方程

(39)

(40)

和气相的N-S方程组成

-μgσ(Ug)U+κgxTg]=

(41)

在方程(39)中,g为重力加速度，ρs为颗粒相的物质密度;pg为气相压力;ms为颗粒相单个颗粒的质量;下表s表示颗粒相相关参数; 下表g表示气相相关参数。 气相和颗粒相的相互作用通过组分间作用力实现

(42)

(43)

其中平衡态满足

(44)

参数r描述非弹性碰撞的能量损失率。

基于离散空间直接建模方法，我们构建由微观和宏观演化方程耦合的颗粒相的UGKS格式[92]。其中微观分布函数的数值演化方程为

(45)

(46)

(47)

对应的，颗粒相宏观量的演化方程为

(48)

(49)

(50)

(51)

(52)

(53)

其中源项为

速度和能量方程中的参数为

颗粒相UGKS格式多尺度性质的核心在于其在微观分布函数演化方程(45)和宏观量的演化方程(48)中的数值通量由动理学方程(39)界面处积分解构建

e-t/τsfs,0(x0-ut,vj-ω1t)

(54)

因而可以在稀薄颗粒流和连续流域分别给出颗粒自由输运和双流体方程对应的通量。气相的演化方程为基于GKS的有限体积格式[24]

Lgw1:

(55)

(56)

Lgw3:

(57)

γTg3=

γTg4=

气相方程的数值通量由基于网格界面积分解的GKS通量给出，其中网格界面积分解为

e-t/τf0,NS(x0-vt,v-gt)

(58)

与颗粒相不同，气相网格界面积分解的初值由N-S方程对应的一阶CE展开f0,NS给出[24]。颗粒相和气相的演化相耦合，构成描述气固离散两相流的多尺度数值格式UGKS-Multiphase[92]。

这里我们给出UGKS-Multiphase对激波驱动的颗粒层实验的数值模拟，实验装置设置如图13所示，颗粒层初始位置为150 mm，在颗粒层下方110 mm，颗粒层上方43 mm和718 mm处有三个压力感应器。初始激波强度为Ma=1.3，气相密度1.2 kg/m3，颗粒相密度2500 kg/m3。计算数值分辨率为Δx=3.75×10-3，Δv=4.4375。我们分别计算了不同颗粒层流域的情形，包括稀薄流域颗粒层厚度2 mm 和连续流域颗粒层厚度20 mm的情况，都能够和实验数据相吻合,见图14。

图13 激波驱动的颗粒层实验装置设置Fig.13 Set up of shock driven particle bed experiment

(a) 稀薄流域颗粒层厚度2 mm

(b) 连续流域颗粒层厚度20 mm图14 激波驱动的颗粒层实验UGKS的数值 模拟结果和实验数据的对比Fig.14 Comparison the pressure profile of pressure driven particle bed

2.4 多尺度光子输运模拟

根据介质的透光性，光子输运可以分为在光性薄介质的稀薄流域、过渡流域和在光性厚介质的扩散流域。传统的隐式蒙卡方法、离散速度方法和扩散方程等方法，都是单尺度数值方法，不能够在全流域适用。根据离散时空直接建模方法，我们发展了基于离散速度的UGKS格式[93-94, 96]和基于统计学粒子的UGKWP 格式[89]。这里我们以纯散射介质内的UGKWP格式为例介绍多尺度光子输运数值方法的建模思想。 纯散射介质光子输运的动理学方程为

(59)

(60)

其中

(61)

宏观守恒量的演化方程为

(62)

UGKWP的核心在于采用网格界面处动理学方程积分解构造数值通量，即

(63)

我们给出采用UGKWP计算的Marshak波和线光源传播问题。对于Marshak波问题，吸收系数设为σ=30/T3，光速c=29.98，物质比热容Cv=0.3。图15给出时间t=0.33,0.66,1.0时刻的光子和物质能量分布。由于Marshak问题的吸收系数与温度相关，低温区为扩散流域，高温区为稀薄流域，UGKWP在不同流域都能够和参考解吻合，同时计算效率高于隐式蒙卡。对于线光源问题，计算参数为σ=1,ε=1,c=1。计算域为[-1.5,1.5]×[-1.5,1.5]，物理空间网格量为201×201。图16给出UGKWP、S8、隐式蒙卡的结果，可以看到UGKWP与参考解吻合，同时没有Sn方法的非物理射线效应。

基于离散空间直接建模方法的UGKS格式还被应用中子输运[97]、双原子气体[101]、多组分气体[102]等领域，DUGKS格式被成功应用于微流动[103-104]、两相流及固液相变[105-107]、混合气体流动[108]、颗粒流动[107]、电渗流动[109]、热辐射[110]等领域，在计算精度和计算效率上相对于传统单尺度方法由体现出明显优势。 DUGKS还被应用于多尺度声子输运的模拟，能够准确模拟从弹道输运流域(ballistic transport regime)到扩散流域(diffusive regime)的声子输运物理[111-112]。在体元尺度渗流和气固颗粒两项流等领域，DUGKS也有具有较高的准确性和计算效率[83, 86]。另一方面，GKS和DUGKS格式也在湍流的直接模拟和湍流建模方面提供了新的思路和手段，有望为湍流建模开辟一条更为有效的研究道路[68, 113-115]。

(a) 辐射温度

(b) 物质温度图15 Marshak问题t=0.33,0.66,1.0时刻 UGKWP的计算结果与IMC方法的比较Fig.15 UGKWP simulation of Marshak problem compared to IMC result at t=0.33,0.66,1.0

2.5 加速机制

在近年的发展中，我们提出了一系列针对UGKS的加速机制，包括速度空间自适应[116]和隐式和多重网格加速方法[71-73]。 这里我们简要介绍UGKS的隐式加速方法。隐式UGKS的宏观迭代控制方程为

(64)

微观迭代方程为

(65)

图16 线光源问题t=1.2的解Fig.16 Line source scattering problem at t=1.2

通过采用多重网格等迭代技术可以快速得到隐式控制方程的解，相对于显式推进的UGKS，计算时间可以提高至少两个数量级。 隐式UGKS在不同流域都快速收敛的关键在于宏观和微观耦合迭代，在连续流域通过宏观迭代驱动微观分布函数快速达到收敛，近年来基于这种想法发展了一些高效的多尺度隐式算法[56, 117]。

3 多尺度建模方法的总结和展望

传统的偏微分方程数值解往往是从偏微分方程出发，通过一定的离散方法得到相应的数值格式，但偏微分方程都是对固定物理尺度建模得到的方程，因而对其直接的离散不能改变其物理本质，因而只能发展出单尺度的方法，比如针对N-S方程的各种算法和直接离散Boltzmann方程的DVM方法。离散时空直接建模思想是回到有限控制体上，把有限大小的控制体当成物理尺度，建立对应于此尺度上的离散演化方程，以及通过界面的数值通量。此格式的构建完全根据物理守恒率，以及在不同网格尺度上的传输过程。格式直接给出了在离散空间上物理量的演化方程和它的演化解。根据网格尺度和分子平均自由程的随地方变化的不同比值，实现了对不同流域的物理描述和直接的数值自适应，即在连续流域收敛到N-S方程数值解，在稀薄流域收敛到动理学方程的解。在中间流域，由于多尺度建模中严格遵从物理守恒率和从非平衡态到平衡态满足熵增的演化规律，方法本身是可靠的。在所有算例中还没有发现过违背物理规律的数值解。基于这种思想，我们还发展了包含气体、等离子体、光子、中子、离散两相流在内的多种输运过程的数值多尺度控制方程，并且在计算精度和效率上相对于传统单尺度格式体现出明显优势。离散时空直接建模思想不仅为物理研究和工程应用提供了有力的数值方法，同时为包含湍流建模在内的建模难题提供了一个新的思路。直接建模思想有深刻的物理基础和独特的建模手段，跳出了传统数值计算的制约，把建立演化方程和得到数值解紧密地结合起来。直接建模为解决多尺度问题和刻画非平衡传输提供了一种重要的解决方案，在未来输运问题研究以及湍流建模的研究中必将发挥重要作用，成为多尺度建模的重要工具。