高超音速稀薄环境中的电磁流控制

《Journal of Fluid Mechanics》：Electromagnetic flow control in hypersonic rarefied environment

　　摘要：统一气体动力学波粒方法（UGKWP）是为部分电离等离子体的多尺度模拟而开发的，它被用于模拟从近似连续态到稀薄条件范围内的半球形区域的电磁流动。据我们所知，这项研究是首次将多尺度等离子体求解器应用于此类问题。在该方法中，中性粒子、离子和电子被视为不同的物种。通过将计算结果与马赫数为4.75的预电离氩气的实验数据进行比较，验证了该方法的数值实现能力，结果显示UGKWP的结果与实验数据非常吻合。进一步在不同克努森数（Knudsen number）下的对比研究表明，稀薄效应减弱了电磁控制的影响。这些发现突显了UGKWP在模拟电磁控制问题中的能力，并强调了稀薄效应对预测流动控制行为的重要作用，从而强调了多尺度建模在等离子体流动应用中的必要性。

1. 引言
磁流体动力学（MHD）热防护技术对于国防系统和深空探索任务具有战略意义。在高速飞行期间，强烈的外部磁场与强冲击波产生的弱电离气体相互作用，从而增加了冲击波的缓冲距离，减少了传递到飞行器表面的热流。这项技术的理论基础可以追溯到20世纪50年代的美苏太空竞赛时期，当时Resler Jr和Sears首次提出了基本的磁气动方程（参考Resler和Sears1958）。他们将高速、高焓流中的弱电离等离子体建模为导电流体，将洛伦兹力（Lorentz force）和焦耳热（Joule heating）项纳入修改后的纳维-斯托克斯方程（Navier–Stokes equations）。在随后的几十年里，MHD热防护的理论、计算和实验方面都取得了显著进展。Bush提供了描述磁场对冲击波缓冲距离影响的解析解（参考Bush1958），而Ziemer和Bush使用马赫数为4.5的球体模型进行了实验验证，标志着早期的重要里程碑。Cambel和Porter进一步研究了霍尔效应（Hall effect）对磁控制效率的负面影响，并使用预电离氩气进行了相应的实验，创建了对数值验证有价值的基准数据集。在这些基础研究之后，计算流体动力学的进步使得能够高保真地模拟MHD流动（参考Poggie和Gaitonde2002；Fujino, Yoshino和Ishikawa2008；Kai, Jun和Weiqiang2017）。最近的数值研究越来越多地纳入了复杂的物理现象，如稀薄气体效应以及等离子体鞘层与壁面的相互作用。例如，Fawley等人使用直接模拟蒙特卡罗方法（DSMC）来模拟过渡态流中的导电性（参考Fawley, Putnam, D’Souza和Borner2022），而Katsurayama等人应用DSMC来模拟稀薄气体效应，将洛伦兹力近似为体力（body force）而不是显式求解带电粒子轨迹（参考Katsurayama, Kawamura, Matsuda和Abe2008）。Jin, Peng和Feng使用扩展的纳维-斯托克斯方法研究了稀薄MHD流动，突出了稀薄气体效应对近连续态下MHD模拟的重要性。Parent等人进一步量化了等离子体鞘层对电子温度的影响（参考Parent, Rajendran, Macheret, Little, Moses, Johnston和Ch Cheatwood2023）。尽管取得了这些进展，大多数数值模拟仍然依赖于MHD方程，但在以高克努森数为特征的高海拔条件下，这些方程的适用性会降低。在这种条件下，基于MHD的模型所固有的连续性假设不再成立（参考Lani2023）。由于需要过于密集的网格和极小的时间步长（这些条件受到平均自由路径尺度（mean free path scale）的限制，基于DSMC的方法在计算上变得不可行（参考Boyd和Schwartzentruber2017）。因此，准确捕捉连续态、过渡态和稀薄态下的流动物理特性需要一种真正多尺度的数值方法。为了解决这些挑战，最近的发展产生了统一气体动力学波粒方法（UGKWP），该方法能够稳健地模拟从稀薄态到连续态的多尺度中性气体流动（参考Zhu, Liu, Zhong和Xu2019；Liu, Zhu和Xu2020）。此后，UGKWP框架系统地扩展到了各种传输现象，包括等离子体（参考Liu和Xu2021；Pu, Liu和Xu2024）、颗粒流动（参考Yang, Shyy和Xu2022，参考Yang, Shyy和Xu2024）、辐射传输（参考Liu, Li, Wang, Song和Xu2023；Yang, Zhu, Liu和Xu2025）以及声子传输（参考Liu, Yang, Zhang, Ji和Xu2025），证明了其捕捉非平衡传输物理现象的能力。最近，UGKWP特别被应用于部分电离等离子体流动（参考Pu和Xu2025）。在连续态下，UGKWP能够准确地再现理想MHD、耗散MHD和多流体模型；当流动过渡到稀薄条件时，该方法无缝地转向基于粒子的描述，有效地捕捉了无碰撞极限下的等离子体行为，并再现了粒子模拟（PIC）动力学。在这项工作中，我们应用UGKWP方法来模拟稀薄环境中的电磁控制问题。本文的其余部分组织如下：第2节简要介绍了用于部分电离等离子体的UGKWP方法；第3节展示了数值测试结果；第4节提供了本研究的主要结论。

2. UGKWP方法
2.1. 控制方程
在本小节中，我们介绍了本研究中使用的物理模型和控制方程。理论上，Bhatnagar–Gross–Krook（BGK）–Maxwell系统可以用来模拟部分电离等离子体（参考Pu和Xu2025）：
\[
\begin{equation}
\begin{aligned}
& \frac{\partial f_\alpha}{\partial t}+\boldsymbol{u}_\alpha \boldsymbol{\cdot }\boldsymbol{\nabla}_{\boldsymbol{x}} f_\alpha + \frac{q_\alpha (\boldsymbol{E}+\boldsymbol{u}_\alpha \times \boldsymbol{B})}{m_\alpha} \boldsymbol{\cdot }\boldsymbol{\nabla}_{\boldsymbol{u}} f_\alpha = Q_{\alpha}, \\
& \frac{\partial \boldsymbol{B}}{\partial t}+\boldsymbol{\nabla}_{\boldsymbol{x}} \times \boldsymbol{E}=0, \quad \frac{\partial \boldsymbol{E}}{\partial t}-c^2 \boldsymbol{\nabla}_{\boldsymbol{x}} \times \boldsymbol{B}=-\frac{1}{\epsilon_0} \boldsymbol{J}, \\
&\boldsymbol{\nabla}_{\boldsymbol{x}} \boldsymbol{\cdot }\boldsymbol{E}=\frac{\rho_c}{\epsilon_0},\quad \boldsymbol{\nabla}_{\boldsymbol{x}} \boldsymbol{\cdot }\boldsymbol{B}=0,
\end{aligned}
\end{equation}
\]
其中
$f_\alpha = f_\alpha (t, \boldsymbol{x}, \boldsymbol{u})$ 是物种$\alpha$（$\alpha =i$表示离子，$\alpha =e$表示电子，$\alpha =n$表示原子）在空间和时间$(\boldsymbol{x},t)$处的微观分布函数，$\boldsymbol{u}$是微观平动速度；$q_\alpha$是物种$\alpha$携带的电荷；$\boldsymbol{E}, \boldsymbol{B}$是电场和磁场强度；$\epsilon_0$是真空的介电常数，$c$是光速；$\boldsymbol{J}$是电流密度；$\rho_c$是电荷密度。这里，$Q_{\alpha}=\sum_{k\neq \alpha}Q_{\alpha k}(f_\alpha, f_k)$是物种$\alpha$与物种$k$之间的碰撞算子，通过Andries–Aoki–Perthame（AAP）模型来近似（参考Andries, Aoki和Perthame2002）：
\[
\begin{equation}
Q_\alpha =\frac{g_\alpha ^M-f_\alpha}{\tau _\alpha},
\end{equation}
\]
其中$g_\alpha ^{M}$是麦克斯韦分布，后碰撞速度$\bar {\boldsymbol{U}}_\alpha$和温度$\bar {T}_\alpha$定义为：
\[
\begin{align}
\bar {\boldsymbol{U}}_\alpha & =\boldsymbol{U}_\alpha +\frac{\tau _\alpha}{m_\alpha} \sum_{k=1}^N 2 \mu _{\alpha k} \chi _{\alpha k} n_k \left (\boldsymbol{U}_k-\boldsymbol{U}_\alpha \right), \\
\bar {T}_\alpha & =T_\alpha -\frac{m_\alpha}{3 k_B}\left (\bar {\boldsymbol{U}}_\alpha -\boldsymbol{U}_\alpha \right)^2+\tau _\alpha \sum_{k=1}^N \frac{4 \mu _{\alpha k} \chi _{\alpha k} n_k}{m_\alpha +m_k}\left (T_k-T_\alpha +\frac{m_k}{3 k_B}\left (\boldsymbol{U}_k-\boldsymbol{U}_\alpha \right)^2\),
\end{align}
\]
其中$\mu _{\alpha k} = m_{\alpha }m_k/(m_{\alpha }+m_k)$是减少了的质量，$N$是系统中的物种数量。物种$\alpha$的平均松弛时间$\tau _\alpha$由$1 / \tau _\alpha =\sum_{k=1}^m \chi _{\alpha k} n_{k}$确定，变硬球（VHS）模型的相互作用系数$\chi _{\alpha k}$为（参考Morse1963）：
\[
\begin{equation}
\chi _{\alpha k}= \frac{(5-2\omega )(7-2\omega )\sqrt{\pi}}{15}\left (\frac{2 k_B T_{\alpha}}{m_\alpha }+\frac{2 k_B T_{k}}{m_k}\right )^{1 / 2}\left (\frac{d_\alpha +d_k}{2}\right)^2,
\end{equation}
\]
温度改变了分子直径，表达式为$ ({d}/{d_{\textit{ref}}})^2 = ( {T}/{T_{\textit{ref}}} )^{({1}/{2}) - \omega },$其中$\omega$是粘度和温度之间的幂系数（参考Bird1976）。后续的数值研究中将给出每种物种的具体直径。对于涉及中性和带电物种的复杂相互作用，将使用实验衍生的拟合曲线。AAP模型的普朗特数（Prandtl number）为单位。未来的工作将扩展多松弛时间、多物种动力学模型，以更准确地捕捉粘度和热导率的同时行为，参考Bisi等人（参考Bisi, Boscheri, Dimarco, Groppi和Martalò2022）以及Li, Zeng和Wu（参考Li, Zeng和Wu2024）的研究。对于高超音速应用，BGK–Maxwell系统可以简化。在许多此类应用中，冲击波层内的电离气体表现出较低的电导率$\sigma ^{\textit{cond}}$（参考Poggie和Gaitonde2002）。这根据欧姆定律导致较小的感应电流密度$\boldsymbol{J}$（$\boldsymbol{J} = \sigma ^{\textit{cond}}(\boldsymbol{E}+\boldsymbol{U}\times \boldsymbol{B})$）。因此，该电流产生的磁场显著弱于外部施加的磁场$\boldsymbol{B}_{\textit{ext}}$。这允许忽略磁场对流动场的畸变，从而可以在本研究中假设一个静态磁场。法拉第定律可以忽略，因为静态磁场不会产生横向电场。此外，双曲散度清除技术实现了电荷中性（参考Munz, Omnes, Schneider, Sonnendrücker和Voss2000）：
\[
\begin{align}
\epsilon_0\frac{\partial \boldsymbol{E}}{\partial t} + c^2\boldsymbol{\nabla }\phi = -\boldsymbol{J},\quad \frac{\partial \phi}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{E} = \frac{\rho _c}{\epsilon_0},
\end{align}
\]
其中$\phi$是人为校正势。为了减少数值刚性，实施了大约10^-4米的德拜长度（Debye length）。因此，物理模型的最终形式为：
\[
\begin{equation}
\begin{aligned}
& \frac{\partial f_\alpha}{\partial t}+\boldsymbol{u}_\alpha \boldsymbol{\cdot }\boldsymbol{\nabla}_{\boldsymbol{x}} f_\alpha + \frac{q_\alpha (\boldsymbol{E}+\boldsymbol{u}_\alpha \times \boldsymbol{B}_{\textit{ext}})}{m_\alpha} \boldsymbol{\cdot }\boldsymbol{\nabla}_{\boldsymbol{u}} f_\alpha = Q_{\alpha}, \\
&\epsilon_0\frac{\partial \boldsymbol{E}}{\partial t} + c^2\boldsymbol{\nabla }\phi = -\boldsymbol{J},\quad \frac{\partial \phi}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{E} = \frac{\rho _c}{\epsilon _0}.
\end{aligned}
\end{equation}
\]
在流体极限$\tau _{\alpha} \rightarrow 0$下，（2.7）转变为多流体方程：
\[
\begin{equation}
\begin{aligned}
& \partial _t \rho _n+\boldsymbol{\nabla}_{\boldsymbol{x}} \boldsymbol{\cdot }\left (\rho _n \boldsymbol{U}_n\right)=0, \\
& \partial _t\left (\rho _n \boldsymbol{U}_n\right )+\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left (\rho _n \boldsymbol{U}_n \boldsymbol{U}_n+p_n \mathbb{I}\right )=R_n, \\
& \partial _t \mathscr{E}_n+\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left (\left (\mathscr{E}_n+p_n\right ) \boldsymbol{U}_n\right )=S_n, \\
& \partial _t \rho _i+\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left (\rho _i \boldsymbol{U}_{\!i}\right )=0, \\
& \partial _t\left (\rho _i \boldsymbol{U}_{\!i}\right )+\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left (\rho _i \boldsymbol{U}_{\!i} \boldsymbol{U}_{\!i}+p_i \mathbb{I}\right )=q_in_i\left (\boldsymbol{E}+\boldsymbol{U}_{\!i} \times \boldsymbol{B}_{\textit{ext}}\right )+R_i, \\
& \partial _t \mathscr{E}_i+\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left (\left (\mathscr{E}_i+p_i\right ) \boldsymbol{U}_{\!i}\right )=q_in_i \boldsymbol{U}_{\!i} \boldsymbol{\cdot }\boldsymbol{E}+S_i, \\
& \partial _t \rho _e+\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left (\rho _e \boldsymbol{U}_e\right )=0, \\
& \partial _t\left (\rho _e \boldsymbol{U}_e\right )+\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left (\rho _e \boldsymbol{U}_e \boldsymbol{U}_e+p_e \mathbb{I}\right )=q_en_e\left (\boldsymbol{E}+\boldsymbol{U}_e \times \boldsymbol{B}_{\textit{ext}}\right )+R_e, \\
& \partial _t \mathscr{E}_e+\boldsymbol{\nabla} _{\boldsymbol{x}} \boldsymbol{\cdot }\left (\left (\mathscr{E}_e+p_e\right ) \boldsymbol{U}_e\right )=q_en_e \boldsymbol{U}_e \boldsymbol{\cdot }\boldsymbol{E}+S_e, \\
&\epsilon _0\frac{\partial \boldsymbol{E}}{\partial t} + c^2\boldsymbol{\nabla }\phi = -\boldsymbol{J},\quad \frac{\partial \phi}{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{E} = \frac{\rho _c}{\epsilon _0},
\end{aligned}
\]
其中$\rho _\alpha, \boldsymbol{U}_\alpha, \mathscr{E}_\alpha$表示物种$\alpha$的质量密度、宏观速度和能量密度，$p_\alpha$是物种$\alpha$的热压。这里，$R_\alpha, S_\alpha$是由于物种间动量和能量传递产生的源项。在强磁场和非相对论性条件下，（2.7）简化为MHD方程。Pu和Xu（参考Pu和Xu2025）提供了这一简化的详细渐近分析。首先，多物种碰撞算子被重写为（2.9）：
$$
Q_\alpha = \frac {g_\alpha ^{eq}-f_\alpha }{\tau _\alpha } + \frac {g_\alpha ^M-g_\alpha ^{eq}}{\tau _\alpha},
$$
其中 $g^{eq}_\alpha$ 表示相同物种内通过碰撞得到的麦克斯韦分布，$g^M_\alpha$ 表示不同物种之间碰撞后的麦克斯韦分布。然后，使用算子分裂方法来求解系统（2.7）和（2.9），时间步长为 $\Delta t$。进一步定义算子（2.10）至（2.13）：
$$
\begin{align}
&\mathcal{L}_1: \frac {\partial f_\alpha }{\partial t}+\boldsymbol{u}_\alpha \boldsymbol{\cdot }\boldsymbol{\nabla} _{\!\boldsymbol{x}} f_\alpha = \frac {g_\alpha ^{eq}-f_\alpha }{\tau _\alpha }, \\
&\mathcal{L}_2: \frac {\partial f_\alpha }{\partial t}+ \frac {q_\alpha (\boldsymbol{E}+\boldsymbol{u}_\alpha \times \boldsymbol{B}_{\textit{ext}})}{m_\alpha } \boldsymbol{\cdot }\boldsymbol{\nabla} _{\!\boldsymbol{u}} f_\alpha = 0,\quad \epsilon _0\frac {\partial \boldsymbol{E}}{\partial t} = -\boldsymbol{J}, \\
&\mathcal{L}_3: \frac {\partial f_\alpha }{\partial t} = \frac {g_\alpha ^M-g_\alpha ^{eq}}{\tau _\alpha }, \\
&\mathcal{L}_4:\epsilon _0\frac {\partial \boldsymbol{E}}{\partial t} + c^2\boldsymbol{\nabla }\phi = 0,\quad \frac {\partial \phi }{\partial t} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{E} = \frac {\rho _c}{\epsilon _0}.
\end{align}
$$
每次更新后，变量表示为：
$$
\begin{align}
&\mathcal{L}_1: f^{n}_\alpha \rightarrow f^{*}_\alpha, \boldsymbol{W}^{n}_\alpha \rightarrow \boldsymbol{W}^{*}_\alpha, \\
&\mathcal{L}_2: f^{*}_\alpha \rightarrow f^{**}_\alpha, \boldsymbol{W}^{*}_\alpha \rightarrow \boldsymbol{W}^{**}_\alpha, \boldsymbol{E}^{n} \rightarrow \boldsymbol{E}^{*}, \\
&\mathcal{L}_3: f^{**}_\alpha \rightarrow f^{n+1}_\alpha, \boldsymbol{W}^{**}_\alpha \rightarrow \boldsymbol{W}^{n+1}_\alpha, \\
&\mathcal{L}_4: \boldsymbol{E}^{*} \rightarrow \boldsymbol{E}^{n+1},\phi ^{n} \rightarrow \phi ^{n+1},
\end{align}
$$
其中 $\boldsymbol{W}_{\alpha } = (\rho _{\alpha }, \rho _{\alpha } \boldsymbol{U}_{\alpha }, \rho _{\alpha } \mathscr{E}_{\alpha })$ 表示物种 $\alpha$ 的质量、动量和能量密度。首先，我们介绍求解 $\mathcal{L}_1$ 的算法。为简化起见，这里省略了物种下标 $\alpha$。不考虑力项的BGK方程表示为（2.15）：
$$
\frac {\partial f}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla} _x f = \frac {g-f}{\tau },
$$
界面处分布函数的时间依赖解可以表示为（2.16）：
$$
f\left (\boldsymbol{x},\boldsymbol{u}, \boldsymbol{\xi }, t\right )=\frac {1}{\tau } \int _{0}^{t} g(\boldsymbol{x}^{\prime }, \boldsymbol{u}, \boldsymbol{\xi },t^{\prime }) e^{-\left (t-t^{\prime }\right ) / \tau }\, \mathrm{d} t^{\prime } + e^{-t / \tau } f_{0}\left (\boldsymbol{x}-\boldsymbol{u} t \right ),
$$
其中 $\tau$ 是局部平均弛豫时间，$f_0$ 是 $t=0$ 时的初始气体分布函数，$g$ 是沿着特征线 $\boldsymbol{x}^{\prime } = \boldsymbol{x} - \boldsymbol{u} t^{\prime }$ 的平衡分布。将平衡分布函数用泰勒级数展开到二阶精度，得到（2.17）：
$$
g^{\prime }=g+ {\boldsymbol{g}}_x \boldsymbol{\cdot }(\boldsymbol{x}^{\prime }-\boldsymbol{x})+g_t(t^{\prime }-t),
$$
将它们代入（2.16），可以得到模拟粒子的数值多尺度演化解（2.18）：
$$
f(\boldsymbol{x}, \boldsymbol{u}, \boldsymbol{\xi }, t)=\left (1-e^{-t / \tau }\right ) g^{+}(\boldsymbol{x}, \boldsymbol{u}, \boldsymbol{\xi }, t )+e^{-t / \tau } f_{0}\left (\boldsymbol{x} - \boldsymbol{u}t\right ),
$$
其中（2.19）描述了BGK方程的解，即时间 $t$ 时的分布函数 $f$ 是初始分布函数 $f_0$ 与平衡状态 $g$ 的泰勒展开的结合。从粒子的角度来看，这意味着粒子在时间区间 $[0,t]$ 内有 $e^{-t/\tau }$ 的概率自由流动，有 $(1 - e^{-t/\tau })$ 的概率与其他粒子碰撞。经过多次碰撞后，粒子分布将达到平衡状态 $g^+$，其自由度可能会降低。基于这种物理图景，UGKWP方法以微观-宏观耦合的方式表示速度分布函数。一部分分布通过平衡分布函数 $g^+$ 理论上捕获，其余部分由随机模拟粒子 $P_k = (m_k, \boldsymbol{x}_k, \boldsymbol{u}_k)$ 表示。这里，$m_k$ 表示模拟粒子 $P_k$ 的质量，对应于相同物种的一簇真实气体粒子，$\boldsymbol{x}_k$ 和 $\boldsymbol{u}_k$ 分别表示模拟粒子 $P_k$ 的位置和速度。根据（2.18），粒子碰撞前的自由流动时间 $t_{\!f}$ 的累积分布函数表示为（2.20）：
$$
F ( t_{\!f} \lt t) = \exp \left ( - t\text{/}\tau \right ),
$$
自由流动时间 $t_{\!f}$ 可以通过 $t_{\!f} = - \tau \ln (\eta )$ 抽样得到，其中 $\eta$ 是服从均匀分布 $\eta \sim U(0,1)$ 的随机变量。对于时间步长 $\Delta t$，满足 $t_{\!f} \geqslant \Delta t$ 的粒子将经历无碰撞的自由流动，而满足 $t_{\!f} \lt \Delta t$ 的粒子将经历碰撞相互作用。UGKWP方法中更新粒子的过程如下：
步骤1. 在时间步长内，让每个粒子 $P_k$ 流动 $\min (\Delta t, t_{\!f,k})$ 的时间。然后，识别并保留无碰撞粒子，同时移除有碰撞的粒子。计算粒子通过细胞界面贡献的自由传输通量，并累加粒子的总守恒量 $\boldsymbol{W}_{\!i}^p$。
步骤2. 更新宏观守恒变量后，从更新后的守恒量 $\boldsymbol{W}_{\!i}$ 计算有碰撞粒子的总守恒量 $\boldsymbol{W}_{\!i}^h$，即 $\boldsymbol{W}_{\!i}^h = \boldsymbol{W}_{\!i} - \boldsymbol{W}_{\!i}^p$。
步骤3. 在时间步长结束时，重建速度分布。计算分析分布 $g^{+,c}$ 并根据更新后的守恒量从分布 $g^{+,f}$ 重新抽样无碰撞粒子，并从累积分布函数 $F(t_{\!f} \lt t) = \exp (-t/\tau )$ 抽样每个粒子 $P_k$ 的自由流动时间 $t_{\!f,k}$。有关该过程的更多详细描述，请参阅 Pu & Xu（参考文献 Pu 和 Xu2025）。在已经概述的过程中，下一段介绍了更新宏观守恒变量的算法。现在使用粒子数值求解微观速度分布的演化。基于有限体积框架，物理细胞 $\varOmega _i$ 上的细胞平均宏观守恒变量通过细胞界面的数值通量进行更新（2.21）：
$$
(\boldsymbol{W}_{\alpha })_i^{n+1} = (\boldsymbol{W}_{\alpha })_i^n - \frac {\Delta t}{|\varOmega _i|}\sum _{s\in \partial \varOmega _i}|l_s|(\mathscr{F}_{\boldsymbol{W}_\alpha })_s,
$$
其中 $l_s\in \partial \varOmega _i$ 是中心为 $\boldsymbol{x}_s$ 且具有外法线向量 $\boldsymbol{n}_{s}$ 的细胞界面；$|l_s|$ 是细胞界面的面积。界面上的数值通量 $(\mathscr{F}_{\boldsymbol{W}_\alpha })_s$ 可以从界面处的分布函数计算得到（2.22）：
$$
(\mathscr{F}_{\boldsymbol{W}_\alpha })_s = \frac {1}{\Delta t}\int _{t^n}^{t^{n+1}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}_sf_\alpha (\boldsymbol{x}_{s},\boldsymbol{u},\boldsymbol{\xi },t)\boldsymbol{\varPsi } \mathrm{d}\boldsymbol{\varXi }\mathrm{d}t,
$$
其中 $\boldsymbol{\varPsi }=(1,\boldsymbol{u}, ({1}/{2})(\boldsymbol{u}^2+\boldsymbol{\xi }^2)$ 是分布函数的守恒矩，$\mathrm{d}\boldsymbol{\varXi }=\mathrm{d}\boldsymbol{u}\,d\boldsymbol{\xi }$ 是相空间中的体积元素。对于细胞界面 $l_s$，在过渡和无碰撞状态下，分布函数 $f_\alpha (\boldsymbol{x}_{s},\boldsymbol{u},\boldsymbol{\xi },t)$ 无法用解析方法表示，因此需要像前一段介绍的那样直接跟踪其演化。这种宏观守恒量的数值通量可以根据（2.16）分为平衡通量和自由流动通量。平衡通量为（2.23）：
$$
(\mathscr{F}_{\boldsymbol{W}})_s^g = \frac {1}{\Delta t}\int _{t^n}^{t^{n+1}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}_s \left [\frac {1}{\tau } \int _{0}^{t} g(\boldsymbol{x}^{\prime }, \boldsymbol{u}, \boldsymbol{\xi },t^{\prime }) e^{-\left (t-t^{\prime }\right ) / \tau } \mathrm{d} t^{\prime }\right ] \boldsymbol{\varPsi } \,\mathrm{d}\boldsymbol{\varXi }\,\mathrm{d}t,
$$
自由流动通量为（2.24）：
$$
(\mathscr{F}_{\boldsymbol{W}})_s^f = \frac {1}{\Delta t}\int _{t^n}^{t^{n+1}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}_s \left [e^{-t / \tau } f_{0}\left (\boldsymbol{x}-\boldsymbol{u} t \right )\right ] \boldsymbol{\varPsi } \,\mathrm{d}\boldsymbol{\varXi }\,\mathrm{d}t.
$$
平衡通量可以直接从宏观守恒变量计算得到。假设在（2.22）中，$\boldsymbol{x}_{s} = \boldsymbol{0}$ 且 $t^{n} = 0$，则细胞界面的平衡状态 $g$ 可以通过守恒约束（2.25）得到：
$$
\int _{}^{}{g\boldsymbol{\varPsi }}\,{\rm d}\boldsymbol{\varXi } = \int _{\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{n} \gt 0}^{}{g^{l}\boldsymbol{\varPsi }}\,{\rm d}\boldsymbol{\varXi } + \int _{\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{n} \lt 0}^{}{g^{r}\boldsymbol{\varPsi }}\,{\rm d}\boldsymbol{\varXi }.
$$
可以获得空间和时间导数（2.26）：
$$
\int _{}^{}{g_x\boldsymbol{\varPsi }}\,{\rm d}\boldsymbol{\varXi } = \boldsymbol{W}_x, \quad \int _{}^{}{g_t\boldsymbol{\varPsi }}\,{\rm d}\boldsymbol{\varXi } = -\!\! \int _{}^{}\boldsymbol{u} \boldsymbol{\cdot }g_x\boldsymbol{\varPsi }\,{\rm d}\boldsymbol{\varXi },
$$
其中 $g^{l}$ 和 $g^{r}$ 分别是根据重构的细胞界面左右侧守恒变量 $\boldsymbol{W}^{l},{\ \boldsymbol{W}}^{r}$ 得到的平衡分布。使用 van Leer 限制器在空间重构中实现二阶精度。在梯度较大的区域，限制器限制梯度以防止非物理振荡。将重构的平衡分布代入平衡通量，得到（2.27）：
$$
(\mathscr{F}_{\boldsymbol{W}})_s^g = \int _{}^{}\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{n}_{s}\left ( C_{1}g_{0} + C_{2}\boldsymbol{u} \boldsymbol{\cdot }{\boldsymbol{g}}_{0x} + C_{3}g_{0t} \right )\varPsi \,{\rm d}\boldsymbol{\varXi },
$$
时间积分系数为（2.28）至（2.30）。接下来，我们考虑自由流动通量 $(\mathscr{F}_{\boldsymbol{W}})_s^f$。如上一小节所述，初始分布部分由解析分布 $g_{a}^{+,c}$ 表示，部分由粒子表示，因此自由流动通量也部分从重构的解析分布 $(\mathscr{F}_{\boldsymbol{W}})_s^{f,w}$ 计算，部分从粒子计算 $(\mathscr{F}_{\boldsymbol{W}})_s^{f,p}$。初始解析分布 $g_{\alpha }^{+,c}$ 的重构过程为（2.31）：
$$
g_{0}^{+,c}\left ( \boldsymbol{x},\boldsymbol{u} \right ) = g_{0}^{+,c} + {\boldsymbol{g}}_{0x}^{+,c} \boldsymbol{\cdot }\boldsymbol{x},
$$
由此得到（2.32）：
$$
(\mathscr{F}_{\boldsymbol{W}})_s^{f,w} = \int _{}^{}\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{n}_s\big ( C_{4}g_{0}^{+,c} + C_{5}\boldsymbol{u} \boldsymbol{\cdot }{g}_{0x}^{+,c} \big )\boldsymbol{\varPsi }\,{\rm d}\boldsymbol{\varXi },
$$
时间积分系数为（2.33）至（2.34）。净粒子通量 $(\mathscr{F}_{\boldsymbol{W}})_s^{f,p}$ 的计算过程为（2.35）：
$$
(\mathscr{F}_{\boldsymbol{W}})_s^{f}= \sum _{k \in P_{\partial \varOmega _{i}^{+}}}^{}\boldsymbol{W}_{P_{k}} - \sum _{k \in P_{\partial \varOmega _{i}^{-}}}^{}\boldsymbol{W}_{P_{k}},
$$
其中 $\boldsymbol{W}_{P_{k}} = ( m_{k},m_{k}\boldsymbol{v}_{k}, ({1}/{2}) m_{k}\boldsymbol{v}_{\boldsymbol{k,}}^{2} ),P_{\partial \varOmega _{i}^{-}}$ 是在时间步长内流出细胞 $\varOmega _{i}$ 的粒子索引集，$P_{\partial \varOmega _{i}^{+}}$ 是流入细胞 $\varOmega _{i}$ 的粒子索引集。最后，守恒变量的更新公式为（2.36）：
$$
\begin{split}
\boldsymbol{W}_{\boldsymbol{i}}^{*} = \boldsymbol{W}_{\boldsymbol{i}}^{\boldsymbol{n}} -&\sum _{s}^{}{\frac {\Delta t}{\left | \varOmega _{i} \right |}\left | l_{s} \right |(\mathscr{F}_{\boldsymbol{W}})_s^{g}} - \sum _{s}^{}{\frac {\Delta t}{\left | \varOmega _{i} \right |}\left | l_{s} \right |(\mathscr{F}_{\boldsymbol{W}})_s^{f,w}} + \frac {1}{\left | \varOmega _{i} \right |}(\mathscr{F}_{\boldsymbol{W}})_s^{f,p}.
\end{split}
$$
通过耦合微观粒子传输和宏观流体动力学的演化，已经解决了宏观守恒变量的演化。接下来，我们求解 $\mathcal{L}_2$。首先隐式求解（2.11）的矩方程。然后使用更新后的 $\boldsymbol{E}^{*}$ 加速粒子。关于矩方程和隐式线性系统的更多细节，请参见 Pu & Xu（参考文献 Pu 和 Xu2025）。接下来，我们求解 $\mathcal{L}_3$；（2.12）的矩方程是一组非线性常微分方程（2.37至2.34）。使用后向欧拉方案离散化并用牛顿迭代求解器求解。最后，在 $\mathcal{L}_4$ 中，通过通量分裂方法求解（2.13）。整个算法流程图如图1所示。图1. UGKWP方法的流程图。与DSMC和PIC方法相比，当前方法避免了与网格大小相关的限制性稳定性条件，即时间步长必须小于平均自由程和平均碰撞时间。为了获得准确的物理分辨率，需要一个大约为$\sqrt{\lambda}$的空间网格尺寸$\Delta x$，其中$\lambda$表示平均自由程（参考文献Guo, Li和Xu 2023）。这种对网格约束的放宽，特别是在$\lambda$较小时的近连续介质状态下，与传统粒子方法相比，显著提高了计算效率。3. 数值研究 在本节中，我们对马赫数为4.75的氩气在半球周围的流动进行了三维模拟，复制了Cambel（参考文献Cambel 1967）报告的实验条件。这里采用的UGKWP方法已经通过大量基准测试和实际稀薄气体及等离子体问题得到了验证（参考文献Zhu, Liu, Zhong和Xu 2019；Liu等人参考文献Liu, Zhu和Xu 2020；Liu和Xu参考文献Liu和Xu 2021；Wei等人参考文献Wei, Cao, Ji和Xu 2023；Long, Wei和Xu参考文献Long, Wei和Xu 2024；Pu和Xu参考文献Pu和Xu 2025；Cao等人参考文献Cao, Wei, Long, Zhong和Xu 2026）。为了进一步验证该代码在高超音速等离子体问题中的准确性，我们首先将结果与实验结果进行比较。然后，我们改变氩原子的克努森数，以观察在不同稀薄程度下控制效果的变化。3.1. 与实验的验证 模拟中，一个高超音速的来流沿$-x$方向撞击半球，如图2(a)所示。半球的半径$r_n = 0.01905$米。如图2(b)所示，模拟了半球的$1/4$部分。第一层的高度为0.00005米，总网格数量为7488个。我们进行了网格独立性研究，使用了五种不同网格分辨率，其第一层高度分别为h = 0.0004、0.0002、0.0001、0.00005和0.000025米。相对于半球体直径（D约0.04米），这些比率分别约为0.01、0.005、0.0025、0.00125和0.000625。每种网格下的氩原子表面压力和热流分布结果在图3中进行了比较。表面压力分布（图a）在所有网格分辨率下几乎保持不变，表明收敛速度很快。相比之下，表面热流（图b）对网格细化更为敏感，随着网格的细化逐渐减小，直到收敛。h = 0.00005米和h = 0.000025米的结果基本相同，证实了网格的独立性。基于此分析，所有后续模拟都采用了h = 0.00005米的网格。图2. 半球上高超音速流的计算域规格：(a)情况设置；(b)网格设置。图3. 不同网格尺寸下的表面压力和热流分布：(a)表面压力；(b)表面热流。在实验中，入流是由等离子体火炬产生的部分电离的氩气流。Cambel（参考文献Cambel 1967）估计电离程度为$\alpha = 0.025$，然后Bisek、Boyd和Poggie（参考文献Bisek, Boyd和Poggie 2010）使用更精细的Saha方程重新估计电离程度为$\alpha = 0.00623$。本研究采用了这个修正后的值。Cambel（参考文献Cambel 1967）认为流动在化学上是冻结的，因此可以将电离程度视为在整个流场中恒定的。详细参数列在表1中。表1. 半球周围稀薄氩气流的模拟参数。在半球中心放置了一个磁偶极子，其磁矩沿x轴方向。其方向与沿停滞线的来流相反。这个偶极子产生的磁场定义为（3.1）\begin{align} \boldsymbol{B}=\frac {B_{\textit{max}}}{2\left (\hat {x}^2+\hat {y}^2+\hat {z}^2\right )^{5 / 2}}\left [\begin{array}{c} 2 \hat {x}^2-\left (\hat {y}^2+\hat {z}^2\right ) \\[3pt] 3 \hat {x} \hat {y} \\[3pt] 3 \hat {x} \hat {z} \end{array}\right ]\!, \end{align}，其中$\hat {x} = x/r_n$，$\hat {y} = y / r_n$，$\hat {z} = z / r_n$代表由半球半径$r_n$归一化的坐标。参数$B_{\textit{max}}$表示在停滞点达到的最大磁场强度。图4显示了(a)面板中的偶极子场磁场线以及(b)面板中的磁场强度大小。在接下来的分析中，我们使用$B_{\textit{max}}$来表征磁场强度。模拟在四个不同的磁场强度下进行：0.223 T、0.316 T、0.387 T和0.447 T。图4. (a)磁场线和(b)磁场强度设置$|\boldsymbol{B}|/|\boldsymbol{B}_{\textit{max}|$。在这种情况下，$\textit{Ar}$ – $\textit{Ar}$、$\textit{Ar}^{+}$ – $\textit{Ar}^{+}$、$e^{-}$ – $e^{-}$的碰撞被建模为VHS碰撞，其中$\omega =0.81$和$d=4.17\times {10^{-10}}$米。对于带电-带电相互作用（${\textit{Ar}^+{-}\textit{Ar}^+$和$e^-{-}e^-$），原则上库仑力是主导的。然而，在这里考虑的弱电离状态下（$n_e / n_{\textit{Ar}} \ll 1$），这种碰撞发生的频率远低于中性粒子相关的碰撞。因此，使用VHS模型来近似这两种相互作用只会引入较小的误差，这对于当前的等离子体条件是可接受的。对于混合物种碰撞$\textit{Ar}$ – $\textit{Ar}^+$、$\textit{Ar}$ – $e^-$、$\textit{Ar}^+$ – $e^-$，截面由Devoto（参考文献Devoto 1973）提供的拟合曲线描述，这些曲线适用于复合散射机制，如极化和电荷交换效应，而这些效应无法通过简化的碰撞模型捕捉到。电子质量被人为增加（$m_{\textit{Ar}^+}/m_{e^-}=10$）以便于计算处理。然而，这种修改降低了电导率。因此，磁场强度也被人为增加，按照相似性参数$S= {\sigma ^{\textit{cond}}|\boldsymbol{B}_0|^2L_0}/{\rho _0 |\boldsymbol{U}_0|}$进行了缩放，更多解释见附录A。在固定的磁场强度$B_{\textit{max}} = 0.223$ T下，对质量比从10到40进行了敏感性实验。如图5所示，Ar的数密度剖面在考虑的范围内只有轻微的变化，支持了这种方法的有效性。为了判断模拟的收敛状态，我们监测了平均关键表面量的时间演变，特别是压力和热流。当这些被监测的量稳定时，认为模拟已经收敛，表明没有显著变化。表2. Ar–Ar+、Ar–e?、Ar+–e?的截面（单位是?2）。图5. 不同质量比下$\textit{Ar}$的数密度剖面比较（$B_{\textit{max}}=0.223T$）。图6. 无外部磁场（上）和有外部磁场（下）时的等高线图，其中第一行是(a) $\textit{Ar}$、(b) $\textit{Ar}^+$、(c) $e^-$的质量密度；第二行是(d) $\textit{Ar}$、(e) $\textit{Ar}^+$、(f) $e^-$的x速度；第三行是(g) $\textit{Ar}$、(h) $\textit{Ar}^+$、(i) $e^-$的温度，在$B_{\textit{max}}=0.447\,\mathrm{T}$时。图7. (a)–(d)电场强度$|\boldsymbol{E}|$和$(e)–(h)电流密度大小$\boldsymbol{J}$和$J_x, J_y, J_z$，(i)–(l)洛伦兹力大小$|\boldsymbol{F}|$和$F_x, F_y, F_z$在$B_{\textit{max}}=0.447\,\mathrm{T}$时的等高线图。上图和下图分别表示有无外部磁场的情况。图8. (a)有无应用磁场$B_{\textit{max}}=0.447\,\mathrm{T}$时，沿停滞线的离子的数密度和归一化x速度。当前公式允许观察原子、离子和电子之间的不同行为，包括由于稀薄效应引起的速度和温度滑移。图6和7显示了等高线分布，上半部分显示了无外部磁场时的模拟，下半部分显示了有背景磁场强度（$B_{\textit{max}}=0.447$ T）时的结果。数据证实了应用磁场可以增强冲击波的停滞距离。图6(b)和6(c)表明，在应用磁场后，离子和电子都显著地被捕获在停滞点附近。图8(a)显示了沿停滞线的离子和电子的相应数密度剖面。这些剖面清楚地表明，在应用磁场后，两种物种的数密度都增加了。这种增加是因为在这个停滞区域，强磁场显著抑制了横向粒子运动。因此，带电粒子更直接地被输送到壁面，提高了局部数密度。相比之下，在肩部区域（肩部区域包括大约$30^\circ$到$60^\circ$的弧度，停滞线作为$0^\circ$参考），洛伦兹力的影响更大，将带电粒子从壁面偏转开，导致与非磁性情况相比数密度降低，如图6(b)和6(c)所示。这些带电粒子与氩原子的碰撞导致停滞区域中氩原子的轻微积累，如图6(a)所示。图6(d)、6(e)和6(f)展示了$\textit{Ar}$、$\textit{Ar}^+$和$e^-$的x速度。在停滞区域，$\textit{Ar}^+$的速度紧密跟随$\textit{Ar}$的速度。在肩部区域，由于磁场的影响，它们的速度由于磁场的作用而分离，电子的速度显著偏离。图6(g)到6(i)展示了$\textit{Ar}$、$\textit{Ar}^+$和$e^-$的x速度。在停滞区域，$\textit{Ar}^+$的速度与$\textit{Ar}$的速度紧密跟随。在冲击波之外，由于磁场的影响，它们的速度分离得更为明显。图6(d)、6(e)和6(f)展示了$\textit{Ar}$、$\textit{Ar}^+$和$e^-$的x速度。在停滞区域内，$\textit{Ar}^+$的速度紧密跟随$\textit{Ar}$的速度。在冲击波之外，由于磁场的影响，它们的速度分离得更为明显。图6(b)和6(c)表明，在应用磁场后，离子和电子都显著地被捕获在停滞点附近。图8(a)显示了沿停滞线的离子和电子的数密度剖面。这些剖面清楚地表明，在应用磁场后，两种物种的数密度都增加了。这种增加是因为在这个停滞区域，强磁场显著抑制了横向粒子运动。因此，带电粒子更直接地被输送到壁面，提高了局部数密度。相比之下，在肩部区域（肩部区域包括大约$30^\circ$到$60^\circ$的弧度，停滞线作为$0^\circ$参考），洛伦兹力的影响更大，将带电粒子从壁面偏转开，导致非磁性情况下的数密度降低，如图6(b)和6(c)所示。这些带电粒子与氩原子的碰撞导致停滞区域中氩原子的轻微积累。图6(d)、6(e)和6(f)展示了$\textit{Ar}$、$\textit{Ar}^+$和$e^-$的x速度。在停滞区域内，$\textit{Ar}^+$的速度紧密跟随$\textit{Ar}$的速度。在冲击波之外，由于磁场的影响，它们的速度分离得更为明显。图6(g)到6(i)展示了$\textit{Ar}$、$\textit{Ar}^+$和$e^-$的x速度。在停滞区域，$\textit{Ar}^+$的速度紧密跟随$\textit{Ar}$的速度。在冲击波之外，由于磁场的影响，它们的速度分离得更为明显。图7(a)展示了有无外部磁场时的电场强度。在无磁场的情况下，电场通常较弱，在壁附近由于鞘层形成而产生局部峰值。在接近壁面的地方，离子密度超过电子密度，产生一个净正电荷，形成一个远离壁面的电场。在冲击波前沿，强压力梯度导致离子和电子之间的电荷分离。应用磁场改变了这种结构，将主要电场峰值移动到肩部区域。这种重新分布可以通过欧姆定律$\boldsymbol{E} \sim \boldsymbol{U} \times \boldsymbol{B}$来解释，该定律表明电场峰值出现在电子漂移速度和磁场的交叉积最大处，如图7(a)所示。图7(e)展示了停滞区域集中的强电流密度。这里的电流密度由$\boldsymbol{J} = q_in_i\boldsymbol{U}_{\!i} + q_en_e\boldsymbol{U}_e$计算得出。图7(i)揭示了该区域附近的强洛伦兹力。洛伦兹力计算为$\boldsymbol{F} = \boldsymbol{J}\times \boldsymbol{B}$。由于洛伦兹力$\boldsymbol{J}\times \boldsymbol{B}$的强度也极大地依赖于电流密度$\boldsymbol{J}$和磁场$\boldsymbol{B}$之间的角度，因此力并不完全与电流密度峰值对齐。图9. 在$B_{\textit{max}}=0.447\,\mathrm{T}$时，(a) $\textit{Ar}$、(c) $\textit{Ar}^+$、(e) $e^{-}$的归一化表面压力和(b) $\textit{Ar}$、(d) $\textit{Ar}^+$、(f) $e^{-}$的表面热流。图8(b)显示了沿停滞线的离子的归一化x分量速度（$U_x$）。磁场的应用导致了额外的减速。这种效应可以归因于磁镜力，其中磁矩$\mu = m v_\perp ^2 / (2|\boldsymbol{B}|$——其中$v_\perp$表示垂直于磁场的速度——在变化缓慢的磁场中大致保持不变（Lee, Kim和MacCormack参考文献Lee, Kim和MacCormack 2015）。当带电粒子接近壁面时，收敛的磁场线将垂直动能转换为平行运动，导致沿停滞线的减速。相应的电子流动重组在图10中展示，其中最初平滑的流线（图a）在磁化作用下变得无序（图b）。图9展示了$\textit{Ar}$、$\textit{Ar}^+$和$e^{-}$的归一化表面压力和热流分布。表面压力和热流的参考值分别为$\rho _{\textit{Ar},\infty } U_{\infty }^2/2$和$\rho _{\textit{Ar},\infty } U_{\infty }^3/2$，其中$\infty$表示自由流值。图9(a)显示，$\textit{Ar}$的表面压力在停滞区域附近增加，在肩部区域减少。对于$\textit{Ar}^+$，也观察到了类似的表面压力趋势。图9(b)显示，$\textit{Ar}$在停滞区域附近的热流减少，这归因于增加的 standoff距离导致温度梯度减小。相比之下，$\textit{Ar}^+$在停滞区域的热流增加，这可能是由于磁场和来自磁场的直接能量输入导致壁面附近的数密度显著增加。$e^{-}$的压力和热流分布与$\textit{Ar}^{+}$相似。鉴于$\textit{Ar}^+$的摩尔分数仅为0.00623，整体表面压力和热流主要由$\textit{Ar}$控制。在半球表面上观察到的停滞压力增加和热流减少与该领域大多数已建立的文献一致（Bush参考文献Bush1958；Poggie和Gaitonde参考文献Poggie and Gaitonde2002；Fujino等人参考文献Fujino, Yoshino and Ishikawa2008）。图10. 电子在壁附近的流线图 (a) 无磁场和 (b) 有磁场 $B_{\textit{max}}=0.447\,\mathrm{T}$。图11. 随磁场强度变化的冲击波standoff距离对于马赫数为4.75的预电离氩气流围绕半球几何体（Cambel参考文献Cambel1967）（实验不确定性为$\pm 10\,\%$）。图11展示了冲击波standoff距离与施加的磁场强度的依赖关系。冲击波前沿位置是使用密度跳跃比来确定的，其定义为（Bisek等人参考文献Bisek, Boyd and Poggie2010）\begin{align} \lim _{M_1\rightarrow \infty } \frac {\rho _2}{\rho _1} = \frac {\gamma + 1}{\gamma - 1}。\end{align} 这里，$M_1$表示上游马赫数，$\gamma$表示比热比，$\rho _1$是上游密度，$\rho _2$是下游密度。对于氩气，$\gamma _{\textit{Ar}} = 5/3$，这个密度比计算为4。我们的模拟结果与实验数据非常吻合。（实验不确定性来源于测量技术的限制。）除了实验结果外，还纳入了Bisek等人（参考文献Bisek, Boyd and Poggie2010）使用三种不同导电模型的模拟数据进行比较。当前公式与Bisek工作的主要区别在于我们考虑了电离气体中的非平衡传输现象，而Bisek的结果是基于MHD方程的。对于围绕半球的超音速流动，当稀薄效应显著时，随着Knudsen数的增加，standoff距离会增加（Long等人参考文献Long, Wei and Xu2024）。3.2. 不同Knudsen数下的控制效应 本小节研究了在不同Knudsen数（0.2, 0.08, 0.044和0.02）下控制参数的影响。Knudsen数由进入的氩气密度决定，对于Kn = 0.2设置为2.18 $\,\times\,\text{10}^{-5}$ kg $\,$ m $^{-3}$；对于Kn = 0.08设置为5.45 $\,\times\,\text{10}^{-5}$ kg $\,$ m $^{-3}$；对于Kn = 0.044设置为1.09 $\,\times\, \text{10}^{-4}$ kg $\,$ m $^{-3}$；对于Kn = 0.02设置为2.18 $\,\times\,\text{10}^{-4}$ kg $\,$ m $^{-3}$。计算几何形状、网格划分和其他详细参数与前一节相同。相似性参数$S$在这些Knudsen数变化中保持不变。根据（A3）和（A5），如（A3）所示，如果氩原子密度增加了一个因子$a$，则参数$S$会相应减少一个因子$a^2$。为了抵消这种效应并保持$S$不变，磁场强度$|\boldsymbol{B}|$必须放大同一个因子$a$。本节中的磁场强度设置为$B_{\textit{max}}=0.447\,\mathrm{T}$。图12. 在$B_{\textit{max}}=0.447\,\mathrm{T}$下，不同Knudsen数（a）Kn = 0.02，（b）Kn = 0.04，（c）Kn = 0.08，（d）Kn = 0.2的氩原子温度等高线图。上下图分别表示没有外部磁场和有外部磁场的情况。图13. 在$B_{\textit{max}}=0.447\,\mathrm{T}$下，不同Knudsen数（a）Kn = 0.02，（b）Kn = 0.04，（c）Kn = 0.08，（d）Kn = 0.2的氩原子的局部Knudsen数等高线图。上下图分别表示没有外部磁场和有外部磁场的情况。图12显示了在不同Knudsen数下的氩原子温度场。随着Knudsen数的减小，冲击波前沿显示出尖锐化的趋势。图13显示了氩原子的局部Knudsen数，以说明氩原子的传输非平衡性。局部Knudsen数定义为$Kn_{local} = {\lambda }/{\rho /|\boldsymbol{\nabla }\rho |}$，其中$\lambda$是局部分子平均自由路径；$\rho$是质量密度。图表显示，随着参考Knudsen数的增加，整个区域的Knudsen数也随之增加。特别是，在较高的Knudsen数Kn = 0.2时，冲击波内部的局部Knudsen数与冲击波外部的大小相当，表明即使在冲击波区域内也存在非平衡状态。图12和图13显示，随着Knudsen数的增加，standoff距离减小。因此，如表3和图14所示，总热流减少也呈现出递减趋势。热流是通过积分整个壁面积计算得出的。这种现象是由于在较高Knudsen数下原子与带电粒子之间的耦合减弱所致。图16和图15分别显示了停滞线沿线的归一化x速度和温度剖面。参考值是自由流速度和温度。如图所示，随着Knudsen数的增加，原子和离子之间的速度和温度差异逐渐增大，反映了越来越偏离平衡状态。这种趋势源于Knudsen数与数密度之间的反比关系：在较高的Knudsen数下，较低的原子密度减少了原子和离子之间的碰撞频率。相反，随着Knudsen数的减小——即原子密度的增加——碰撞变得更加频繁，使得速度和温度趋向平衡。表3. 随Knudsen数变化的热流减少量。图14. 随Knudsen数的变化热流减少百分比。图15. 在（a）Kn = 0.02，（b）Kn = 0.04，（c）Kn = 0.08，（d）Kn = 0.2时，停滞线上氩原子和氩离子的归一化温度剖面。温度差$\Delta T = T_{\textit{Ar}^+} - T_{\textit{Ar}}$也以灰色点显示。图16. 在（a）Kn = 0.02，（b）Kn = 0.04，（c）Kn = 0.08，（d）Kn = 0.2时，停滞线上氩原子和氩离子的归一化x速度剖面。温度差$\Delta U_x = U_{\textit{Ar}^+,x} - U_{\textit{Ar},x}$也以灰色点显示。图17. （a）随着Knudsen数的变化，停滞线上原子和离子之间的最大归一化速度差（面板a）和最大归一化温度差（面板b）。图17展示了随着Knudsen数的增加，原子和离子之间的最大速度差（面板a）和最大温度差（面板b）。这两种差值都大致呈线性增长。这一线性趋势直接对应于图14中观察到的表面热流随Knudsen数增加而减少的现象，证实了较高的稀薄程度导致带电粒子和中性粒子之间的平衡偏离更大。图18进一步展示了有和没有施加磁场时氩原子的x速度。磁场显著减慢了氩原子的速度，但随着Knudsen数的增加，这种减速的幅度减小。与之前的讨论一致，在较高的Knudsen数下，带电粒子对中性粒子的拖拽作用减弱，从而减少了减速效果。这种机制解释了为什么较高的Knudsen条件会减弱电磁控制的有效性。总之，增加Knudsen数增强了带电粒子和中性粒子之间的解耦，从而削弱了电磁控制的效力。图18. 有和没有施加磁场时，停滞线上氩原子的归一化x速度剖面（a）Kn = 0.02，（b）Kn = 0.04，（c）Kn = 0.08，（d）Kn = 0.2。速度差以灰色点表示。4. 结论 本研究将UGKWP方法应用于围绕半球的部分电离等离子体流动的多尺度模拟，涵盖了从接近连续态到高度稀薄的状态。通过对预电离氩气流的实验数据验证，确认了数值实现的准确性和鲁棒性。在不同Knudsen数下控制效应的比较研究表明，稀薄效应可以显著改变电磁流动控制的预测。未来的工作将扩展该框架，包括额外的物理现象，例如壁鞘动力学和在中等到高磁雷诺数下的 comprehensive 磁流体动力学效应。这些先进模拟的可靠性将继续依赖于UGKWP求解器所提供的鲁棒多尺度流场解决方案。总的来说，本研究为推进复杂等离子体流动现象的模拟奠定了坚实基础。致谢 第一位作者想要感谢W. Long、J. Cao、H. Liu和S. Qin在这个项目上的宝贵讨论。资金 本项工作得到了中国国家重点研发计划（项目编号2022YFA1004500）、国家自然科学基金（92371107）和香港研究资助委员会（16208324）的支持。利益声明 作者们报告没有利益冲突。附录A. 人造电子质量效应 根据Drude的导电模型（Yoshino, Fujino & Ishikawa参考文献Yoshino, Fujino and Ishikawa2010），\begin{align} \sigma ^{\textit{cond}} = \frac {n_e e^2}{m_e\sum _{s\neq e}\nu _{e,s}}, \end{align} 其中$\nu _{e,s}$是物种$s$与电子之间的碰撞频率，如（Chen等人参考文献Chen1984）\begin{align} \nu _{e,s} = n_s \sigma _{e,s}^{\textit{col}} \bar {v}_{th}, \end{align} 给出，其中$n_s$是物种$s$的数密度；$\sigma _{e,s}^{\textit{col}}$是物种$s$与电子之间的截面；$\bar {v}_{th}$是物种$s$与电子之间的平均相对速度，大约为电子的平均速度$\bar {v}_{th}=\sqrt { {8k_B T_e}/{\pi m_e}}$。因此，导电率表示为\begin{equation} \sigma ^{\textit{cond}} = \frac {n_e e^2}{m_e\sum _{s\neq e}n_s \sigma _{e,s}^{\textit{col}}\sqrt {\frac {8k_{\!B} T_e}{\pi m_e}}。\end{equation} 在当前工作中，$\sigma _{e,s}^{\textit{col}}$由拟合曲线给出，$n_e, T_e$和$n_s$由流场计算得出，因此它们的表达式中不包含$m_e$，那么我们有\begin{align} \sigma ^{\textit{cond}} \propto \frac {1}{\sqrt {m_e}}, \end{align} 这意味着随着$m_e$的增加，导电率也会减小。standoff距离与无量纲磁相互作用参数成正比（Poggie & Gaitonde参考文献Poggie and Gaitonde2002），\begin{align} S = \frac {\sigma ^{\textit{cond}}_0 |\boldsymbol{B}_0|^2 L_0}{\rho _0 |\boldsymbol{U}|_0}, \end{align} 其中下标$0$表示特征量。为了使$S$保持不变，当导电率减小时，$|\boldsymbol{B}_0|^2$必须增加。考虑到实际的氩离子（$\textit{Ar}^+$）与电子（$e^-$）的质量比大约为72811，而本文中使用的比率为10，因此电导率相应减少了大约85倍，所以需要将$|\boldsymbol{B}_0|^2$增加85倍以保持所需的相似性。附录B.跨物种碰撞项的离散化（2.38式）的离散形式为（B1）：
$$
\begin{align}
&\frac {\boldsymbol{U}_\alpha ^{n+1} - \boldsymbol{U}_\alpha ^{**}}{\Delta t} = \sum _{k=1}^N \frac {2 m_k \chi _{\alpha k} n_k}{m_\alpha + m_k} (\boldsymbol{U}_k^{n+1} - \boldsymbol{U}_\alpha ^{n+1}),\\
&\frac {E_\alpha ^{n+1} - E_\alpha ^{**}}{\Delta t} = \sum _{k=1}^N \frac {2m_{k}\chi _{\alpha k} n_k}{m_{\alpha }+m_{k}} \boldsymbol{U}_{\alpha }^{n+1}\big (\boldsymbol{U}_k^{n+1}-\boldsymbol{U}_\alpha ^{n+1}\big )\\
&+\sum _{k=1}^N \frac {4 m_k \chi _{\alpha k} n_k}{(m_\alpha + m_k)^2} \left [ m_k \left(E_k^{n+1} - \frac {1}{2} \big(\boldsymbol{U}_k^{n+1}\big)^2 \right) - m_\alpha \left(E_\alpha ^{n+1} - \frac {1}{2}\big(\boldsymbol{U}_\alpha ^{n+1}\big)^2 \right) + \frac {m_k}{2} \big(\boldsymbol{U}_k^{n+1} - \boldsymbol{U}_\alpha ^{n+1}\big)^2 \right ]\!.
\end{align}
$$
速度和能量的残差雅可比矩阵为（B2）：
$$
\begin{align}
& \frac {\partial R_{U_\alpha }}{\partial \boldsymbol{U}_k} = \delta _{\alpha k}\left(\frac {1}{\Delta t} + \sum _{k=1}^N \frac {2 m_k \chi _{\alpha k} n_k}{m_\alpha + m_k} \right) - (1-\delta _{\alpha k}) \left(\frac {2 m_k \chi _{\alpha k} n_k}{m_\alpha + m_k}\right)\!,\\
& \frac {\partial R_{E_\alpha }}{\partial E_k} = \delta _{\alpha k}\left(\frac {1}{\Delta t} + \sum _{k=1}^N \frac {4 m_k m_\alpha \chi _{\alpha k} n_k}{(m_\alpha + m_k)^2}\right) - (1-\delta _{\alpha k}) \frac {4 m_k^2 \chi _{\alpha k} n_k}{(m_\alpha + m_k)^2}, \\
&\frac {\partial R_{E_\alpha }}{\partial \boldsymbol{U}_k} = \delta _{\alpha k}\sum _{k=1}^N \frac {4 m_k \chi _{\alpha k} n_k}{(m_\alpha + m_k)^2} \left ( -m_\alpha \boldsymbol{U}_\alpha + m_k (\boldsymbol{U}_k - \boldsymbol{U}_\alpha ) \right )\\
&+\quad +(1-\delta _{\alpha k})\left(\frac {4 m_k^2 \chi _{\alpha k} n_k}{(m_\alpha + m_k)^2} \left ( \boldsymbol{U}_k - m_k (\boldsymbol{U}_k - \boldsymbol{U}_\alpha ) \right )\!, \end{align}
$$
其中（B3）为：
$$
\begin{align}
\delta _{\alpha k} = \left \{ \begin{array}{cc}
1 & \text{如果 } \alpha = k,\\
0 & \text{否则}.
\end{array} \right.
\end{align}
$$
残差表示为（B4）：
$$
\begin{align}
&R_{U_\alpha } = \frac {\boldsymbol{U}_\alpha ^{n+1} - \boldsymbol{U}_\alpha ^{**}}{\Delta t} - \sum _{k=1}^N \frac {2 m_k \chi _{\alpha k} n_k}{m_\alpha + m_k} \big(\boldsymbol{U}_k^{n+1} - \boldsymbol{U}_\alpha ^{n+1}\big),\\
&R_{E_\alpha } = \frac {E_\alpha ^{n+1} - E_\alpha ^{**}}{\Delta t} - \sum _{k=1}^N \frac {2m_{k}\chi _{\alpha k} n_k}{m_{\alpha }+m_{k}} \boldsymbol{U}_{\alpha }^{n+1} \big (\boldsymbol{U}_k^{n+1}-\boldsymbol{U}_\alpha ^{n+1}\big )\\
&+\sum _{k=1}^N \frac {4 m_k \chi _{\alpha k} n_k}{(m_\alpha + m_k)^2} \left [ m_k \left(E_k^{n+1} - \frac {1}{2}\big(\boldsymbol{U}_k^{n+1}\big)^2\right) - m_\alpha \left(E_\alpha ^{n+1} - \frac {1}{2}\big(\boldsymbol{U}_\alpha ^{n+1}\big)^2\right) + \frac {m_k}{2} \big(\boldsymbol{U}_k^{n+1} - \boldsymbol{U}_\alpha ^{n+1}\big)^2 \right ]\!.
\end{align}
$$
使用的迭代方法是牛顿-拉夫森法。当残差的相对变化量低于预定义的容忍度 $e_{\textit{tol}}$ 时，迭代过程达到收敛。具体来说，如果第 $k$ 次迭代时的残差表示为 $R^k$，前一次迭代 $k-1$ 时的残差表示为 $R^{k-1}$，则当 $|R^k - R^{k-1}| \lt e_{\textit{tol}}$ 时，迭代视为收敛。