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

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

Info

Publication number
CN111967149A
CN111967149A CN202010765269.4A CN202010765269A CN111967149A CN 111967149 A CN111967149 A CN 111967149A CN 202010765269 A CN202010765269 A CN 202010765269A CN 111967149 A CN111967149 A CN 111967149A
Authority
CN
China
Prior art keywords
particle
solving
interpolation
particles
algorithm
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010765269.4A
Other languages
English (en)
Other versions
CN111967149B (zh
Inventor
蔡文劲
金晓林
黄桃
李斌
杨中海
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN202010765269.4A priority Critical patent/CN111967149B/zh
Publication of CN111967149A publication Critical patent/CN111967149A/zh
Application granted granted Critical
Publication of CN111967149B publication Critical patent/CN111967149B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)

Abstract

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

Description

一种用于粒子模拟算法的粒子运动半插值求解方法
技术领域
本发明属于粒子模拟领域，具体涉及在采用粒子模拟方法模拟物理问题时，一种用于粒子模拟算法的粒子运动半插值求解方法，采用每几个时间步长内有一个时间步长采用插值求解而其余时间步长直接求解的方式。
背景技术
在使用粒子模拟方法对物理问题进行模拟计算时，求解主要流程如图1所示，其中Δt为时间步长，ρg为网格格点上的电荷源密度，
Figure BDA0002614329550000011
为网格格点上的电流源密度，
Figure BDA0002614329550000012
分别为网格格点上的电场值和磁场值，
Figure BDA0002614329550000013
为粒子所受到的力，
Figure BDA0002614329550000014
分别为粒子的运动速度和粒子在空间上的位置。
在每个时间步长内，求解过程可以分为四个步骤：先通过牛顿洛伦兹运动方程求解出当前时刻粒子的位置和运动速度；然后求解空间上每个网格格点处的电荷源密度和电流源密度；再使用麦克斯韦方程组结合已经求出的电流源密度值迭代计算出空间上每个网格的电磁场值；最后利用插值的方式求解出粒子运动时所受到的力；接着继续使用牛顿洛伦兹运动方程求解下一时刻粒子的位置和运动速度，进行下一个时间步长内的求解。如此反复迭代，直至达到预先设置的迭代步数或计算结果收敛为止。
在整个求解过程中，粒子运动求解是最占计算时间的求解部分之一，其包括了求解出粒子运动时所受到的力以及求解粒子的位置和运动速度两个步骤。为保证求解的准确度和提高算法的并行性，最常使用的是权重分配法计算粒子所受到的电磁场值和Boris推动粒子运动算法。Boris算法的主要思路是将粒子运动的求解过程分解为半加速-旋转-半加速三个步骤，实现方式如下：
1、使用网格点权重分配的方式求解粒子当前位置的电磁场值；
2、使用以下公式更新半个时间步长内的粒子速度
Figure BDA0002614329550000015
Figure BDA0002614329550000016
其中c为光速，
Figure BDA0002614329550000017
为当前时刻的粒子速度，γn为相对论因子，u为相对论速度，Δt为时间步长，e为粒子电荷量，me为粒子质量，E为粒子所受到的电场值。
3、进行旋转求解。以三维为例使用以下公式更新u
Figure BDA0002614329550000021
Figure BDA0002614329550000022
Figure BDA0002614329550000023
Figure BDA0002614329550000024
Figure BDA0002614329550000025
Figure BDA0002614329550000026
Figure BDA0002614329550000027
Figure BDA0002614329550000028
其中Bx、By、Bz为粒子所受到的磁场分量。
4、使用下面公式再次半加速迭代求解出最终粒子的速度u
Figure BDA0002614329550000029
5、根据速度求解粒子新的位置。
由上述算法可以看出：在计算粒子所受到的力时，网格点权重分配算法需要读取存储网格点电磁场值的内存，而粒子模拟庞大的网格数目和粒子数目会使得需要访问的内存数目剧增；同时Boris推动粒子运动算法需要计算许多的中间变量，庞大的粒子数目会导致中间变量的计算次数更加频繁，进而增加计算时间。综合以上原因，在粒子模拟求解过程中常规的粒子运动的求解计算负担大、计算效率低下。
发明内容
针对上述存在问题或不足，为解决粒子模拟算法在上述粒子运动求解过程中效率低的问题，本发明提供了一种用于粒子模拟算法的粒子运动半插值求解方法。
一种用于粒子模拟算法的粒子运动半插值求解方法，具体技术步骤如下：
步骤1、设粒子在单位时间步长内的运动最多可以跨越1个网格；开辟3个数组VxNt、VyNt、VzNt来存储当前时刻粒子x、y、z三个坐标方向上的速度分量，6个数组VxNt-1、VyNt-1、VzNt-1、VxNt-2、VyNt-2和VzNt-2分别对应存储粒子在前两个时刻三个坐标方向上的速度分量，将其数值均初始化为粒子初始时刻的值。
Figure BDA0002614329550000031
和
Figure BDA0002614329550000032
分别代表当前时刻和前两个时刻粒子的运动速度。
步骤2、先计算Nt％N的值，然后进行判定；Nt表示当前时刻为第Nt个时间步长，N表示每N个时间步长进行一次粒子运动半插值求解，N≥2，％为取余函数。
若Nt％N的值不为0，使用网格点权重分配的方式求解出粒子运动时所受到的力，接着使用Boris算法求解下一时刻粒子的运动速度。
若Nt％N＝0，采用粒子运动半插值求解方法求解。使用前两个时刻和当前时刻的速度值来插值计算出下一时刻的粒子的速度值。
进一步的，所述插值方式为牛顿插值公式、拉格朗日插值公式或/和其它插值公式。
步骤3、利用步骤2求解出的下一时刻速度值更新该粒子的位置，并以此时间步长循环更新直至预设时间。
进一步的，本发明适用于一维、二维及三维粒子模拟算法。
本发明通过采用每几个时间步长内有一个时间步长采用插值求解而其余时间步长直接求解的方式，避免了常规粒子运动求解算法中庞大的网格数目和粒子数目导致的需要访问的内存数目过多和中间变量的计算次数过于频繁等问题，降低了粒子运动求解部分的计算时间，提高了整体的计算效率。本发明求解方法适合于粒子模拟方法的电磁、静电和静磁模型，可与其它步骤耦合在一起，形成粒子模拟算法完整的求解流程。
附图说明
图1为现有粒子模拟算法流程示意图；
图2为本发明粒子模拟算法流程示意图；
图3为实施例20周期的折叠波导测试实例立体示意图；
图4为实施例20周期的折叠波导测试实例截面示意图；
图5为采用粒子运动半插值计算和不采用该方法计算测试实例的输出信号对比示意图；
图6为输出信号对比局部放大示意图；
图7为采用粒子运动半插值计算和不采用该方法计算测试实例所需的时间对比示意图。
具体实施方式
下面结合附图和实施例对本发明作进一步详细说明。
本实施例三维粒子模拟算法为例，测试了一个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个时间步长进行一次粒子运动半插值求解；
如图2所示：
步骤1、开辟9个大小为np×nt＝22×11507＝253154的数组VxNt-1、VyNt-1、VzNt-1、VxNt-2、VyNt-2、VzNt-2、VxNt、VyNt和VzNt，前6个数组存储粒子在前两个时刻的三个速度分量，后3个数组存储当前时刻粒子三个速度分量，将其数值都初始化为粒子初始时刻的值。
步骤2、假设当前时刻为第N个时间步长，先计算N％2的值，其中％为取余函数。
若Nt％N的值不为0，使用网格点权重分配的方式求解出粒子运动时所受到的力和使用Boris推动粒子运动算法更新粒子运动速度。
若Nt％N＝0，采用粒子运动半插值求解方法求解，这里以牛顿插值公式为例更新粒子运动速度：
VxNt+1＝3VxNt-3VxNt-1+VxNt-2 (12)
VyNt+1＝3VyNt-3VyNt-1+VyNt-2 (13)
VzNt+1＝3VzNt-3VzNt-1+VzNt-2 (14)
其中VxNt+1、VyNt+1、VzNt+1为下一时刻的粒子速度，VxNt、VyNt、VzNt为当前时刻的粒子速度，VxNt-1、VyNt-1、VzNt-1、VxNt-2、VyNt-2和VzNt-2分别为粒子在前两个时刻的速度。
步骤3、利用步骤2求解出的速度值更新粒子的位置。
将上述步骤与电磁模型的其它流程耦合，从n＝0迭代直至达到n＝11507。
图5和图6显示了该测试实施例采用粒子运动半插值计算和全部时间步长内进行完整求解两种方法下求得的输出信号的电压幅值，可以看出二者完全吻合；图7则展示了采用本发明粒子运动半插值算法和常规粒子运动算法的整体粒子模拟求解所需的计算时间，以及粒子运动部分总的求解时间，可以看出相较于全部时间步长内进行完整求解的方式，粒子运动半插值求解方法无论是整体求解时间还是单独粒子运动部分的求解时间都更小。
综上可见，本发明通过使用粒子运动半插值求解方法，避免了全部时间步长内进行完整的粒子运动求解，提高了粒子模拟方法在模拟物理问题时的计算效率。

