利用拟随机数对稀薄气体福克 - 普朗克粒子方法进行方差降低

《Mathematics and Computers in Simulation》：Variance reduction in the Fokker–Planck particle method for rarefied gases using quasi-random numbers

　　

# 摘要

Fokker–Planck（FP）粒子方法通过将直接模拟蒙特卡洛（DSMC）算法中常用的二元碰撞替换为漂移–扩散过程，从而加速稀薄气体模拟。与所有粒子方法一样，FP方法本质上是随机的，这会导致宏观量出现统计波动，并需要大量粒子数才能获得准确结果。

本文研究了拟随机数的使用，该方法能更均匀地对分布进行采样，从而降低方差。然而，标准拟蒙特卡洛方法中，每条轨迹的随机输入均由单个高维点生成，这在空间离散化设置中会失去有效性，因为粒子输运和边界相互作用会导致粒子子集不断变化。

为解决这一问题，我们将FP方法与阵列随机拟蒙特卡洛（Array-RQMC）技术相结合，该技术在每个时间步上局部重构低偏差结构。我们将所提出的FP–Array-RQMC方法与已建立的方差缩减技术在均匀和非均匀问题中进行了比较。该方法在粒子数足够大时实现了更优的收敛率和更小的估计器误差。

# 引言

对于稀薄气体，其平均自由程λ远大于特征长度L，即克努森数Kn=λ/L不可忽略。此类气体可能显著偏离热平衡，因此无法由Navier–Stokes–Fourier方程等连续介质模型准确描述[1]。相反，需要一种考虑粒子相互作用的动力学描述。对于稀薄气体，这些相互作用可建模为二元碰撞，从而导出Boltzmann方程[2]。

由于该动力学描述导致高维相空间，其数值近似在计算上极具挑战性[3]。因此，粒子方法直接模拟蒙特卡洛（DSMC）被广泛使用，该方法通过对粒子集合交替进行自由飞行和二元碰撞步骤来推进[4]。然而，在低克努森数区域，碰撞频率很高，使得显式处理二元碰撞在计算上代价高昂。

Fokker–Planck（FP）粒子方法旨在通过将显式二元碰撞替换为漂移–扩散过程来加速模拟。粒子沿随机轨迹演化，且仅通过速度矩耦合，因此计算复杂度不再随碰撞频率变化。漂移和扩散系数的选取方式使得Boltzmann方程在保留其关键物理特性的前提下得到近似[5]、[6]、[7]。

DSMC和FP方法都是本质随机的，因此从有限粒子集合估计的宏观量会表现出统计波动。一个自然的目标是在不增加粒子数的情况下降低这些波动。针对动力学方程和相关蒙特卡洛方法，已考虑了各种方差缩减技术[8]、[9]、[10]。FP粒子方法以随机微分方程形式表述，为对随机输入应用拟随机数提供了自然框架。

与相同样本量的伪随机采样相比，拟随机数通过形成低偏差点集更均匀地采样基础分布，从而提高了矩估计的精度[11]、[12]。在随机过程中，跨时间步保持这种低偏差结构通常通过将每条轨迹的随机输入从单个高维点生成来实现，但在具有空间离散的粒子方法中，由于粒子输运和边界条件导致粒子子集的组成随时间变化，更新是在粒子子集上执行的，这使得该方法变得无效。阵列随机拟蒙特卡洛（Array-RQMC）技术通过在每个时间步上通过粒子状态与拟随机采样的耦合来局部重构低偏差结构，解决了这一问题[13]。Array-RQMC方法的优点已在多个应用中得到验证，包括期权定价[14]和随机化学反应网络[15]。密切相关的策略也已应用于空间均匀的Boltzmann及相关动力学方程[16]、[17]、[18]，以及空间非均匀扩散问题[19]、[20]。然而，据我们所知，Array-RQMC此前尚未应用于因粒子输运和扩散边界相互作用导致在每个单元格的粒子子集上局部执行的基于网格的粒子方法中。

