一种用于粒子模拟算法的粒子运动半插值求解方法

摘要

本发明属于粒子模拟领域，具体为一种用于粒子模拟算法的粒子运动半插值求解方法。本发明采用每几个时间步长内有一个时间步长采用插值求解而其余时间步长直接求解的方式，避免了常规粒子运动求解算法中庞大的网格数目和粒子数目导致的需要访问的内存数目过多和中间变量的计算次数过于频繁等问题，降低了粒子运动求解部分的计算时间，提高了整体的计算效率。本发明求解方法适合于粒子模拟方法的电磁、静电和静磁模型，可与其它步骤耦合在一起，形成粒子模拟算法完整的求解流程。

权利要求

1.一种用于粒子模拟算法的粒子运动半插值求解方法，其特征在于，包括以下步骤：步骤1、设粒子在单位时间步长内的运动最多可以跨越1个网格；开辟3个数组VxNt、VyNt、VzNt来存储当前时刻粒子x、y、z三个坐标方向上的速度分量，6个数组VxNt‑1、VyNt‑1、VzNt‑1、VxNt‑2、VyNt‑2和VzNt‑2分别对应存储粒子在前两个时刻三个坐标方向上的速度分量，将其数值均初始化为粒子初始时刻的值； 和 分别代表当前时刻和前两个时刻粒子的运动速度；

步骤2、先计算Nt％N的值，然后进行判定；Nt表示当前时刻为第Nt个时间步长，N表示每N个时间步长进行一次粒子运动半插值求解，N≥2，％为取余函数；

若Nt％N的值不为0，使用网格点权重分配的方式求解出粒子运动时所受到的力，接着使用Boris算法求解下一时刻粒子的运动速度；所述网格点权重分配的方式指在计算粒子所受到的力时，网格点权重分配算法需要读取存储网格点电磁场值的内存；

若Nt％N＝0，采用粒子运动半插值求解方法求解，使用前两个时刻和当前时刻的速度值来插值计算出下一时刻的粒子的速度值；

步骤3、利用步骤2求解出的下一时刻速度值更新该粒子的位置，并以此时间步长循环更新直至预设时间。

2.如权利要求1所述用于粒子模拟算法的粒子运动半插值求解方法，其特征在于：所述插值方式为牛顿插值公式或/和拉格朗日插值公式。

3.如权利要求1所述用于粒子模拟算法的粒子运动半插值求解方法，其特征在于：适用于一维、二维及三维粒子模拟算法。

说明书

一种用于粒子模拟算法的粒子运动半插值求解方法

技术领域

[0001] 本发明属于粒子模拟领域，具体涉及在采用粒子模拟方法模拟物理问题时，一种用于粒子模拟算法的粒子运动半插值求解方法，采用每几个时间步长内有一个时间步长采用插值求解而其余时间步长直接求解的方式。

背景技术

[0002] 在使用粒子模拟方法对物理问题进行模拟计算时，求解主要流程如图1所示，其中Δt为时间步长，ρg为网格格点上的电荷源密度， 为网格格点上的电流源密度， 分别为网格格点上的电场值和磁场值， 为粒子所受到的力， 分别为粒子的运动速度和粒子在空间上的位置。

[0003] 在每个时间步长内，求解过程可以分为四个步骤：先通过牛顿洛伦兹运动方程求解出当前时刻粒子的位置和运动速度；然后求解空间上每个网格格点处的电荷源密度和电流源密度；再使用麦克斯韦方程组结合已经求出的电流源密度值迭代计算出空间上每个网格的电磁场值；最后利用插值的方式求解出粒子运动时所受到的力；接着继续使用牛顿洛伦兹运动方程求解下一时刻粒子的位置和运动速度，进行下一个时间步长内的求解。如此反复迭代，直至达到预先设置的迭代步数或计算结果收敛为止。

[0004] 在整个求解过程中，粒子运动求解是最占计算时间的求解部分之一，其包括了求解出粒子运动时所受到的力以及求解粒子的位置和运动速度两个步骤。为保证求解的准确度和提高算法的并行性，最常使用的是权重分配法计算粒子所受到的电磁场值和Boris推动粒子运动算法。Boris算法的主要思路是将粒子运动的求解过程分解为半加速‑旋转‑半加速三个步骤，实现方式如下：

[0005] 1、使用网格点权重分配的方式求解粒子当前位置的电磁场值；

[0006] 2、使用以下公式更新半个时间步长内的粒子速度

[0007]

[0008]

[0009] 其中c为光速， 为当前时刻的粒子速度，γn为相对论因子，u为相对论速度，Δt为时间步长，e为粒子电荷量，me为粒子质量，E为粒子所受到的电场值。

[0010] 3、进行旋转求解。以三维为例使用以下公式更新u

[0011]

[0012]

[0013]

[0014]

[0015]

[0016]

[0017]

[0018]

[0019] 其中Bx、By、Bz为粒子所受到的磁场分量。

[0020] 4、使用下面公式再次半加速迭代求解出最终粒子的速度u

[0021]

[0022] 5、根据速度求解粒子新的位置。

[0023] 由上述算法可以看出：在计算粒子所受到的力时，网格点权重分配算法需要读取存储网格点电磁场值的内存，而粒子模拟庞大的网格数目和粒子数目会使得需要访问的内存数目剧增；同时Boris推动粒子运动算法需要计算许多的中间变量，庞大的粒子数目会导致中间变量的计算次数更加频繁，进而增加计算时间。综合以上原因，在粒子模拟求解过程中常规的粒子运动的求解计算负担大、计算效率低下。

发明内容

[0024] 针对上述存在问题或不足，为解决粒子模拟算法在上述粒子运动求解过程中效率低的问题，本发明提供了一种用于粒子模拟算法的粒子运动半插值求解方法。

