稀薄气体动理论中介观尺度数值模拟的加速方法

杨伟奇; YANG Weiqi

doi: 10.11887/j.cn.202402007

杨伟奇

国防科技大学空天科学学院, 湖南长沙 410073

基金项目: 国家自然科学基金资助项目（12302382，U1730247）；湖南省自然科学基金青年基金资助项目（2022JJ40542）

详细信息

作者简介

杨伟奇（1989—），男，安徽宿州人，讲师，博士，E-mail: yangweiqi@nudt.edu.cn

中图分类号: O354

文献标识码: A

文章编号: 1001-2486(2024)02-070-09

Acceleration scheme for the mesoscopic numerical simulation in the kinetic theory of rarefied gas

YANG Weiqi

College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073 , China

摘要

针对滑移流、早期过渡流区域采用离散速度法(discrete velocity method, DVM)求解Boltzmann方程收敛速度极慢、计算资源消耗大的难题，提出全流场耦合介观/宏观方程加速方法。在介观层面基于有限差分的DVM求解Boltzmann方程，在宏观层面基于有限体积的压力耦合方程组半隐式方法求解矩方程，充分利用纳维-斯托克斯-傅里叶/R26矩方程在低克努森数下的快速收敛特性，对介观方程进行加速。基于高阶Hermite多项式重构分布函数，完成宏观方程与介观方程之间的数据传输。仿真结果表明：全流场耦合介观/宏观方程的加速方法在滑移流、早期过渡流区域具有显著加速性能，最大降低了95.28%的计算时长；但对中、大克努森数流域，加速性能大幅度下降。

关键词

稀薄气体 / Boltzmann方程 / 离散速度法 / R26矩方法 / 加速方法

Abstract

To overcome the difficulty that DVM(discrete velocity method) for solving the Boltzmann equation has extremely slow convergence speed and high computational resource consumption in the slip flow and early transition flow regimes, an acceleration scheme which coupling the mesoscopic/macroscopic equations in the full flow area was proposed. Boltzmann model equation could be solved based on the DVM using finite difference method at the mesoscopic level, and the moment equations could be solved based on the semi-implicit method for pressure linked equations using finite volume method at the macroscopic level. Fast convergence characteristic of the Navier-Stokes-Fourier/R26 moment equations in the low Knudsen number regime was fully utilized to accelerate mesoscopic equation. The distribute function could be reconstructed from the high-order Hermite polynomial function, so that the data transfer between the macroscopic and the mesoscopic equations could be done. Simulation results indicate that: the acceleration scheme, coupling the mesoscopic/macroscopic equations in the full flow area, shows great acceleration performance in the slip and early transition regime, which is able to save up to 95.28% computational time cost. However, the acceleration performance decreases significantly in the middle and large Knudsen number regime.

Keywords

rarefied gas / Boltzmann equation / discrete velocity method / R26 moment method / acceleration scheme

1 Boltzmann方程和速度分布函数的矩 1.1 Maxwell漫反射边界条件 1.2 速度分布函数的高阶矩 2 拓展流体力学方程 3 全流场耦合DVM/R26加速方法 4 算例与讨论 4.1 加速度驱动Poiseuille流动 4.2 顶盖驱动方腔流动 5 结论

经典纳维-斯托克斯-傅里叶（Navier-Stokes-Fourier，NSF）方程组在连续流与近连续流（Kn＜10^-3）区域能够较好地描述气体的输运特性。当Kn增大时，如滑移流区域（10^-3≤Kn≤10^-1），气体分子的非平衡效应开始变得显著，此时NSF方程组不再准确^[1]。传统的做法是将NSF方程组结合速度滑移、温度跳跃边界条件，用以近似描述气体在滑移流的非平衡效应^[2]。此外，以Burnett/超Burnett方程以及R13/R26矩方程为代表的拓展流体力学方程，也主要适用于滑移流与早期过渡流^[3-6]（Kn≤1.0）。当Kn继续增大至过渡流（10^-1＜Kn≤10）与自由分子流（Kn＞10）时，必须引入分子动理论与Boltzmann方程，在介观层间描述气体分子之间的输运以及碰撞特性。

介观Boltzmann方程是7维积分-微分方程，同时方程中含有复杂的非线性碰撞项，精确求解全维Boltzmann方程难以实现。目前Boltzmann 方程的求解方法主要分为两大类：①以直接模拟蒙特卡罗（direct simulation Monte Carlo，DSMC）为代表的概率密度方法^[7]；②以离散速度法（discrete velocity method，DVM）为代表的直接求解Boltzmann法^[8]。大量研究^[9-15]表明，经典DSMC与DVM在近连续流、滑移流与早期过渡流对计算资源消耗巨大，迭代收敛速度极慢。

研究人员在近20年致力于降低Boltzmann方程的计算资源消耗，并提高计算效率^[12-13]。首先，研究人员提出简化的Boltzmann模型方程，如BGK^[16]（Bhatnagar-Gross-Krook）、ES-BGK^[17]，Shakhov模型方程^[18]，将Boltzmann方程中粒子之间的非线性碰撞项简化为分布函数随时间的松弛来降低计算量级。同时，研究人员提出多种Boltzmann模型方程的加速算法，Li等提出气体动理论统一算法（gas-kinetic unified algorithm，GKUA）来解决高超声速来流条件下再入飞行器的数值模拟问题^[12]，Xu等和Guo等分别提出统一气体动理学格式^[13]（unified gas-kinetic scheme，UGKS）和离散统一气体动理学格式^[14]（discrete unified gas-kinetic scheme，DUGKS）等渐进保持迭代格式，将Boltzmann方程数值求解的适用范围推广到了近连续流。Wang等^[15]对比经典DVM与DUGKS在不同Kn下的计算效率和计算时长时认为：经典DVM在较大Kn时计算效率较高，而在近连续流与早期过渡流时，DUGKS收敛速度较DVM更快。进一步地，Zhu等提出隐格式UGKS^[19-20]与DUGKS^[21]，相对原显格式具有显著的加速效果，并可以适用于全流域的稳态问题求解。

