CN108446429A - 一种应用于pic静电模型的粒子受力有限元求解算法 - Google Patents

一种应用于pic静电模型的粒子受力有限元求解算法 Download PDF

Info

Publication number
CN108446429A
CN108446429A CN201810112819.5A CN201810112819A CN108446429A CN 108446429 A CN108446429 A CN 108446429A CN 201810112819 A CN201810112819 A CN 201810112819A CN 108446429 A CN108446429 A CN 108446429A
Authority
CN
China
Prior art keywords
particle
grid
pic
electric charge
formula
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
CN201810112819.5A
Other languages
English (en)
Other versions
CN108446429B (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 CN201810112819.5A priority Critical patent/CN108446429B/zh
Publication of CN108446429A publication Critical patent/CN108446429A/zh
Application granted granted Critical
Publication of CN108446429B publication Critical patent/CN108446429B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

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)

Abstract

本发明属于粒子模拟PIC的数值模拟领域，具体涉及一种应用于PIC静电模型的粒子受力有限元求解算法。本发明使用完全非结构化网格，该网格能够更好的拟合模型边界的形状，使得在复杂边界情况下PIC静电模型的粒子受力求解具有更高的计算精度；将用于求解无源电磁场分布、热分析、力学分析等无粒子源问题的FEM方法结合到典型的PIC方法中，在保持典型PIC方法的计算简单、快速的优良特性的同时，利用FEM得到更高的有限元计算精度；由于FEM方法既可以很好地匹配复杂边界，又可以根据模拟需要使用非均匀网格，且不受数值稳定性条件的限制，因此可以在保持计算精度的条件下，优化空间网格和时间步长，从而大幅地提高模拟效率。

Description

