CN108280309A - 一种应用于pic静电模型的电位有限元求解算法 - Google Patents
一种应用于pic静电模型的电位有限元求解算法 Download PDFInfo
- Publication number
- CN108280309A CN108280309A CN201810114106.2A CN201810114106A CN108280309A CN 108280309 A CN108280309 A CN 108280309A CN 201810114106 A CN201810114106 A CN 201810114106A CN 108280309 A CN108280309 A CN 108280309A
- Authority
- CN
- China
- Prior art keywords
- pic
- grid
- current potential
- electric charge
- static electric
- 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
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
技术领域
本发明属于粒子模拟(Particle-in-cell,简写为PIC)的数值模拟领域,具体涉及一种应用于PIC静电模型的电位有限元求解算法。
背景技术
PIC方法是一种被广泛应用于带电粒子与电磁场相互作用物理问题中的数值模拟方法,它通过跟踪大量带电粒子在外加及自洽电磁场中的运动并统计平均而得到宏观特性及运动规律。经过几十年的发展,PIC模拟方法已经成为研究带电粒子与电磁场相互作用物理问题的一种强有力的数值手段,广泛应用于带电粒子与电磁场相互作用所涉及的许多领域,如磁约束聚变等离子体、惯性约束聚变等离子体、核爆、空间等离子体、人造等离子体(包括电子枪、离子源等)、电推进、自由电子激光以及电真空器件等。
PIC方法按照求解电磁场方程形式的不同而分为静电模型、电磁模型和静磁模型,其中静电模型主要适用于静电分离为主要物理矛盾的带电粒子与时变静电场相互作用问题,如电推进系统中的离子引出过程、朗缪尔振荡、电子枪和收集极中电子的运动轨迹演化过程等。
静电模型求解的核心步骤如下:
1、电位求解,即通过求解静电场满足的离散泊松方程,得出所有网格点上的电位;
2、粒子受力求解,即通过相关网格点上的电位值得到网格内的电位分布,并求解其负梯度得出粒子所在位置处电场,然后求解受力;
3、推动粒子运动,即通过求解离散粒子运动方程,更新粒子的动量及位置等运动信息;
4、电荷分配,即根据粒子所在的位置求得其对周围网格点电荷的贡献,然后将所有粒子对网格点上的电荷贡献累加得到网格点上的电荷密度;
不断循环如上过程,直到计算结果收敛或人为设置的时间为止。
其中步骤1电位求解是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所示,可以看出此网格划分是将一个结构化的六面体又进一步划分成五个四面体,整体还是结构化网格。
IFE方法相比于FD方法能在一定程度上解决PIC静电模型电位求解中的复杂边界问题,但是,依然存在如下缺陷:
1、IFE方法虽然使用侵入式非结构化网格,但是整体还是基于结构化网格的划分,由于结构化网格的边界规整,对于不规则模型的边界匹配度不高,产生的误差较大;
2、从图1中可以看出,IFE方法使用的四面体网格的角度为直角,而根据有限元基本理论,四面体越接近正四面体,网格质量越好,计算精度越高,而嵌入式有限元网格与正四面体相差甚远,因此采用嵌入式有限元网格在同等网格规模的情况下,计算精度远远低于正四面体网格。另外,网格质量越好,最终形成的系数矩阵条件数越好,求解速度越快。
发明内容
针对上述存在的问题和不足,为解决FD和IFE方法在电位求解中边界匹配度不高、求解精度不高的问题,本发明提供了一种应用于PIC静电模型的电位有限元(FEM)求解方法。具体技术方案如下:
步骤1、电位求解。
采用全局非结构化网格,其三维网格划分实例如图2所示。
PIC静电模型中电位的求解是二阶偏微分方程定义的边值问题。
首先将该边值问题等价变换为如下变分问题:
其中:
其中α,β,γ是于与区域物理性质有关的已知参量,f是源或激励函数。然后将计算区域体V离散为M个四面体,在每个网格单元内,未知函数Φ表示为:
Φe(x,y,z)=ae+bex+cey+dez (3)
其中上标e代表某一网格单元。将Φ在四个顶点处的值带入(3)式,可由克莱姆法则解得系数ae,be,ce,de,并将其带回(3)式,整理可得:
其中下标j表示e网格单元中的j号顶点,网格单元的插值函数为:
然后将(4)式带入(2)式,得:
分别取Fe对每个顶点的偏导,并写成矩阵形式:
将(7)式给出的单元方程对所有单元求和,然后强加驻点条件,可得到方程组,并结合边界条件,解之即得每个网格顶点上的电位。
步骤2、粒子受力求解,通过相关网格点上的电位值得到网格内的电位分布,并求解其负梯度得出粒子所在位置处电场,然后求解受力。
步骤3、推动粒子运动
通过求解离散运动方程,更新粒子的动量及位置等运动信息;
步骤4、电荷分配
根据粒子所在的位置求得其对周围网格点电荷的贡献,然后将所有粒子对网格点上的电荷贡献累加得到网格点上的电荷密度;
步骤2至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算法。
求解电位满足的Possion方程:
式中,ρ为电荷密度,ε0为真空介电常数。
变换为等价变分问题为:
因为求解的是Possion方程,所以αx=1,αy=1,αz=1,β=0,则可得:
然后将三维计算区域体V离散为M个四面体,网格划分情况如图4所示。在每个网格单元内,未知函数Φ为:
Φe(x,y,z)=ae+bex+cey+dez (11)
将Φ在四个顶点处的值带入(11)式,并由克莱姆法则可以解得
其中Ve为四面体体积,进而得:
同样地也可以解出be,ce,de。并将解出的ae,be,ce,de带入(11)式,合并同类项,可得
其中下标j表示e网格单元中的j号顶点,网格单元的插值函数为:
然后将(14)式带入(10)式,得:
取Fe对每个顶点的偏导:
分别取Fe对每个顶点的偏导,并对F应用驻点条件,则针对网格单元e可构成线性方程组:
其中,
最后,将M个如式(18)这样的线性方程组组合成一个线性方程组,求解方程组即得每个网格顶点上的电位。
步骤2、粒子受力求解,通过相关网格点上的电位值得到网格内的电位分布,并求解其负梯度得出粒子所在位置处电场,然后求解受力;
步骤3、推动离子运动
通过求解离散运动方程,更新离子的动量及位置等运动信息;
步骤4、电荷分配
根据离子所在的位置求得其对周围网格点电荷的贡献,然后将所有离子对网格点上的电荷贡献累加得到网格点上的电荷密度;
步骤2至4的求解可采用结构化、浸入式非结构化以及完全非结构化等网格。
循环步骤1至4,直至达到收敛条件或模拟终止条件,最后进行数值诊断。采用本发明中算法对这一实例进行PIC静电模拟的结果如图3所示。
Claims (2)
1.一种应用于PIC静电模型的电位有限元求解算法,具体如下:
步骤1、电位求解,采用全局非结构化网格;
PIC静电模型中电位的求解是二阶偏微分方程定义的边值问题,首先将该边值问题等价变换为如下变分问题:
其中:
其中α,β,γ是于与区域物理性质有关的已知参量,f是源或激励函数;然后将计算区域体V离散为M个四面体,在每个网格单元内,未知函数Φ表示为:
Φe(x,y,z)=ae+bex+cey+dez (3)
其中上标e代表某一网格单元,将Φ在四个顶点处的值带入(3)式,可由克莱姆法则解得系数ae,be,ce,de,并将其带回(3)式,整理可得:
其中下标j表示e网格单元中的j号顶点,网格单元的插值函数为:
然后将(4)式带入(2)式,得:
分别取Fe对每个顶点的偏导,并写成矩阵形式:
将(7)式给出的单元方程对所有单元求和,然后强加驻点条件,可得到方程组,将得到的方程组组合成一个矩阵方程,并结合边界条件,解之即得每个网格顶点上的电位;
步骤2、粒子受力求解,通过相关网格点上的电位值得到网格内的电位分布,并求解其负梯度得出粒子所在位置处电场,然后求解受力;
步骤3、推动粒子运动,通过求解离散运动方程,更新粒子的动量及位置等运动信息;
步骤4、电荷分配,根据粒子所在的位置求得其对周围网格点电荷的贡献,然后将所有粒子对网格点上的电荷贡献累加得到网格点上的电荷密度;
循环步骤1至4,直至达到收敛条件或模拟终止条件,最后进行数值诊断。
2.如权利要求1所述应用于PIC静电模型的电位有限元求解算法,其特征在于:所述步骤2至4的求解可采用结构化、浸入式非结构化或完全非结构化网格。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810114106.2A CN108280309B (zh) | 2018-02-05 | 2018-02-05 | 一种应用于pic静电模型的电位有限元求解方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810114106.2A CN108280309B (zh) | 2018-02-05 | 2018-02-05 | 一种应用于pic静电模型的电位有限元求解方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108280309A true CN108280309A (zh) | 2018-07-13 |
CN108280309B CN108280309B (zh) | 2021-12-03 |
Family
ID=62807679
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810114106.2A Active CN108280309B (zh) | 2018-02-05 | 2018-02-05 | 一种应用于pic静电模型的电位有限元求解方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108280309B (zh) |
Cited By (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111967149A (zh) * | 2020-08-03 | 2020-11-20 | 电子科技大学 | 一种用于粒子模拟算法的粒子运动半插值求解方法 |
CN114722670A (zh) * | 2022-04-02 | 2022-07-08 | 电子科技大学 | 一种用于二维静电粒子模型的电势分布有限差分求解方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR100674204B1 (ko) * | 2005-03-16 | 2007-01-24 | 주식회사 브이엠티 | 스퍼터 이온펌프의 배기 방법 및 그에 대한 구조 |
CN103226638A (zh) * | 2013-04-24 | 2013-07-31 | 兰州空间技术物理研究所 | 一种电推进器产生的等离子体分布特性数值模拟预估方法 |
CN107577639A (zh) * | 2017-08-19 | 2018-01-12 | 电子科技大学 | 一种应用于ecr离子源数值模拟的mpm混合算法 |
-
2018
- 2018-02-05 CN CN201810114106.2A patent/CN108280309B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
KR100674204B1 (ko) * | 2005-03-16 | 2007-01-24 | 주식회사 브이엠티 | 스퍼터 이온펌프의 배기 방법 및 그에 대한 구조 |
CN103226638A (zh) * | 2013-04-24 | 2013-07-31 | 兰州空间技术物理研究所 | 一种电推进器产生的等离子体分布特性数值模拟预估方法 |
CN107577639A (zh) * | 2017-08-19 | 2018-01-12 | 电子科技大学 | 一种应用于ecr离子源数值模拟的mpm混合算法 |
Non-Patent Citations (2)
Title |
---|
刘畅 等: "基于PIC方法的离子发动机光学系统粒子模拟", 《北京航空航天大学学报》 * |
王二蒙 等: "基于IFE-PIC和MCC方法的离子推进器加速栅极孔壁腐蚀机理仿真", 《高电压技术》 * |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111967149A (zh) * | 2020-08-03 | 2020-11-20 | 电子科技大学 | 一种用于粒子模拟算法的粒子运动半插值求解方法 |
CN111967149B (zh) * | 2020-08-03 | 2022-11-04 | 电子科技大学 | 一种用于粒子模拟算法的粒子运动半插值求解方法 |
CN114722670A (zh) * | 2022-04-02 | 2022-07-08 | 电子科技大学 | 一种用于二维静电粒子模型的电势分布有限差分求解方法 |
CN114722670B (zh) * | 2022-04-02 | 2023-08-04 | 电子科技大学 | 一种用于二维静电粒子模型的电势分布有限差分求解方法 |
Also Published As
Publication number | Publication date |
---|---|
CN108280309B (zh) | 2021-12-03 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Remelli et al. | Meshsdf: Differentiable iso-surface extraction | |
Lim et al. | Surface reconstruction techniques: a review | |
CN108416107A (zh) | 一种应用于pic的推动粒子运动有限元算法 | |
CN107577639B (zh) | 一种应用于ecr离子源数值模拟的mpm混合模型模拟方法 | |
CN111259599B (zh) | 一种用于分析托卡马克中高能量粒子约束性能的测试粒子模拟方法 | |
CN113158527B (zh) | 一种基于隐式fvfd计算频域电磁场的方法 | |
CN108280309A (zh) | 一种应用于pic静电模型的电位有限元求解算法 | |
Koomullil et al. | Flow simulation using generalized static and dynamic grids | |
CN108460188A (zh) | 一种应用于pic静电模型的电荷分配有限元fem求解算法 | |
Bai et al. | An implicit particle-in-cell model based on anisotropic immersed-finite-element method | |
CN108446429A (zh) | 一种应用于pic静电模型的粒子受力有限元求解算法 | |
CN116776695A (zh) | 基于超高阶有限元技术的一维电磁计算方法、系统及设备 | |
Tskhakaya et al. | PIC/MC code BIT1 for plasma simulations on HPC | |
CN111339688A (zh) | 基于大数据并行算法求解火箭仿真模型时域方程的方法 | |
CN111968113B (zh) | 一种基于最优传输映射的脑影像二维卷积深度学习方法 | |
Pfeiffer et al. | Hyperbolic divergence cleaning, the electrostatic limit, and potential boundary conditions for particle-in-cell codes | |
Zhang et al. | Two-dimensional finite element mesh generation algorithm for electromagnetic field calculation | |
Kumar et al. | A multi-GPU based accurate algorithm for simulations of gas-liquid flows | |
CN110163981A (zh) | 一种基于运动相似性的引导发丝提取方法 | |
Spirkin | A three-dimensional particle-in-cell methodology on unstructured Voronoi grids with applications to plasma microdevices | |
Vshivkova et al. | 3D numerical modeling of the plasma beam expansion using the MHD-kinetic approach | |
CN117217130B (zh) | 一种逐层推进的壁面距离确定方法、装置以及介质 | |
Akalin et al. | An accelerated BEM formulation for the forward problem solution of ESI using realistic head models | |
Jin et al. | Analytically integratable zero-restlength springs for capturing dynamic modes unrepresented by quasistatic neural networks | |
Duarte et al. | Mean curvature flow and applications |
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 |