一种圆柱坐标系下电流电荷分配方法及系统

1.本发明属于高功率微波技术领域，具体涉及一种圆柱坐标系下电流电荷分配方法及系统。


背景技术：

2.全电磁粒子模拟(particle in cell，pic)是模拟具有非线性和强耦合性的真空电子器件的重要方法之一，由于其对基本物理过程忽略较少，又被称为第一性原理方法。
3.pic方法中一个关键环节是带电粒子与电磁场的耦合，其中包含电荷分配和电流分配。电荷分配主要有线性分配方法和体积加权分配方法，两者在直角坐标系下是等价的。在柱坐标系下，由于网格是非均匀的，这导致线性分配方法和体积加权法不等价，且线性分配方法越靠近轴误差越大，体积加权法考虑到网格形状因素，因此与线性分配相比误差也更小。常用的电流分配方法有两种，一种是线性插值，另一种是电荷守恒方法，其中线性插值方法会造成电荷不守恒，而电荷守恒方法则不会。
4.目前国内三维柱坐标系下全电磁粒子模拟软件如unipic和chipic在电荷分配上都采用线性分配方法，在三维柱坐标系中这种分配方法的误差在越靠近轴时越大。电流分配方法采用线性插值方法，不满足电流连续性条件。


技术实现要素：