近期，Su等提出了介观方程的加速方法，基于宏观方程对小Kn问题进行加速收敛^[22]。Yuan等借鉴DUGKS思想，引入独特的宏观参数预测技术，不但能够满足全流域稳态问题的求解要求，还能够满足全Kn下快速收敛要求^[23-25]。目前宏观参数的预测大多基于NS方程组，本文针对R13/R26矩方法，为介观尺度加速方法提供新的思路。前期，作者首次将宏观层面的矩方程与介观Boltzmann方程耦合起来，最大降低了99%的计算内存，提高了R26矩方程组的计算精度^[8]。同时，作者成功将耦合介观/宏观方程应用到热驱稀薄流动的数值模拟中^[9]。壁面处耦合DVM/R26矩方程组有效地降低了计算资源消耗，具有较大的应用潜质，但对于计算效率与计算时长消耗的优化能力有限。因此，本文针对壁面处耦合算法的局限性，提出全流场耦合迭代介观/宏观方程，并研究耦合方法的计算效率与适用范围。

1 Boltzmann方程和速度分布函数的矩

介观尺度Boltzmann-Shakhov模型方程^[7]可以写为：

\frac{\partial f}{\partial t} + C \cdot \frac{\partial f}{\partial x} + a \cdot \frac{\partial f}{\partial C} = - \frac{1}{τ} (f - f^{e q})

(1)

式中：C=（C_x，C_y，C_z）为气体分子离散速度；x=（x，y，z）、t分别为空间和时间； f=f（t，x，C）为气体分子速度分布函数；τ为平均松弛时间；a为外力加速度；f^eq为Shakhov平衡态速度分布函数，可由Maxwell平衡态速度分布函数加上热流修正项得到，如式（2）所示：

\begin{matrix} f^{e q} = \frac{ρ}{(2 π R T)^{3 / 2}} e x p [- \frac{(C - u)^{2}}{2 R T}] \\ \{1 + (1 - P r) \frac{(C - u) \cdot q}{5 p R T} [\frac{(C - u)^{2}}{R T} - 5]\} \end{matrix}

(2)

式中：ρ、u、T、q分别代表宏观参数中密度、速度矢量、温度与热流矢量；Pr为普朗特数，对于单原子分子普朗特常量Pr=2/3；R为理想气体常数；克努森数Kn定义为分子平均自由程λ₀与参考长度L₀的比值，如式（3）所示：

K n = \frac{λ_{0}}{L_{0}}

(3)

分子平均自由程λ₀与参考黏性μ₀的关系可以写为：

λ_{0} = \frac{μ_{0}}{p_{0}} \sqrt{\frac{π R T_{0}}{2}}

(4)

式中，p₀、T₀分别为参考压力与参考温度，服从理想气体状态方程p₀=ρ₀RT₀。基于离散速度法求解Shakhov模型方程的迭代格式^[15，26-27]如式（5）所示：

(\frac{1}{Δ t^{n}} + C \cdot \nabla + \frac{1}{τ^{n}}) Δ f^{n} = \frac{1}{τ^{n}} (f^{e q, n} - f^{n}) - C \cdot \nabla f^{n}

(5)

式中，

\nabla

代表差分格式，Δfⁿ=fⁿ⁺¹-fⁿ，上标n与n+1分别代表第n与n+1次迭代次数。在隐格式离散速度法中，首先忽略速度分布函数的时间导数项，考虑式（5）等号左边为时间隐格式，空间离散格式可以采用二阶迎风格式，同时控制式（5）等号右边为时间显格式，任意方向上分布函数的二阶迎风格式可以写为：