Claims (3)

1.一种用于粒子模拟算法的粒子运动半插值求解方法，其特征在于，包括以下步骤：
步骤1、设粒子在单位时间步长内的运动最多可以跨越1个网格；开辟3个数组VxNt、VyNt、VzNt来存储当前时刻粒子x、y、z三个坐标方向上的速度分量，6个数组VxNt-1、VyNt-1、VzNt-1、VxNt-2、VyNt-2和VzNt-2分别对应存储粒子在前两个时刻三个坐标方向上的速度分量，将其数值均初始化为粒子初始时刻的值；
Figure FDA0002614329540000011
和
Figure FDA0002614329540000012
分别代表当前时刻和前两个时刻粒子的运动速度。
步骤2、先计算Nt％N的值，然后进行判定；Nt表示当前时刻为第Nt个时间步长，N表示每N个时间步长进行一次粒子运动半插值求解，N≥2，％为取余函数；
若Nt％N的值不为0，使用网格点权重分配的方式求解出粒子运动时所受到的力，接着使用Boris算法求解下一时刻粒子的运动速度；
若Nt％N＝0，采用粒子运动半插值求解方法求解，使用前两个时刻和当前时刻的速度值来插值计算出下一时刻的粒子的速度值。
步骤3、利用步骤2求解出的下一时刻速度值更新该粒子的位置，并以此时间步长循环更新直至预设时间。
2.如权利要求1所述用于粒子模拟算法的粒子运动半插值求解方法，其特征在于：所述插值方式为牛顿插值公式或/和拉格朗日插值公式。
3.如权利要求1所述用于粒子模拟算法的粒子运动半插值求解方法，其特征在于：适用于一维、二维及三维粒子模拟算法。
CN202010765269.4A 2020-08-03 2020-08-03 一种用于粒子模拟算法的粒子运动半插值求解方法 Active CN111967149B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010765269.4A CN111967149B (zh) 2020-08-03 2020-08-03 一种用于粒子模拟算法的粒子运动半插值求解方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010765269.4A CN111967149B (zh) 2020-08-03 2020-08-03 一种用于粒子模拟算法的粒子运动半插值求解方法

