耦合宏观/介观方法的振荡热驱稀薄流动

杨伟奇; 钮耀斌; YANG Weiqi; NIU Yaobin

耦合宏观/介观方法的振荡热驱稀薄流动

doi: 10.11887/j.issn.1001-2486.23090015

杨伟奇，钮耀斌

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

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

详细信息

作者简介

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

通讯作者

钮耀斌（1983—），男，浙江嘉兴人，副研究员，博士，硕士生导师，E-mail: nyb2007@nudt.edu.cn

中图分类号: O354

文献标识码: A

文章编号: 1001-2486(2025)06-168-10

Hybrid macroscopic/mesoscopic scheme for oscillatory thermally driven rarefied gas flows

YANG Weiqi ， NIU Yaobin

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

摘要

对二维方腔中顶盖温度振荡驱动稀薄流动进行研究，分析克努森数Kn与顶盖温度振荡频率对流场参数的影响规律。在靠近壁面处基于介观方法求解Shakhov模型方程，中心流域基于宏观方法减少计算资源消耗。宏观/介观虚拟边界利用高阶Hermite多项式重构分布函数，以封闭数值迭代格式。仿真结果表明：耦合方法相较全流场介观尺度方法，对方腔垂直中线处温度的预测结果相符，最大计算误差为0.23%，计算内存消耗降低约69.91%。耦合方法能够捕获振荡热驱稀薄流动在大Kn流域的非线性现象，此时水平方向速度分布不再服从周期振荡余弦函数规律，上升时间远大于下降时间。黏性穿透层厚度与受扰区域随着Kn的增大而增大，随着St的增大而减小。

关键词

稀薄气流 / 非平衡热驱流动 / 离散速度法 / 矩方法 / 耦合方法

Abstract

Thermally-induced oscillatory rarefied gas flow inside a two-dimensional rectangular cavity was investigated. The effects of the Knudsen numbers and the oscillation frequency of lid temperature on the flow parameters were analyzed. The Shakhov model equation was solved numerically based on the mesoscopic approach in the near-wall region, and the macroscopic approach was adopted in the bulk flow region to reduce the computational cost. To close the numerical iteration procedure, the velocity distribution functions, served as the pseudo boundary between macroscopic and mesoscopic methods, were reconstructed using the high-order Hermite polynomials. Numerical simulations demonstrate that the temperature profile at the central vertical of the cavity predicted by the hybrid method is in good agreement with results from the mesoscopic method, with maximum error 0.23%. Besides, the computational memory cost can be saved up to about 69.91%. The hybrid approach is able to capture the nonlinear phenomenon in the thermally-induced oscillatory rarefied gas flow under high Kn numbers, where the horizontal velocity no longer obeys the law of periodic oscillating cosine function, and the rise time of the horizontal velocity is much longer than the fall time. The thickness of the viscous penetration layer and the disturbed region increases as the Kn number increases, and decreases as the St number increases.

Keywords

rarefied gas flows / thermally-induced non-equilibrium flows / discrete velocity method / moment method / hybrid method

1 问题描述 2 数值方法 2.1 Boltzmann-Shakhov模型方程与DVM 2.2 拓展流体力学方程组与R26矩方法 2.3 耦合宏观/介观方法 3 算例与讨论 3.1 耦合宏观/介观方法计算精度对比 3.2 非平衡效应对流场宏观参数的影响 3.3 顶盖温度振荡频率对流场宏观参数的影响 4 结论

随着精密制造领域微机电系统（micro electro mechanical systems，MEMS）在近十年的高速发展，微型器件广泛应用在国防科技与航空航天领域^[1-2]。为了保证元器件的高分辨率与高敏感度，其内部运动部件通常工作在非平衡稀薄流状态，气体对运动部件表面产生显著的黏滞阻尼效应^[2-3]。在微尺度或稀薄流动下，MEMS的壁面温度梯度与非均匀温度场驱动气体分子流动，如热蠕流动^[4]、Knudsen泵^[5]等热驱流动，温度的变化将大幅降低测量精度^[6]，导致测量误差较大，甚至出现错误结果^[7]。因此，精确高效地描述其流动机理，对预测设计MEMS器件热力学性能具有重要价值。

在MEMS中，描述气体稀薄程度的克努森数Kn一般定义为：分子平均自由程λ与特征尺度L₀的比值^[8]。根据Kn的不同，流体可以划分为连续流域（Kn≤0.001）、滑移流域（0.001＜Kn≤0.1）、过渡流域（10^-1＜Kn≤10）与自由分子流域（Kn＞10）。在连续流域，基于非定常纳维-斯托克斯-傅里叶（Navier-Stokes-Fourier，NSF）方程的振荡流动受到了广泛的研究^[9-12]，如平板间振荡气流流动传热研究^[9-10]；顶盖振荡驱动方腔流动中，不同雷诺数与振动频率对于涡系位置与流线的影响研究等^[11-12]。微纳尺度MEMS器件一般处于滑移流域与早期过渡流域（10^-1＜Kn≤1），气体的非平衡特性可以通过高阶宏观方法或介观方法进行描述。