{\frac{\partial f^{n + 1}}{\partial x}|}_{x = x_{i}} = \{\begin{matrix} \frac{3 f_{i}^{n + 1} - 4 f_{i - 1}^{n + 1} + f_{i - 2}^{n + 1}}{2 Δ x}, C_{x} ⩾ 0 \\ \frac{- f_{i + 2}^{n + 1} + 4 f_{i + 1}^{n + 1} - 3 f_{i}^{n + 1}}{2 Δ x}, C_{x} < 0 \end{matrix}

(6)

其中：x_i为x方向上的网格点坐标，将式（6）代入式（5）并简化可得：

(\frac{1}{τ} + \frac{3 C_{x}}{2 Δ x}) f_{i}^{n + 1} = \{\begin{matrix} \frac{1}{τ} f_{i}^{e q, n} - g_{x} \frac{\partial f_{i}^{n}}{\partial C_{x}} + \\ C_{x} \frac{4 f_{i - 1}^{n + 1} - f_{i - 2}^{n + 1}}{2 Δ x}, C_{x} ⩾ 0 \\ \frac{1}{τ} f_{i}^{e q, n} - g_{x} \frac{\partial f_{i}^{n}}{\partial C_{x}} + \\ C_{x} \frac{f_{i + 2}^{n + 1} - 4 f_{i + 1}^{n + 1}}{2 Δ x}, C_{x} < 0 \end{matrix}

(7)

其中，g_x为流体驱动加速度，式（7）为基于经典DVM求解Boltzmann模型方程隐格式迭代格式。

1.1 Maxwell漫反射边界条件

Maxwell镜面反射/漫反射壁面边界条件认为：当气体分子碰撞壁面后，比例为α的气体分子通过漫反射的方式返回流场；比例为1-α的气体分子通过镜面反射的方式返回到流场。同时，在分子碰撞前后流场保持质量守恒。α=1和α=0分别代表全漫反射条件与全镜面反射情况，壁面处分布函数的边界条件f^w为：

f^{w} = \{\begin{matrix} α f_{M}^{w} + (1 - α) f (- C_{i} n_{i}), C_{i} n_{i} ⩾ 0 \\ f (C_{i} n_{i}), C_{i} n_{i} < 0 \end{matrix}

(8)

其中：C_i为气体分子离散速度C的张量表达形式；n_i为垂直于壁面的单位张量，指向气体；i为张量运算角标；f^w_M为Maxwell速度分布函数，如式（9）所示：

f_{M}^{w} = \frac{ρ^{w}}{\sqrt[3]{2 π R T^{w}}} e x p [- \frac{{(C_{i} - u_{i}^{w})}^{2}}{2 R T^{w}}]

(9)

式中，ρ^w为碰撞后的气体密度，u^w_i、T^w分别为壁面处速度与温度。

1.2 速度分布函数的高阶矩

当得到全流场气体分子速度分布函数f=f（t，x，C）时，气体密度ρ、速度u、温度T、应力张量σ_ij、热流q等宏观量可以通过分布函数对离散速度求矩获得，如式（10）所示：

\{\begin{matrix} ρ = \int f d C \\ u_{i} = \frac{1}{ρ} \int C_{i} f d C \\ T = \frac{2}{3 ρ} \int c^{2} d C \\ σ_{i j} = \int c_{i} c_{j} f d C - δ_{i j} \frac{ρ T}{2} \\ q_{i} = \frac{1}{2} \int c^{2} c_{i} f d C \end{matrix}

(10)

其中，

c_{i} = C_{i} - u_{i} ， c^{2} = \sum {(C_{i} - u_{i})}^{2}

，δ_ij为狄拉克函数，下角标i，j为张量计算角标。进一步可以得到R13与R26矩方程组中宏观参数的表达形式，其中R13矩方程组对应的高阶矩表达式为：

\{\begin{matrix} m_{i j k} = \int c_{< i} c_{j} c_{k >} f d C \\ R_{i j} = \int c_{< i} c_{j >} c^{2} f d C - 7 R T σ_{i j} \\ Δ = \int c^{4} f d C - 15 p R T \end{matrix}

(11)

其中，下角标“＜”“＞”为高阶矩张量运算符号，其运算法则可参考文献^[2，6]，在此不做赘述。

R26矩方程组对应的高阶矩表达式为：

\{\begin{matrix} ϕ_{i j k l} = \int c_{< i} c_{j} c_{k} c_{l >} f d C \\ ψ_{i j k} = \int c_{< i} c_{j} c_{k >} c^{2} f d C - 9 R T m_{i j k} \\ Ω_{i} = \int c^{4} c_{i} f d C - 28 R T q_{i} \end{matrix}

(12)

本文的无量纲格式为：

v_{0} = \sqrt{2 R T_{0}} ， \bar{u} = \frac{u}{v_{0}} ， \bar{x} = \frac{x}{L_{0}} ， \bar{t} = \frac{v_{0} t}{L_{0}} ， \bar{a} = \frac{L_{0} a}{v_{0}^{2}} ， \bar{C} = \frac{C}{v_{0}} ， \bar{T} = \frac{T}{T_{0}} ， \bar{τ} = \frac{v_{0} τ}{L_{0}} ， \bar{f} = \frac{v_{0}^{3} f}{ρ_{0}} ， \bar{ρ} = \frac{ρ}{ρ_{0}} ， \bar{p} = \frac{p}{p_{0}} ， \bar{μ} = \frac{μ}{μ_{0}}

。其中，v₀为参考速度，

\bar{u}

、

\bar{x}

等分别代表对应参数的无量纲量。

高阶矩的无量纲格式为：

{\bar{σ}}_{i j} = \frac{σ_{i j}}{ρ_{0} v_{0}^{2}} ， {\bar{q}}_{i} = \frac{q_{i}}{ρ_{0} v_{0}^{3}} ， {\bar{m}}_{i j k} = \frac{m_{i j k}}{ρ_{0} v_{0}^{3}} ， {\bar{ϕ}}_{i j k l} = \frac{ϕ_{i j k l}}{ρ_{0} v_{0}^{4}} ， {\bar{R}}_{i j} = \frac{R_{i j}}{ρ_{0} v_{0}^{4}} ， {\bar{ψ}}_{i j k} = \frac{ψ_{i j k}}{ρ_{0} v_{0}^{5}} ， {\bar{Ω}}_{i} = \frac{Ω_{i}}{ρ_{0} v_{0}^{5}} ， \bar{Δ} = \frac{Δ}{ρ_{0} v_{0}^{4}} 。

2 拓展流体力学方程

经典NSF方程的控制方程可以认为是质量、能量与动量三大守恒方程，如式（13）所示：

\{\begin{matrix} \frac{\partial ρ}{\partial t} + \frac{\partial ρ u_{i}}{\partial x_{i}} = 0 \\ \frac{\partial ρ u_{i}}{\partial t} + \frac{\partial ρ u_{i} u_{j}}{\partial x_{j}} + \frac{\partial σ_{i j}}{\partial x_{j}} = - \frac{\partial p}{\partial x_{i}} + ρ a_{i} \\ \frac{\partial ρ T}{\partial t} + \frac{\partial ρ u_{i} T}{\partial x_{i}} + \frac{2 \partial p_{i}}{3 R \partial q_{i}} = - \frac{2}{3 R} (p \frac{\partial u_{i}}{\partial x_{i}} + σ_{i j} \frac{\partial u_{j}}{\partial x_{i}}) \end{matrix}

(13)

在三大守恒方程中，切应力σ_ij和热流q_i可以通过Fourier导热定理得到，也可以通过Boltzmann方程推导其控制方程，如式（14）、式（15）所示：

\begin{matrix} \frac{\partial σ_{i j}}{\partial t} + \frac{\partial u_{k} σ_{i j}}{\partial x_{k}} + \frac{\partial m_{i j k}}{\partial x_{k}} \\ = - \frac{p}{μ} σ_{i j} - 2 p \frac{\partial u_{< i}}{\partial x_{j >}} - \frac{4}{5} \frac{\partial q_{< i}}{\partial x_{j >}} - 2 σ_{k < i} \frac{\partial u_{j >}}{\partial x_{k}} \end{matrix}

(14)

\begin{matrix} \frac{\partial q_{i}}{\partial t} + \frac{\partial u_{j} q_{i}}{\partial x_{j}} + \frac{1}{2} \frac{\partial R_{i j}}{\partial x_{j}} \\ = - \frac{2}{3} \frac{p}{μ} q_{i} - \frac{5}{2} p R \frac{\partial T}{\partial x_{i}} - \frac{7 σ_{i k} R}{2} \frac{\partial T}{\partial x_{k}} - R T \frac{\partial σ_{i k}}{\partial x_{k}} + \\ \frac{σ_{i j}}{ρ} (\frac{\partial p}{\partial x_{j}} + \frac{\partial σ_{j k}}{\partial x_{k}}) - \frac{2}{5} \times (\frac{7}{2} q_{k} \frac{\partial u_{i}}{\partial x_{k}} + q_{k} \frac{\partial u_{k}}{\partial x_{i}} + \\ q_{i} \frac{\partial u_{k}}{\partial x_{k}}) - \frac{1}{6} \frac{\partial Δ}{\partial x_{i}} - m_{i j k} \frac{\partial u_{j}}{\partial x_{k}} \end{matrix}

(15)

结合式（13）~（15）得到R13矩控制方程，由于高阶矩m_ijk、R_ij、Δ未知，为了封闭R13矩方程组，Struchtrup等基于Boltzmann方程推导其截断形式^[2]，如式（16）~（18）所示：

m_{i j k} = - 2 μ \frac{\partial (σ_{< i j} / ρ)}{\partial x_{k >}} - \frac{8 μ}{5 p} q_{< i} \frac{\partial u_{j}}{\partial x_{k >}}

(16)

\begin{matrix} R_{i j} = - \frac{24}{5} μ \frac{\partial (q_{< i j} / ρ)}{\partial x_{j >}} - \frac{24 μ}{5 p} q_{< i} \frac{\partial R T}{\partial x_{j >}} - \frac{24 μ}{7 ρ} \times \\ (σ_{k < i} \frac{\partial u_{j >}}{\partial x_{k}} + σ_{k < i} \frac{\partial u_{k}}{\partial x_{j >}} - \frac{2}{3} σ_{i j} \frac{\partial u_{k}}{\partial x_{k}}) - \frac{4}{7} \frac{σ_{k < i} σ_{j > k}}{ρ} \end{matrix}

(17)

\begin{matrix} Δ = - 12 μ \frac{\partial (q_{k} / ρ)}{\partial x_{k}} - 20 \frac{μ}{p} q_{k} \frac{\partial R T}{\partial x_{k}} - \\ 12 \frac{μ}{ρ} σ_{i j} \frac{\partial u_{i}}{\partial x_{j}} - \frac{σ_{i j} σ_{i j}}{ρ} \end{matrix}

(18)

对比经典NSF方程，Struchtrup等将经典流体力学控制方程中的宏观参数由ρ、u_i、T拓展至ρ、u_i、T、σ_ij、q_i。进一步，Gu等将R13矩方程的控制方程拓展到R26矩方程，其高阶矩的控制方程如式（19）所示：

\{\begin{matrix} \frac{\partial m_{i j k}}{\partial t} + \frac{\partial u_{l} m_{i j k}}{\partial x_{l}} + \frac{\partial ϕ_{i j k l}}{\partial x_{l}} \\ = - \frac{3}{2} \frac{p}{μ} m_{i j k} - 3 p \frac{\partial (σ_{< i j} / ρ)}{\partial x_{k >}} + M_{i j k} \\ \frac{\partial R_{i j}}{\partial t} + \frac{\partial u_{k} R_{i j}}{\partial x_{k}} + \frac{\partial ψ_{i j k}}{\partial x_{k}} \\ = - \frac{7}{6} \frac{p}{μ} R_{i j} - \frac{28}{5} p \frac{\partial (q_{< i} / ρ)}{\partial x_{j >}} + R_{i j} \\ \frac{\partial Δ}{\partial t} + \frac{\partial Δ u_{i}}{\partial x_{i}} + \frac{\partial Ω_{i}}{\partial x_{i}} = - \frac{3}{2} \frac{p}{μ} Δ - 8 p \frac{\partial (q_{i} / ρ)}{\partial x_{i}} + χ \end{matrix}

(19)

R26控制方程中源项

M_{i j k}

、

R_{i j}

、

χ

表达式如式（20）~（22）所示：

M_{i j k} = - 3 σ_{< i j} \frac{\partial R T}{\partial x_{k >}} + \frac{3 σ_{< i j}}{ρ} (\frac{\partial p}{\partial x_{k >}} + \frac{\partial σ_{k > m}}{\partial x_{m}}) - \frac{12}{5} q_{< i} \frac{\partial u_{j}}{\partial x_{k >}} - 3 m_{m < i j} \frac{\partial u_{k >}}{\partial x_{m}}

(20)

\begin{matrix} R_{i j} = - \frac{2}{3} \frac{p}{μ} \frac{σ_{k < i} σ_{j > k}}{ρ} + (\frac{8}{3} R T σ_{i j} - \frac{2}{7} R_{i j} \frac{\partial u_{k}}{\partial x_{k}} - \\ \frac{4}{7} (7 R T σ_{k < i} + R_{k < i}) (\frac{\partial u_{j >}}{\partial x_{k}} + \frac{\partial u_{k}}{\partial x_{j >}}) - \frac{14}{15} Δ \frac{\partial u_{< i}}{\partial x_{j >}} - \\ 2 R_{k < i} \frac{\partial u_{j >}}{\partial x_{k}} + \frac{28}{5} \frac{q_{< i}}{ρ} \frac{\partial u_{j > k}}{\partial x_{k}} + 2 \frac{m_{i j k}}{ρ} \frac{\partial σ_{k l}}{\partial x_{l}} + \\ \frac{28}{5} R T q_{< i} (\frac{\partial p}{p \partial x_{j >}} - 2 \frac{\partial T}{T \partial x_{j >}}) - 2 ϕ_{i j k l} \frac{\partial u_{k}}{\partial x_{l}} + \\ m_{i j k} (\frac{2}{ρ} \frac{\partial p}{\partial x_{k}} - 9 \frac{\partial R T}{\partial x_{k}}) + \frac{14}{3} \frac{σ_{i j}}{ρ} (\frac{\partial q_{m}}{\partial x_{m}} + σ_{m l} \frac{\partial u_{m}}{\partial x_{l}}) \end{matrix}

(21)

\begin{matrix} χ = - \frac{2}{3} \frac{p}{μ} \frac{σ_{k j} σ_{j k}}{ρ} - \frac{4}{3} Δ \frac{\partial u_{k}}{\partial x_{k}} - 4 (2 R T σ_{k l} + R_{k l}) \frac{\partial u_{k}}{\partial x_{l}} + \\ 8 \frac{q_{k}}{ρ} (\frac{\partial σ_{k l}}{\partial x_{l}}) + R T q_{k} (8 \frac{\partial p}{p \partial x_{k}} - 28 \frac{\partial T}{T \partial x_{k}}) \end{matrix}

(22)

Gu等给出了R26矩方程以及封闭方程中高阶矩φ_ijkl、ψ_ijk、Ω_i的表达式^[6]。至此，式（13）~（15）、式（19），结合封闭方程式（20）~（22），构成R26矩方程组。

3 全流场耦合DVM/R26加速方法

基于Boltzmann方程的离散速度法在滑移流与早期过渡流区域计算精度高但收敛速度极慢；相反宏观矩方程组在小Kn区域计算精度低，但收敛速度快，占用内存少。因此，采用全流场耦合介观/宏观方程的形式，可以有效地结合宏观R26矩方程组计算效率高和介观Boltzmann方程计算精度高的优点。

作者前期首次提出了壁面处耦合DVM/R26矩方程算法，最大降低了99%的计算内存^[8]。在此基础上，本文提出全流场耦合DVM/R26加速方法。在加速方法中，宏观R26矩方程组通过Hermite多项式重构分布函数f，并传输给介观Boltzmann方程，分布函数的重构方程为：

\begin{matrix} f = f^{e q} \sum_{n = 0}^{\infty} \frac{1}{n!} a_{A}^{(n)} H_{A}^{(n)} \\ = f^{e q} (a^{(0)} H^{(0)} + a_{i}^{(1)} H_{i}^{(1)} + \frac{1}{2!} a_{i j}^{(2)} H_{i j}^{(2)} + \\ \frac{1}{3!} a_{i j k}^{(3)} H_{i j k}^{(3)} + \dots) \end{matrix}

(23)

其中，

H_{A}^{（ n ）}

与

a_{A}^{（ n ）}

分别为Hermite函数与对应的系数。基于R13矩方程组得到的宏观参数来重构分布函数f^（3），可以写为：

\begin{matrix} f^{(3)} = f_{M} \{1 + \frac{c_{i} q_{i}}{p R T} (\frac{c^{2}}{5 R T} - 1) + \frac{m_{i j k} c_{i} c_{j} c_{k}}{6 p (R T)^{2}} + \\ \frac{σ_{i j} c_{i} c_{j}}{2 p R T} + \frac{R_{i j} c_{i} c_{j}}{4 p (R T)^{2}} (\frac{c^{2}}{7 R T} - 1) + \\ \frac{Δ}{8 p R T} [\frac{c^{4}}{15 (R T)^{2}} - \frac{2 c^{2}}{3 R T} + 1]\} \end{matrix}

(24)

基于R26矩方程组得到的宏观参数来重构分布函数f^（5），可以写为：

\begin{matrix} f^{(5)} = f_{M} \{1 + \frac{σ_{i j} c_{i} c_{j}}{2 p R T} + \frac{c_{i} q_{i}}{p R T} (\frac{c^{2}}{5 R T} - 1) + \\ \frac{ψ_{i j k} c_{i} c_{j} c_{k}}{12 p (R T)^{3}} (\frac{c^{2}}{9 R T} - 1) + \frac{m_{i j k} c_{i} c_{j} c_{k}}{6 p (R T)^{2}} + \\ \frac{ϕ_{i j k} c_{i} c_{j} c_{k} c_{l}}{24 p (R T)^{3}} + \frac{R_{i j} c_{i} c_{j}}{4 p (R T)^{2}} (\frac{c^{2}}{7 R T} - 1) + \\ \frac{c_{i} Ω_{i}}{40 p (R T)^{2}} [\frac{c^{4}}{7 (R T)^{2}} - \frac{2 c^{2}}{R T} + 5] + \\ \frac{Δ}{8 p R T} [\frac{c^{4}}{15 (R T)^{2}} - \frac{2 c^{2}}{3 R T} + 1]\} \end{matrix}

(25)

至此，Hermite多项式成为连接介观Boltzmann模型方程与宏观R26矩方程数据传输的桥梁，全流场耦合迭代介观/宏观加速方法的计算流程图如图1所示。

图1全流场耦合介观/宏观方程加速方法的计算流程图

Fig.1Calculation flow chart of the acceleration scheme for coupled mesoscopic/macroscopic in the full flow area

具体可以写为：

1）首先基于DVM与Maxwell镜面反射/漫反射壁面边界条件，在全流场求解Boltzmann模型方程，并将壁面处的宏观量作为矩方程的壁面边界条件；

2）评估宏观量（ρ，u_i，T，σ_ij，q_i，m_ijk，R_ij，Δ）是否满足收敛条件；

3）在宏观方程虚拟边界处，将DVM计算出的全场宏观量（ρ，u_i，T，σ_ij，q_i，m_ijk，R_ij，Δ）输入矩方程中，作为矩方程的边界条件；

4）求解宏观方程，例如NSF方程、R13以及R26矩方程，设置矩方程收敛条件，直至宏观方程得到稳态收敛解；

5）通过宏观方程的解计算Hermite多项式

H_{A}^{（ n ）}

以及相应的系数

a_{A}^{（ n ）}

；

6）利用Hermite多项式展开，如式（25），重构分布函数；

7）给介观方程提供边界条件，回到步骤（Ⅰ）然后迭代求解介观Boltzmann方程，直到宏观参数满足收敛条件。

经典DVM从任意流场初值进行数值迭代，本文提出的全流场耦合DVM/矩方程的加速思路在于：充分利用宏观方程在低Kn下的加速作用，由宏观方程计算流场稳态解，在此基础上，由介观Boltzmann方程对稳态解进行优化，宏观方程为介观方程的求解提供了较为准确的初始值。

4 算例与讨论

针对本文提出的全流场耦合DVM/R26加速迭代方法，选取加速度驱动Poiseuille流、顶盖驱动方腔流典型算例进行验证。对于本文中的算例，NSF、R26矩方法、DVM与耦合DVM/R26加速算法共用同一套计算网格，其中宏观方程的求解基于有限体积法，DVM基于有限差分法，宏观方法与介观方法在整个计算域耦合。DVM计算结果通过了空间网格和速度网格无关性测试，可认为是标准解^[1]。模型中气体分子介质为氩气，黏性服从Sutherland法则，如式（26）所示，进而对比加速方法与经典DVM计算时长与计算精度。

μ = μ_{0} {(\frac{T}{T_{0}})}^{1.5} \frac{T_{0} + S}{T + S}

(26)

式中，参考温度与参考黏性分别为T₀=273 K、μ₀=21.25×10^-6 Pa·s，气体常数R与Sutherland常数为R=208 J/（kgK）、S=144 K。Maxwell壁面边界条件采用漫反射边界条件α=1。

4.1 加速度驱动Poiseuille流动

全流场耦合DVM/R26加速方法首先在加速度驱动Poiseuille流动中进行验证，上、下壁面分别位于y=-L₀/2与y=L₀/2，参考长度与参考温度分别为L₀=10^-5 m、T₀=273 K，外界驱动力作用于x正向，无量纲大小为

{\bar{a}}_{i}

=L₀a_i/（2RT₀）=0.01，在空间上沿y方向分布101个非均匀网格点，无限长速度空间被截断为

[- 6 \sqrt{2 R T_{0}} ， {6 \sqrt{2 R T_{0}}]}^{3}

，并基于Newton-Cotes积分方法将其离散为64×64×24个非均匀网格点。

当来流Kn分别为0.01，0.1，0.5，1.0时，分别采用宏观NSF、R26矩方法、介观DVM与全流场耦合DVM/R26加速方法进行迭代求解，对比计算结果如图2所示，计算时长如表1所示。

图2加速度驱动Poiseuille流动中速度

{\bar{u}}_{x}

沿垂直中心线的分布

Fig.2

{\bar{u}}_{x}

-velocity along the vertical centerline in the acceleration-driven Poiseuille flow

表1加速度驱动Poiseuille流动计算时间消耗对比

Tab.1 Comparison of calculation time consumption of the acceleration-driven Poiseuille flow

从图2（a）可以看出：当Kn=0.01时，宏观NSF、R26矩方法、介观DVM与耦合加速方法均能够准确地捕获流体的非平衡效应，最大无量纲速度

{\bar{u}}_{x}

=0.229位于管道中心线处。随着Kn升高到0.1、0.5，如图2（b）~（c）所示，宏观NSF方程组对管道中心处最大速度的计算误差分别为8.20%和27.53%，R26矩方法对最大速度的计算误差分别为1.10%和6.43%。当Kn=1.0时，如图2（d）所示，宏观NSF方程组失效，R26矩方法计算精度大幅度降低。而全流场耦合DVM/R26加速方法由于最终切换为介观DVM，因此最终计算结果与DVM计算结果一致。

从表1可以看出，DVM在Kn=0.01、Kn=0.1时分别需要142.20 min、19.18 min才能收敛，而宏观方法大大缩短了稳态求解的计算时长；全流场耦合DVM/R26加速方法收敛需要5.4 min、4.72 min。这是由于全流场耦合DVM/R26加速方法将宏观方程较为准确的稳态解作为DVM迭代初始解，因而大大降低了计算时间。反观大Kn情况，如Kn=0.5、Kn=1.0时，由于R26矩方法需要更小的松弛因子才能收敛，故迭代时长会大幅度增加，全流场耦合DVM/R26加速方法的加速效果随之下降。

4.2 顶盖驱动方腔流动

其次，全流场耦合DVM/R26加速方法的性能在2D顶盖驱动方腔流动中进行进一步验证，方腔的边长L=10^-5 m，顶盖以速度u₀=10 m/s进行移动，方腔内部参考温度与参考黏性分别设置为T₀=273 K，μ₀=21.25×10^-6 Pa·s，初始压力P₀通过Kn反算得到，壁面温度设置为T^w=T₀，流动示意图如图3所示。在物理空间内，方腔被划分为101×101非均匀网格，在靠近壁面处进行加密。在速度空间内，将连续速度空间截断为

[- 6 \sqrt{2 R T_{0}} ， {6 \sqrt{2 R T_{0}}]}^{3}

，并基于Newton-Cotes积分方法将其离散为64×64×24个非均匀网格点。

图3顶盖驱动方腔流动示意图

Fig.3Schematic diagram of the lid-driven cavity flow

当来流Kn分别为0.01，0.1，0.5，1.0时，分别采用宏观NSF、R26矩方法、介观DVM与全流场耦合DVM/R26加速方法进行迭代求解，对比计算结果如图4所示，计算时长如表2所示。

图4顶盖驱动方腔流中速度u_x沿垂直中心线的分布

Fig.4u_x-velocity along the vertical centerline in the lid-driven cavity flow

从图4（a）可以看出：在滑移流Kn=0.01时，所有算法对速度分布的计算结果一致，均能够较为准确地模拟气体的非平衡流动状态。随着Kn升高到0.1、0.5，如图4（b）~（c）所示，宏观NSF方程组对方腔中心线处速度的计算误差分别为12.59%和47.56%，R26矩方法对最大速度的计算误差分别为6.74%和8.02%。当Kn=1.0时，如图2（d）所示，NSF方程组失效，R26矩方法无法获得稳态解。在所有算例中，全流场耦合DVM/R26加速方法均能够恢复到DVM计算结果，加速方法相较宏观方法大大提高了计算精度。

表2顶盖驱动方腔流计算时间消耗对比

Tab.2 Comparison of calculation time consumption of the lid-driven cavity flow

从表2可以看出：DVM在低Kn下计算时长较长，这是由于此时分子平均自由程较小时，分子之间的碰撞与迁移传播速度较慢，因此需要更多的迭代步数才能获得全流场稳态解。宏观方程虽然计算精度较低，但低Kn下计算时长远小于介观DVM。全流场耦合DVM/R26加速方法在滑移区与早期过渡区，如Kn=0.01、Kn=0.1时，具有显著加速作用，最大降低95.28%的计算时长，同时保持了DVM的计算精度。但对于中、大Kn情况，由于宏观方程的松弛因子需要很小才能满足数值迭代的稳定性，因此全流场耦合DVM/R26加速方法的加速能力大幅度下降，耦合DVM/R26方法应尽早切换至DVM求解器以快速获得稳态解。

5 结论

本文针对稀薄气体在滑移流、早期过渡流的数值模拟，提出了全流场耦合介观/宏观方程加速迭代方法。仿真结果表明：一方面，随着Kn的增加，宏观方程的计算精度逐渐降低，同时在大Kn下，R26矩方法易发散，而耦合加速方法最终切换为介观DVM，因此对任意Kn都能维持DVM计算精度。另一方面，求解Boltzmann方程的离散速度法在近连续流、滑移流、早期过渡流收敛速度极慢，耦合加速方法充分利用了宏观方程在低Kn下快速收敛的特点，最大降低了95.28%的计算时长，但耦合加速方法对于大Kn情况加速作用不明显。本文所提方法对Kn适用范围的下限受二阶迎风格式DVM的限制，无法模拟连续流、近连续流问题。对Kn适用范围的上限受R26矩方法稳定性限制，Kn≤1.0。因此，本文所提出的加速方法仅适用于滑移流、早期过渡流的数值模拟问题研究。

致谢

英国达斯勃利实验室（STFC Daresbury Lab）顾晓军教授在矩方法、耦合方法上提供了帮助和指导，谨致谢意！

图1全流场耦合介观/宏观方程加速方法的计算流程图

Fig.1Calculation flow chart of the acceleration scheme for coupled mesoscopic/macroscopic in the full flow area

下载: 全尺寸图片

图2加速度驱动Poiseuille流动中速度

{\bar{u}}_{x}

沿垂直中心线的分布

Fig.2

{\bar{u}}_{x}

-velocity along the vertical centerline in the acceleration-driven Poiseuille flow

下载: 全尺寸图片

图3顶盖驱动方腔流动示意图

Fig.3Schematic diagram of the lid-driven cavity flow

下载: 全尺寸图片

图4顶盖驱动方腔流中速度u_x沿垂直中心线的分布

Fig.4u_x-velocity along the vertical centerline in the lid-driven cavity flow

下载: 全尺寸图片

表1加速度驱动Poiseuille流动计算时间消耗对比

下载: 全尺寸图片

表2顶盖驱动方腔流计算时间消耗对比

下载: 全尺寸图片

图1全流场耦合介观/宏观方程加速方法的计算流程图

Fig.1Calculation flow chart of the acceleration scheme for coupled mesoscopic/macroscopic in the full flow area

图2加速度驱动Poiseuille流动中速度

{\bar{u}}_{x}

沿垂直中心线的分布

Fig.2

{\bar{u}}_{x}

-velocity along the vertical centerline in the acceleration-driven Poiseuille flow

图3顶盖驱动方腔流动示意图

Fig.3Schematic diagram of the lid-driven cavity flow

图4顶盖驱动方腔流中速度u_x沿垂直中心线的分布

Fig.4u_x-velocity along the vertical centerline in the lid-driven cavity flow

表1加速度驱动Poiseuille流动计算时间消耗对比

表2顶盖驱动方腔流计算时间消耗对比

图1全流场耦合介观/宏观方程加速方法的计算流程图

Fig.1Calculation flow chart of the acceleration scheme for coupled mesoscopic/macroscopic in the full flow area

图2加速度驱动Poiseuille流动中速度

{\bar{u}}_{x}

沿垂直中心线的分布

Fig.2

{\bar{u}}_{x}

-velocity along the vertical centerline in the acceleration-driven Poiseuille flow

图3顶盖驱动方腔流动示意图

Fig.3Schematic diagram of the lid-driven cavity flow

图4顶盖驱动方腔流中速度u_x沿垂直中心线的分布

Fig.4u_x-velocity along the vertical centerline in the lid-driven cavity flow

表1加速度驱动Poiseuille流动计算时间消耗对比

表2顶盖驱动方腔流计算时间消耗对比

SU W, WANG P, ZHANG Y H,et al. A high-order hybridizable discontinuous Galerkin method with fast convergence to steady-state solutions of the gas kinetic equation[J]. Journal of Computational Physics,2019,376:973-991.

STRUCHTRUP H. Macroscopic transport equations for rarefied gas flows[M]. Berlin, Heidelberg: Springer-Verlag Berlin Heidelberg,2005:145-160.

SHAVALIYEV M S. Super-Burnett corrections to the stress tensor and the heat flux in a gas of Maxwellian molecules[J]. Journal of Applied Mathematics and Mechanics,1993,57(3):573-576.

GU X J, BARBER R W, EMERSON D R. How far can 13 moments go in modeling microscale gas phenomena?[J]. Nanoscale and Microscale Thermophysical Engineering,2007,11(1/2):85-97.

STRUCHTRUP H, TORRILHON M. Regularization of Grad′s 13 moment equations:derivation and linear analysis[J]. Physics of Fluids,2003,15(9):2668-2680.

GU X J, EMERSON D R. A high-order moment approach for capturing non-equilibrium phenomena in the transition regime[J]. Journal of Fluid Mechanics,2009,636:177-216.

BIRD G A. Molecular gas dynamics and the direct simulation of gas flows[M]. Oxford: Clarendon Press,1994.

YANG W Q, GU X J, WU L,et al. A hybrid approach to couple the discrete velocity method and method of moments for rarefied gas flows[J]. Journal of Computational Physics,2020,410:109397.

YANG W Q, GU X J, EMERSON D R,et al. Modelling thermally induced non-equilibrium gas flows by coupling kinetic and extended thermodynamic methods[J]. Entropy,2019,21(8):816.

YANG W Q, TANG S, YANG H. Analysis of the moment method and the discrete velocity method in modeling non-equilibrium rarefied gas flows.a comparative study[J]. Applied Sciences,2019,9(13):2733.

YANG W Q, GU X J, EMERSON D,et al. Comparative study of the discrete velocity and the moment method for rarefied gas flows[C]//Proceedings of 31st International Symposium on Rarefied GAS Dynamics: RGD31,2019.

LI Z H, ZHANG H X. Gas-kinetic numerical studies of three-dimensional complex flows on spacecraft re-entry[J]. Journal of Computational Physics,2009,228(4):1116-1138.

XU K, HUANG J C. A unified gas-kinetic scheme for continuum and rarefied flows[J]. Journal of Computational Physics,2010,229(20):7747-7764.

GUO Z L, WANG R J, XU K. Discrete unified gas kinetic scheme for all Knudsen number flows:Ⅱ.thermal compressible case[J]. Physical Review E,2015,91(3):033313.

WANG P, HO M T, WU L,et al. A comparative study of discrete velocity methods for low-speed rarefied gas flows[J]. Computers & Fluids,2018,161:33-46.

BHATNAGAR P L, GROSS E P, KROOK M. A model for collision processes in gases.Ⅰ.small amplitude processes in charged and neutral one-component systems[J]. Physical Review,1954,94(3):511-525.

JR HOLWAY L H. New statistical models for kinetic theory:methods of construction[J]. Physics of Fluids,1966,9(9):1658-1673.

SHAKHOV E M. Generalization of the Krook kinetic relaxation equation[J]. Fluid Dynamics,1968,3(5):95-96.

ZHU Y J, ZHONG C W, XU K. Implicit unified gas-kinetic scheme for steady state solutions in all flow regimes[J]. Journal of Computational Physics,2016,315:16-38.

ZHU Y J, ZHONG C W, XU K. An implicit unified gas-kinetic scheme for unsteady flow in all Knudsen regimes[J]. Journal of Computational Physics,2019,386:190-217.

PAN D X, ZHONG C W, ZHUO C S. An implicit discrete unified gas-kinetic scheme for simulations of steady flow in all flow regimes[J]. Communications in Computational Physics,2019,25(5):1469-1495.

SU W, ZHU L H, WANG P,et al. Can we find steady-state solutions to multiscale rarefied gas flows within dozens of iterations?[J]. Journal of Computational Physics,2020,407:109245.

YUAN R F, LIU S, ZHONG C W. A novel multiscale discrete velocity method for model kinetic equations[J]. Communications in Nonlinear Science and Numerical Simulation,2021,92:105473.

YUAN R F, LIU S, ZHONG C W. A multi-prediction implicit scheme for steady state solutions of gas flow in all flow regimes[J]. Communications in Nonlinear Science and Numerical Simulation,2021,92:105470.

YUAN R F, ZHONG C W. A conservative implicit scheme for steady state solutions of diatomic gas flow in all flow regimes[J]. Computer Physics Communications,2020,247:106972.

YANG J Y, HUANG J C. Rarefied flow computations using nonlinear model Boltzmann equations[J]. Journal of Computational Physics,1995,120(2):323-339.

TITAREV V A. Conservative numerical methods for model kinetic equations[J]. Computers & Fluids,2007,36(9):1446-1459.

Home

About Journal

Guide for Authors

Editorial Board

Publication Statement

Open Access Statement

Contact

Journal Subscription

Rss

AI assistant

Chinese

1 Boltzmann方程和速度分布函数的矩

2 拓展流体力学方程

3 全流场耦合DVM/R26加速方法

4 算例与讨论

5 结论