5.本发明所要解决的技术问题在于针对上述现有技术中的不足，提供一种圆柱坐标系下电流电荷分配方法及系统，通过适当调整角向分割个数，从而减少空间网格剖分个数，能够合理发挥器件的轴对称性，并大大降低计算量。
6.本发明采用以下技术方案：
7.一种圆柱坐标系下电流电荷分配方法，包括以下步骤：
8.s1、获得粒子的位置信息如果粒子不在轴上网格内，根据非轴上电荷分配计算出权重因子α，β，γ，根据非轴上电荷分配公式将带电粒子的电量分配到网格节点上；
9.如果粒子在轴上网格内，当r≤δr/2时，根据轴上电荷分配计算出权重因子α，β，γ；当r>δr/2时，根据非轴上电荷分配计算出权重因子α，β，γ，根据轴上电荷分配公式将带电粒子的电量分配到网格节点上；
10.s2、判断粒子运动是否在轴上网格内，如果运动不在轴上网格内，分别计算出粒子起始位置和末位置的权重因子α1，β1，γ1，α2，β2，γ2，计算δα，δβ，δγ，进行非轴上网格电流分配；
11.如果运动在轴上网格内，判断粒子运动是否穿过r=δr/2临界线，如果穿过就将位移从临界线处分为两段，即第一段和第二段计算第一段的权重因子α1，β1，γ1，α2，β2，γ2及对应的δα，δβ，δ
γ，进行轴上网格电流分配，计算出第二段的权重因子α3，β3，γ3，α4，β4，γ4及对应的δα，δβ，δγ，进行轴上网格电流分配；
12.如果没穿过r＝δr/2临界线，分别计算出粒子起始位置和末位置的权重因子α1，β1，γ1，α2，β2，γ2及对应的δα，δβ，δγ，进行轴上网格电流分配。
13.具体的，步骤s1中，非轴上网格电荷分配具体为：
14.s101、在圆柱坐标系中，将电荷q作为体积的粒子云，粒子云的中心坐标为δz、δr和分别为粒子所在网格z方向、r方向和角向的剖分步长，r
p
为粒子所在位置r方向坐标；
15.s102、计算对应粒子云与以网格八个节点为中心的立体单元的重叠体积δv
i,j,k
,δv
i+1,j,k
,δv
i,j+1,k
,δv
i+1,j+1,k
,δv
i,j,k+1
,δv
i+1,j,k+1
,δv
i,j+1,k+1
,δv
i+1,j+1,k+1
；
16.s103、根据步骤s102得到的重叠体积计算其与粒子云体积的比值；
17.s104、将步骤s103得到的体积的比值作为权重因子，将权重因子与电荷q相乘得到网格节点上所分配的电荷量。
18.进一步的，步骤s104中，网格节点上所分配的电荷量具体为：
19.q
i,j,k
＝(1
‑
α)(1
‑
β)(1
‑
γ)q
20.q
i+1,j,k
＝α(1
‑
β)(1
‑
γ)q
21.q
i,j+1,k
＝(1
‑
α)β(1
‑
γ)q
22.q
i+1,j+1,k
＝αβ(1
‑
γ)q
23.q
i,j,k+1
＝(1
‑
α)(1
‑
β)γq
24.q
i+1,j,k+1
＝α(1
‑
β)γq
25.q
i,j+1,k+1
＝(1
‑
α)βγq
26.q
i+1,j+1,k+1
＝αβγq
27.其中，α为z方向权重因子，β为r方向权重因子，γ为方向权重因子，q为带电粒子所带电荷量，q
i,j,k
为第(i,j,k)网格节点上分配到的电荷量。
28.具体的，步骤s1中，轴上网格电荷分配具体为：
29.s105、以r＝δr/2为临界线求出粒子云的体积，如果粒子所在位置满足r
p
≤δr/2，则将粒子看作体积为的粒子云；如果粒子所在位置满足r
p
＞δr/2，则将粒子看作体积为的粒子云；
30.s106、计算粒子云与以网格六个节点为中心的立体单元的重叠体积δv
i,0,k
,δv
i+1,0,k
,δv
i,1,k
,δv
i+1,1,k
,δv
i,1,k+1
,δv
i+1,1,k+1
；
31.s107、计算步骤s106中计算出的体积与粒子云体积的比值；
32.s108、让步骤s107中得到的体积比值作为权重因子，将权重因子与电荷q相乘得到网格节点上所分配的电荷量。
33.进一步的，步骤s108中，网格节点上所分配的电荷量具体为：
34.q
i,0,k
＝(1
‑
α)(1
‑
β)q
35.q
i+1,0,k
＝α(1
‑
β)q
36.q
i,1,k
＝(1
‑
α)β(1
‑
γ)q
37.q
i+1,1,k
＝αβ(1
‑
γ)q
38.q
i,1,k+1
＝(1
‑
α)βγq
39.q
i+1,1,k+1
＝αβγq
40.其中，α为z方向权重因子，β为r方向权重因子，γ为方向权重因子，q为带电粒子所带电荷量，q
i,0,k
为第(i,0,k)网格节点上分配到的电荷量。
41.具体的，步骤s2中，设δt时间内，电荷q在该网格内从运动到所形成的电流，非轴上网格的电流分配方法具体为：
42.s201、计算出粒子在处的权重因子α1，β1，γ1，以及粒子在处的权重因子α2，β2，γ2，采用柱坐标系下非轴上网格电荷分配方法分别将两个位置的粒子q分配到网格节点上，即和
43.s202、求出电荷变化量δq
i,j,k
，δq
i+1,j,k
，δq
i,j+1,k
，δq
i+1,j+1,k
，δq
i,j,k+1
，δq
i+1,j,k+1
，δq
i,j+1,k+1
，δq
i+1,j+1,k+1
；
44.s203、根据电流连续性方程列出12个电流和电荷量变化量之间的关系方程；
45.s204、求解步骤s203中的方程组，得到特解，根据特解得到网格边上12个电流的计算方法，将计算的电荷分配权重因子带入特解中得到非轴上网格边上的电流大小。
46.进一步的，步骤s204中，特解具体为：
[0047][0048][0049][0050][0051][0052][0053][0054][0055][0056][0057][0058]
[0059]
其中，δα为粒子末位置与起始位置z方向权重因子差值，δβ为粒子末位置与起始位置r方向权重因子差值，δγ为粒子末位置与起始位置方向权重因子差值，为粒子末位置与起始位置z方向权重因子平均值，为粒子末位置与起始位置r方向权重因子平均值，为粒子末位置与起始位置方向权重因子平均值，i
z,j+1/2,k,m
为z方向第(j,k,m)网格边上的电流值。
[0060]
具体的，步骤s2中，设δt时间内，电荷q在网格内从运动到形成电流，轴上网格的电流分配如下：
[0061]
s205、当粒子运动的位移穿过r＝δr/2临界线，将位移分为和两段，然后两段分别进行电流分配；
[0062]
s206、求出电荷变化量δq
i,0,k
,δq
i+1,0,k
,δq
i,1,k
,δq
i+1,1,k
,δq
i,1,k+1
,δq
i+1,1,k+1
。对于穿过临界线的情况，要分别求出两段位移的电荷变化量；
[0063]
s207、根据电流连续性方程列出9个电流和电荷量变化量之间的方程关系；
[0064]
s208、计算步骤s207的方程组，得到特解，根据特解得到网格边上9个电流的计算方法，将计算的电荷分配权重因子带入特解中得到非轴上网格边上的电流大小。
[0065]
进一步的，步骤s208中，特解具体为：
[0066][0067][0068][0069][0070][0071][0072][0073][0074][0075]
其中，δα为粒子末位置与起始位置z方向权重因子差值，δβ为粒子末位置与起始位置r方向权重因子差值，δγ为粒子末位置与起始位置方向权重因子差值，为粒子末
位置与起始位置z方向权重因子平均值，为粒子末位置与起始位置r方向权重因子平均值，为粒子末位置与起始位置方向权重因子平均值，为z方向第(i,0,k)网格边上的电流值。
[0076]
本发明的另一技术方案是，一种圆柱坐标系下电流电荷分配系统，包括：
[0077]
电荷分配模块，获得粒子的位置信息如果粒子不在轴上网格内，根据非轴上电荷分配计算出权重因子α，β，γ，根据非轴上电荷分配公式将带电粒子的电量分配到网格节点上；
[0078]
如果粒子在轴上网格内，当r≤δr/2时，根据轴上电荷分配计算出权重因子α，β，γ；当r>δr/2时，根据非轴上电荷分配计算出权重因子α，β，γ，根据轴上电荷分配公式将带电粒子的电量分配到网格节点上；
[0079]
电流分配模块，判断粒子运动是否在轴上网格内，如果运动不在轴上网格内，分别计算出粒子起始位置和末位置的权重因子α1，β1，γ1，α2，β2，γ2，计算δα，δβ，δγ，进行非轴上网格电流分配；
[0080]
如果运动在轴上网格内，判断粒子运动是否穿过r＝δr/2临界线，如果穿过就将位移从临界线处分为两段，即第一段和第二段计算第一段的权重因子α1，β1，γ1，α2，β2，γ2及对应的δα，δβ，δγ，进行轴上网格电流分配，计算出第二段的权重因子α3，β3，γ3，α4，β4，γ4及对应的δα，δβ，δγ，进行轴上网格电流分配；
[0081]
如果没穿过r＝δr/2临界线，分别计算出粒子起始位置和末位置的权重因子α1，β1，γ1，α2，β2，γ2及对应的δα，δβ，δγ，进行轴上网格电流分配。
[0082]
与现有技术相比，本发明至少具有以下有益效果：
[0083]
本发明一种圆柱坐标系下电流电荷分配方法，为了解决粒子模拟中电磁场与带电粒子耦合的问题，耦合问题包含有电荷分配和电流分配。由于yee网格(如图1)中电流密度定义在网格边中心，要求出电流密度就首先要得到电流，由电流连续性方程可知网格边上的电流计算方法与该网格边两端节点上电荷变化量有关，由于带电粒子不一定恰好处于网格节点上，因此需要采用适当的分配方法，在满足电荷守恒的前提下将电荷量分配到网格节点上，在满足电流连续性方程的条件下计算出网格边上的电流，进而可计算出电流密度。
[0084]
进一步的，在三维柱坐标系中，由于轴上网格和非轴上网格形状不同，因而电荷分配到网格节点的计算方式也不同，因此需要将非轴上网格单独讨论。
[0085]
进一步的，由电流连续性方程可知网格边上的电流计算方法与该网格边两端节点上电荷变化量有关，将电荷量分配到网格节点上是为了计算出该网格边上的电流。
[0086]
进一步的，在三维柱坐标系中，由于轴上网格和非轴上网格形状不同，因而电荷分
配到网格节点的计算方式也不同，因此需要将轴上网格单独讨论。
[0087]
进一步的，由电流连续性方程可知网格边上的电流计算方法与该网格边两端节点上电荷变化量有关，将电荷量分配到网格节点上是为了计算出该网格边上的电流。
[0088]
进一步的，电磁场依靠电流源进行推进，这就要求计算出到网格边上的电流密度，并将其作为电磁场推进的源。因此首先计算出网格上的电流，再除以相应的面积可很容易得到电流密度。
[0089]
进一步的，根据电流连续性方程列出的方程组有无穷多个解，只要找出一组满足方程的特解即可，可将villasenor等人提出的电流分配方法作为特解。该特解具有电流连续的优点，且充分利用了电荷权重中的权重因子，减少了计算量。
[0090]
进一步的，电磁场依靠电流源进行推进，这就要求计算出到网格边上的电流密度，并将其作为电磁场推进的源。因此首先计算出网格上的电流，再除以相应的面积可很容易得到电流密度。
[0091]
进一步的，特解具有电流连续的优点，且充分利用了电荷权重中的权重因子，减少了计算量。
[0092]
综上所述，本发明采用的电荷分配方法解决了三维柱坐标系下电荷分配问题，具有电荷守恒的优点，尤其对于轴上网格电荷分配；同时本发明提出的电流分配方法，解决三维柱坐标系下电流分配问题，具有电流连续性的优点。
[0093]
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
[0094]
图1为yee网格及场量定义图；
[0095]
图2为粒子云投影图，其中，(a)为z
‑
r投影，(b)为投影，(c)为投影；
[0096]
图3为三维柱坐标系下电荷守恒电流分配图；
[0097]
图4为轴上网格的电流分布图；
[0098]
图5为粒子模拟流程图；
[0099]
图6为电荷分配的流程图；
[0100]
图7为电流分配的流程图。
具体实施方式
[0101]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0102]
应当理解，当在本说明书和所附权利要求书中使用时，术语“包括”和“包含”指示所描述特征、整体、步骤、操作、元素和/或组件的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素、组件和/或其集合的存在或添加。
[0103]
还应当理解，在本发明说明书中所使用的术语仅仅是出于描述特定实施例的目的
而并不意在限制本发明。如在本发明说明书和所附权利要求书中所使用的那样，除非上下文清楚地指明其它情况，否则单数形式的“一”、“一个”及“该”意在包括复数形式。
[0104]
还应当进一步理解，在本发明说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合。
[0105]
在附图中示出了根据本发明公开实施例的各种结构示意图。这些图并非是按比例绘制的，其中为了清楚表达的目的，放大了某些细节，并且可能省略了某些细节。图中所示出的各种区域、层的形状及它们之间的相对大小、位置关系仅是示例性的，实际中可能由于制造公差或技术限制而有所偏差，并且本领域技术人员根据实际所需可以另外设计具有不同形状、大小、相对位置的区域/层。
[0106]
本发明提供了一种圆柱坐标系下电流电荷分配方法，推导并实现三维圆柱坐标系下cic电荷权重分配方法，包括非轴上网格电荷分配和轴上网格电荷分配；推导并实现三维圆柱坐标系下电荷守恒的电流分配方法，包括非轴上网格电流分配和轴上网格电流分配；在处理轴上网格时，以半径向步长为临界线，将电流进行分段分配。
[0107]
请参阅图1，磁场定义在网格面中心，电荷密度定义在网格顶点，通过求解newton
‑
lorentz方程得到带电粒子的运动状态与位移：
[0108][0109][0110]
对其进行中心差分离散得：
[0111][0112][0113]
由于带电粒子的位置不一定恰好在网格边或面上，因此在推进粒子时需要将电磁场插值到粒子位置处，并且将电荷以及电荷运动形成的电流分配到网格顶点和网格边上。
[0114]
本发明一种圆柱坐标系下电流电荷分配方法，包括以下步骤：
[0115]
s1、柱坐标系下cic电荷分配，包括非轴上网格电荷分配和轴上网格电荷分配；
[0116]
电荷分配方法一般有pic和cic。在三维情况下，pic实际上相当于线性插值或者分配方法，是利用三维网格单元八个顶点上的值通过线性插值得到网格内任一点处的相应值，或者将网格内任一点处的值按照与线性插值相同的权重分配到网格顶点上去的方法。cic把宏粒子看作是均匀的、网格单元大小的粒子云，根据粒子云与以相邻节点为中心的对偶网格交叉的体积占粒子云体积的比重来分配电荷。
[0117]
非轴上网格电荷分配具体为：
[0118]
s101、在圆柱坐标系中，把电荷q看作体积为的粒子云(δz、δr和分别为粒子所在网格z方向、r方向和角向的剖分步长，r
p
为粒子所在位置r方向坐标)，粒子云的中心坐标为粒子云分别在z
‑
r，坐标面上的投影如图2所示。
[0119]
s102、计算出该粒子云与以网格八个节点为中心的立体单元的重叠体积δv
i,j,k
,
δv
i+1,j,k
,δv
i,j+1,k
,δv
i+1,j+1,k
,δv
i,j,k+1
,δv
i+1,j,k+1
,δv
i,j+1,k+1
,δv
i+1,j+1,k+1
；
[0120]
s103、计算步骤s102得到的体积与粒子云体积的比值；
[0121][0122][0123][0124][0125][0126][0127][0128][0129][0130]
s104、让步骤s103中得到的体积比值作为权重因子，该权重因子与电荷q相乘得到网格节点上所分配的电荷量。
[0131]
q
i,j,k
＝(1
‑
α)(1
‑
β)(1
‑
γ)q
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(14)
[0132]
q
i+1,j,k
＝α(1
‑
β)(1
‑
γ)q
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0133]
q
i,j+1,k
＝(1
‑
α)β(1
‑
γ)q
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0134]
q
i+1,j+1,k
＝αβ(1
‑
γ)q
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(17)
[0135]
q
i,j,k+1
＝(1
‑
α)(1
‑
β)γq
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(18)
[0136]
q
i+1,j,k+1
＝α(1
‑
β)γq
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(19)
[0137]
q
i,j+1,k+1
＝(1
‑
α)βγq
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(20)
[0138]
q
i+1,j+1,k+1
＝αβγq
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(21)
[0139]
轴上网格电荷分配具体为：
[0140]
s105、以r＝δr/2为临界线求出粒子云的体积；
[0141]
s1051、如果粒子所在位置满足r
p
≤δr/2，则将粒子看作体积为的粒子云；
[0142]
s1052、如果粒子所在位置满足r
p
＞δr/2，则将粒子看作体积为的
粒子云。
[0143]
δz和δr分别为粒子所在网格z方向和r方向上的剖分步长，r
p
为粒子所在位置r方向坐标,粒子云的中心坐标为
[0144]
s106、计算出该粒子云与以网格六个节点为中心的立体单元的重叠体积δv
i,0,k
,δv
i+1,0,k
,δv
i,1,k
,δv
i+1,1,k
,δv
i,1,k+1
,δv
i+1,1,k+1
；
[0145]
(1)当时
[0146][0147][0148][0149][0150][0151][0152]
(2)当时
[0153][0154][0155][0156][0157][0158][0159]
s107、计算步骤s106中计算出的体积与粒子云体积的比值；
[0160][0161]
[0162][0163][0164][0165][0166][0167]
s108、让步骤s107中得到的体积比值作为权重因子，该权重因子与电荷q相乘得到网格节点上所分配的电荷量。
[0168]
q
i,0,k
＝(1
‑
α)(1
‑
β)q
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(29)
[0169]
q
i+1,0,k
＝α(1
‑
β)q
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(30)
[0170]
q
i,1,k
＝(1
‑
α)β(1
‑
γ)q
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(31)
[0171]
q
i+1,1,k
＝αβ(1
‑
γ)q
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(32)
[0172]
q
i,1,k+1
＝(1
‑
α)βγq
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(33)
[0173]
q
i+1,1,k+1
＝αβγq
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(34)
[0174]
s2、柱坐标系下的电流分配
[0175]
电荷的运动会产生电流，而根据yee网格和离散maxwell方程，这样得到的电流需要分配到网格线上。常用的分配方法有两种，一种是线性插值，另一种是电荷守恒方法。本发明中采用了电荷守恒的方法。
[0176]
请参阅图3，设δt时间内，电荷q在该网格内从运动到所形成的电流，非轴上网格的电流分配方法具体为：
[0177]
s201、计算出粒子在处的权重因子α1，β1，γ1，以及粒子在处的权重因子α2，β2，γ2，采用柱坐标系下cic电荷分配方法分别将两个位置的粒子q分配到网格节点上，即和和
[0178]
s202、求出电荷变化量δq
i,j,k
，δq
i+1,j,k
，δq
i,j+1,k
，δq
i+1,j+1,k
，δq
i,j,k+1
，δq
i+1,j,k+1
，δq
i,j+1,k+1
，δq
i+1,j+1,k+1
；
[0179]
s203、根据电流连续性方程列出12个电流和电荷量变化量之间的关系方程。
[0203][0204]
设δt时间内，电荷q在该网格内从运动到所形成的电流如图4所示，轴上网格的电流分配方法具体如下：
[0205]
s205、当r
p
≤δr/2和r
p
＞δr/2时权重因子计算方式不同，因此如果粒子运动的位移穿过了r＝δr/2临界线，就将该位移分为和两段，然后两段分别进行电流分配。下面就粒子是否穿过临界线进行讨论：
[0206]
s2051、粒子位移没有穿过临界线
[0207]
计算出粒子在处的权重因子α1，β1，γ1，以及粒子在处的权重因子α2，β2，γ2，采用柱坐标系下cic电荷分配方法分别将两个位置的粒子q分配到网格节点上，即和和
[0208]
s2052、粒子位移穿过临界线
[0209]
计算出粒子在处的权重因子α1,β1,γ1，以及粒子在处的权重因子α2,β2,γ2，采用柱坐标系下cic电荷分配方法分别将两个位置的粒子q分配到网格节点上得到和同理求出第二段位移的α1，β1，γ1和α2，β2，γ2以及和
[0210]
s206、求出电荷变化量δq
i,0,k
,δq
i+1,0,k
,δq
i,1,k
,δq
i+1,1,k
,δq
i,1,k+1
,δq
i+1,1,k+1
。对于穿过临界线的情况，要分别求出两段位移的电荷变化量；
[0211]
s207、根据电流连续性方程列出9个电流和电荷量变化量之间的方程关系；
[0212][0213][0214][0215]
[0216][0217][0218]
s208、计算步骤s207的方程组；
[0219]
轴上网格与非轴上网格所列的电流连续性方程组的不同之处是角向电流少了两个，z方向电流少了一个，但其他电流与网格节点上电荷变化量之间的关系保持一样，因此这些非轴上网格的电流的解作为轴上网格电流的解，进而可求出轴上网格电流分配方法的一个特解：
[0220][0221][0222][0223][0224][0225][0226][0227][0228][0229]
其中，
[0230]
δα＝α2‑
α1，δβ＝β2‑
β1，δγ＝γ2‑
γ1[0231][0232]
本发明再一个实施例中，提供一种圆柱坐标系下电流电荷分配系统，该系统能够用于实现上述圆柱坐标系下电流电荷分配方法，具体的，该圆柱坐标系下电流电荷分配系统包括电荷分配模块以及电流分配模块。
[0233]
其中，电荷分配模块，获得粒子的位置信息如果粒子不在轴上网格内，根据非轴上电荷分配计算出权重因子α，β，γ，根据非轴上电荷分配公式将带电粒子的电量分配到网格节点上；
[0234]
如果粒子在轴上网格内，当r≤δr/2时，根据轴上电荷分配计算出权重因子α，β，γ；当r>δr/2时，根据非轴上电荷分配计算出权重因子α，β，γ，根据轴上电荷分配公式将
带电粒子的电量分配到网格节点上；
[0235]
电流分配模块，判断粒子运动是否在轴上网格内，如果运动不在轴上网格内，分别计算出粒子起始位置和末位置的权重因子α1，β1，γ1，α2，β2，γ2，计算δα，δβ，δγ，进行非轴上网格电流分配；
[0236]
如果运动在轴上网格内，判断粒子运动是否穿过r＝δr/2临界线，如果穿过就将位移从临界线处分为两段，即第一段和第二段计算第一段的权重因子α1，β1，γ1，α2，β2，γ2及对应的δα，δβ，δγ，进行轴上网格电流分配，计算出第二段的权重因子α3，β3，γ3，α4，β4，γ4及对应的δα，δβ，δγ，进行轴上网格电流分配；
[0237]
如果没穿过r＝δr/2临界线，分别计算出粒子起始位置和末位置的权重因子α1，β1，γ1，α2，β2，γ2及对应的δα，δβ，δγ，进行轴上网格电流分配。
[0238]
本发明再一个实施例中，提供了一种终端设备，该终端设备包括处理器以及存储器，所述存储器用于存储计算机程序，所述计算机程序包括程序指令，所述处理器用于执行所述计算机存储介质存储的程序指令。处理器可能是中央处理单元(central processing unit，cpu)，还可以是其他通用处理器、数字信号处理器(digital signal processor、dsp)、专用集成电路(application specific integrated circuit，asic)、现成可编程门阵列(field
‑
programmable gatearray，fpga)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等，其是终端的计算核心以及控制核心，其适于实现一条或一条以上指令，具体适于加载并执行一条或一条以上指令从而实现相应方法流程或相应功能；本发明实施例所述的处理器可以用于圆柱坐标系下电流电荷分配方法的操作，包括：
[0239]
获得粒子的位置信息如果粒子不在轴上网格内，根据非轴上电荷分配计算出权重因子α，β，γ，根据非轴上电荷分配公式将带电粒子的电量分配到网格节点上；如果粒子在轴上网格内，当r≤δr/2时，根据轴上电荷分配计算出权重因子α，β，γ；当r>δr/2时，根据非轴上电荷分配计算出权重因子α，β，γ，根据轴上电荷分配公式将带电粒子的电量分配到网格节点上；
[0240]
判断粒子运动是否在轴上网格内，如果运动不在轴上网格内，分别计算出粒子起始位置和末位置的权重因子α1，β1，γ1，α2，β2，γ2，计算δα，δβ，δγ，进行非轴上网格电流分配；如果运动在轴上网格内，判断粒子运动是否穿过r＝δr/2临界线，如果穿过就将位移从临界线处分为两段，即第一段和第二段计算第一段的权重因子α1，β1，γ1，α2，β2，γ2及对应的δα，δβ，δγ，进行轴上网格电流分配，计算出第二段的权重因子α3，β3，γ3，α4，β4，γ4及对应的δα，δβ，δγ，应的δα，δβ，δγ，进行轴上网格电流分配；如果没穿过r＝δr/2临界线，分别计算出粒子起始位置和末位置的权重因子α1，β1，γ1，α2，β2，γ2及对应的δα，δβ，δγ，进行轴上网格电流分配。
[0241]
本发明再一个实施例中，本发明还提供了一种存储介质，具体为计算机可读存储介质(memory)，所述计算机可读存储介质是终端设备中的记忆设备，用于存放程序和数据。可以理解的是，此处的计算机可读存储介质既可以包括终端设备中的内置存储介质，当然也可以包括终端设备所支持的扩展存储介质。计算机可读存储介质提供存储空间，该存储空间存储了终端的操作系统。并且，在该存储空间中还存放了适于被处理器加载并执行的一条或一条以上的指令，这些指令可以是一个或一个以上的计算机程序(包括程序代码)。需要说明的是，此处的计算机可读存储介质可以是高速ram存储器，也可以是非不稳定的存储器(non
‑
volatile memory)，例如至少一个磁盘存储器。
[0242]
可由处理器加载并执行计算机可读存储介质中存放的一条或一条以上指令，以实现上述实施例中有关圆柱坐标系下电流电荷分配方法的相应步骤；计算机可读存储介质中的一条或一条以上指令由处理器加载并执行如下步骤：
[0243]
获得粒子的位置信息如果粒子不在轴上网格内，根据非轴上电荷分配计算出权重因子α，β，γ，根据非轴上电荷分配公式将带电粒子的电量分配到网格节点上；如果粒子在轴上网格内，当r≤δr/2时，根据轴上电荷分配计算出权重因子α，β，γ；当r>δr/2时，根据非轴上电荷分配计算出权重因子α，β，γ，根据轴上电荷分配公式将带电粒子的电量分配到网格节点上；
[0244]
判断粒子运动是否在轴上网格内，如果运动不在轴上网格内，分别计算出粒子起始位置和末位置的权重因子α1，β1，γ1，α2，β2，γ2，计算δα，δβ，δγ，进行非轴上网格电流分配；如果运动在轴上网格内，判断粒子运动是否穿过r＝δr/2临界线，如果穿过就将位移从临界线处分为两段，即第一段和第二段计算第一段的权重因子α1，β1，γ1，α2，β2，γ2及对应的δα，δβ，δγ，进行轴上网格电流分配，计算出第二段的权重因子α3，β3，γ3，α4，β4，γ4及对应的δα，δβ，δγ，应的δα，δβ，δγ，进行轴上网格电流分配；如果没穿过r＝δr/2临界线，分别计算出粒子起始位置和末位置的权重因子α1，β1，γ1，α2，β2，γ2及对应的δα，δβ，δγ，进行轴上网格电流分配。
[0245]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。通常在此处附图中的描述和所示的本发明实施例的组件可以通过各种不同的配置来布置和设计。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0246]
实施例
[0247]
请参阅图5，采用本发明方法应用于真空电子器件的模拟仿真中，具体调用步骤如下：
[0248]
第一步：按照圆柱坐标系进行网格剖分；
[0249]
第二步：设置初始时刻粒子的初始位置和速度；
[0250]
第三步：设置初始时刻的电场和磁场；
[0251]
第四步：将电场和磁场插值到粒子所在位置，进而计算出带电粒子的受力；
[0252]
第五步：求解newton
‑
lorentz方程得到粒子在第四步中求出的受力下的位移和速度；
[0253]
第六步：判断粒子是否碰到器件边界，如果碰到边界就吸收掉该粒子；
[0254]
第七步：将新时刻的粒子分配到网格节点上，将粒子运动形成的电流分配到网格边上；
[0255]
第八步：求解maxwell方程组，粒子运动形成的电流作为推进场的源得到下一时刻的电磁场；
[0256]
第九步：场边界截断条件，吸收到达计算区域边界的场，使之尽量不产生反射。
[0257]
解决三维柱坐标系下粒子模拟中电磁场与带电粒子耦合问题，其中包括电荷分配、电流分配以及电磁场插值等，且本发明具有电荷守恒和电流连续性等优点。
[0258]
请参阅图6，电荷分配算法实现步骤如下：
[0259]
第一步：获得粒子的位置信息如果粒子在轴上网格内则执行第二步，否则执行第四步；
[0260]
第二步：如果r≤δr/2，根据轴上电荷分配计算出权重因子α，β，γ；如果r>δr/2，根据非轴上电荷分配计算出权重因子α，β，γ；
[0261]
第三步：根据轴上电荷分配公式将带电粒子的电量分配到网格节点上，执行第六步；
[0262]
第四步：根据非轴上电荷分配计算出权重因子α，β，γ；
[0263]
第五步：根据非轴上电荷分配公式将带电粒子的电量分配到网格节点上；
[0264]
第六步：结束。
[0265]
请参阅图7，电流分配算法实现步骤如下：
[0266]
第一步：判断粒子运动是否在轴上网格内，如果该运动不在轴上网格内，就分别计算出粒子起始位置和末位置的权重因子α1，β1，γ1，α2，β2，γ2，然后执行第五步，否则执行第二步；
[0267]
第二步：判断粒子运动是否穿过了r＝δr/2临界线，如果穿过了就将位移从该临界线处分为两段，即第一段和第二段
[0268]
第三步：计算出第一段的权重因子α1，β1，γ1，α2，β2，γ2；
[0269]
第四步：计算出δα，δβ，δγ，执行第七步，如果两个位移段都计算完就执行第八步，否则计算出第二段的权重因子α3，β3，γ3，α4，β4，γ4，执行第四步；
[0270]
第五步：计算出δα，δβ，δγ，
[0271]
第六步：进行非轴上网格电流分配，执行第八步。
[0272]
第七步：进行轴上网格电流分配；
[0273]
第八步：结束
[0274]
综上所述，本发明一种圆柱坐标系下电流电荷分配方法及系统，解决了三维柱坐
标系下电磁场与带电粒子耦合问题，提出的三维柱坐标系下电荷分配方法和电流分配方法具有电荷守恒和电流连续性等优点，能够解决粒子模拟中场与粒子耦合问题。
[0275]
本领域内的技术人员应明白，本技术的实施例可提供为方法、系统、或计算机程序产品。因此，本技术可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本技术可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、cd
‑
rom、光学存储器等)上实施的计算机程序产品的形式。
[0276]
本技术是参照根据本技术实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
[0277]
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
[0278]
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
[0279]
以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。