在宏观尺度，纳维-斯托克斯（Navier-Stokes，NS）方程组与速度滑移、温度跳跃边界条件可以描述滑移流域的非平衡特性^[13-15]。Emerson等^[9]基于NS方程组，得到了圆柱间振荡Couette流动速度剖面与切应力的解析解。Guo等^[16]与Tang等^[17]结合修正的壁面边界条件或修正的分子平均自由程拓展了NS方程组对Kn的适用范围，但这种方法仅对简单壁面几何较为高效。为了进一步拓展宏观方程的适用范围，基于气体动理学方程CE（Chapman-Enskog）展开的Burnett^[18]（二阶展开）与超Burnett方程组^[19]（三阶展开），与基于Hermite多项式展开^[20]的Grad矩方法^[21]得到了大量的研究。Struchtrup等^[22-23]与Gu等^[24]分别对Grad矩方法进行规整化处理，得到R13/R26矩方法。Tatsios等^[25]从半幅Hermite多项式出发，改进了矩方法的壁面边界条件。规整R13/R26矩方法与半幅矩方法成功应用到低频振荡稀薄气体流动^[23]、Stokes第二问题^[24]、平板振荡Couette流动中^[24-25]，仿真结果表明线性R26矩方法在早期过渡流域与气体动理论结果符合较好。矩方法在滑移流域与早期过渡流域计算效率高，但无法捕获大Kn流域非平衡效应。

在介观尺度，振荡稀薄流动中的气体黏滞阻尼效应一般通过Boltzmann方程来进行描述^[8]。其中，Park等^[26]、Bahukudumbi^[27]与Sharipov等^[10]研究了全Kn流域下任意振荡频率对一维Couette流动的影响，并给出了近连续流域（Kn≤0.1）与自由分子流域（Kn＞10）的解析解。Wu等基于快速谱方法求解线性Boltzmann方程，分别研究了顶盖振荡驱动方腔稀薄流^[13]与二维振荡方腔流中声传播问题^[28]，给出了高频振荡下的解析解。Wu等指出：二维方腔中振荡顶盖处的黏滞阻力小于一维振荡Couette流动中的黏滞阻力，同时黏滞阻力与靠近顶盖处的声速受方腔长宽比影响较大。在过渡流域，振荡稀薄流动很难得到解析解，一般采用直接模拟蒙特卡罗^[8]（direct simulation Monte Carlo，DSMC）或离散速度法^[29]（discrete velocity method，DVM）等数值方法求解。Frangi等基于BGK（Bhatnagar-Gross-Krook）模型在宽Kn范围内预测了MEMS的黏滞阻尼，仿真结果与实验结果符合较好^[30]。Emerson等基于NS方程和DSMC，在全流域内预测了圆柱间振荡Couette流动，分析了Kn对速度与切应力剖面的影响^[9]。Tsimpoukis等针对平板间非线性振荡稀薄流，对比了基于DSMC方法与DVM-BGK方法的轴向速度与流量计算结果，最大误差在10%以内^[31]。Wang等使用Shakhov模型方程，基于离散统一气体动理学格式（discrete unified gas kinetic scheme，DUGKS）分别研究了二维^[32]与三维^[33]顶盖振荡驱动方腔流中Kn、长宽比、振动频率对黏滞阻尼的影响，对比了线性与非线性振动条件下流体的传热特性。

Yang等提出了耦合宏观/介观数值迭代方法，并成功求解管道流、方腔流等定常来流的稳态问题^[34-35]，而非定常条件下周期性振荡热驱稀薄流动尚无相关研究。耦合宏观/介观方法可以看为宏观尺度拓展流体力学与介观尺度分子动理论之间的桥梁，能够平衡稀薄流动模拟中的计算资源消耗与计算精度要求。

1 问题描述

考虑2D方腔中温度振荡驱动稀薄流动，其中坐标原点O位于方腔左下角，方腔的长和高分别为L=10^-5 m和H=10^-5 m，如图1所示。左/右壁面与底边壁面温度T₀设为参考温度，顶盖温度T_w为：

T_{w} = T_{0} + d T c o s (ω t)

(1)

式中，t为时间，dT=0.3T₀与ω分别为温度振荡的振幅与频率。论文重点开展温度周期性振荡驱动稀薄流动，其中特征尺度的无量纲参数为：

\{\begin{matrix} A = \frac{L}{H} \\ S t = \frac{ω H}{v_{0}} \\ K n = \frac{λ}{H} \end{matrix}

(2)

其中，A为方腔的长宽比，St为Strouhal数，v₀=

\sqrt{2 R T_{0}}

为参考速度，R=208 J/（kg·K）为理想气体常数，λ为分子平均自由程。

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

(3)

式中，压力p₀服从理想状态方程p₀=ρRT₀，ρ为密度，黏性μ服从Sutherland黏性定律：

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

(4)

式中，参考黏性μ₀与参考温度T₀分别为μ₀=21.25×10^-6 Pa·s、T₀=273 K，氩气分子对应的Sutherland常量为S=144 K。

图1方腔中温度振荡驱动稀薄流动示意图

Fig.1Schematic diagram of the thermally-induced oscillatory rarefied gas flows in a rectangular cavity

2 数值方法

2.1 Boltzmann-Shakhov模型方程与DVM

介观尺度基于Shakhov模型方程^[36]描述速度分布函数演化规律，如式（5）所示。

\frac{\partial f}{\partial t} + C_{i} \frac{\partial f}{\partial x_{i}} = - \frac{1}{τ} (f - f^{e q})

(5)

式中，f=f（t，x_i，C_i）为气体分子速度分布函数（velocity distribution function，VDF），x_i与C_i分别为物理空间坐标与分子离散速度，τ=μ/p为松弛时间，p为压力，f^eq为Shakhov模型中平衡态分布函数：

\begin{matrix} f^{e q} = f_{s} = \frac{ρ}{(2 π R T)^{3 / 2}} e x p [- \frac{{(C_{i} - u_{i})}^{2}}{2 R T}] \cdot \\ \{1 + (1 - P r) \frac{(C_{i} - u_{i}) q_{i}}{5 p R T} [\frac{{(C_{j} - u_{j})}^{2}}{R T} - 5]\} \end{matrix}

(6)

其中，u_i、T与q_i分别为流体的速度、温度与热流，Pr为普朗特数。

Shakhov模型方程中，VDF在速度空间内离散为非均匀速度点，基于Newton-Cotes数值积分方法求解VDF。在物理空间内，用基于二阶迎风格式的有限差分方法迭代Shakhov模型方程的迁移项与碰撞项^[34]。对于非稳态问题，Shakhov模型方程中VDF随时间的演化方程可以写为：

(7)

其中：Δt为时间步长，通过CFL（Courant-Friedrichs-Lewy）条件计算。

结合Maxwell镜面反射/漫反射壁面边界条件，得到全流场速度分布函数f=f（t，x，C）。其中气体密度ρ，速度u，温度T，应力张量σ_ij，热流q，高阶矩m_ijk、R_ij、Δ、φ_ijkl、ψ_ijk、Ω_i等流场宏观参数与无量纲格式如参考文献^[34]所示。

2.2 拓展流体力学方程组与R26矩方法

直接给出拓展流体力学中R26矩方程组的控制方程，如式（8）~（15）所示，其推导过程参考文献^[37]，在此不做赘述。

\frac{\partial ρ}{\partial t} + \frac{\partial ρ u_{i}}{\partial x_{i}} = 0

(8)

\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}

(9)

\frac{\partial ρ T}{\partial t} + \frac{\partial ρ u_{i} T}{\partial x_{i}} = - \frac{2}{3 R} \frac{\partial q_{i}}{\partial x_{i}} - \frac{2}{3 R} p \frac{\partial u_{i}}{\partial x_{i}} - \frac{2}{3 R} σ_{i j} \frac{\partial u_{j}}{\partial x_{i}}

(10)

\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} + Σ_{i j} - 2 p \frac{\partial u_{< i}}{\partial x_{j >}}

(11)

\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} + Q_{i} - \frac{5}{2} p R \frac{\partial T}{\partial x_{i}}

(12)

(13)

(14)

(15)

其中，Σ_ij、Q_i、

为源项，见文献^[37]。R26矩方法将传统流体力学中的自由变量{ρ，u_i，T}拓展为{ρ，u_i，T，σ_ij，q_i，m_ijk，R_ij，Δ}，引入的高阶矩信息用以描述流体的非平衡效应。基于Hermite多项式展开，Struchtrup等^[22]与Gu等^[37]分别为R13/R26矩方法推导了高阶矩表达式，用于封闭矩方程组。同时，Gu等^[38]与Torrilhon等^[39]给出了R13/R26矩方法的壁面边界条件。

2.3 耦合宏观/介观方法

基于介观尺度方法模拟温度振荡驱动稀薄流动时，将消耗巨大的计算资源，而基于宏观矩方法模拟稀薄流动时难以捕获大Kn处的非平衡效应。因此，提出耦合宏观/介观快速迭代方法，平衡稀薄流动模拟中的计算精度与计算代价。

方腔中温度振荡驱动稀薄流动，计算域可以分解为两个部分：靠近壁面处基于介观尺度描述气体的非平衡效应，介观尺度计算域厚度为l，远离壁面位置处基于宏观尺度矩方法简化计算量，如图2所示。由于计算域的分解，矩方法与DVM的壁面边界条件均不再封闭，在图2虚线处为通过Hermite多项式重构分布函数来构建的宏观/介观方程的虚拟边界。

图2基于耦合宏观/介观方法的温度振荡驱动稀薄流动示意图

Fig.2Schematic diagram of the thermally-induced oscillatory rarefied gas flows based on the hybrid macroscopic/mesoscopic algorithm

在靠近壁面处，基于DVM与Maxwell漫反射壁面边界条件计算得到VDF，VDF对离散速度求矩可以作为R26矩方程组的虚拟边界条件。反之，当计算得到宏观方程的高阶矩时，可以通过Hermite多项式重构分布函数，Hermite多项式的展开形式如式（16）所示：

(16)

式中，

H_{A}^{（ n ）}

与

a_{A}^{（ n ）}

分别为Hermite函数与对应的系数。

理论上认为，Hermite多项式的无穷维度展开可以精确恢复到非平衡态分布函数。但由于宏观方程维度的限制，只能通过截断的Hermite多项式重构分布函数，对于R26矩方程组，最多能够恢复至5阶精度，f^（5）表达式如式（17）所示。

\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 l} 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{ψ_{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{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}

(17)

其中，f_M为平衡态分布函数。

基于宏观方程重构的分布函数可以为介观尺度Shakhov模型方程提供虚拟边界条件。为保证物理空间离散的二阶精度，虚拟边界层的厚度至少为两层。耦合宏观/介观快速迭代算法数据传输方式与计算域划分方法如图3所示。通过选取合适的壁面介观动理论计算域的厚度，耦合宏观/介观尺度方法可以应用在振荡热驱稀薄流动的数值模拟中。

图3耦合宏观/介观方法计算域与数据传输方法

Fig.3Computational sub-domain and data transfer of the hybrid macroscopic/mesoscopic algorithm

3 算例与讨论

针对方腔中振荡热驱稀薄流动，考虑模型中气体分子介质为氩气，所有壁面条件采用Maxwell漫反射壁面边界条件。方腔长宽比A=1，在物理空间上沿x、y方向各分布201个非均匀网格点。速度空间被截断为

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

，并基于Newton-Cotes积分方法将其离散为64×64×24个非均匀网格点，物理空间与速度空间均经过了网格无关性测试。当来流Kn分别为0.01，0.1，1.0，10时，耦合宏观/介观尺度方法中介观尺度计算域的厚度l与参考长度H的比值l/H分别取为：0.04，0.1，0.2，0.5。前期工作表明：耦合宏观/介观快速迭代方法能够满足Kn≤1时的计算精度与计算效率要求。当Kn=10时，气体的非平衡效应大，矩方程计算精度不足，故选取l/H=0.5，基于全流场DVM捕获精确解。

3.1 耦合宏观/介观方法计算精度对比

当A=1，St=2，Kn=0.1时，对比耦合宏观/介观方法与全流场DVM的计算精度与计算代价。一个温度振荡周期内，方腔沿垂直中线处两种方法计算的无量纲温度

\bar{T}

分布对比如图4所示。

图4一个振荡周期内方腔垂直中线处温度

\bar{T}

（x=L/2）变化

Fig.4Dynamic temperature profiles

\bar{T}

(x=L/2) at the central vertical of the cavity within an oscillation period

图4中实线为耦合宏观/介观方法计算结果，符号为全流场DVM计算结果，将一个完整振荡周期分解为0，0.25π，0.5π，0.75π，分别对比周期内特征时刻计算精度。仿真结果选取开始时刻为全流场达到周期性的“稳态过程”（即下一个振荡周期的流场参数与本周期流场参数完全相同）。仿真表明耦合宏观/介观方法与全流场DVM在全振荡周期内计算结果一致，能够较好地捕获非平衡条件下振荡热驱方腔流动。当靠近壁面边界处的介观尺度计算域为10%时，沿垂直中线温度分布的最大计算误差为0.23%。同时，全流场DVM计算内存消耗约为14.82 GB，耦合宏观/介观方法计算内存消耗约为4.46 GB，计算内存消耗减少约69.91%，进入周期性的“稳态过程”计算时长减少约70.83%。后文仿真均基于耦合宏观/介观方法计算，其中当Kn=10时，选择l/H=0.5，耦合方法自动切换为全流场DVM以保证计算精度。

3.2 非平衡效应对流场宏观参数的影响

由于方腔顶盖温度周期性振荡，温度梯度在非平衡条件下驱动流场运动。考虑稀薄效应对流场参数的影响，选取（x，y）=（0.9L，0.9H）特征点位置，对比当St=2，A=1，Kn取0.01，0.1，1.0，10时，特征点处的流场参数变化，其无量纲温度、密度、x方向与y方向的速度变化如图5所示。

图5方腔特征点处宏观参数随时间变化

Fig.5Time evolution of the macroscopic parameters at the characteristic point of the cavity

选取无量纲时间3π≤

\bar{t}

≤5π分析流场参数，此时流场均进入周期性“稳态过程”。图5（a）表明随着Kn的增大，温度变化随热驱振荡的相位滞后逐渐增大，同时靠近顶盖处的温度跳跃逐渐增大，不同Kn下最大温度响应的幅值偏差为1.53%。因为振荡温度驱动的本质为分子与壁面以及分子之间的碰撞产生的能量转移，因此碰撞的响应存在延迟，且方腔靠近下壁面处的气体延迟较上层气体更大。图5（d）表明滑移流域与近连续流域，在温度梯度驱动作用下最大y向速度

{\bar{u}}_{y}

振动幅值分别为0.035与0.022；当Kn取1，10时，

{\bar{u}}_{y}

振动幅值分别为0.013与0.012。这是由于在近连续流域与早期过渡流域，分子间碰撞频率高，能量传导快，温度梯度产生的驱动力相较高稀薄流域更大，因此速度幅值也更大。由图5可以发现Kn对密度与速度振荡幅值的影响远大于其对温度振荡幅值的影响。密度与速度振荡幅值随Kn变化最大偏差分别为5.8%与87.25%。

在滑移流域与早期过渡流域（Kn取0.01，0.1），流场流动特性接近，宏观参数随时间变化为周期振荡余弦函数。当Kn取1，10时，平均自由程较大，气体的非平衡效应显著，放大了问题的非线性特性。此时x方向速度

{\bar{u}}_{x}

不再服从余弦函数振荡规律，上升时间远大于下降时间，同时上升沿存在两个上升台阶，如图5（c）所示。

进一步对比不同Kn条件下，在一个振荡周期内的0，0.25π，0.5π，0.75π时刻，方腔垂直中线处速度

{\bar{u}}_{y}

变化曲线，结果如图6所示。

仿真结果表明，随着Kn的增大，最大

{\bar{u}}_{y}

产生的位置由靠近壁面处逐渐过渡到方腔中心线位置，受扰运动区域增大，黏性穿透层厚度也随着Kn的增大而增大。这是因为对于硬球模型，黏性穿透层厚度δ_u约为：

δ_{u} = \sqrt{\frac{K n \cdot H^{2}}{S t \sqrt{π}}}

(18)

图6方腔垂直中线处速度

{\bar{u}}_{y}

（x=L/2）变化曲线

Fig.6Variation curve of velocity

{\bar{u}}_{y}

(x=L/2) at the central vertical of the cavity

当St一定而Kn增大时，顶盖振荡热驱的作用距离与分子平均自由程成正相关，如图6所示。同时，从图6（a）~（d）可以看出，从滑移流域到过渡流域，

{\bar{u}}_{y}

的最大值从0.041（Kn=0.01）先增加到0.062（Kn=0.1），然后随着Kn的增大减小到0.031（Kn=10）。图5、图6表明Kn大于1之后，流场参数变化对Kn的继续增加不敏感。

在t取0，0.25π与t取0.5π，0.75π时刻，

{\bar{u}}_{y}

变化曲线并非严格对称，且随着Kn的增大，不对称性加剧。这是由于顶盖升温过程中，速度的驱动力在高温顶盖向其余三个低温壁面的传热传质过程中产生；在顶盖降温过程中，速度的驱动力在三个高温壁面向低温顶盖传热传质过程中产生。前半周期与后半周期的驱动力大小不同，当非平衡效应显著时，加剧了前、后半周期的非对称现象。

3.3 顶盖温度振荡频率对流场宏观参数的影响

当顶盖温度振荡频率不同时，分析振荡频率对流场参数的影响规律。当A=1，Kn=0.1，St为2，4，10，30，70时，方腔垂直中线处无量纲y向速度、温度、压力与热流分布如图7所示，定义y向平均速度

|{\bar{u}}_{y}|

为：

|{\bar{u}}_{y}| (x = L / 2) = \frac{1}{A} \int_{0}^{A} {\bar{u}}_{y} (x = L / 2) d y

(19)

由图7（b）可以看出，随着顶盖温度振荡频率的增大，顶盖处温度跳跃增大。具体地，St取2，4，10，30，70时，相对顶盖温度的温度跳跃分别为：0.115，0.129，0.156，0.169，0.173。热流沿方腔垂直中线的变化规律与温度变化规律相近，如图7（d）所示。同时可以发现，对于低频振荡条件（St≤10），顶盖温度变化能够传导至底面处，最大y向速度

|{\bar{u}}_{y}|

发生位置靠近方腔中心位置。但高频振荡条件下（St≥30），方腔中下部（y/H＜0.5）温度几乎不随顶盖温度振荡频率变化而变化，此时方腔中下部的y向速度亦趋于0，最大

|{\bar{u}}_{y}|

发生位置为靠近顶盖处，如图7（a）~（b）所示。这是由于当Kn一定时，黏性穿透层厚度δ_u随着St的增大而减小，如式（18）所示，当St≥30时，受扰运动区域显著减小。

图7不同St下垂直中线处（x=L/2）宏观参数的包络线

Fig.7Envelope curves of macroscopic parameters along the vertical centerline where x=L/2 under different values of St

St取2，4，10，30，70对应的y向平均速度

|{\bar{u}}_{y}|

分别为：0.041，0.029，0.021，0.008，0.003。显然，随着St的增加，顶盖温度变化频率远大于分子平均碰撞频率，在一个变化周期内，扰动能量没有充分的时间通过分子间碰撞传播到更远的区域，振荡热驱产生的y向平均速度

|{\bar{u}}_{y}|

逐渐降低。进一步分析当A=1，Kn取0.01，0.1，1.0，10时，y向平均速度

|{\bar{u}}_{y}|

随St与Kn的变化曲线，结果如图8所示。

图8表明：从滑移流域到早期过渡流域，在相同Kn下，流场宏观参数的响应延迟随着St的增大而增大（一个振荡周期内，分子碰撞次数减少，扰动传播能力减小，纵向延迟增大），因此y向平均速度

|{\bar{u}}_{y}|

随着St的增大而减小。但当Kn≥1时，y向平均速度

|{\bar{u}}_{y}|

随着St的增大先增大后减小，后趋于平稳。当St≥40时，Kn越大分子平均自由程越大，当分子间碰撞次数一定时，扰动传播距离远，因此

|{\bar{u}}_{y}|

随着Kn的增大而增大。

对比A=1，Kn=0.1，St取2，4，10，30，70时，方腔顶盖处切应力

{\bar{P}}_{x y}

分布，结果如图9所示。定义y向平均切应力

|{\bar{P}}_{x y}|

为：

|{\bar{P}}_{x y}| (y = H) = \frac{1}{A} \int_{0}^{A} {\bar{P}}_{x y} (y = H) d x

(20)

图8不同Kn下平均速度

|{\bar{u}}_{y}|

随St的变化曲线

Fig.8Variation curve of average velocity

|{\bar{u}}_{y}|

with St under different Kn

图9不同St下方腔顶盖处切应力

{\bar{P}}_{x y}

包络线

Fig.9Envelope curves of shear stress

{\bar{P}}_{x y}

at the lid of the cavity for different values of St

从图9可以看出：切应力在方腔的左上与右上方达到极值，最大值约为0.024，随着St的增大，顶盖处的切应力逐渐减小。在顶盖温度高频振荡条件下（St≥40时），方腔中心线处的切应力

{\bar{P}}_{x y}

约为零，同时切应力曲线沿方腔中心线对称分布。St取2，4，10，30，70的顶盖平均切应力

|{\bar{P}}_{x y}|

分别为：0.003 0，0.002 7，0.001 9，0.001 1，0.000 8。当Kn取0.01，0.1，1.0，10时，顶盖平均切应力

|{\bar{P}}_{x y}|

随St与Kn的变化曲线如图10所示。

与图8结论相似，在大振荡频率条件下（St≥15），

|{\bar{P}}_{x y}|

随着Kn的增大而增大，随着St的增大而减小。黏性穿透层厚度与扰动传播距离亦随着St的增大而减小，随着Kn的增大而增大。当St=70时，Kn取0.01，0.1，1.0，10对应的顶盖平均切应力

|{\bar{P}}_{x y}|

分别为：4.38×10^-4，7.88×10^-4，8.39×10^-4，8.42×10^-4。

图10不同Kn下平均切应力

|{\bar{P}}_{x y}|

随St的变化曲线

Fig.10Variation curve of average shear stress

|{\bar{P}}_{x y}|

with St under different Kn

4 结论

基于耦合宏观/介观方法研究了方腔中顶盖振荡热驱稀薄流动特性，得到结论如下：

1）首次将耦合宏观/介观方法的适用性推广到非定常振荡热驱稀薄流动。在模拟顶盖振荡热驱稀薄流动时，基于耦合宏观/介观方法与全流场DVM，沿垂直中线温度分布的最大计算误差为0.23%，计算内存消耗减少约69.91%。

2）耦合方法能够精确地捕获不同Kn与St下流场的非平衡效应。在大Kn流域，振荡热驱问题呈现显著非线性现象，x方向速度

{\bar{u}}_{x}

不再服从周期振荡余弦函数规律，上升时间远大于下降时间，同时上升沿存在两个上升台阶。密度与速度振荡幅值随Kn变化最大偏差分别为5.8%与87.25%。

3）顶盖位置处的温度跳跃现象随着振荡频率St的增大而增大。黏性穿透层厚度δ_u与受扰区域随着Kn的增大而增大，随着St的增大而减小，平均速度

|{\bar{u}}_{y}|

与平均切应力

|{\bar{P}}_{x y}|

与δ_u变化规律基本相符。当St≤10时，顶盖振荡温度产生的扰动能够传导至方腔底边，但当St≥30时，方腔中下部的y向速度与热流趋于0。

致谢

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

图1方腔中温度振荡驱动稀薄流动示意图

Fig.1Schematic diagram of the thermally-induced oscillatory rarefied gas flows in a rectangular cavity

下载: 全尺寸图片

图2基于耦合宏观/介观方法的温度振荡驱动稀薄流动示意图

Fig.2Schematic diagram of the thermally-induced oscillatory rarefied gas flows based on the hybrid macroscopic/mesoscopic algorithm

下载: 全尺寸图片

图3耦合宏观/介观方法计算域与数据传输方法

Fig.3Computational sub-domain and data transfer of the hybrid macroscopic/mesoscopic algorithm

下载: 全尺寸图片

图4一个振荡周期内方腔垂直中线处温度

\bar{T}

（x=L/2）变化

Fig.4Dynamic temperature profiles

\bar{T}

(x=L/2) at the central vertical of the cavity within an oscillation period

下载: 全尺寸图片

图5方腔特征点处宏观参数随时间变化

Fig.5Time evolution of the macroscopic parameters at the characteristic point of the cavity

下载: 全尺寸图片

图6方腔垂直中线处速度

{\bar{u}}_{y}

（x=L/2）变化曲线

Fig.6Variation curve of velocity

{\bar{u}}_{y}

(x=L/2) at the central vertical of the cavity

下载: 全尺寸图片

图7不同St下垂直中线处（x=L/2）宏观参数的包络线

Fig.7Envelope curves of macroscopic parameters along the vertical centerline where x=L/2 under different values of St

下载: 全尺寸图片

图8不同Kn下平均速度

|{\bar{u}}_{y}|

随St的变化曲线

Fig.8Variation curve of average velocity

|{\bar{u}}_{y}|

with St under different Kn

下载: 全尺寸图片

图9不同St下方腔顶盖处切应力

{\bar{P}}_{x y}

包络线

Fig.9Envelope curves of shear stress

{\bar{P}}_{x y}

at the lid of the cavity for different values of St

下载: 全尺寸图片

图10不同Kn下平均切应力

|{\bar{P}}_{x y}|

随St的变化曲线

Fig.10Variation curve of average shear stress

|{\bar{P}}_{x y}|

with St under different Kn

下载: 全尺寸图片

图1方腔中温度振荡驱动稀薄流动示意图

Fig.1Schematic diagram of the thermally-induced oscillatory rarefied gas flows in a rectangular cavity

图2基于耦合宏观/介观方法的温度振荡驱动稀薄流动示意图

Fig.2Schematic diagram of the thermally-induced oscillatory rarefied gas flows based on the hybrid macroscopic/mesoscopic algorithm

图3耦合宏观/介观方法计算域与数据传输方法

Fig.3Computational sub-domain and data transfer of the hybrid macroscopic/mesoscopic algorithm

图4一个振荡周期内方腔垂直中线处温度

\bar{T}

（x=L/2）变化

Fig.4Dynamic temperature profiles

\bar{T}

(x=L/2) at the central vertical of the cavity within an oscillation period

图5方腔特征点处宏观参数随时间变化

Fig.5Time evolution of the macroscopic parameters at the characteristic point of the cavity

图6方腔垂直中线处速度

{\bar{u}}_{y}

（x=L/2）变化曲线

Fig.6Variation curve of velocity

{\bar{u}}_{y}

(x=L/2) at the central vertical of the cavity

图7不同St下垂直中线处（x=L/2）宏观参数的包络线

Fig.7Envelope curves of macroscopic parameters along the vertical centerline where x=L/2 under different values of St

图8不同Kn下平均速度

|{\bar{u}}_{y}|

随St的变化曲线

Fig.8Variation curve of average velocity

|{\bar{u}}_{y}|

with St under different Kn

图9不同St下方腔顶盖处切应力

{\bar{P}}_{x y}

包络线

Fig.9Envelope curves of shear stress

{\bar{P}}_{x y}

at the lid of the cavity for different values of St

图10不同Kn下平均切应力

|{\bar{P}}_{x y}|

随St的变化曲线

Fig.10Variation curve of average shear stress

|{\bar{P}}_{x y}|

with St under different Kn

引用提醒

图(10)

引用本文

杨伟奇, 钮耀斌. 耦合宏观/介观方法的振荡热驱稀薄流动[J]. 国防科技大学学报, 2025, 47(6): 168-177.

复制

YANG W Q, NIU Y B. Hybrid macroscopic/mesoscopic scheme for oscillatory thermally driven rarefied gas flows[J]. Journal of National University of Defense Technology, 2025, 47(6): 168-177.

Copy

计量

图1方腔中温度振荡驱动稀薄流动示意图

Fig.1Schematic diagram of the thermally-induced oscillatory rarefied gas flows in a rectangular cavity

图2基于耦合宏观/介观方法的温度振荡驱动稀薄流动示意图

Fig.2Schematic diagram of the thermally-induced oscillatory rarefied gas flows based on the hybrid macroscopic/mesoscopic algorithm

图3耦合宏观/介观方法计算域与数据传输方法

Fig.3Computational sub-domain and data transfer of the hybrid macroscopic/mesoscopic algorithm

图4一个振荡周期内方腔垂直中线处温度

\bar{T}

（x=L/2）变化

Fig.4Dynamic temperature profiles

\bar{T}

(x=L/2) at the central vertical of the cavity within an oscillation period

图5方腔特征点处宏观参数随时间变化

Fig.5Time evolution of the macroscopic parameters at the characteristic point of the cavity

图6方腔垂直中线处速度

{\bar{u}}_{y}

（x=L/2）变化曲线

Fig.6Variation curve of velocity

{\bar{u}}_{y}

(x=L/2) at the central vertical of the cavity

图7不同St下垂直中线处（x=L/2）宏观参数的包络线

Fig.7Envelope curves of macroscopic parameters along the vertical centerline where x=L/2 under different values of St

图8不同Kn下平均速度

|{\bar{u}}_{y}|

随St的变化曲线

Fig.8Variation curve of average velocity

|{\bar{u}}_{y}|

with St under different Kn

图9不同St下方腔顶盖处切应力

{\bar{P}}_{x y}

包络线

Fig.9Envelope curves of shear stress

{\bar{P}}_{x y}

at the lid of the cavity for different values of St

图10不同Kn下平均切应力

|{\bar{P}}_{x y}|

随St的变化曲线

Fig.10Variation curve of average shear stress

|{\bar{P}}_{x y}|

with St under different Kn

苑伟政, 常洪龙, 谢建兵. MEMS集成设计与制造技术进展[J]. 机械工程学报,2023,59(19):176-186.YUAN W Z, CHANG H L, XIE J B. Advances in MEMS integrated design and manufacturing technology[J]. Journal of Mechanical Engineering,2023,59(19):176-186.(in Chinese)

张旭辉. 航天微系统模块关键共性支撑技术研究与实现[D]. 西安: 西北大学,2022:2-32.ZHANG X H. Research and implementation of key common support technology for aerospace microsystem module[D]. Xi′an: Northwest University,2022:2-32.(in Chinese)

ZENGERLE T, ABABNEH A, SEIDEL H. Generalized damping model for MEMS oscillators from molecular to viscous flow regime[J]. Eng,2022,3(1):124-141.

HATTORI M. Poiseuille and thermal transpiration flows of a dense gas between two parallel plates[J]. Journal of Fluid Mechanics,2023,962: A20.

ZHAO L L, WANG X W, ZHANG Z J. Numerical investigation into the compression characteristics of a multi-stage Knudsen pump with rectangular channels[J]. The European Physical Journal Plus,2023,138(3):302.

王燕. MEMS流量传感器温度补偿方法研究[D]. 济南: 山东大学,2022.WANG Y. Study of temperature compensation methods for MEMS flow sensors[D]. Jinan: Shandong University,2022.(in Chinese)

EJEIAN F, AZADI S, RAZMJOU A,et al. Design and applications of MEMS flow sensors:a review[J]. Sensors and Actuators A: Physical,2019,295:483-502.

BIRD G A. Molecular gas dynamics and the direct simulation of gas flows[M]. New York: Oxford University Press,1994:15-113.

EMERSON D R, GU X J, STEFANOV S K,et al. Nonplanar oscillatory shear flow:from the continuum to the free-molecular regime[J]. Physics of Fluids,2007,19(10):107105.

SHARIPOV F, KALEMPA D. Oscillatory Couette flow at arbitrary oscillation frequency over the whole range of the Knudsen number[J]. Microfluidics and Nanofluidics,2008,4(5):363-374.

DUCK P W. Oscillatory flow inside a square cavity[J]. Journal of Fluid Mechanics,1982,122:215-234.

BIRINGEN S, DANABASOGLU G. Oscillatory flow with heat transfer in a square cavity[J]. Physics of Fluids A: Fluid Dynamics,1989,1(11):1796-1812.

WU L, REESE J M, ZHANG Y H. Oscillatory rarefied gas flow inside rectangular cavities[J]. Journal of Fluid Mechanics,2014,748:350-367.

LOBASOV A S, MINAKOV A V, RUDYAK V Y. The investigation of the velocity slip and the temperature jump effect on the heat transfer characteristics in a microchannel[J]. Case Studies in Thermal Engineering,2022,31:101791.

USAFZAI W K, ALY E H, ALSHOMRANI A S,et al. Multiple solutions for nanofluids flow and heat transfer in porous medium with velocity slip and temperature jump[J]. International Communications in Heat and Mass Transfer,2022,131:105831.

GUO Z L, SHI B C, ZHENG C G. An extended Navier-Stokes formulation for gas flows in the Knudsen layer near a wall[J]. Europhysics Letters,2007,80(2):24001.

TANG G H, ZHANG Y H, GU X J,et al. Lattice Boltzmann modelling Knudsen layer effect in non-equilibrium flows[J]. Europhysics Letters,2008,83(4):40008.

BOBYLEV A V. On some properties of generalized Burnett equations of hydrodynamics[J]. Journal of Statistical Physics,2022,188(1):6.

JADHAV R S, YADAV U, AGRAWAL A. OBurnett equations:thermodynamically consistent continuum theory beyond the Navier-Stokes regime[J]. Journal of Heat and Mass Transfer,2023,145(6):060801.

KHAN N, KHAN M I, HUSAIN S,et al. Certain results on generalized Humbert-Hermite polynomials[J]. Research in Mathematics,2023,10(1):2215577.

ARTALE V, CONFORTO F, MARTALÒ G,et al. Shock structure solutions of Grad 13-moment equations for binary gas mixtures[J]. Wave Motion,2022,115:103055.

STRUCHTRUP H, FREZZOTTI A. Twenty-six moment equations for the Enskog-Vlasov equation[J]. Journal of Fluid Mechanics,2022,940: A40.

STRUCHTRUP H,ÖTTINGER H C. Thermodynamically admissible 13-moment equations[J]. Physics of Fluids,2022,34:017105.

GU X J, EMERSON D R. Modeling oscillatory flows in the transition regime using a high-order moment method[J]. Microfluidics and Nanofluidics,2011,10(2):389-401.

TATSIOS G, TSIMPOUKIS A, VALOUGEORGIS D. The half-range moment method in harmonically oscillating rarefied gas flows[J]. Fluids,2021,6(1):17.

PARK J H, BAHUKUDUMBI P, BESKOK A. Rarefaction effects on shear driven oscillatory gas flows:a direct simulation Monte Carlo study in the entire Knudsen regime[J]. Physics of Fluids,2004,16(2):317-330.

BAHUKUDUMBI P, PARK J H, BESKOK A. A unified engineering model for steady and quasi-steady shear-driven gas microflows[J]. Microscale Thermophysical Engineering,2003,7(4):291-315.

WU L. Sound propagation through a rarefied gas in rectangular channels[J]. Physical Review E,2016,94:053110.

ZHANG Q D, ZHUO C S, MU J L,et al. A multiscale discrete velocity method for diatomic molecular gas[J]. Physics of Fluids,2023,35(7):076118.

FRANGI A, FREZZOTTI A, LORENZANI S. On the application of the BGK kinetic model to the analysis of gas-structure interactions in MEMS[J]. Computers & Structures,2007,85(11/12/13/14):810-817.

TSIMPOUKIS A, VASILEIADIS N, TATSIOS G,et al. Nonlinear oscillatory fully-developed rarefied gas flow in plane geometry[J]. Physics of Fluids,2019,31(6):067108.

WANG P, ZHU L H, SU W,et al. Nonlinear oscillatory rarefied gas flow inside a rectangular cavity[J]. Physical Review E,2018,97:043103.

WANG P, SU W, ZHUANG Y H. Oscillatory rarefied gas flow inside a three dimensional rectangular cavity[J]. Physics of Fluids,2018,30(10):102002.

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.

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

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.

GU X J, EMERSON D R. A computational strategy for the regularized 13 moment equations with enhanced wall-boundary conditions[J]. Journal of Computational Physics,2007,225(1):263-283.

TORRILHON M, STRUCHTRUP H. Boundary conditions for regularized 13-moment-equations for micro-channel-flows[J]. Journal of Computational Physics,2008,227(3):1982-2011.

首页

期刊介绍

投稿指南

编委会

出版声明

开放获取声明

联系我们

期刊订阅

Rss

AI检索

English

1 问题描述

2 数值方法

3 算例与讨论

4 结论