DSMC方法中边界处粒子初始化的建模：泊松分布的贡献

《Journal of Computational Physics》：Modeling of particle initiation at the boundary in the DSMC method: the Poisson distribution contribution

　　

• 标准DSMC边界条件会导致静止气体产生高达6%的系统性误差。• 当粒子通量在边界上采用泊松分布时，该误差消失。• 研究结果对于DSMC模拟低速流动至关重要。

引言

直接模拟蒙特卡洛（DSMC）方法[1]于20世纪60年代初发展，用于模拟稀薄气体流动。此后数十年间，该方法得到显著改进，降低了计算复杂度，提高了模拟结果的精度[[2], [3], [4], [5], [6], [7], [8], [9]]，并将其应用范围从航天器再入和轨道空气动力学[10,11]扩展到真空技术[[12], [13], [14], [15]]和MEMS器件[[16], [17], [18], [19]]。目前，DSMC方法是求解稀薄气体动力学问题的主要数值工具[[20], [21], [22]]。

为了获得精确解，正确设置边界条件至关重要。最初，该方法是为求解涉及超音速流动的问题而设计的。对于这类流动，通常使用流场边界条件或真空边界条件。对于流场边界条件，假定已知流动密度、速度和温度，并基于这些数据确定进入计算区域的分子通量。真空边界条件用于下游边界，假定无分子进入计算区域。

流场边界条件的建模通常基于面方法或体积方法。在面方法中，利用已知的流动参数确定通过边界表面的粒子通量[[23], [24], [25]]。这是最常见的方法，实际上在Bird的专著[1]中有描述。另一种替代方法基于使用计算区域外的虚拟体积储层[23,[25], [26], [27], [28], [29], [30], [31]]。在该体积储层中，模型粒子被随机分布，然后在时间步长内按照其速度运动。进入计算域的粒子用于模拟，而留在储层单元中的粒子被删除。这种方法通常称为体积生成储层方法[28,29]。该方法还用于高效耦合DSMC和连续区域的混合方法中使用缓冲单元[[32], [33], [34], [35], [36]]，并且可以通过同时使用缓冲单元和储层使其更加复杂[37]。

对于亚音速流动，建模情况更为复杂，并且开发了一系列特殊的基于压力的边界条件[38]，如粒子通量守恒方案[39]、基于特征理论的方案[17,18]和零阶外推方案[16,40,41]。在亚音速流动中，整体流动速度可能远小于热速度，粒子的随机热运动可能显著影响模拟结果，导致大的统计波动。为了克服DSMC方法在低速流动中的这一缺陷，提出了各种方法，如信息保持方法[37,[42], [43], [44]]、低扩散方法[45]、线性化Bhatnagar-Gross-Krook（BGK）近似修正[46]和方差缩减方法[[47], [48], [49]]。

近年来，研究气体涨落对DSMC模拟结果的影响以及DSMC方法调查气体涨落的能力的兴趣日益增长[28,[50], [51], [52], [53], [54], [55]]。研究表明，DSMC的随机性方面——通常被认为是不需要的数值噪声——对于布朗运动和热力学涨落等统计物理问题是一个有用的特性。

在该领域的开创性工作中，Hadjiconstantinou等人[54]利用平衡统计力学的结果，特别是单元中粒子数泊松涨落的假设，开发了估计DSMC计算中气体宏量统计误差的方法论。研究表明，为了正确估计统计误差，必须考虑相关性。还表明，在单元中粒子数呈泊松涨落条件下，离开单元的粒子分布也呈泊松分布。

在该领域的另一项重要工作中，Tysanner和Garcia[28]分析了泊松分布对正确制定边界条件的效果。考虑了面方法和体积方法来向计算域注入粒子。模拟了一个处于平衡状态的气体区域。为了量化动量和密度涨落之间的相关性，作者采用了偏差速度估计器——样本平均方法（SAM）[28]。为了避免使用体积方法时产生的非物理相关性，提出根据具有相应均值值的泊松分布来建模虚拟储层中的粒子数。研究表明，使用标准面方法[1]会导致非物理相关性。然而，采用泊松分布的面方法可以消除这些非物理相关性。此外，需要回顾的是，先前研究[22,24]已指出使用泊松分布来采样每个时间步进入计算区域的随机整数粒子数的必要性。