一种应用于PIC静电模型的粒子受力有限元求解算法
技术领域
本发明属于粒子模拟(Particle-in-cell，简写为PIC)的数值模拟领域，具体涉及一种应用于PIC静电模型的粒子受力有限元求解算法。
背景技术
PIC方法是一种被广泛应用于带电粒子与电磁场相互作用物理问题中的数值模拟方法，它通过跟踪大量带电粒子在外加及自洽电磁场中的运动并统计平均而得到宏观特性及运动规律。经过几十年的发展，PIC模拟方法已经成为研究带电粒子与电磁场相互作用物理问题的一种强有力的数值手段，广泛应用于带电粒子与电磁场相互作用所涉及的许多领域，如磁约束聚变等离子体、惯性约束聚变等离子体、核爆、空间等离子体、人造等离子体(包括电子枪、离子源等)、电推进、自由电子激光以及电真空器件等。
PIC方法按照求解电磁场方程形式的不同而分为静电模型、电磁模型和静磁模型，其中静电模型主要适用于静电分离为主要物理矛盾的带电粒子与时变静电场相互作用问题，如电推进系统中的离子引出过程、朗缪尔振荡、电子枪和收集极中电子的运动轨迹演化过程等。
静电模型求解的核心步骤如下：
1、电位求解，即通过求解静电场满足的离散泊松方程，得出所有网格点上的电位；
2、粒子受力求解，即通过相关网格点上的电位值得到网格内的电位分布，并求解其负梯度得出粒子所在位置处电场，然后求解受力；
3、推动粒子运动，即通过求解离散粒子运动方程，更新粒子的动量及位置等运动信息；
4、电荷分配，即根据粒子所在的位置求得其对周围网格点电荷的贡献，然后将所有粒子对网格点上的电荷贡献累加得到网格点上的电荷密度；
不断循环如上过程，直到计算结果收敛或人为设置的时间为止。
其中步骤2粒子受力求解是PIC静电模型必不可少的核心步骤之一，该步骤的精确且高效的求解对于PIC静电模型的整体求解精度和效率的控制十分重要。截止目前为止，PIC静电模型中粒子受力求解主要有两种方法，分别是有限差分(FD)方法和嵌入式有限元(IFE)方法。
FD方法：在PIC静电模型粒子受力求解的应用中，FD方法是通过采用结构化网格离散求解区域的，因此当求解粒子受力时，需要由粒子所属结构化网格上的电位通过采用插值的方法得出粒子所在位置处的电场，进而求解受力。
FD方法完全基于结构化网格，因此判断粒子所属网格、由网格上电位插值得到粒子所在位置处的电场等算法形式简单、易于理解，但是在PIC静电模型的应用中存在如下缺点：
1、FD方法采用的是由正交线划分而成的结构化网格，对于复杂曲形边界的拟合较差，从而使得复杂曲形边界附近的粒子受力求解的精度较低；
2、由于FD方法对网格尺寸均匀性的要求比较高，因此受限于模拟系统中细小物理结构的限制，必须在全部求解区域划分足够小的网格才能满足计算精度要求，从而使得总网格数巨大，而模拟粒子数目正比于总网格数，这就导致粒子受力求解的计算量十分庞大；
3、FD方法严格受限于数值稳定性条件的限制，即在对PIC静电模型的数值模拟中，如果空间网格尺寸很小，则时间步长也会随之取的很小，这会进一步增加对粒子受力时间循环求解的FD数值模拟负担。
针对FD方法在PIC静电模型粒子受力求解应用中出现的不足，Kafafy和Wang在2003年，提出了可应用于PIC静电模型粒子受力求解中的IFE方法。
IFE方法：在PIC静电模型粒子受力求解的应用中，IFE方法通过采用侵入式非结构化网格离散求解区域。侵入式非结构化网格划分情况如图1所示，可以看出此网格划分相当于有2重网格，其中第1重网格是结构化网格，第2重网格是将第1重网格中的每一个结构化网格进一步划分成五个四面体的侵入式非结构化网格。IFE方法在PIC静电模型的应用中，电位求解采用第2重侵入式非结构化网格，而粒子受力求解是在第1重结构化网格内进行的，因此IFE方法并没有克服FD方法在PIC静电模型粒子受力求解应用中的缺点。
发明内容
针对上述存在的问题和不足，为解决FD和IFE方法在粒子受力求解中对边界匹配度不高、求解精度不高以及数值模拟负担大的问题，本发明提供了一种应用于PIC静电模型的粒子受力有限元FEM求解方法。具体技术方案如下：
步骤1、电位求解。
采用全局非结构化网格，其三维网格划分实例如图2所示。
利用该非结构化网格离散求解区域和静电泊松方程，然后采用FEM的方法求解泊松方程来得到网格点上的电位。
步骤2、粒子受力求解。
采用基于非结构化网格的插值算法。
首先确定粒子所属的非结构化网格，然后根据粒子所在网格，由步骤1的计算结果，得到该四面体网格节点上的电位，再根据(1)式求解粒子所在位置处的电场。
其中E代表电场矢量，Φ代表单元内电位分布函数，上标e代表某一网格单元编号。将Φ在四个顶点处的值及顶点坐标带入(2)式：
Φe(x,y,z)＝ae+bex+cey+dez (2)
其中，x，y，z表示单元内任意位置全局坐标，可由克莱姆法则解得系数ae,be,ce,de，并将其带回(2)式，整理可得：
其中下标j表示e网格单元中的第j号顶点，为网格单元的插值函数，表示为：
其中V为网格单元的体积，将(3)式带入(1)式可得：
其中为坐标方向矢量，(5)式即为粒子所在位置处的电场矢量求解公式，进而采用公式(6)求解粒子受力。
F＝qE (6)
其中q为粒子所带电量。
步骤3、推动粒子运动
通过求解离散运动方程，更新粒子的动量及位置等运动信息；
步骤4、电荷分配
根据粒子所在的位置求得其对周围网格点电荷的贡献，然后将所有粒子对网格点上的电荷贡献累加得到网格点上的电荷密度；
步骤3至4的求解可采用结构化、浸入式非结构化或完全非结构化网格。
循环步骤1至4，直至达到收敛条件或模拟终止条件，最后进行数值诊断。
本发明适用于二维及三维结构，适用于二维时，网格划分从四面体网格变为三角形网格。
相对于PIC静电模型粒子受力求解的FD方法和IFE方法，本发明的有益效果体现在：
1、使用完全非结构化网格，该网格能够更好的拟合模型边界的形状，使得在复杂边界情况下PIC静电模型的粒子受力求解具有更高的计算精度；
2、将用于求解无源电磁场分布、热分析、力学分析等无粒子源问题的FEM方法结合到典型的PIC方法中，在保持典型PIC方法的计算简单、快速的优良特性的同时，利用FEM得到更高的有限元计算精度；
3、由于FEM方法既可以很好地匹配复杂边界，又可以根据模拟需要使用非均匀网格，且不受数值稳定性条件的限制，因此可以在保持计算精度的条件下，优化空间网格和时间步长，从而大幅地提高模拟效率。
附图说明
图1为PIC静电模型求解的IFE网格示意图；
图2为PIC静电模型求解的FEM网格示意图；
图3为七孔双栅离子光学系统的PIC静电模型计算实例示意图；
图4为七孔双栅离子光学系统的PIC静电模型计算实例网格划分示意图。
具体实施方式
下面通过实施实例对本发明作进一步详细说明。
以离子推进器七孔双栅离子光学系统为例，其示意图如图3所示。采用本发明中算法对这一实例进行PIC静电模拟的具体实施步骤如下：
步骤1、电位求解。
采用全局非结构化网格，利用该非结构化网格离散求解区域和静电泊松方程，然后采用FEM的方法求解泊松方程来得到网格点上的电位。
步骤2、离子受力求解。
采用基于非结构化网格的插值算法。
首先确定离子所属的非结构化网格，然后根据离子所在网格，由步骤1的计算结果，得到该四面体网格节点上的电位，再根据(7)式求解离子所在位置处的电场。
其中上标e代表某一网格单元。将Φ在四个顶点处的值及顶点坐标带入(8)式：
Φe(x,y,z)＝ae+bex+cey+dez (8)
可由克莱姆法则解得ae,be,ce,de，并将其带回(8)式，整理可得：
其中下标j表示e网格单元中的j号顶点，网格单元的插值函数为：
将(9)式带入(7)式可得：
(11)式即为离子所在位置处的电场求解公式，进而采用公式(12)求解离子受力。
F＝qE (12)
步骤3、推动离子运动
通过求解离散运动方程，更新离子的动量及位置等运动信息；
步骤4、电荷分配
根据离子所在的位置求得其对周围网格点电荷的贡献，然后将所有离子对网格点上的电荷贡献累加得到网格点上的电荷密度；
步骤3至4的求解可采用结构化、浸入式非结构化或完全非结构化网格。
循环步骤1至4，直至达到收敛条件或模拟终止条件，最后进行数值诊断。采用本发明中算法对这一实例进行PIC静电模拟的结果如图3所示。

Claims (3)

1.一种应用于PIC静电模型的粒子受力有限元求解方法，具体如下：
步骤1、电位求解；采用全局非结构化网格，利用该非结构化网格离散求解区域和静电泊松方程，然后采用FEM的方法求解泊松方程来得到网格点上的电位；
步骤2、粒子受力求解，采用基于非结构化网格的插值算法；
首先确定粒子所属的非结构化网格，然后根据粒子所在网格，由步骤1的计算结果，得到该四面体网格节点上的电位，再根据(1)式求解粒子所在位置处的电场；
Ee＝-▽Φe (1)
其中E代表电场矢量，Φ代表单元内电位分布函数，上标e代表某一网格单元编号；将Φ在四个顶点处的值及顶点坐标带入(2)式：
Φe(x,y,z)＝ae+bex+cey+dez (2)
其中，x，y，z表示单元内任意位置全局坐标，可由克莱姆法则解得系数ae,be,ce,de，并将其带回(2)式，整理可得：
其中下标j表示e网格单元中的第j号顶点，为网格单元的插值函数，表示为：
其中V为网格单元的体积，将(3)式带入(1)式可得：
其中为坐标方向矢量，(5)式即为粒子所在位置处的电场矢量求解公式，进而采用公式(6)求解粒子受力；
F＝qE (6)
其中q为粒子所带电量；
步骤3、推动粒子运动，通过求解离散运动方程，更新粒子的动量及位置等运动信息；
步骤4、电荷分配，根据粒子所在的位置求得其对周围网格点电荷的贡献，然后将所有粒子对网格点上的电荷贡献累加得到网格点上的电荷密度；
循环步骤1至4，直至达到收敛条件或模拟终止条件，最后进行数值诊断。
2.如权利要求1所述应用于PIC静电模型的粒子受力有限元求解方法，其特征在于：本发明适用于二维及三维结构，适用于二维时，网格划分从四面体网格变为三角形网格。
3.如权利要求1所述应用于PIC静电模型的粒子受力有限元求解方法，其特征在于：步骤3至4的求解采用结构化、浸入式非结构化或完全非结构化网格。
CN201810112819.5A 2018-02-05 2018-02-05 一种应用于pic静电模型的粒子受力有限元求解算法 Active CN108446429B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810112819.5A CN108446429B (zh) 2018-02-05 2018-02-05 一种应用于pic静电模型的粒子受力有限元求解算法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810112819.5A CN108446429B (zh) 2018-02-05 2018-02-05 一种应用于pic静电模型的粒子受力有限元求解算法

Publications (2)

Publication Number Publication Date
CN108446429A true CN108446429A (zh) 2018-08-24
CN108446429B CN108446429B (zh) 2021-07-06

Family

ID=63191775

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810112819.5A Active CN108446429B (zh) 2018-02-05 2018-02-05 一种应用于pic静电模型的粒子受力有限元求解算法

Country Status (1)

Country Link
CN (1) CN108446429B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113987898A (zh) * 2021-10-27 2022-01-28 电子科技大学 一种放电室中气体放电的粒子模拟加速方法

Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101063731A (zh) * 2006-04-30 2007-10-31 深圳大学 分布式反馈抑制半导体光放大器
CN101240412A (zh) * 2007-02-05 2008-08-13 富士通株式会社 磁控溅射设备的设计支持方法、系统和程序
US20130040092A1 (en) * 2011-08-08 2013-02-14 David A. Moore Carpet waste composite product and method for making same
CN103412988A (zh) * 2013-08-01 2013-11-27 电子科技大学 基于相移降阶模型周期结构的三维电磁场仿真模拟方法
CN104239623A (zh) * 2014-09-05 2014-12-24 兰州空间技术物理研究所 一种基于多时标粒子推动获得卫星表面电位的方法
CN104281740A (zh) * 2014-09-05 2015-01-14 兰州空间技术物理研究所 一种基于非均匀网格划分获得卫星表面电位的方法
CN105678002A (zh) * 2016-01-12 2016-06-15 中国科学技术大学 等离子体粒子-场自洽系统长期大规模高保真模拟方法
CN106067140A (zh) * 2016-05-31 2016-11-02 武汉大学 一种社会网络事件检测的混合指标量子群智能方法
US20160377814A1 (en) * 2015-06-29 2016-12-29 Coriant Advanced Technology, LLC Optimized 2x2 3db multi-mode interference coupler
CN106855957A (zh) * 2015-12-09 2017-06-16 四川大学 基于相似日和最小二乘支持向量机的工厂母线负荷预测
CN107526860A (zh) * 2017-03-31 2017-12-29 福州大学 基于电场能建模技术的vlsi标准单元布局方法
CN107577639A (zh) * 2017-08-19 2018-01-12 电子科技大学 一种应用于ecr离子源数值模拟的mpm混合算法

Patent Citations (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101063731A (zh) * 2006-04-30 2007-10-31 深圳大学 分布式反馈抑制半导体光放大器
CN101240412A (zh) * 2007-02-05 2008-08-13 富士通株式会社 磁控溅射设备的设计支持方法、系统和程序
US20130040092A1 (en) * 2011-08-08 2013-02-14 David A. Moore Carpet waste composite product and method for making same
CN103412988A (zh) * 2013-08-01 2013-11-27 电子科技大学 基于相移降阶模型周期结构的三维电磁场仿真模拟方法
CN104239623A (zh) * 2014-09-05 2014-12-24 兰州空间技术物理研究所 一种基于多时标粒子推动获得卫星表面电位的方法
CN104281740A (zh) * 2014-09-05 2015-01-14 兰州空间技术物理研究所 一种基于非均匀网格划分获得卫星表面电位的方法
US20160377814A1 (en) * 2015-06-29 2016-12-29 Coriant Advanced Technology, LLC Optimized 2x2 3db multi-mode interference coupler
CN106855957A (zh) * 2015-12-09 2017-06-16 四川大学 基于相似日和最小二乘支持向量机的工厂母线负荷预测
CN105678002A (zh) * 2016-01-12 2016-06-15 中国科学技术大学 等离子体粒子-场自洽系统长期大规模高保真模拟方法
CN106067140A (zh) * 2016-05-31 2016-11-02 武汉大学 一种社会网络事件检测的混合指标量子群智能方法
CN107526860A (zh) * 2017-03-31 2017-12-29 福州大学 基于电场能建模技术的vlsi标准单元布局方法
CN107577639A (zh) * 2017-08-19 2018-01-12 电子科技大学 一种应用于ecr离子源数值模拟的mpm混合算法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
黄桃 等: "包含碰撞效应的离子引出三维FEM-PIC模拟研究", 《第十八届全国等离子体科学技术会议摘要集》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113987898A (zh) * 2021-10-27 2022-01-28 电子科技大学 一种放电室中气体放电的粒子模拟加速方法
CN113987898B (zh) * 2021-10-27 2024-05-17 电子科技大学 一种放电室中气体放电的粒子模拟加速方法

Also Published As

Publication number Publication date
CN108446429B (zh) 2021-07-06

Similar Documents

Publication Publication Date Title
CN108416107A (zh) 一种应用于pic的推动粒子运动有限元算法
Groth et al. Global three‐dimensional MHD simulation of a space weather event: CME formation, interplanetary propagation, and interaction with the magnetosphere
CN111259599B (zh) 一种用于分析托卡马克中高能量粒子约束性能的测试粒子模拟方法
Bayyuk et al. A simulation technique for 2-D unsteady inviscid flows around arbitrarily moving and deforming bodies of arbitrary geometry
Marchand PTetra, a tool to simulate low orbit satellite–plasma interaction
CN106446432B (zh) 一种求解材料大变形的最优输运无网格方法
CN107577639B (zh) 一种应用于ecr离子源数值模拟的mpm混合模型模拟方法
CN102819647B (zh) 一种非均质材料随机微观结构有限元建模方法
JP5255714B2 (ja) 三次元の流体シミュレーション方法
CN108145975B (zh) 一种三维运动物体的磁场正演系统和方法
CN108460188A (zh) 一种应用于pic静电模型的电荷分配有限元fem求解算法
CN108280309A (zh) 一种应用于pic静电模型的电位有限元求解算法
JP6157941B2 (ja) プラズマシミュレーション方法及びプラズマシミュレーションプログラム
CN108446429A (zh) 一种应用于pic静电模型的粒子受力有限元求解算法
Frenod et al. An exponential integrator for a highly oscillatory Vlasov equation
Shen et al. An asynchronous and parallel time-marching method: Application to three-dimensional MHD simulation of solar wind
Chadwick et al. Recursive grain remapping scheme for phase‐field models of additive manufacturing
CN106227982A (zh) 一种电磁继电器静态特性计算方法及装置
Pfeiffer et al. Hyperbolic divergence cleaning, the electrostatic limit, and potential boundary conditions for particle-in-cell codes
Back et al. An axisymmetric PIC code based on isogeometric analysis
Cavenago Extraction layer models for negative ion sources
CN108268697A (zh) 一种高效率电推进羽流等离子体并行仿真方法
Spirkin A three-dimensional particle-in-cell methodology on unstructured Voronoi grids with applications to plasma microdevices
CN110163981A (zh) 一种基于运动相似性的引导发丝提取方法
Gan et al. Research on nondestructive measurement solving algorithm based on grid slice volume

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