Publications (2)

Publication Number Publication Date
CN111967149A true CN111967149A (zh) 2020-11-20
CN111967149B CN111967149B (zh) 2022-11-04

Family

ID=73363615

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010765269.4A Active CN111967149B (zh) 2020-08-03 2020-08-03 一种用于粒子模拟算法的粒子运动半插值求解方法

Country Status (1)

Country Link
CN (1) CN111967149B (zh)

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2000348014A (ja) * 1999-06-01 2000-12-15 Nec Corp 仮想粒子の追跡計算方法、プログラムを記録した記録媒体および計算機
CN107633123A (zh) * 2017-09-13 2018-01-26 浙江工业大学 一种用于光滑粒子流体动力学模拟出血及处理加速的方法
CN107908913A (zh) * 2017-12-22 2018-04-13 中国海洋大学 基于并行计算机的地球动力数模算法
CN108280309A (zh) * 2018-02-05 2018-07-13 电子科技大学 一种应用于pic静电模型的电位有限元求解算法
CN108416107A (zh) * 2018-02-05 2018-08-17 电子科技大学 一种应用于pic的推动粒子运动有限元算法
CN110442919A (zh) * 2019-07-12 2019-11-12 西安空间无线电技术研究所 一种基于gpu架构的微波部件微放电数值模拟方法
CN111241742A (zh) * 2019-12-27 2020-06-05 西安交通大学 一种多相流计算方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2000348014A (ja) * 1999-06-01 2000-12-15 Nec Corp 仮想粒子の追跡計算方法、プログラムを記録した記録媒体および計算機
CN107633123A (zh) * 2017-09-13 2018-01-26 浙江工业大学 一种用于光滑粒子流体动力学模拟出血及处理加速的方法
CN107908913A (zh) * 2017-12-22 2018-04-13 中国海洋大学 基于并行计算机的地球动力数模算法
CN108280309A (zh) * 2018-02-05 2018-07-13 电子科技大学 一种应用于pic静电模型的电位有限元求解算法
CN108416107A (zh) * 2018-02-05 2018-08-17 电子科技大学 一种应用于pic的推动粒子运动有限元算法
CN110442919A (zh) * 2019-07-12 2019-11-12 西安空间无线电技术研究所 一种基于gpu架构的微波部件微放电数值模拟方法
CN111241742A (zh) * 2019-12-27 2020-06-05 西安交通大学 一种多相流计算方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
F. MACKAY: "Divergence-free magnetic field interpolation and charged", 《HTTPS://AGUPUBS.ONLINELIBRARY.WILEY.COM/DOI/PDF/10.1029/2005JA011382》 *
李斌: "改进Kriging方法在螺旋桨测量中的应用", 《中国机械工程》 *
郑飞腾: "空心阴极虚火花放电初始电离过程的PIC/MCC模拟研究", 《中国优秀博硕士学位论文全文数据库(硕士)基础科学辑》 *

