CN108460188B - 一种应用于pic静电模型的电荷分配有限元fem求解算法 - Google Patents

一种应用于pic静电模型的电荷分配有限元fem求解算法 Download PDF

Info

Publication number
CN108460188B
CN108460188B CN201810114097.7A CN201810114097A CN108460188B CN 108460188 B CN108460188 B CN 108460188B CN 201810114097 A CN201810114097 A CN 201810114097A CN 108460188 B CN108460188 B CN 108460188B
Authority
CN
China
Prior art keywords
solving
grid
pic
particles
fem
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.)
Active
Application number
CN201810114097.7A
Other languages
English (en)
Other versions
CN108460188A (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 CN201810114097.7A priority Critical patent/CN108460188B/zh
Publication of CN108460188A publication Critical patent/CN108460188A/zh
Application granted granted Critical
Publication of CN108460188B publication Critical patent/CN108460188B/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
    • 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静电模型的电荷分配有限元FEM求解算法。本发明使用完全非结构化网格，以更好的拟合模型边界的形状，使得在复杂边界情况下PIC静电模型的电荷分配具有更高的计算精度；将用于求解无源电磁场分布、热分析、力学分析等无粒子源问题的有限元FEM求解方法结合到典型的PIC方法中，在保持典型PIC方法的计算简单、快速的优良特性的同时，利用FEM得到更高的有限元计算精度；由于FEM方法既可以很好地匹配复杂边界，又可以根据模拟需要使用非均匀网格，且不受数值稳定性条件的限制，因此可以在保持计算精度的条件下，优化空间网格和时间步长，从而大幅地提高模拟效率。

Description

一种应用于PIC静电模型的电荷分配有限元FEM求解算法
技术领域
本发明属于粒子模拟(Particle-in-cell，简写为PIC)的数值模拟领域，具体涉及一种应用于PIC静电模型的电荷分配有限元FEM求解算法。
背景技术
PIC方法是一种被广泛应用于带电粒子与电磁场相互作用物理问题中的数值模拟方法，它通过跟踪大量带电粒子在外加及自洽电磁场中的运动并统计平均而得到宏观特性及运动规律。经过几十年的发展，PIC模拟方法已经成为研究带电粒子与电磁场相互作用物理问题的一种强有力的数值手段，广泛应用于带电粒子与电磁场相互作用所涉及的许多领域，如磁约束聚变等离子体、惯性约束聚变等离子体、核爆、空间等离子体、人造等离子体(包括电子枪、离子源等)、电推进、自由电子激光以及电真空器件等。
PIC方法按照求解电磁场方程形式的不同而分为静电模型、电磁模型和静磁模型，其中静电模型主要适用于静电分离为主要物理矛盾的带电粒子与时变静电场相互作用问题，如电推进系统中的离子引出过程、朗缪尔振荡、电子枪和收集极中电子的运动轨迹演化过程等。
静电模型求解的核心步骤如下：
1、电位求解，即通过求解静电场满足的离散泊松方程，得出所有网格点上的电位；
2、粒子受力求解，即通过相关网格点上的电位值得到网格内的电位分布，并求解其负梯度得出粒子所在位置处电场，然后求解受力；
3、推动粒子运动，即通过求解离散粒子运动方程，更新粒子的动量及位置等运动信息；
4、电荷分配，即根据粒子所在的位置求得其对周围网格点电荷的贡献，然后将所有粒子对网格点上的电荷贡献累加得到网格点上的电荷密度；
不断循环如上过程，直到计算结果收敛或人为设置的时间为止。
其中步骤4电荷分配是PIC静电模型必不可少的核心步骤之一，该步骤的精确且高效的求解对于PIC静电模型的整体求解精度和效率的控制十分重要。截止目前为止，PIC静电模型中电荷分配主要有两种方法，分别是有限差分(FD)方法和嵌入式有限元(IFE)方法。
FD方法：在PIC静电模型电荷分配的应用中，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、粒子受力求解，通过相关网格点上的电位值得到网格内的电位分布，并求解其负梯度得出粒子所在位置处电场，然后求解受力；
步骤3、推动粒子运动
通过求解离散运动方程，更新粒子的动量及位置等运动信息；
步骤4、电荷分配
采用完全非结构化网格。
对于每个粒子的单步运算，粒子从前一时刻位置P1经过一个时间步长运动到当前时刻所在位置P2，选取P1和P2连线的中点P3，然后定位P3所在网格E1。将粒子所带电量Q除以E1的体积，即可得到粒子在E1内的空间电荷密度ρe。
将网格E1在四个顶点处的坐标值和电位值带入(1)式：
Φe(x,y,z)＝ae+bex+cey+dez (1)
其中，x，y，z分别表示单元内任意位置全局坐标，Φ为单元内的电位分布函数，上标e代表网格单元编号，可由克莱姆法则解得系数ae,be,ce,de，并将其带回(1)式，整理可得：
Figure BDA0001570156110000031
其中下标j表示e网格单元中的第j号顶点，
Figure BDA0001570156110000032
为网格单元的插值函数，表示为：
Figure BDA0001570156110000033
其中V为网格单元的体积。
对于E1的第j个顶点，利用(3)式乘以ρe，即可得到该空间电荷密度分配在该顶点上的值
Figure BDA0001570156110000034
Figure BDA0001570156110000035
将所有粒子在当前时间步长中分配到同一网格顶点的空间电荷密度叠加，从而得到步骤1所需的空间电荷密度。
步骤2至3的求解可采用结构化、浸入式非结构化或完全非结构化网格。
循环步骤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、粒子受力求解，通过相关网格点上的电位值得到网格内的电位分布，并求解其负梯度得出粒子所在位置处电场，然后求解受力；
步骤3、推动离子运动
通过求解离散运动方程，更新离子的动量及位置等运动信息；
步骤4、电荷分配。
采用完全非结构化网格。
对于每个离子的单步运算，离子从前一时刻位置P1经过一个时间步长运动到当前时刻所在位置P2，选取P1和P2连线的中点P3，然后定位P3所在网格E1。将离子所带电量Q除以E1的体积，即可得到离子在E1内的空间电荷密度ρe。
将网格E1在四个顶点处的值带入(5)式：
Φe(x,y,z)＝ae+bex+cey+dez (5)
可由克莱姆法则解得ae,be,ce,de，并将其带回(5)式，整理可得：
Figure BDA0001570156110000041
其中下标j表示e网格单元中的j号顶点，网格单元的插值函数为：
Figure BDA0001570156110000042
对于E1的第j个顶点，利用(7)式乘以ρe，即可得到该空间电荷密度分配在该顶点上的值
Figure BDA0001570156110000043
Figure BDA0001570156110000051
将所有离子在当前时间步长中分配到同一网格顶点的空间电荷密度叠加，从而得到步骤1所需的空间电荷密度。
步骤2至3的求解可采用结构化、浸入式非结构化以及完全非结构化等网格。
循环步骤1至4，直至达到收敛条件或模拟终止条件，最后进行数值诊断。采用本发明中算法对这一实例进行PIC静电模拟的结果如图3所示。

Claims (2)

1.一种应用于PIC静电模型的电荷分配有限元FEM求解方法，具体如下：
步骤1、电位求解；
采用全局非结构化网格，利用该非结构化网格离散求解区域和静电泊松方程，然后采用FEM的方法求解泊松方程来得到网格点上的电位；
步骤2、粒子受力求解，通过相关网格点上的电位值得到网格内的电位分布，并求解其负梯度得出粒子所在位置处电场，然后求解受力；
步骤3、推动粒子运动；通过求解离散运动方程，更新粒子的动量及位置运动信息；
步骤4、电荷分配，采用完全非结构化网格；
对于每个粒子的单步运算，粒子从前一时刻位置P1经过一个时间步长运动到当前时刻所在位置P2，选取P1和P2连线的中点P3，然后定位P3所在网格E1；将粒子所带电量Q除以E1的体积，即可得到粒子在E1内的空间电荷密度ρe；
将网格E1在四个顶点处的坐标值和电位值带入(1)式：
Φe(x,y,z)＝ae+bex+cey+dez (1)
其中，x，y，z分别表示单元内任意位置全局坐标，Φ为单元内的电位分布函数，上标e代表网格单元编号，由克莱姆法则解得系数ae,be,ce,de，并将其带回(1)式，整理可得：
Figure FDA0002972728350000011
其中下标j表示e网格单元中的第j号顶点，
Figure FDA0002972728350000012
为网格单元的插值函数，表示为：
Figure FDA0002972728350000013
其中V为网格单元的体积；
对于E1的第j个顶点，利用(3)式乘以ρe，即可得到该空间电荷密度分配在该顶点上的值
Figure FDA0002972728350000014
Figure FDA0002972728350000015
将所有粒子在当前时间步长中分配到同一网格顶点的空间电荷密度叠加，从而得到步骤1所需的空间电荷密度；
循环步骤1至4，直至达到收敛条件或模拟终止条件，最后进行数值诊断。
2.如权利要求1所述应用于PIC静电模型的电荷分配有限元FEM求解方法，其特征在于：所述步骤2至3的求解采用结构化、浸入式非结构化或完全非结构化网格。
CN201810114097.7A 2018-02-05 2018-02-05 一种应用于pic静电模型的电荷分配有限元fem求解算法 Active CN108460188B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810114097.7A CN108460188B (zh) 2018-02-05 2018-02-05 一种应用于pic静电模型的电荷分配有限元fem求解算法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810114097.7A CN108460188B (zh) 2018-02-05 2018-02-05 一种应用于pic静电模型的电荷分配有限元fem求解算法

Publications (2)

Publication Number Publication Date
CN108460188A CN108460188A (zh) 2018-08-28
CN108460188B true CN108460188B (zh) 2021-06-01

Family

ID=63239723

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810114097.7A Active CN108460188B (zh) 2018-02-05 2018-02-05 一种应用于pic静电模型的电荷分配有限元fem求解算法

Country Status (1)

Country Link
CN (1) CN108460188B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111581876B (zh) * 2020-05-13 2023-07-07 电子科技大学 一种用于圆柱坐标系粒子模拟的粒子源求解方法
CN113705969B (zh) * 2021-07-26 2022-12-09 西安交通大学 一种圆柱坐标系下电流电荷分配方法及系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101276378A (zh) * 2007-03-26 2008-10-01 中国航天科技集团公司第五研究院第五一○研究所 卫星太阳阵表面充电数值模拟预估的方法
CN104281740A (zh) * 2014-09-05 2015-01-14 兰州空间技术物理研究所 一种基于非均匀网格划分获得卫星表面电位的方法
CN106159296A (zh) * 2016-04-28 2016-11-23 江苏科技大学 一种固体氧化物燃料电池复合电极性质的计算方法
CN106354925A (zh) * 2016-08-26 2017-01-25 东北电力大学 一种基于电晕放电空间电势分布的模拟方法
CN107577639A (zh) * 2017-08-19 2018-01-12 电子科技大学 一种应用于ecr离子源数值模拟的mpm混合算法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103440163B (zh) * 2013-09-09 2016-06-08 中国科学院近代物理研究所 使用gpu并行实现的基于pic模型的加速器仿真方法
CN107273625A (zh) * 2017-06-22 2017-10-20 西南科技大学 一种三维不可压缩非定常n‑s方程有限元数值求解方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101276378A (zh) * 2007-03-26 2008-10-01 中国航天科技集团公司第五研究院第五一○研究所 卫星太阳阵表面充电数值模拟预估的方法
CN104281740A (zh) * 2014-09-05 2015-01-14 兰州空间技术物理研究所 一种基于非均匀网格划分获得卫星表面电位的方法
CN106159296A (zh) * 2016-04-28 2016-11-23 江苏科技大学 一种固体氧化物燃料电池复合电极性质的计算方法
CN106354925A (zh) * 2016-08-26 2017-01-25 东北电力大学 一种基于电晕放电空间电势分布的模拟方法
CN107577639A (zh) * 2017-08-19 2018-01-12 电子科技大学 一种应用于ecr离子源数值模拟的mpm混合算法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Development of a parallelized 3D electrostatic PIC-FEM code and its applications;J-S WU 等;《Computer Physics Communications》;20070731;第177卷(第1-2期);98-101 *
应用非结构性网格之平行化三维PIC-FEM程式的研究与发展;许国贤;《https://ir.nctu.edu.tw/handle/11536/81369》;20051231;1-235 *

Also Published As

Publication number Publication date
CN108460188A (zh) 2018-08-28

Similar Documents

Publication Publication Date Title
CN108416107B (zh) 一种应用于pic的推动粒子运动有限元算法
CN111259599B (zh) 一种用于分析托卡马克中高能量粒子约束性能的测试粒子模拟方法
Feng et al. Validation of the 3D AMR SIP–CESE solar wind model for four Carrington rotations
CN107577639B (zh) 一种应用于ecr离子源数值模拟的mpm混合模型模拟方法
Hara et al. Test cases for grid-based direct kinetic modeling of plasma flows
CN108460188B (zh) 一种应用于pic静电模型的电荷分配有限元fem求解算法
Wang et al. Electron–ion coupling in mesothermal plasma beam emission: Full particle PIC simulations
CN111709150A (zh) 任意形状磁体磁场空间分布仿真模拟方法
CN108280309B (zh) 一种应用于pic静电模型的电位有限元求解方法
CN107944113B (zh) 一种计算三维高速平动目标电磁散射场的方法
Chu et al. AN IMMERSED-FINITE-ELEMENT PARTICLE-IN-CELL SIMULATION TOOL FOR PLASMA SURFACE INTERACTION.
CN108446429B (zh) 一种应用于pic静电模型的粒子受力有限元求解算法
Bai et al. An implicit particle-in-cell model based on anisotropic immersed-finite-element method
Frenod et al. An exponential integrator for a highly oscillatory Vlasov equation
Zagórski et al. 2-D fluid model for discharge analysis of the RF-driven prototype ion source for ITER NBI (SPIDER)
CN108983290B (zh) 一种三维横向各向同性介质中旅行时确定方法及系统
Kovaleski Calculation of the ion extraction boundary of a plasma ion source
Cui et al. Development of a parallel multi-dimensional grid-based Vlasov solver for plasma plume simulation
Spirkin A three-dimensional particle-in-cell methodology on unstructured Voronoi grids with applications to plasma microdevices
Mayushan et al. On improving the particle-in-cell simulations accuracy for sources of charged particles
Xue et al. An Efficient iterative MoM-PO Hybrid Method for Calculating Scattered Fields of the Multiscale and Multiphysics Scatterers
Abe et al. Analysis of H− extraction in the Linac4 negative ion source by 2.5 D particle simulation
Wang Immersed finite element particle-in-cell modeling of surface charging in rarefied plasmas
CN110955955B (zh) 一种基于对数拟合的尘埃等离子体数值计算提速方法
Prudencio et al. Parallel 3d finite element numerical modelling of dc electron guns

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