本文工作中，我们将Array-RQMC集成到FP粒子方法中。对于均匀设置，我们专门化了标准Array-RQMC论证，该方法为每个标记粒子产生无偏估计器和正确的路径律。然后我们将路径律论证扩展到具有预设分单元时间相关系数和扩散边界条件的空间离散化设置。我们随后考虑了均匀和非均匀测试案例，并将所提出的方法与伪随机采样和其他已建立的方差缩减技术进行了比较。数值实验表明，该方法在粒子数足够大时实现了更优的收敛率和更小的估计器误差。

第2节回顾了FP方程和相应的粒子算法。第3节介绍了拟随机数和Array-RQMC方法。第4节展示了它们如何集成到FP方法中。第5节报告了数值结果，第6节总结了主要发现。

# 章节摘要

## Fokker–Planck方法

稀薄气体原则上可通过追踪每个独立粒子来建模，但由于其数量极其庞大，这通常不可行。因此，气体通过粒子位置和速度上的分布函数f(x,c,t)=ρ(x,t)f?(c|x,t)进行统计描述，其中ρ为质量密度，f?为位置x∈?3、时间t∈??处的速度c∈?3的条件概率密度函数。

在实际应用中，我们主要关注气体的宏观量。

## 阵列随机拟蒙特卡洛方法

设Y∈??为具有概率密度函数p_Y(y)的随机向量。我们通常对可测函数g:??→?的期望值μ=E[g(Y)]=∫?? g(y)p_Y(y)dy感兴趣。

在许多情况下，由于密度p_Y(y)仅隐式可得或被积函数g(y)p_Y(y)过于复杂而无法以闭式积分，期望值μ无法解析求解。因此，使用蒙特卡洛（MC）方法通过对Y的N次独立实现{Y^(i)}_(i=1)^N来近似μ，得到估计器。

## 将Array-RQMC集成到Fokker–Planck方法中

FP方法通过追踪一组代表粒子{(X^(i),M^(i))}_(i=1)^N的轨迹来模拟稀薄气体，从而提供分布函数f(x,c,t)的离散近似。粒子速度根据Langevin方程（公式(13))演化，粒子位置通过自由飞行在离散时间步中演化（公式(14)）。由于漂移和扩散系数依赖于宏观量，这些量必须在局部评估。

## 数值结果

本节展示了数值实验，以评估FP–Array-RQMC方法与替代采样策略在估计宏观量的误差和收敛性方面的表现。我们首先介绍用于比较的采样策略，然后介绍本节中使用的误差度量。共考虑了四个测试案例。第一个是向预设储层温度弛豫的均匀情况，其中漂移和扩散。

## 结论

我们通过Array-RQMC技术将拟随机数集成到FP方法中，该技术在每个时间步上通过粒子状态和随机更新增量的耦合来局部重构低偏差结构。这一表述使得拟随机数在具有粒子输运和边界相互作用的空间离散化设置中能够得到有效使用。对于预设的分单元时间相关系数和扩散边界条件，我们证明了每个。

# CRediT作者贡献声明

Lukas Netterdon：撰写—原稿，软件，方法，概念化。Veronica Montanaro：撰写—审阅与编辑，方法，概念化。Hossein Gorji：撰写—审阅与编辑，监督，方法，概念化。Manuel Torrilhon：撰写—审阅与编辑，监督，方法，概念化。

# 利益冲突声明

作者声明，他们没有任何已知的竞争性经济利益或个人关系可能影响本文所报告的工作。

# 致谢

本研究由德国研究基金会（DFG）优先项目SPP 2410"流体力学中的双曲平衡定律：复杂性、尺度和随机性（CoScaRa）"资助，项目编号525660607"稀薄多原子气体随机粒子模拟的模型级联"。

Lukas Netterdon | Veronica Montanaro | Hossein Gorji | Manuel Torrilhon