在更广泛的背景下，偏离泊松统计可以被视为不正确建模引起的非物理行为的指标[28,53]。这些效应通常很微妙（例如，SAM速度与正确值的偏差比声速小三个数量级[28]）。然而，如果物理条件极端或系统演化对微小扰动高度敏感，它们可能变得显著。

本工作的发起源于最近观察到的效应，即需要出乎意料的大量模型粒子来获得低速流动开放式边界的高精度解，即亚音速蒸发到半空间[[56], [57], [58], [59], [60], [61]]，以及我们先前关于DSMC方法误差的研究[8,62,63]。低速流动的极端情况是静止气体。利用一维模型问题中静止气体的例子，讨论了DSMC方法中向计算域引入粒子的各种方法。详细研究了DSMC方法参数（时间步长、模拟粒子数、单元大小）和向计算域引入粒子的方法（有无泊松分布）对计算气体宏量系统性误差的影响。虽然泊松分布对体积储层的重要性已经确立[28]，但标准Bird面通量粒子输入方法[1]尚未针对这类系统性误差进行分析。在本工作中，我们证明广泛使用的标准方法实际上是非物理的，并且在低速流动中会产生显著的偏差。

DSMC方法

采用DSMC方法的传统方案[1]。物理域被划分为计算单元，模拟气体状态在离散时间步中变化。建模过程分为两部分：首先是粒子在空间中的无碰撞运动，其次是粒子间碰撞的建模。每个分子的状态由其空间位置和速度矢量确定。粒子间的碰撞使用"无时间计数器"模拟。

问题描述

考虑了静止气体具有不同边界条件的一维模型问题。指定了一个充满静止气体的区域，其密度为n?，温度为T?。使用硬球模型模拟粒子碰撞过程[1]。区域的长度指定在5λ < L < 500λ范围内，其中λ = 1/(n?σ2)为平均自由程。对于典型计算，设置单元大小Δx = 0.25λ，时间步长Δt = 0.2t?，其中t? = λ/u。

标准方法的计算误差

图2展示了使用标准恒定粒子输入方法（见第3.2节）进行开放式边界计算时，不同模型粒子数5 < Nλ < 200条件下的密度、温度、压力和与x轴平行的速度分量剖面（域大小L = 10λ）。对于模型粒子数Nλ = 10，最大密度值超过精确值1.1%，而温度比精确值低0.07%。同时，注意到一个显著的增加...

标准方法的二维计算

上述结果基于一维计算。了解进行多维计算时结果将如何变化是有趣的。作为示例，我们考虑相同问题（静止气体）的二维解。计算矩形域长度L = 10λ，宽度H可变（图18a）。在计算域的侧面（y = 0和y = H处），采用镜面或漫反射...

泊松方法在亚音速蒸发模拟中的应用

作为一个实际应用，可以考虑慢速蒸发过程中水向半空间蒸发时温度跳跃的确定问题[61]。作为示例，我们将分析文献[61]表1中实验数据N 11，其中液面温度T_S = 306.277 K，p_S = 5071 Pa，气相温度T_∞ = 305.94 K，质量通量J_∞ = 42.731 g/(m2·s)。

该问题在一维表述中求解。问题描述先前已在[56,57]中描述。一个域...

结论

使用DSMC方法求解了一个静止气体的一维问题，以分析向计算域引入粒子的不同方法对所得解精度的影响。本工作的主要结果是，在每个时间步向计算域引入几乎恒定（无波动）粒子数的标准方法，会导致宏量（密度...）的系统性误差。

未引用的参考文献

[64,65]

CRediT作者贡献声明

A.A. Morozov：撰写——审阅与编辑，撰写——原稿，可视化，验证，监督，方法论，调查，形式化分析，数据整理，概念化。M.Yu. Plotnikov：撰写——审阅与编辑，撰写——原稿，可视化，验证，方法论，调查，形式化分析，数据整理。E.Ya. Gatapova：撰写——审阅与编辑，监督，项目管理，资金获取。

利益冲突声明

作者声明，没有已知的竞争性财务利益或个人关系可能影响本文报告的工作。

致谢

本研究由俄罗斯科学基金会支持（项目编号25-19-00959，https://rscf.ru/en/project/25-19-00959/）。

A.A. Morozov|M.Yu. Plotnikov|E.Ya. Gatapova