保持局域高斯定律的本征PIC方法：一种用于等离子体模拟的自洽场粒子更新方案

《Computer Physics Communications》：Intrinsic local Gauss’s law preserving PIC method: A self-consistent field-particle update scheme for plasma simulations

　　赵中华|徐振利|殷倩|周盛高
中国香港香港理工大学应用数学系

摘要
高斯定律是实现物理上精确的粒子在单元格模拟中的关键要求，因为违反该定律往往会在长期模拟中导致严重错误。本文提出了一种新方法，能够在无需辅助场修正或特殊电流注入技术的情况下，为Vlasov-Ampère系统内在地满足高斯定律。在该方法中，通过将粒子运动沿计算网格的每个维度分解为多个子步骤，从而实现对电场的局部更新，并使其与粒子运动保持一致。为进一步获得无旋电场，还开发了一种局部更新方案，在满足高斯定律的前提下最小化电场自由能。该方法并非直接在网格上求解泊松方程或安培方程，而是以拉格朗日-粒子形式隐式满足安培定律。由此产生的局部算法，其计算成本与每次松弛迭代中的局部网格操作次数成线性关系，可灵活与粒子运动的各类时间离散化方法结合使用。理论分析表明，该方法确实能够精确保持离散形式的高斯定律。通过对经典基准问题的数值测试，包括朗道阻尼、双流不稳定性以及迪奥克特龙不稳定性，验证了该方法的优势。此外，由于该方法的局部特性，有望成为大规模等离子体近邻通信并行模拟中的有力工具。

引言
Vlasov方程是一种描述带电粒子在电磁场作用下演化过程的动力学方程[1]。它与玻尔兹曼方程[2]略有不同，后者忽略了短程碰撞，而更强调由自洽场主导的相互作用。在完整的电磁学框架中，Vlasov方程与完整的麦克斯韦方程组相耦合，其中电场和磁场可通过粒子运动产生的电荷密度和电流来求解。在静电近似下，这一框架可简化为两种等效形式[3][4]：即通过泊松方程计算电场的Vlasov-Poisson系统，以及基于局部电荷守恒方程从安培方程推导电场的Vlasov-Ampère系统，这两种形式均忽略了磁效应。由于只有对于理想化系统（如线性波）才存在解析解，因此需要借助数值方法来研究由Vlasov方程描述的现象。总体而言，这些数值方法可分为两大类：欧拉方法，即直接求解Vlasov方程；以及基于粒子的方法，通过离散的宏观粒子来近似描述粒子分布。欧拉方法包括间断伽辽金法[5]、有限体积法[6]、半拉格朗日法[7][8]、矩量法[9]以及其他方法[10][11]，适用于高精度的相空间分析；而属于基于粒子方法的粒子在单元格模拟[12][13][14][15]，则因其计算可扩展性而在高维、大规模模拟中占据主导地位。有一些研究通过相空间重映射来降低粒子在单元格模拟中的数值噪声，以提高精度[16][17]。此外，还有许多尝试致力于证明将Vlasov方程与泊松方程耦合时的粒子方法的收敛性[18][19]。

在粒子在单元格模拟中保持关键的数学和物理性质，对于确保数值结果的准确性与可靠性至关重要，尤其是在长期模拟中。例如，那些能够保持跨尺度一致性的渐近保持格式[20][21]，已被证明在德拜长度接近准中性极限时解决多尺度耦合问题方面不可或缺。此外，在连续系统中，总能量应理论上保持守恒，它由电场能量和动能共同构成。但有限网格的不稳定性可能引发非物理性的加热或冷却现象，这促使人们大力发展能量守恒格式[22][23][24][25][26][27]。同时，动量守恒方法[28][29]也同样重要。在各种守恒要求中，质量守恒和高斯定律是实现稳健粒子在单元格模拟的基础，相关研究已对此进行了深入探讨。正如第2.3节将展示的，局部质量守恒的设计本身有助于维持高斯定律。Anderson等人[30]指出，未能遵循高斯定律会在粒子在单元格模拟中带来灾难性后果。