[0025] 一种用于粒子模拟算法的粒子运动半插值求解方法，具体技术步骤如下：

[0026] 步骤1、设粒子在单位时间步长内的运动最多可以跨越1个网格；开辟3个数组VxNt、VyNt、VzNt来存储当前时刻粒子x、y、z三个坐标方向上的速度分量，6个数组VxNt‑1、VyNt‑1、VzNt‑1、VxNt‑2、VyNt‑2和VzNt‑2分别对应存储粒子在前两个时刻三个坐标方向上的速度分量，将其数值均初始化为粒子初始时刻的值。 和 分别代表当前时刻和前两个时刻粒子的运动速度。

[0027] 步骤2、先计算Nt％N的值，然后进行判定；Nt表示当前时刻为第Nt个时间步长，N表示每N个时间步长进行一次粒子运动半插值求解，N≥2，％为取余函数。

[0028] 若Nt％N的值不为0，使用网格点权重分配的方式求解出粒子运动时所受到的力，接着使用Boris算法求解下一时刻粒子的运动速度。

[0029] 若Nt％N＝0，采用粒子运动半插值求解方法求解。使用前两个时刻和当前时刻的速度值来插值计算出下一时刻的粒子的速度值。

[0030] 进一步的，所述插值方式为牛顿插值公式、拉格朗日插值公式或/和其它插值公式。

[0031] 步骤3、利用步骤2求解出的下一时刻速度值更新该粒子的位置，并以此时间步长循环更新直至预设时间。

[0032] 进一步的，本发明适用于一维、二维及三维粒子模拟算法。

[0033] 本发明通过采用每几个时间步长内有一个时间步长采用插值求解而其余时间步长直接求解的方式，避免了常规粒子运动求解算法中庞大的网格数目和粒子数目导致的需要访问的内存数目过多和中间变量的计算次数过于频繁等问题，降低了粒子运动求解部分的计算时间，提高了整体的计算效率。本发明求解方法适合于粒子模拟方法的电磁、静电和静磁模型，可与其它步骤耦合在一起，形成粒子模拟算法完整的求解流程。

附图说明

[0034] 图1为现有粒子模拟算法流程示意图；

[0035] 图2为本发明粒子模拟算法流程示意图；

[0036] 图3为实施例20周期的折叠波导测试实例立体示意图；

[0037] 图4为实施例20周期的折叠波导测试实例截面示意图；

[0038] 图5为采用粒子运动半插值计算和不采用该方法计算测试实例的输出信号对比示意图；

[0039] 图6为输出信号对比局部放大示意图；

[0040] 图7为采用粒子运动半插值计算和不采用该方法计算测试实例所需的时间对比示意图。

具体实施方式

[0041] 下面结合附图和实施例对本发明作进一步详细说明。

[0042] 本实施例三维粒子模拟算法为例，测试了一个20周期的折叠波导模型，如图3、图4所示。单周期长度p＝0.8mm，模拟参数为：空间步长Δx＝Δz＝0.05mm，Δy＝0.045mm，Δt＝9.04340828e‑5ns，模拟时间步数nt＝11507，输入信号频率f＝139GHz，平均功率为0.5W，采用0.5mm*0.09mm的带状电子注，每个时间步长内发射电子数目np＝22，发射电压V＝15750，发射电流I＝0.12A。每2个时间步长进行一次粒子运动半插值求解；

[0043] 如图2所示：

[0044] 步骤1、开辟9个大小为np×nt＝22×11507＝253154的数组VxNt‑1、VyNt‑1、VzNt‑1、VxNt‑2、VyNt‑2、VzNt‑2、VxNt、VyNt和VzNt，前6个数组存储粒子在前两个时刻的三个速度分量，后3个数组存储当前时刻粒子三个速度分量，将其数值都初始化为粒子初始时刻的值。

[0045] 步骤2、假设当前时刻为第N个时间步长，先计算N％2的值，其中％为取余函数。

[0046] 若Nt％N的值不为0，使用网格点权重分配的方式求解出粒子运动时所受到的力和使用Boris推动粒子运动算法更新粒子运动速度。

[0047] 若Nt％N＝0，采用粒子运动半插值求解方法求解，这里以牛顿插值公式为例更新粒子运动速度：

[0048] VxNt+1＝3VxNt‑3VxNt‑1+VxNt‑2   (12)

[0049] VyNt+1＝3VyNt‑3VyNt‑1+VyNt‑2   (13)

[0050] VzNt+1＝3VzNt‑3VzNt‑1+VzNt‑2   (14)

[0051] 其中VxNt+1、VyNt+1、VzNt+1为下一时刻的粒子速度，VxNt、VyNt、VzNt为当前时刻的粒子速度，VxNt‑1、VyNt‑1、VzNt‑1、VxNt‑2、VyNt‑2和VzNt‑2分别为粒子在前两个时刻的速度。

[0052] 步骤3、利用步骤2求解出的速度值更新粒子的位置。

[0053] 将上述步骤与电磁模型的其它流程耦合，从n＝0迭代直至达到n＝11507。

[0054] 图5和图6显示了该测试实施例采用粒子运动半插值计算和全部时间步长内进行完整求解两种方法下求得的输出信号的电压幅值，可以看出二者完全吻合；图7则展示了采用本发明粒子运动半插值算法和常规粒子运动算法的整体粒子模拟求解所需的计算时间，以及粒子运动部分总的求解时间，可以看出相较于全部时间步长内进行完整求解的方式，粒子运动半插值求解方法无论是整体求解时间还是单独粒子运动部分的求解时间都更小。

[0055] 综上可见，本发明通过使用粒子运动半插值求解方法，避免了全部时间步长内进行完整的粒子运动求解，提高了粒子模拟方法在模拟物理问题时的计算效率。