Also Published As

Publication number Publication date
CN111967149B (zh) 2022-11-04

Similar Documents

Publication Publication Date Title
US7191161B1 (en) Method for constructing composite response surfaces by combining neural networks with polynominal interpolation or estimation techniques
US8204925B2 (en) Controlling or analyzing a process by solving a system of linear equations in real-time
JP7133894B2 (ja) データに基づくインタラクティブ３ｄエクスペリエンス
JP2008541203A (ja) 相互作用粒子系のための固定時間ステップダイナミックスソルバー
CN104881510A (zh) 一种直升机旋翼/尾桨气动干扰数值仿真方法
CN104679958A (zh) 基于弹簧模型的球b样条编针织物形变仿真的方法
CN110414053B (zh) 一种快速确定部件微放电阈值的时域数值模拟方法
CN108416107B (zh) 一种应用于pic的推动粒子运动有限元算法
CN110610050A (zh) 基于改进径向基函数变形算法的翼型气动减阻方法
CN115618498B (zh) 一种飞行器跨流域流场的预测方法、装置、设备及介质
CN111709150A (zh) 任意形状磁体磁场空间分布仿真模拟方法
Yu et al. SPACE code for beam-plasma interaction
CN111967149B (zh) 一种用于粒子模拟算法的粒子运动半插值求解方法
CN113821915B (zh) 一种轴对称电子光学系统的快速计算方法
CN115688212A (zh) 一种基于物质点法的软体机器人仿真方法
CN113325350B (zh) 一种基于离散网格的高性能梯度线圈设计方法
CN113987898B (zh) 一种放电室中气体放电的粒子模拟加速方法
CN109118561B (zh) 一种基于位置的层次化动态模拟方法
CN105550424A (zh) 一种基于rbf网格变形插值序列的筛选方法
CN108268697A (zh) 一种高效率电推进羽流等离子体并行仿真方法
Coco et al. 3-D finite-element analysis of TWT grid electron guns
Stueber Discharge chamber primary electron modeling activities in 3-dimension
Yiyu et al. A FPGA implementation of the two-dimensional Digital Huygens' Model
CN116720409B (zh) 一种运动时变色散介质目标的电磁散射场计算方法
Gibbons et al. Flexible three-dimensional modeling of electric thrusters in vacuum chambers

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant