CONSTRUCRION AND ANALYSIS OF A NEW COMPUTABLE MODEL FOR BOLTZMANN EQUATION
-
摘要: 临近空间飞行器因各部件尺寸差异较大, 在高空高马赫数条件下飞行会出现多流区共存的多尺度复杂非平衡流动现象, 流场中的气体分子速度分布函数与当地位置、流场分子速度、气体密度、流动速度、温度、热流矢量、应力张量等相关. 通过分析玻尔兹曼方程的一阶查普曼−恩斯科近似解, 构造了一种同时考虑热流矢量和应力张量影响、满足玻尔兹曼方程高阶碰撞矩的跨流域统一可计算模型方程, 并在数学上分析了其守恒性、H定理等基本属性, 证明了新模型方程与玻尔兹曼方程的相容性, 给出了新模型与现有模型如沙克霍夫(Shakhov)模型等的递进关系, 基于碰撞动力学确定了各流域统一气体分子碰撞松弛参数表达式. 在气体动理论统一算法中采用新建模型及现有模型模拟了一维激波结构、二维近空间飞行环境平板和多体圆柱干扰流动, 并与直接模拟蒙特卡洛方法对比分析, 结果表明在流场中粘性效应显著的区域新建模型能更好地捕捉激波位置, 尤其是在激波内部新模型模拟的宏观参数分布与直接模拟蒙特卡洛方法结果符合更好, 验证了新模型的有效性和可靠性, 同时说明在非平衡显著的流动区域碰撞松弛模型受多参数共同作用的影响.
-
关键词:
- 玻尔兹曼方程 /
- 气体动理论统一算法 /
- 查普曼−恩斯科近似解 /
- 考虑热流矢量与应力张量影响的碰撞模型 /
- 气体分子碰撞松弛参数
Abstract: Due to large differences of geometric scale between the components of near space vehicles, the multi-scale complex non-equilibrium flow phenomenon will appear in many flow field regions when vehicles flying at high Mach number and high altitudes. In those regions the gas molecular velocity distribution functions are related to the local molecular velocities and macroscopic parameters, such as velocities, temperatures, heat flux vectors and stress tensors. By analyzing the first-order Chapman−Enskog approximate solution of Boltzmann equation, a new computable collision relaxation model is constructed, which considers the influence of heat flux vector and stress tensor, and satisfies the high-order collision moments of Boltzmann equation. The basic properties such as conservation law and H theorem are analyzed mathematically. The compatibility between the new model equation and Boltzmann equation is proved. The relationships between the new model and old models such as Shakhov and Belyi models are given. The expression of the collision relaxation parameter is determined by using molecular collision dynamics. As examples, one dimensional shock profiles and two dimensional flows around a flat plate and two side-by-side cylinders in near space environments are simulated by gas kinetic unified algorithm with different models. By comparing with results of DSMC method, it shows that in one dimensional problems the results of Shakhov model with heat flux is better than the new model because of smaller 1D shear stress, but in two dimension the new model can capture the position of shock wave better than the other two models due to higher dimensional shear stresses leading to more distinct viscosity effect. Especially for macro parameters in shock waves, the results of the new model accord with those of DSMC better. Then the validity and reliability of the new model is verified. These results illuminate that the collision relaxation model is influenced by multi-parameters together in the flow field when the non-equilibrium effects are quite distinct. -
引言
近年来, 随着对临近空间战略意义认识的加深, 各航空航天大国发展了众多临近空间飞行器. 为满足长航时、大载荷、高超声速巡航等需求, 这类临近空间飞行器通常整体尺寸较大, 各部件间差异显著, 飞行过程会出现多物理多场耦合、多流区共存的多尺度复杂流动现象, 传统的数值模拟方法如基于宏观连续介质假设的纳维−斯托克斯(N-S)方程、模拟微观分子运动与碰撞的直接模拟蒙特卡洛(DSMC)方法难以满足多尺度一体化模拟需求, 需要研究发展新的模拟手段[1-4]. 本文从介观气体分子速度分布函数满足的方程即玻尔兹曼方程[5-7]出发开展多尺度流动一体化模拟研究.
玻尔兹曼方程作为气体动理学理论的基本方程, 能描述各流域气体流动输运现象, 但它是一个高度非线性的积分−微分方程, 其理论解和直接数值求解均难以实现[7]. 为利用玻尔兹曼方程能统一描述各流域气体流动的特性, 众多学者提出了系列替代求解技术, 其中应用最广泛的是模型方程方法, 即用一个相对简单的表达式替换玻尔兹曼方程碰撞积分项[8], 再对简化后的偏微分方程进行求解.
模型方程作为玻尔兹曼方程的近似, 必须满足玻尔兹曼方程的基本属性, 如总和不变量性质、H定理、正确的输运系数等. 第一个也是最著名的模型方程是巴特纳−格罗斯−克鲁克模型[8](BGK模型), 但通过对其进行查普曼−恩斯科(Chapman−Enskog)展开分析得到的普朗特数(Pr)为1, 而单原子气体的Pr约为2/3[9]. 为克服这一缺陷, 国内外学者提出了一系列模型方程, 包括椭球统计−BGK (ES-BGK)模型[10-11]、沙克霍夫(Shakhov)模型[12-13]、刘(Liu)模型[14]、贝利(Belyi)模型[15]等, 这些模型方程均假定其碰撞频率与分子速度无关. 此外, 也有学者发展了考虑碰撞频率与分子速度相关性的模型方程[16]以及适于多原子气体、考虑热力学非平衡的模型方程[17].
大量数学分析表明[16, 18-20], 上述模型方程均能满足玻尔兹曼方程的基本属性, 且与原始玻尔兹曼方程的碰撞项具有类似的松弛特性, 能保持与纳维−斯托克斯方程、直接模拟蒙特卡洛方法的收敛一致, 因结构简单、贴近跨流域流体物理变化特点、可计算技术针对性强而得到广泛应用[21], 如泰塔雷夫(Titarev)建立了物理空间和速度空间均采用非结构网格的内斯韦塔依(Nesvetay)求解器[22-23], 针对返回式空间飞行器等高马赫数流动与直接模拟蒙特卡洛方法进行了系列对比分析. 李志辉等[24-31]从研究描述各流域均适用的气体分子速度分布函数方程出发, 吸收计算数学指数型积分求解原理, 发展了气体分子运动论离散速度坐标法, 研制经改进的高斯−埃尔米特积分方法等系列离散速度数值积分技术, 消除气体分子速度分布函数对速度空间的连续依赖性, 将玻尔兹曼模型方程化为在各个离散速度坐标点处具有非线性源项的双曲型守恒方程, 基于计算流体力学数值求解技术, 构造直接求解速度分布函数演化更新的气体动理论数值格式和大规模并行计算策略, 建立了气体动理论统一算法(GKUA), 开展了自由分子流到连续流复杂高超声速气动力/热问题大规模并行计算及其在航天工程中的应用研究, 同时对跨流域多体干扰流动问题、跨流域非定常流动问题进行了研究, 获得较好的成功. 徐昆等[32-37]利用上述离散速度坐标法与BGK类型格式相耦合, 发展了统一气体动理学格式(UGKS), 使之适于模拟稀薄流到连续流气体流动问题.
另一方面, 从流动反映的物理现象来看, 当流场中出现非平衡现象时, 某点的速度分布函数已经不再是简单的麦克斯韦平衡态分布, 而是呈现各种形态, 如双峰分布、偏心分布等[38]. 这些分布与当地的气体分子速度、流场宏观速度、温度以及热流矢量、应力张量等宏观输运性质相关, 这从玻尔兹曼方程的查普曼−恩斯科分析解[5, 39]可以看出. 因此在对玻尔兹曼方程碰撞积分进行物理分析、构造碰撞项的模型方程时应同时考虑这些因素, 而常用的玻尔兹曼碰撞积分模型如沙克霍夫模型没有考虑应力张量的影响, ES-BGK模型属于各向异性的类麦克斯韦分布, 没有考虑热流矢量的影响. 当模型方程不能全面反映这些因素时会导致模拟结果出现系统误差, 例如在稀薄过渡流区流动的克努森层内黏性效应占主导, 应考虑应力张量的影响, 沙克霍夫模型与直接模拟蒙特卡洛方法模拟结果存在较大差别[40].
为此, 本文将从玻尔兹曼方程出发, 基于查普曼−恩斯科渐近分析解[9]的数学推导, 构造一种同时考虑热流矢量和应力张量影响、满足玻尔兹曼方程高阶碰撞矩的综合型可计算碰撞模型, 并在数学上分析其基本属性, 建立与玻尔兹曼方程的联系, 给出新模型与现有模型方程的关系, 同时基于碰撞动力学确定该可计算模型碰撞松弛参数表达式. 最后在气体动理论统一算法框架下, 使用该新建玻尔兹曼方程可计算模型, 通过模拟一维激波结构、二维近空间飞行环境平板和多体圆柱干扰流动来验证新模型的有效性和可靠性.
1. 模型方程的构造
玻尔兹曼方程是描述气体分子速度分布函数
$f\left( {{\boldsymbol{V}}({V_x},{V_y},{V_{\textit{z}}});{\boldsymbol{r}}(x,y,{\textit{z}});t} \right)$ 在位置空间与速度空间随时间演化更新的控制方程, 其一般形式[39]如下$$ \frac{{\partial f}}{{\partial t}} + \sum\limits_i {{V_i}\frac{{\partial f}}{{\partial {r_i}}}} = S\left( f \right) = {\left( {\frac{{\partial f}}{{\partial t}}} \right)_{\rm{col}}} $$ (1) 其中,
$i = x$ ,$y$ ,${\textit{z}}$ ,$ {V_x} $ ,$ {V_y} $ ,$ {V_z} $ 为气体分子速度分量;$ x $ ,$ y $ ,${\textit{z}}$ 为空间位置坐标,$ t $ 为时间;$S\left( f \right)$ ,$ {\left( {\partial f/\partial t} \right)_{\rm{col} }} $ 为碰撞积分项.在确定速度分布函数
$f$ 后, 可以通过对其取矩得到宏观流动参数, 如密度$ \rho $ , 宏观速度分量$ {U_x} $ ,$ {U_y} $ ,$ {U_z} $ , 温度$ T $ , 应力张量分量$ {P_{ij}} $ ; 热流矢量分量$ {q_x} $ ,$ {q_y} $ ,${q_{\textit{z}}}$ 分别为$$ \rho = mn = \iiint_{{{\mathbb{R}}^3}} {mf{\rm{d}} W} $$ (2) $$ \rho {U_i} = mn{U_i} = \iiint_{{{\mathbb{R}}^3}} {m{V_i}f{\rm{d}} W} $$ (3) $$ \rho {e_{\rm{tr}}} = \frac{3}{2}n{{{k}}_{\rm{B}}}T = \iiint_{{{\mathbb{R}}^3}} {\frac{m}{2}{{\left| {\boldsymbol{C}} \right|}^2}f{\rm{d}} W} $$ (4) $$ {P_{ij}} = \iiint_{{{\mathbb{R}}^3}} {m{C_i}{C_j}f{\rm{d}} W} $$ (5) $$ {q_i} = \iiint_{{{\mathbb{R}}^3}} {\frac{m}{2}{C_i}{{\left| {\boldsymbol{C}} \right|}^2}f{\rm{d}} W} $$ (6) 其中,
$ {{\mathbb{R}}^3} $ 为三维实数空间,${{\rm{d}}} W = {{\rm{d}}} {V_x}{{\rm{d}}} {V_y}{{\rm{d}}} {V_{\textit{z}}}$ (下同),$ m $ 为分子质量,$ n $ 为分子数密度,${{{{k}}}_{\rm{B}}}$ 为玻尔兹曼常数;$i$ ,$j = x$ ,$y$ ,${\textit{z}}$ ;$ {\boldsymbol{C}} $ 为分子热运动速度,$ {C_i} = {V_i} - {U_i} $ ,$ {\left| {\boldsymbol{C}} \right|^2} = {C^2} = C_x^2 + C_y^2 + C_{\textit{z}}^2 $ ,${e_{\rm{tr}}}$ 为单位质量动量. 此外流场压力$p$ , 黏性切应力张量分量$ {\tau _{ij}} $ 分别为$$ p = n{{{k}}_{\rm{B}}}T \quad\qquad$$ (7) $$ {\tau _{ij}} = - \left( {{P_{ij}} - {\delta _{ij}}p} \right) $$ (8) 这里, 符号
$ {\delta _{ij}} $ 表示当$i = j$ 时为1, 当$i \ne j$ 时为0.气体动理论统一算法所用玻尔兹曼方程可计算模型[24-31]可以写为
$$ \frac{{\partial f}}{{\partial t}} + \sum\limits_i {{V_i}\frac{{\partial f}}{{\partial {r_i}}}} = {S_{\rm{m}} }\left( f \right) = \nu \left( {{f^{\rm{N}} } - f} \right) $$ (9) 这里,
$ \nu $ 为气体分子碰撞松弛参数,$ {f^{\rm{N}} } $ 为当地平衡态速度分布函数. 常用的沙克霍夫模型[12]和贝利模型[15]的表达式分别为$$ f_{\rm{S}} ^{\rm{N}} = {f_{\rm{M}} }\left[ {1 + \left( {1 - Pr} \right)\frac{{m\displaystyle\sum\limits_i {{C_i}{q_i}} }}{{p{{{k}}_{\rm{B}}}T}}\left( {\frac{{m{C^2}}}{{5{{{k}}_{\rm{B}}}T}} - 1} \right)} \right] $$ (10) $$ f_{\rm{VVB}}^{\rm{N}} = {f_{\rm{M}} }\left[ {1 + \frac{{Pr - 1}}{{Pr}}\left( { - m\sum\limits_{i,j} {\frac{{{\tau _{ij}}{C_i}{C_j}}}{{2p{{{k}}_{\rm{B}}}T}}} } \right)} \right]\quad\;\;\; $$ (11) 这里
${f_{\rm{M}} }$ 为麦克斯韦速度分布函数, 即$$ {f_{\rm{M}} } = \frac{n}{{{{\left( {2{{\text{π}} }{{{k}}_{\rm{B}}}T/m} \right)}^{3/2}}}}{{\rm{e}} ^{ - \frac{{m{C^2}}}{{2{{{k}}_{{\rm{B}}}}T}}}} $$ (12) 任何一个模型方程都必须满足玻尔兹曼方程碰撞积分的如下3个基本属性[39].
(1)守恒性(质量、动量、能量守恒), 即
$$ \iiint_{{{\mathbb{R}}^3}} {m{S_{\rm{m}} }\left( f \right){\rm{d}} W} = 0 \;\;$$ (13) $$ \iiint_{{{\mathbb{R}}^3}} {m{V_i}{S_{\rm{m}} }\left( f \right){\rm{d}} W} = 0 $$ (14) $$ \iiint_{{{\mathbb{R}}^3}} {\frac{m}{2}{V^2}{S_{\rm{m}} }\left( f \right){\rm{d}} W} = 0 $$ (15) 其中
${V^2} = V_x^2 + V_y^2 + V_{\textit{z}}^2$ .(2) H定理, 即
$$ H\left( t \right) = \iiint_{{{\mathbb{R}}^3}} {f\ln f{\rm{d}} W} $$ (16) $$ \frac{\partial }{{\partial t}}H\left( t \right) \leqslant 0 \qquad\qquad$$ (17) (3)在平衡态时, 有
${S_{\rm{m}} }\left( {{f_{\rm{M}} }} \right) = 0$ .此外, 对于单原子气体, 模型方程的查普曼−恩斯科分析得到的普朗特数
$ Pr $ 应约为2/3.为研究构造新的玻尔兹曼方程可计算模型, 首先分析其一阶查普曼−恩斯科分析解[9], 其可以写为
$$ f = {f_{\rm{M}} }\left[ {1 + \varPhi \left( {\boldsymbol{C}} \right)} \right]\qquad\qquad\qquad\qquad $$ (18) $$ \begin{split} \varPhi \left( {\boldsymbol{C}} \right) = &\left[ {\frac{m}{{p{{{k}}_{\rm{B}}}T}}\left( {\frac{{m{C^2}}}{{5{{{k}}_{\rm{B}}}T}} - 1} \right)\sum\limits_i {{C_i}{q_i}} } \right] - \\& { m\sum\limits_{i,j} {\frac{{{\tau _{ij}}{C_i}{C_j}}}{{2p{{{k}}_{\rm{B}}}T}}} } \end{split} $$ (19) 显然,
$\varPhi \left( {\boldsymbol{C}} \right)$ 的第一项与常用的沙克霍夫模型[12]中第二项除系数外表达式一致, 而$\varPhi \left( {\boldsymbol{C}} \right)$ 的第二项与贝利模型[15]的第二项除系数外表达式也一致.从玻尔兹曼方程的查普曼−恩斯科分析解可以看出, 流场中气体的速度分布函数与局部宏观流动参数如速度、热流、应力等相关, 而沙克霍夫模型没有考虑应力的影响, 贝利模型没有考虑热流的影响. 当流场中局部区域的温度、速度梯度较大即热流和应力影响显著时, 将出现较为严重的非平衡效应, 仅考虑单一因素难以准确模拟复杂非平衡现象. 因此, 可以依据玻尔兹曼方程的查普曼−恩斯科分析解以及常用模型方程的表达式, 构造如下同时考虑热流和应力影响的新型玻尔兹曼方程可计算模型
$$ \frac{{\partial f}}{{\partial t}} + \sum\limits_i {{V_i}\frac{{\partial f}}{{\partial {r_i}}}} = {\nu ^{\rm{N}}}\left( {f_{\rm{mix}}^{\rm{N}} - f} \right) $$ (20) $$ f_{\rm{mix}}^{\rm{N}} = {f_{\rm{M}}}\left[ {1 + A \cdot {\varPhi _{\rm{A}} }\left( {\boldsymbol{C}} \right) + B \cdot {\varPhi _{\rm{B}} }\left( {\boldsymbol{C}} \right)} \right] $$ (21) $$ {\varPhi _{\rm{A}} }\left( {\boldsymbol{C}} \right) = \frac{m}{{p{{{k}}_{\rm{B}}}T}}\left( {\frac{{m{C^2}}}{{5{{{k}}_{\rm{B}}}T}} - 1} \right)\sum\limits_i {{C_i}{q_i}} $$ (22) $$ {\varPhi _{\rm{B}} }\left( {\boldsymbol{C}} \right) = - m\sum\limits_{i,j} {\frac{{{\tau _{ij}}{C_i}{C_j}}}{{2p{{{k}}_{\rm{B}}}T}}} $$ (23) 这里, A和B是与分子速度无关的常量,
$ {\nu ^{\rm{N}}} $ 为碰撞松弛时间的倒数. 显然, 当$ A = 1 - Pr $ ,$ B = 0 $ 时对应沙克霍夫模型[12],$ A = 0 $ ,$ B = \left( {Pr - 1} \right)/Pr $ 时对应贝利模型[15].2. 新型玻尔兹曼模型方程的数学分析
在构造玻尔兹曼方程的一个新模型后, 需要对其性质进行分析, 验证其是否满足玻尔兹曼方程碰撞积分项的基本属性.
2.1 守恒性
由于
$$ \iiint_{{{\mathbb{R}}^3}} {mf_{\rm{mix}}^{\rm{N}}{\rm{d}} W} = \rho $$ (24) $$ \iiint_{{{\mathbb{R}}^3}} {m{V_i}f_{\rm{mix}}^{\rm{N}}{\rm{d}} W} = mn{U_i} = \rho {U_i} $$ (25) $$ \iiint_{{{\mathbb{R}}^3}} {\frac{m}{2}{V^2}f_{\rm{mix}}^{\rm{N}}{\rm{d}} W} = \frac{{3n{{{k}}_{\rm{B}}}T + mn{U^2}}}{2} $$ (26) $$ \iiint_{{{\mathbb{R}}^3}} {m{C_i}{C_j}f_{\rm{mix}}^{\rm{N}}{\rm{d}} W} = p{\delta _{ij}} - B{\tau _{ij}} $$ (27) $$ \iiint_{{{\mathbb{R}}^3}} {m{C_i}{C^2}f_{\rm{mix}}^{\rm{N}}{\rm{d}} W} = 2A{q_i} $$ (28) 其中
${U^2} = U_x^2 + U_y^2 + U_{\textit{z}}^2$ . 据此对式(20)按$ m $ ,$ m{V_i} $ ,$ m{V^2}/2 $ 取矩, 可得$$ \frac{{\partial \rho }}{{\partial t}} + \sum\limits_i {\frac{\partial }{{\partial {r_i}}}\left( {\rho {U_i}} \right)} = 0 $$ (29) $$ \frac{{\partial \left( {\rho {U_i}} \right)}}{{\partial t}} + \sum\limits_j {\frac{\partial }{{\partial {r_j}}}\left( {\rho {U_i}{U_j}} \right)} = - \frac{{\partial p}}{{\partial {r_i}}} + \sum\limits_j {\frac{{\partial {\tau _{ij}}}}{{\partial {r_j}}}} $$ (30) $$ \begin{split} & \frac{\partial }{{\partial t}}\rho \left( {{e_{\rm{tr}}} + \frac{1}{2}{U^2}} \right) + \sum\limits_j {\frac{\partial }{{\partial {r_j}}}\left[ {\rho {U_j}\left( {h + \frac{1}{2}{U^2}} \right)} \right]}= \\&\qquad { \sum\limits_j {\frac{\partial }{{\partial {r_j}}}\left[ {\sum\limits_k {\left( {{\tau _{jk}} - {U_k}} \right)} - {q_j}} \right]} } \end{split} $$ (31) 其中
$h = {e_{\rm{tr}}} + p/\rho $ . 上述3个方程分别对应质量、动量和能量守恒方程. 由此可知新模型方程满足玻尔兹曼方程碰撞积分项的守恒性.在平衡态时热流和应力张量均为零, 显然此时当地麦克斯韦分布是模型方程的解.
2.2 查普曼−恩斯科分析
为开展查普曼−恩斯科分析, 首先将新模型方程的无量纲形式写成如下形式
$$ \xi \left( {\frac{{\partial \hat f}}{{\partial \hat t}} + \sum\limits_i {{{\hat V}_i}\frac{{\partial \hat f}}{{\partial {{\hat r}_i}}}} } \right) = {\hat \nu ^{\rm{N}}}\left( {\hat f_{\rm{mix}}^{\rm{N}} - \hat f} \right) $$ (32) 其中
$\xi $ 为与来流克努森数(Kn)同量级的小量, 顶标$ \wedge $ 表示无量纲量. 按查普曼−恩斯科分析的流程[9], 定义$$ \hat f = {\hat f_{\rm{M}} }\left( {1 + \xi {{\hat{\phi}}_1}} \right) $$ (33) $$ f = {f_{\rm{M}} }\left( {1 + {\phi _1}} \right) $$ (34) 经过系列数学运算并回归有量纲量后, 可得
$$ \begin{split} & {\varphi _1} = A \cdot {\varPhi _{\rm{A}} }\left( {\boldsymbol{C}} \right) + B \cdot {\varPhi _{\rm{B}} }\left( {\boldsymbol{C}} \right) - \\&\qquad \frac{1}{{{\nu ^{\rm{N}} }}}\left[ \sum\limits_j {{C_j}\left( {\frac{{m{C^2}}}{{2{{{k}}_{\rm{B}}}T}} - \frac{5}{2}} \right)\frac{{\partial T}}{{T\partial {r_j}}}}+ \right.\\&\qquad\left. \frac{m}{{{{{k}}_{\rm{B}}}T}}\sum\limits_{i,j} {\left( {{C_i}{C_j} - \frac{{{C^2}}}{3}{\delta _{ij}}} \right)\frac{{\partial {U_i}}}{{\partial {r_j}}}} \right] \end{split} $$ (35) 由于
$$ \begin{split} & {\tau _{ij}} = - \left( {\iiint_{{{\mathbb{R}}^3}} {m{C_i}{C_j}f{\rm{d}} W} - {\delta _{ij}}p} \right) =\\&\qquad - \iiint_{{{\mathbb{R}}^3}} {m{C_i}{C_j}{\phi _1}{\rm{d}} W} = B{\tau _{ij}} + \\&\qquad \frac{{n{{{k}}_{\rm{B}}}T}}{{{\nu ^{\rm{N}}}}}\left( {\frac{{\partial {U_i}}}{{\partial {r_j}}} + \frac{{\partial {U_j}}}{{\partial {r_i}}} - \frac{2}{3}{\delta _{ij}}\sum\limits_k {\frac{{\partial {U_k}}}{{\partial {r_k}}}} } \right) \end{split}$$ (36) $$ \begin{split} & {q_i} = \iiint_{{{\mathbb{R}}^3}} {\frac{m}{2}{C_i}{C^2}f{\rm{d}} W} = \iiint_{{{\mathbb{R}}^3}} {\frac{m}{2}{C_i}{C^2}{\phi _1}{\rm{d}} W} = \\&\qquad A{q_i} - \frac{{n{{{k}}_{\rm{B}}}T}}{{{\nu ^{\rm{N}} }}}\frac{{5{{{k}}_{\rm{B}}}}}{{2m}}\frac{{\partial T}}{{\partial {r_i}}} \end{split} $$ (37) 则
$$ {\tau _{ij}} = \frac{{n{{{k}}_{\rm{B}}}T}}{{\left( {1 - B} \right){\nu ^{\rm{N}} }}}\left( {\frac{{\partial {U_i}}}{{\partial {r_j}}} + \frac{{\partial {U_j}}}{{\partial {r_i}}} - \frac{2}{3}{\delta _{ij}}\sum\limits_k {\frac{{\partial {U_k}}}{{\partial {r_k}}}} } \right) $$ (38) $$ {q_i} = - \frac{{n{{{k}}_{\rm{B}}}T}}{{{\nu ^{\rm{N}} }\left( {1 - A} \right)}}\frac{{5{{{k}}_{\rm{B}}}}}{{2m}}\frac{{\partial T}}{{\partial {r_i}}} $$ (39) 进而黏性系数
$\mu $ 、热传导系数$\kappa $ 、普朗特数$ Pr $ 为$$ \mu = \frac{1}{{1 - B}}\frac{p}{{{\nu ^{\rm{N}} }}} $$ (40) $$ \kappa = \frac{5}{2}\frac{{{{{k}}_{\rm{B}}}}}{m}\frac{p}{{{\nu ^{\rm{N}} }\left( {1 - A} \right)}} = \frac{5}{2}\frac{{{{{k}}_{\rm{B}}}}}{m}\frac{{1 - B}}{{1 - A}}\mu $$ (41) $$ Pr = \frac{{1 - A}}{{1 - B}} $$ (42) 可知对于沙克霍夫模型[12],
$ A = 1 - Pr $ ,$ B = 0 $ ,${\nu ^{\rm{S}}} = p/\mu $ ; 对于贝利模型[15],$ A = 0 $ ,$ B = \left( {Pr - 1} \right)/Pr $ ,${\nu ^{\rm{VVB}}} = Pr \cdot p/\mu $ , 这与前述一致. 通过选择合适的A, B可以得到正确的普朗特数Pr.2.3 高阶矩分析
为了进一步验证新模型方程在数学上与玻尔兹曼方程的相容性, 这里从麦克斯韦分子满足的玻尔兹曼方程高阶矩进行分析. 由于麦克斯韦分子的特殊性, 即其碰撞积分与相对速度无关, 可以得到对应的高阶矩.
定义碰撞积分项的矩为
$$ \Delta Q = \iiint_{{{\mathbb{R}}^3}} {Q{{\left( {\frac{{\partial f}}{{\partial t}}} \right)}_{\rm{col}}}{\rm{d}} W} $$ (43) 这里
$Q$ 为分子速度的某种组合. 当$Q = {V_i}{V_j}$ 和${V_i}{V^2}$ 时有[7]$$ \Delta {{Q_1}}= \Delta \left( {{V_i}{V_j}} \right)= \frac{p}{\mu }\frac{{{\tau _{ij}}}}{m} $$ (44) $$ \Delta{{Q_2}} = \Delta \left( {{V_i}{V^2}} \right) = \frac{{2p}}{\mu }\left( {\sum\limits_j {\frac{{{U_j}{\tau _{ij}}}}{m}} - \frac{{Pr{q_i}}}{m}} \right) $$ (45) 对于新模型方程,
$Q$ 代入碰撞积分矩, 有$$ \Delta {{Q_1}} = {\nu ^{\rm{N}}}\left( {1 - B} \right)\frac{{{\tau _{ij}}}}{m} $$ (46) $$ \Delta {{Q_2}} = {\nu ^{\rm{N}}}\left[ {2\left( {A - 1} \right)\frac{{{q_i}}}{m} + 2\left( {1 - B} \right)\sum\limits_j {\frac{{{U_j}{\tau _{ij}}}}{m}} } \right] $$ (47) 对比可知
$$ {\nu ^{\rm{N}}} = \frac{{Pr}}{{1 - A}}\frac{p}{\mu } = \frac{1}{{1 - B}}\frac{p}{\mu } = {\nu _{\rm{coe}}}\frac{p}{\mu } $$ (48) 该结论与上一节一致, 这进一步证明了新模型方程与玻尔兹曼方程的相容性, 与沙克霍夫模型[12]和贝利模型[15]的一致性.
2.4 H定理的证明
定义[39]
$$ H\left( t \right) = \iiint_{{{\mathbb{R}}^3}} {f\ln f{\rm{d}} W} $$ (49) 为简单起见, 这里仅考虑
$$ \frac{{\partial f}}{{\partial t}} = {\nu ^{\rm{N}}}\left( {f_{\rm{mix}}^{\rm{N}} - f} \right) $$ (50) 的情况. 则
$$ \begin{split} & \frac{\partial }{{\partial t}}H\left( t \right) = \iiint_{{{\mathbb{R}}^3}} {\frac{{\partial f}}{{\partial t}}\ln f{\rm{d}} W} + \iiint_{{{\mathbb{R}}^3}} {f\frac{{\partial \ln f}}{{\partial t}}{\rm{d}} W} = \\&\qquad \iiint_{{{\mathbb{R}}^3}} {{\nu ^{\rm{N}} }\left( {f_{\rm{mix}}^{\rm{N}} - f} \right)\ln f{\rm{d}} W} \end{split} $$ (51) 即
$$ \begin{split} & \frac{\partial }{{\partial t}}H\left( t \right) = {\nu ^{\rm{N}} }\iiint_{{{\mathbb{R}}^3}} {\left( {f_{\rm{mix}}^{\rm{N}} - f} \right)\ln \frac{f}{{f_{\rm{mix}}^{\rm{N}}}}{\rm{d}} W} +\\&\qquad {\nu ^{\rm{N}}}\iiint_{{{\mathbb{R}}^3}} {\left( {f_{\rm{mix}}^{\rm{N}} - f} \right) \left[ {\ln {f_{\rm{M}} } + \ln \left({1 + \varPhi _{\rm{mix}}^{\rm{N}}} \right)} \right]{\rm{d}} W} \end{split} $$ (52) 其中
$\varPhi _{\rm{mix}}^{\rm{N}} = \left( {1 - \dfrac{{Pr}}{{{\nu _{\rm{coe}}}}}} \right){\varPhi _{\rm{A}} }\left( {\boldsymbol{C}} \right) + \left( {1 - \dfrac{1}{{{\nu _{\rm{coe}}}}}} \right){\varPhi _{\rm{B}} }\left( {\boldsymbol{C}} \right)$ .由于当
$ - 1 < x \leqslant 1$ 时有$$ \ln \left( {1 + x} \right) = x - \frac{{{x^2}}}{2} + \frac{{{x^3}}}{3} - \cdots + {\left( { - 1} \right)^{l + 1}}\frac{{{x^l}}}{l} $$ (53) 且从查普曼−恩斯科分析的过程可知分布函数
$f$ 相对${f_{\rm{M}} }$ 的扰动为一小量, 因此有$$ \begin{split} & \iiint_{{{\mathbb{R}}^3}} {f_{\rm{mix}}^{\rm{N}}\ln \left( {1 + \varPhi _{\rm{mix}}^{\rm{N}}} \right){\rm{d}} W} = \iiint_{{{\mathbb{R}}^3}} {f_{\rm{mix}}^{\rm{N}}\varPhi _{\rm{mix}}^{\rm{N}}{\rm{d}} W} = \\&\qquad \frac{{2{{\left( {1 - Pr/{\nu _{\rm{coe}}}} \right)}^2}{q^2}}}{{5p{{{k}}_{\rm{B}}}T \cdot {{{k}}_{\rm{B}}}T/m}} + \frac{{{{\left( {1 - 1/{\nu _{\rm{coe}}}} \right)}^2}\displaystyle\sum\limits_{i,j} {\tau _{ij}^2} }}{{2p{{{k}}_{\rm{B}}}T}} \end{split} $$ (54) 和
$$ \begin{split} & \iiint_{{{\mathbb{R}}^3}} {f\ln \left( {1 + \varPhi _{\rm{mix}}^{\rm{N}}} \right){\rm{d}} W} = \iiint_{{{\mathbb{R}}^3}} {f\varPhi _{\rm{mix}}^{\rm{N}}{\rm{d}} W} = \\&\qquad \left( {1 - \frac{{Pr}}{{{\nu _{\rm{coe}}}}}} \right)\frac{{m{q^2}}}{{p{{{k}}_{\rm{B}}}T}}\frac{2}{{5{{{k}}_{\rm{B}}}T}} - \left( {1 - \frac{1}{{{\nu _{\rm{coe}}}}}} \right)\frac{{\displaystyle\sum\limits_{i,j} {\tau _{ij}^2} }}{{2p{{{k}}_{\rm{B}}}T}} \end{split} $$ (55) 合并后, 有
$$ \begin{split} & \frac{\partial }{{\partial t}}H\left( t \right) = {\nu ^{\rm{N}}}\iiint_{{{\mathbb{R}}^3}} {\left( {f_{\rm{mix}}^{\rm{N}} - f} \right)\ln \frac{f}{{f_{\rm{mix}}^{\rm{N}}}}{\rm{d}} W} + \\&\qquad \frac{{2m{q^2}{\nu ^{\rm{N}}}}}{{5p{{{k}}_{\rm{B}}}T \cdot {{{k}}_{\rm{B}}}T}}\left[ {{{\left( {1 - \frac{{Pr}}{{{\nu _{\rm{coe}}}}}} \right)}^2} - \left( {1 - \frac{{Pr}}{{{\nu _{\rm{coe}}}}}} \right)} \right] +\\&\qquad \frac{{{\nu ^{\rm{N}} }}}{{2p{{{k}}_{\rm{B}}}T}}\left[ {{{\left( {1 - \frac{1}{{{\nu _{\rm{coe}}}}}} \right)}^2} + \left( {1 - \frac{1}{{{\nu _{\rm{coe}}}}}} \right)} \right]\sum\limits_{i,j} {\tau _{ij}^2} \end{split} $$ (56) 由于
$\left( {f_{\rm{mix}}^{\rm{N}} - f} \right)$ 与$\ln \left( {f/f_{\rm{mix}}^{\rm{N}}} \right)$ 反号, 有$$ {\nu ^{\rm{N}}}\iiint_{{{\mathbb{R}}^3}} {\left( {f_{\rm{mix}}^{\rm{N}} - f} \right)\ln \frac{f}{{f_{\rm{mix}}^{\rm{N}}}}{\rm{d}} W} \leqslant 0 $$ (57) 当
$f_{\rm{mix}}^{\rm{N}} = f$ 时等号成立. 对于式(56)后两项, 当$$ {\left( {1 - \frac{{Pr}}{{{\nu _{\rm{coe}}}}}} \right)^2} - \left( {1 - \frac{{Pr}}{{{\nu _{\rm{coe}}}}}} \right) \leqslant 0 $$ (58) 且
$$ {\left( {1 - \frac{1}{{{\nu _{\rm{coe}}}}}} \right)^2} + \left( {1 - \frac{1}{{{\nu _{\rm{coe}}}}}} \right) \leqslant 0 $$ (59) 时, 有
$\dfrac{\partial }{{\partial t}}H\left( t \right) \leqslant 0$ , 得到$$ \max \left( {Pr,1/2} \right) \leqslant {\nu _{\rm{coe}}} \leqslant 1 $$ (60) 即当上述不等式成立, 且Kn较小时, 新模型方程满足H定理. 从上面数学分析过程来看,
${\nu _{\rm{coe}}}$ 在其取值范围内可以任意选取, 但在实际物理过程中即分布函数的碰撞松弛中${\nu _{\rm{coe}}}$ 需要具有一定的物理意义.2.5 碰撞松弛参数的确定
在气体动理学理论中, 根据碰撞间隔理论[39]可以推导出碰撞时间间隔
$\eta $ 、黏性系数$\mu $ 、压力$p$ 之间的关系为$$ \mu = p\eta $$ (61) 显然, 若将
$\eta $ 直接视为玻尔兹曼模型方程中气体分子碰撞松弛参数$\nu $ 的倒数, 则上式与沙克霍夫模型一致, 对应于${\nu _{\rm{coe}}} $ = 1. 但上述碰撞间隔没有考虑到分子碰撞后的速度残留[39], 因此需要对其进行修正, 通常将其取为扰动驰豫时间, 即$$ \eta = \frac{{4\lambda }}{{{{\text{π}} }\bar C}}{, }\;\;\;\bar C = \sqrt {\frac{{8{{{k}}_{\rm{B}}}T}}{{{{\text{π}} }m}}} $$ (62) 表征为分布函数松弛至其初始值的
$1/{{\rm{e}}} $ 所历经的时间(${{\rm{e}}} $ 为自然对数的底), 这里$\lambda $ 为平均自由程.此外在直接模拟蒙特卡洛方法中, 对于变径硬球(VHS)模型有[41]
$$ \lambda = \frac{{2\left( {5 - 2\omega } \right)\left( {7 - 2\omega } \right)\mu }}{{15\sqrt 2 {{\text{π}} }n{{\left( {m{{{k}}_{\rm{B}}}T/{{\text{π}} }} \right)}^{1/2}}}} $$ (63) 其中,
$\omega $ 为分子的变径硬球模型参数, 则$$ \eta = \frac{{4\lambda }}{{{{\text{π}} }\bar C}} = \frac{{2\left( {5 - 2\omega } \right)\left( {7 - 2\omega } \right)\mu }}{{15{{\text{π}} }p}} $$ (64) 由此可得
$$ {\nu _{\rm{coe}}} = \frac{{15{{\text{π}} }}}{{2\left( {5 - 2\omega } \right)\left( {7 - 2\omega } \right)}} $$ (65) 上述
${\nu _{\rm{coe}}}$ 是通过微观分子碰撞间隔理论得到的, 在气体动理学介观理论中, 为了表征分布函数$f$ 松弛至$f_{\rm{mix}}^{\rm{N}}$ 的$1/{{\rm{e}}} $ 的物理意义, 同时使${\nu _{\rm{coe}}}$ 满足式(60), 本文针对新模型方程取$$ {\nu _{\rm{coe}}} = \frac{{15{{\text{π}} }}}{{2\left( {5 - 2\omega } \right)\left( {7 - 2\omega } \right)}}\left( {1 - \frac{1}{{{\rm{e}}} }} \right) $$ (66) 对于单原子气体氩气(Ar),
$\omega $ = 0.81, 此时${\nu _{\rm{coe}}} $ = 0.819; 对于空气,$\omega $ = 0.77,${\nu _{\rm{coe}}} $ = 0.788 4.3. 数值分析
3.1 一维激波结构模拟分析
为验证本文新建模型的有效性, 首先模拟了一维激波结构. 设置激波马赫数为
$M{a_{\rm{s}} } $ = 8, 上游条件为: 数密度${n_1} $ = 1 × 1023 m−3, 温度${T_1} $ = 200 K, 速度${V_1} $ = 2105.738 m/s; 下游条件:${n_2} $ = 3.821 × 1020 m−3,${T_2} $ = 4174.414 K,${V_2} $ = 551.111 m/s; 气体为氩气; 计算中激波上下游取当地平衡分布, 空间网格间距取上游气体分子平均自由程的0.05倍.图1给出了采用沙克霍夫模型、本文新建模型、贝利模型模拟的激波轮廓线与直接模拟蒙特卡洛方法结果的对比, 分别用“Shakhov” “new model” “VVB”和“DSMC”表示, 其中直接模拟蒙特卡洛方法结果由贝德(Bird)的一维可视化直接模拟(DS1V)代码[41]计算得到,
${\lambda _1}$ 表示激波上游气体子平均自由程, 并将密度轮廓$\left( {\rho - {\rho _1}} \right)/\left( {{\rho _2} - {\rho _1}} \right) $ = 0.5置于$x/{\lambda _1} $ = 0处. 图2给出了3种模型模拟的一维激波内部无量纲热流和切应力变化曲线. 从图1可以看出本文新建模型能较好地反映$M{a_{\rm{s}} } $ = 8的激波轮廓, 验证了其有效性; 对于几种模型而言, 沙克霍夫模型与直接模拟蒙特卡洛方法结果符合最好, 本文模型次之; 差别主要在激波上游靠近$x/{\lambda _1} $ = 0区域, 从图2可以看出该区域无量纲切应力${\tau _{xx}}$ 变化很小, 而无量纲热流${q_x}$ 变化较大、且各模型结果区别也较大, 仅考虑了热流项的沙克霍夫模型与直接模拟蒙特卡洛方法结果最为接近, 说明在一维激波内部宏观切应力对内部参数分布影响很小, 这可能是由于一维问题的宏观切应力仅出现在x方向且无量纲值相对热流较小, 黏性效应不显著. 此外模型方程模拟的温度轮廓在激波上游均出现了提前抬起的现象, 这可能是碰撞频率与分子速度无关造成的.![]()
图 1 不同模型模拟的激波轮廓与直接模拟蒙特卡洛方法对比Figure 1. Comparison of shock profiles simulated by GKUA with different models and DSMC![]()
图 2 不同模型模拟的无量纲热流和切应力对比Figure 2. Comparison of non-dimensional heat fluxes and shear stresses simulated by GKUA with different models3.2 平板绕流模拟分析
为研究新模型方程模拟飞行器再入时在近空间环境黏性效应急剧增加的非平衡流场中黏性应力和热流影响较大的区域流动问题, 首先考察平板绕流. 通常对于这类问题, 在克努森层内越靠近平板前缘, 温度和速度梯度越大, 非平衡效应越显著, 对应的热流和黏性效应影响也越发严重.
来流条件: 马赫数
$M{a_\infty } $ = 5, 来流气体为氩气, 来流温度和物面温度均为200 K, 数密度为${n_\infty } $ = 1.175 × 1020 m−3. 平板长度为1 m、宽0.04 m. 计算中, 为了揭示高马赫数附面层内流动与传热机制, 直接模拟蒙特卡洛方法的背景网格间距取3.62 mm, 统一算法贴近物面第一排网格间距为2.5 mm, 均远小于来流分子平均自由程10 mm. 统一算法模拟中除了所用新、旧碰撞模型外其他计算条件完全一样, 物面采用与直接模拟蒙特卡洛方法一致的无穿透麦克斯韦散射分布, 来流边界为来流无量纲麦克斯韦平衡分布, 出口边界根据当地分子运动方向按特征边界条件处理. 图3给出了本文采用新模型方程模拟的温度分布云图(用GKUA_NewModel表示的下半流场)与二维可视化直接模拟(DS2V)软件[42]模拟结果(用DSMC表示的上半流场)的对比情况. 从图中可以看出, 两个结果几乎完全一致, 而直接模拟蒙特卡洛方法尤其在尾部低速流场存在明显的统计波动. 这说明新模型模拟的流场结果是可信的.![]()
图 3 本文模拟的温度云图与直接模拟蒙特卡洛方法结果对比Figure 3. Comparison of temperature contours simulated by GKUA with new model and DSMC图4给出了在
$x $ = 0.1 m截面上压力和温度沿y方向不同碰撞模型的统一算法模拟结果与直接模拟蒙特卡洛方法的对比, 其中本文建立的新模型用“new model”表示, 沙克霍夫模型[12]用“Shakhov”表示, 贝利模型[15]用“VVB”表示. 从图中可以看出, 整体上使用新模型的统一算法结果与直接模拟蒙特卡洛方法模拟值相容, 局部附面层区域$y$ < 0.15 m时, 相比沙克霍夫模型, 新模型统一算法结果与直接模拟蒙特卡洛方法模拟值较为接近; 而且还可看出对此黏性应力影响极为严重的附面层流动, 使用贝利模型的统一算法结果与直接模拟蒙特卡洛方法模拟值符合更好; 在0.15 m$ < y < $ 0.225 m时新模型的统一算法结果与直接模拟蒙特卡洛方法模拟几乎完全重合, 而沙克霍夫模型和贝利模型的统一算法计算值与直接模拟蒙特卡洛方法结果存在一定偏差. 从流场结果来看, 对于平板流动在靠近头部的区域, 稀薄效应显著、非平衡现象突出, 尤其是靠近物面的克努森层内, 因黏性效应导致剪切层内应力张量的影响远大于热流的影响, 故而使用贝利模型的统一算法结果能更好地与直接模拟蒙特卡洛方法模拟值吻合; 在远离附面层的流场内部, 应力张量与热流影响相当, 均会对当地气体分子速度分布函数产生影响, 故而同时考虑了热流与应力张量的新模型统一算法计算值能与直接模拟蒙特卡洛方法几乎重合; 在有限的计算流场外缘($y$ > 0.225 m)流动处于近平衡态, 直接模拟蒙特卡洛方法依靠麦克斯韦平衡态分布随机采样模拟分子的运动与碰撞, 而统一算法基于求解当地气体分子速度分布函数的演化更新实时捕捉宏观流动量变化分布, 使得不同模型得到的统一算法计算值较为一致, 均与直接模拟蒙特卡洛方法模拟存在一定的偏差, 反映出进一步发展求解玻尔兹曼模型速度分布函数方程的气体动理论统一算法有其独到新颖性与应用价值.![]()
图 4$x $ = 0.1 m截面上压力和温度沿y方向不同碰撞模型的统一算法计算值与直接模拟蒙特卡洛方法模拟结果的对比Figure 4. Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of$x $ = 0.1 m图5和图6分别给出了在
$x $ = 0.5 m和0.8 m截面上压力和温度沿y方向不同碰撞模型的统一算法计算结果与直接模拟蒙特卡洛方法模拟值对比. 随着x增大剪切层范围也扩大, 在这些区域内流体黏性的影响较大, 因而考虑了应力张量的碰撞模型, 如使用3种碰撞模型得到的统一算法计算分布变化趋势彼此一致, 尤其是贝利模型与本文新提出的模型方程统一算法计算分布近乎完全重合, 能较好模拟这些区域的局部流动, 整体上都与直接模拟蒙特卡洛方法模拟结果吻合. 而当y增大到剪切层外缘时, 热流的影响开始显现加剧的某些区域, 使用沙克霍夫模型得到的统一算法结果更吻合直接模拟蒙特卡洛方法模拟值. 同样在当y进一步增大到远离附面层临近计算区域外缘, 流动梯度仍然较大的近平衡态局部流场, 几种模型得到的统一算法计算值均与直接模拟蒙特卡洛方法模拟值存在一定的系统误差, 究其实质极可能缘于直接模拟蒙特卡洛方法完全使用麦克斯韦平衡态分布随机统计模拟且存在统计涨落, 导致了两种方法体制的完全不一样, 对于流场梯度较大的局部流动模拟, 彼此结果间存在一定的计算偏差.![]()
图 5$x $ = 0.5 m截面上压力和温度沿y方向不同碰撞模型统一算法计算结果与直接模拟蒙特卡洛方法模拟值比较Figure 5. Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of$x $ = 0.5 m![]()
图 6$x $ = 0.8 m截面上压力和温度沿y方向不同碰撞模型统一算法模拟结果与直接模拟蒙特卡洛方法的对比Figure 6. Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of$x $ = 0.8 m从图4 ~ 图6可以看出, 由于本文建立的玻尔兹曼方程新模型同时兼顾了流场中应力和热流对气体分子速度分布函数的影响, 整体而言能较好地模拟全域流动. 从统一算法计算
$x $ = 0.1 m, 0.5 m和0.8 m压力和温度变化来看, 物面压力从1.5 Pa降至0.9 Pa和0.8 Pa, 温度从380 K降至320 K和300 K, 相对于壁温200 K, 物面温度跳跃逐渐降低, 说明高马赫数绕流平板前缘强压缩流动在近空间飞行环境稀薄效应影响下逐渐变得稀疏, 沿平板物面压力、温度不断减小, 过渡到平板尾部低压常规流.图7给出了驻点线上不同碰撞模型得到的统一算法计算温度分布与直接模拟蒙特卡洛方法模拟值对比. 平板头部压缩激波层是非平衡效应最为严重的区域, 温度梯度和速度梯度均较大. 从图中可以看出, 在近空间飞行环境激波内部, 本文建立的玻尔兹曼方程新模型统一算法计算结果始终与直接模拟蒙特卡洛方法模拟值符合较好. 这进一步说明了需要同时考虑流场中黏性应力和热流对气体分子速度分布函数演化更新的影响. 同时也看出基于玻尔兹曼方程可计算建模统一算法计算结果在驻点激波前部流动趋于近平衡态时温度高于直接模拟蒙特卡洛方法模拟值, 这与一维激波结构中温度轮廓的变化规律类似.
![]()
图 7 驻点线上不同碰撞模型统一算法计算的温度分布与直接模拟蒙特卡洛方法模拟值比较Figure 7. Comparison of temperature simulated by GKUA with different collision models and DSMC along the stagnation line3.3 双圆柱干扰流动模拟分析
为验证新建模型统一算法对激波干扰的模拟能力, 这里考察双圆柱干扰流动. 来流条件和计算条件及边界条件与上节一致, 圆柱直径为1 m, 两圆柱圆心间距为4 m. 在该条件下上下两个圆柱的头部脱体激波会相互作用, 形成复杂的激波干扰流动.
图8给出了分别使用3种模型方程的统一算法计算温度分布云图与直接模拟蒙特卡洛方法模拟值的对比. 两圆柱脱体激波相互作用后在
$x = - $ 0.19 m的位置形成一个正激波类马赫杆, 跨过马赫杆后温度和压力增加、速度降低. 而马赫杆两端形成的斜激波仅对圆柱尾流产生影响, 将导致单边圆柱法向力很小. 从云图可以看出, 使用贝利模型统一算法得到核心区温度分布与直接模拟蒙特卡洛方法模拟值最为接近, 但马赫杆位置稍向后, 而沙克霍夫模型差别最大, 对马赫杆位置而言本文新建模型统一算法与直接模拟蒙特卡洛方法符合最好. 表1给出了3种碰撞模型统一算法及直接模拟蒙特卡洛方法模拟得到的单边圆柱轴向力系数及法向力系数, 以及不同模型所得轴向力系数与直接模拟蒙特卡洛方法模拟结果的偏差. 从表中可以看出, 本文新建立的碰撞模型统一算法计算结果与直接模拟蒙特卡洛方法模拟值最大偏差在4.6%以内, 能较好地模拟气动力特性, 精度优于常用的沙克霍夫模型.![]()
图 8 3种模型方程计算得到的温度分布云图与直接模拟蒙特卡洛方法结果的对比Figure 8. Comparison of temperature contours simulated by GKUA with three collision models and DSMC表 1 不同模型/方法所得轴向力系数及其偏差与法向力系数Table 1. The axial and normal force coefficients by different models and DSMC method and errors between axial force coefficientsModel/Method Ca Error/% Cn DSMC 1.287 0 — 3.82 × 10−4 New model 1.3458 4.57 5.76 × 10−3 Shakhov model 1.3760 6.92 8.25 × 10−3 Belyi model 1.3261 3.04 3.90 × 10−3 图9给出了对称面(即
$y $ = 0)上马赫杆附近沿x方向不同碰撞模型统一算法计算得到的压力、温度和速度分布与直接模拟蒙特卡洛方法模拟值的对比情况. 本文新建碰撞模型较好地求解了二维绕流的激波位置及激波前后宏观参数的变化趋势, 明显优于沙克霍夫模型, 在激波内部模拟结果优于贝利模型. 这是因为激波前后温度和速度梯度较大, 但$y $ = 0截面上仅存在x方向热流${q_x}$ , 而x和y方向的切应力${\tau _{xx}}$ 和${\tau _{yy}}$ 同时存在, 速度梯度对分布函数的影响更大, 使得本文新建碰撞模型能更好地模拟激波位置和内部参数变化. 这与模拟一维激波结构所得结论不一致, 可能是由于随着维度的增加多方向宏观切应力导致黏性效应更加显著, 其对分子速度分布函数的影响增大, 与热流共同作用影响激波位置及其内部参数的分布.![]()
图 9$y $ = 0 m截面上压力、温度和速度沿$ x $ 方向不同碰撞模型模拟结果与直接模拟蒙特卡洛方法的对比Figure 9. Comparison of pressure, temperature and velocity simulated by GKUA with different collision models and DSMC along the direction of x on the section of$y $ = 0 m4. 结论
本文通过分析玻尔兹曼方程的一阶查普曼−恩斯科分析解, 结合沙克霍夫模型和贝利模型构造过程及数学推导与物理分析, 建立了一种同时考虑热流矢量和应力张量影响、满足玻尔兹曼方程高阶碰撞矩的综合可计算碰撞模型, 并在数学上分析了其守恒性、H定理等基本属性, 证明了新模型方程与玻尔兹曼方程的相容性, 给出了新模型与现有模型的关系, 通过碰撞动力学确定了碰撞松弛参数统一表达式. 经过模拟一维激波结构、二维近空间飞行环境平板和多体圆柱干扰流动, 并与直接模拟蒙特卡洛方法及不同碰撞模型方程的统一算法结果对比分析, 得出以下结论:
(1)在稀薄效应显著、非平衡现象突出的流动区域, 热流矢量和应力张量对当地气体分子速度分布函数的影响较大, 在构造碰撞松弛模型时需要同时考虑两者的影响;
(2)在激波捕捉方面, 流场中黏性效应显著的区域沙克霍夫模型结果与直接模拟蒙特卡洛方法模拟值差别较大, 本文新建模型与贝利模型结果接近, 在激波内部本文模型与直接模拟蒙特卡洛方法符合更好, 而在黏性效应较小时仅考虑热流时结果与直接模拟蒙特卡洛方法符合更好, 进一步说需要考虑多参数对分子速度分布函数的综合影响;
(3)本文新建模型的碰撞松弛参数中引入了与自然对数的底
$ {{\rm{e}}} $ 相关的常数作为自由参数, 间接表征了分布函数松弛的物理意义, 虽然能满足模型方程的H定理, 但在流动的某些区域仍然与直接模拟蒙特卡洛方法存在差距. 若能构造一个与当地非平衡参数相关联的碰撞松弛参数, 将气体分子速度分布函数的松弛过程与当地流动参数的梯度耦合, 应该能更好地模拟再入跨流域非平衡流动过程, 这也是本文下一步拟开展的工作. -
![]()
图 1 不同模型模拟的激波轮廓与直接模拟蒙特卡洛方法对比
Figure 1. Comparison of shock profiles simulated by GKUA with different models and DSMC
![]()
图 2 不同模型模拟的无量纲热流和切应力对比
Figure 2. Comparison of non-dimensional heat fluxes and shear stresses simulated by GKUA with different models
![]()
图 3 本文模拟的温度云图与直接模拟蒙特卡洛方法结果对比
Figure 3. Comparison of temperature contours simulated by GKUA with new model and DSMC
![]()
图 4
$x $ = 0.1 m截面上压力和温度沿y方向不同碰撞模型的统一算法计算值与直接模拟蒙特卡洛方法模拟结果的对比Figure 4. Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of
$x $ = 0.1 m![]()
图 5
$x $ = 0.5 m截面上压力和温度沿y方向不同碰撞模型统一算法计算结果与直接模拟蒙特卡洛方法模拟值比较Figure 5. Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of
$x $ = 0.5 m![]()
图 6
$x $ = 0.8 m截面上压力和温度沿y方向不同碰撞模型统一算法模拟结果与直接模拟蒙特卡洛方法的对比Figure 6. Comparison of pressure and temperature simulated by GKUA with different collision models and DSMC along the direction of y on the section of
$x $ = 0.8 m![]()
图 7 驻点线上不同碰撞模型统一算法计算的温度分布与直接模拟蒙特卡洛方法模拟值比较
Figure 7. Comparison of temperature simulated by GKUA with different collision models and DSMC along the stagnation line
![]()
图 8 3种模型方程计算得到的温度分布云图与直接模拟蒙特卡洛方法结果的对比
Figure 8. Comparison of temperature contours simulated by GKUA with three collision models and DSMC
![]()
图 9
$y $ = 0 m截面上压力、温度和速度沿$ x $ 方向不同碰撞模型模拟结果与直接模拟蒙特卡洛方法的对比Figure 9. Comparison of pressure, temperature and velocity simulated by GKUA with different collision models and DSMC along the direction of x on the section of
$y $ = 0 m表 1 不同模型/方法所得轴向力系数及其偏差与法向力系数
Table 1 The axial and normal force coefficients by different models and DSMC method and errors between axial force coefficients
Model/Method Ca Error/% Cn DSMC 1.287 0 — 3.82 × 10−4 New model 1.3458 4.57 5.76 × 10−3 Shakhov model 1.3760 6.92 8.25 × 10−3 Belyi model 1.3261 3.04 3.90 × 10−3 
下载: 导出CSV
-
[1] Schouler M, Prevereaud Y, Mieussens L. Survey of flight and numerical data of hypersonic rarefied flows encountered in earth orbit and atmospheric reentry. Progress in Aerospace Sciences, 2020, 118: 100638 doi: 10.1016/j.paerosci.2020.100638
[2] Xu SS, Wu ZN, Li Q, et al. Hybrid continuum /DSMC computation of rocket mode lightcraft flow in near space with high temperature and rarefaction effect. Computers and Fluids, 2009, 38: 1394-1404 doi: 10.1016/j.compfluid.2008.01.024
[3] 周恒, 张涵信. 空气动力学的新问题. 中国科学: 物理学 力学 天文学, 2015(45): 104709 (Zhou Heng, Zhang Hanxin. New problems of aerodynamics. Scientia Sinica:Physica,Mechanica &Astronomica, 2015(45): 104709 (in Chinese) [4] 方方, 田园, 赵攀等. 空间返回航天器气动外形设计与需求分析. 空气动力学学报, 2018, 36(5): 816-825 (Fang Fang, Tian Yuan, Zhao Pan, et al. Aerodynamics shape designs and requirement analysis of re-entry spacecraft. Acta Aerodynamics Sinica, 2018, 36(5): 816-825 (in Chinese) doi: 10.7638/kqdlxxb-2017.0027 [5] Struchtrup H. Macroscopic Transport Equations for Rarefied Gas Flows. Berlin: Springer, 2005
[6] Cercignani C. Rarefied Gas Dynamics: from Basic Concepts to Actual Calculations. Cambridge: Cambridge University Press, 2000
[7] 沈青. 稀薄气体动力学. 北京: 国防工业出版社, 2003 (Shen Qing. Rarefied Gas Dynamics. Beijing: National Defense Industry Press, 2003 (in Chinese))
[8] Bhatnagar PL, Gross EP, Krook M. A model collision processes in gases: I. Small amplitude processes is charged and neutral one-component system. Physical Review, 1954, 94(3): 511-525 doi: 10.1103/PhysRev.94.511
[9] Vincenti WG, Kruger CH. Introduction to Physical Gas Dynamics. New York: Wiley, 1965
[10] Holway LH. New statistical models for kinetic theory: methods of construction. Physics of Fluids, 1966, 9(9): 1658-1673 doi: 10.1063/1.1761920
[11] Brull S, Schneider J. On the ellipsoidal statistical model for polyatomic gases. Continuum Mechanics and Thermodynamics, 2009, 20: 489-508 doi: 10.1007/s00161-009-0095-3
[12] Shakhov EM. Generalization of the Krook kinetic relaxation equation. Fluid Dynamics, 1968, 3(5): 95-96
[13] Shakhov EM. Kinetic model equations and numerical results//14th International Symposium on Rarefied Gas Dynamics, Tokyo, 1984: 137-148
[14] Liu GJ. A method for constructing a model form for the Boltzmann equation. Physics of Fluids A, 1990, 2(2): 277-280 doi: 10.1063/1.857777
[15] Belyi VV. Derivation of model kinetic equation. Europhysics Letters, 2015, 111: 40011 doi: 10.1209/0295-5075/111/40011
[16] Zheng YS, Struchtrup H. Ellipsoidal statistical BGK model with velocity-dependent collision frequency. Physics of Fluids, 2015, 17: 127103
[17] Christos T, Gian PG, Dimitris V, et al. Effect of vibrational degrees of freedom on the heat transfer in polyatomic gases confined between parallel plates. International Journal of Heat and Mass Transfer, 2016, 102: 162-173 doi: 10.1016/j.ijheatmasstransfer.2016.06.010
[18] Brull S, Schneider J. A new approach for the ellipsoidal statistical model. Continuum Mechanics and Thermodynamics, 2008, 20: 63-74 doi: 10.1007/s00161-008-0068-y
[19] Belyi VV. On the model kinetic description of plasma and a Boltzmann gas of hard spheres. Journal of Statistical Mechanics: Theory and Experiment, 2009, 6: 06001
[20] Li ZH, Fang M, Jiang XY, et al. Convergence proof of the DSMC method and the gas kinetic unified algorithm for the Boltzmann equation. Science China, Physics, Mechanics & Astronomy, 2013, 56(2): 404-417
[21] Frolova A, Titarev V. Recent progress on supercomputer modelling of high-speed rarefied gas flows using kinetic equations. Supercomputing Frontiers and Innovations, 2018, 5(3): 117-121
[22] Titarev V. Application of model kinetic equations to hypersonic rarefied gas flows. Computers and Fluids, 2018, 169: 62-70 doi: 10.1016/j.compfluid.2017.06.019
[23] Titarev V, Frolova AA, Rykov VA. Comparison of the Shakhov kinetic equation and DSMC method as applied to space vehicle aerothermodynamics. Journal of Computational and Applied Mathematics, 2020, 364: 112354 doi: 10.1016/j.cam.2019.112354
[24] Li ZH, Zhang HX. Study on gas kinetic unified algorithm for flows from rarefied transition to continuum. Journal of Computational Physics, 2004, 193: 708-738 doi: 10.1016/j.jcp.2003.08.022
[25] Li ZH, Peng AP, Zhang HX, et al. Rarefied gas flow simulations using high-order gas kinetic unified algorithms for Boltzmann model equations. Progress in Aerospace Sciences, 2015, 74: 81-113 doi: 10.1016/j.paerosci.2014.12.002
[26] Peng AP, Li ZH, Wu JL, et al. Implicit gas kinetic unified algorithm based on multi-block docking grid for multi-body reentry flows covering all flow regimes. Journal of Computational Physics, 2016, 327: 919-942 doi: 10.1016/j.jcp.2016.09.050
[27] Li ZH, Hu WQ, Peng AP, et al. Gas kinetic unified algorithm for plane external force-driven flows covering all flow regimes by modeling of Boltzmann equation. International Journal for Numerical Methods in Fluids, 2020, 92: 922-949 doi: 10.1002/fld.4812
[28] 李志辉, 蒋新宇, 吴俊林等. 转动非平衡玻尔兹曼模型方程统一算法与全流域绕流计算应用. 力学学报, 2014, 46(3): 336-351 (Li Zhihui, Jiang Xinyu, Wu Junlin, et al. Gas kinetic unified algorithm for Boltzmann model equation in rotational non-equilibrium and its application to the whole range flow regimes. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(3): 336-351 (in Chinese) doi: 10.6052/0459-1879-13-246 [29] 彭傲平, 李志辉, 吴俊林等. 基于多块对接网格的隐式气体运动论统一算法应用研究. 力学学报, 2016, 48(1): 95-101 (Peng Aoping, Li Zhihui, Wu Junlin, et al. An application of implicit gas kinetic unified algorithm based on multi block patched grid. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(1): 95-101 (in Chinese) doi: 10.6052/0459-1879-14-279 [30] 彭傲平, 李志辉, 吴俊林等. 含振动能激发Boltzmann模型方程气体动理论统一算法验证与分析. 物理学报, 2017, 66(20): 204703 (Peng Aoping, Li Zhihui, Wu Junlin, et al. Validation and analysis of gas kinetic unified algorithm for solving Boltzmann model equation with vibrational energy excitation. Acta Physica Sinica, 2017, 66(20): 204703 (in Chinese) doi: 10.7498/aps.66.204703 [31] 李志辉, 梁杰, 李中华等. 跨流域空气动力学模拟方法与返回舱再入气动研究. 空气动力学学报, 2018, 36(5): 826-847 (Li Zhihui, Liang Jie, Li Zhonghua, et al. Simulation methods of aerodynamics covering various flow regimes and applications to aerodynamic characteristics of re-entry spacecraft module. Acta Aerodynamics Sinica, 2018, 36(5): 826-847 (in Chinese) doi: 10.7638/kqdlxxb-2018.0121 [32] Xu K, Huang JC. A unified gas kinetic scheme for continuum and rarefied flows. Journal of Computational Physics, 2010, 229: 7747-7764 doi: 10.1016/j.jcp.2010.06.032
[33] Chen SZ, Xu K, Lee CB, et al. A unified gas kinetic scheme with moving mesh and velocity space adaptation. Journal of Computational Physics, 2012, 231: 6643-6664 doi: 10.1016/j.jcp.2012.05.019
[34] Huang JC, Xu K, Yu PB. A unified gas kinetic scheme for continuum and rarefied flows II: multi-dimensional cases. Communications in Computational Physics, 2012, 12(3): 662-690 doi: 10.4208/cicp.030511.220911a
[35] Huang JC, Xu K, Yu P. A unified gas kinetic scheme for continuum and rarefied flows III: microflow simulations. Communications in Computational Physics, 2013, 14(5): 1147-1173 doi: 10.4208/cicp.190912.080213a
[36] Liu C, Xu K, Sun QH, et al. A unified gas kinetic scheme for continuum and rarefied flows IV: full Boltzmann and model equations. Journal of Computational Physics, 2016, 314: 305-340 doi: 10.1016/j.jcp.2016.03.014
[37] Wang Z, Yan H, Li QB, et al. Unified gas kinetic scheme for diatomic molecular flow with translational, rotational, and vibrational modes. Journal of Computational Physics, 2017, 350: 237-259 doi: 10.1016/j.jcp.2017.08.045
[38] Baranger C, Claudel J, Hérouard N, et al. Locally refined discrete velocity grids for stationary rarefied flow simulations. Journal of Computational Physics, 2014, 257: 572-593 doi: 10.1016/j.jcp.2013.10.014
[39] Chapman S, Cowling TG. The Mathematical Theory of Non-uniform Gases. Cambridge: Cambridge University Press, 1970
[40] Kudryavtsev AN, Shershnev AA. A numerical method for simulation of microflows by solving directly kinetic equations with WENO schemes. Journal of Scientific Computing, 2013, 57: 42-73 doi: 10.1007/s10915-013-9694-z
[41] Bird GA. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Oxford: Claredon Press, 1994
[42] Bird GA. The DSMC Method. Create Space Independent Publishing Platform, 2013
计量
- 文章访问数: 1556
- HTML全文浏览量: 515
- PDF下载量: 77