本研究提出了一种内在的、局部保持高斯定律的粒子在单元格模拟框架（GP-PIC），该框架可与其他守恒机制相结合，从而为Vlasov-Ampère系统实现全面的守恒性质。第2.3节对高斯定律保持方面的最新进展进行了系统梳理，这些进展可分为两类：一类是场修正方法[31][32][33][34]，通过对电场进行后处理来消除高斯定律的不一致性；另一类是电荷守恒方法[35][36][37][38][39]，通过精心设计的电流密度来强制满足离散的连续性方程。在电荷守恒方法中，符合旋度条件的FETD-PIC格式[40][41]通过在非结构化网格上使用Whitney基函数来满足高斯定律。相比之下，我们的GP-PIC则通过具有线性复杂度的局部显式松弛步骤来实现相同目标，它在时间离散化方面更具灵活性，也能更便捷地整合其他守恒定律。本研究提出的这种保持高斯定律的框架与现有传统方法有本质区别。该算法通过统一的场-粒子推进方案，同时更新电场和粒子位置，同时严格保持离散形式的高斯定律，精度可达机器精度水平。具体而言，首先根据“之字形”方案[38]的思想，将粒子运动分解为平行于坐标轴的多个子步骤。在每个方向上的子步骤中，根据该轴上粒子运动所引起的电荷通量，对电场进行局部更新，以满足离散高斯定律。随后，通过局部无旋松弛机制将这些局部电场更新结果传播到整个区域，从而最小化电场自由能。与最初为全电磁粒子在单元格模拟开发的“之字形”方案不同，我们的方法专为静电（Vlasov-Ampère）场景设计，通过自由能最小化机制内在地同时满足高斯定律和无旋条件。整个过程通过其对场传播和粒子动力学的自洽处理，避免了对外部泊松或安培求解器的依赖。这类局部算法最早被提出用于蒙特卡罗模拟[42][43][44]，以及静电系统的分子动力学模拟[45][46][47][48]。后来，它又被推广用于求解连续静电模型，如泊松-玻尔兹曼方程[49][50][51]和泊松-能斯特-普朗克方程[52][53]。此外，这些研究还证明了局部无旋松弛算法具有线性复杂度。其在保持高斯定律方面的独特优势以及内在的局部特性，使得该方法在需要同时保持多种性质（质量、动量、能量）的等离子体模拟中具有显著优势。

本文的其余结构如下：第2节概述了Vlasov型模型和粒子在单元格模拟方法，并回顾了那些用于保持高斯定律的电场处理的相关工作。第3节详细介绍了所提出的局部保持高斯定律的粒子在单元格模拟方法。第4节进行了若干数值实验。最后，第5节总结了研究结论。

片段摘录
Vlasov-Ampère方程
常用于描述惯性约束聚变、天体物理学等领域 的Vlasov方程如下：
?fs(x,v,t)?t+v·?xfs(x,v,t)+F(x,t)·?vfs(x,v,t)=0
其中，fs(x, v, t)表示在特定位置x、速度v和时间t处，种类为s的粒子分布。F(x, t)为外部力。通常，力场包括静电力和洛伦兹力，即F(x,t)=qs(E+v∧B)/ms，其中qs为粒子所带电荷，ms为质量，E为电场强度，B为磁场强度。

离散化
在粒子在单元格模拟中，对于方程(2.9)和(2.10)的经典时间离散化方法称为“跳跃法”，即位置在时间步长的整数倍处定义，即xp(mΔt)；而速度则在半整数倍处定义：vpm+1/2=vp((m+1/2)Δt)。它们的更新方式如下：
xm+1=xm+Δtvm+12,
vm+32=vm+12?ΔtEm+1(xp)，
其中Em+1(xp)=E(xpm+1)是利用网格数据插值得到的电场值。在半时间步长处的初始速度为……

数值结果
为验证所提出的GP-PIC方法的有效性，我们进行了系列2D2V数值实验，涉及朗道阻尼、双流不稳定性以及迪奥克特龙不稳定性等场景。在这些模拟中，计算域均为Ωx=[0,L]×[0,L]，并采用周期性边界条件。为满足周期性边界条件所要求的电中性条件，背景中均匀分布着不可移动的离子，即仅有……

结论
在等离子体模拟中严格遵循高斯定律的重要性不言而喻。本研究提出了一种用于静电极限下粒子在单元格模拟中电场更新的高斯定律保持方法。该方法通过将粒子运动沿计算网格的每个维度分解为多个子步骤，根据粒子运动所引起的电通量来局部更新电场，从而实现能够满足高斯定律的模拟方法。

CRediT作者贡献说明
赵中华：写作——审阅与编辑、指导、方法论、概念设计。
徐振利：写作——审阅与编辑、指导、方法论、概念设计。
殷倩：写作——初稿撰写、可视化、软件应用、方法论、研究实施、数据整理、概念设计。
周盛高：写作——审阅与编辑、验证、指导、方法论、概念设计。

利益冲突声明
作者声明不存在任何可能影响本文研究的已知财务利益或个人关系。