电磁粒子模拟(Electromagnetic Particle-in-Cell, EMPIC)方法^[1-8]被广泛应用于很多等离子体相关现象和应用中^[9-14],如雷达和通信系统、聚变等离子体加热、射频加速器、电介质击穿,以及高功率微波、毫米波和太赫兹波等. EMPIC方法的数学模型包含描述电磁现象的Maxwell方程组和描述带电粒子分布的Vlasov方程,其中国际单位制下的Maxwell方程组为

(1.1a) $ \nabla \times {B} = \mu_0 {J} + \frac{1}{c^2}\frac{\partial {E} }{\partial t}, \\ $

(1.1b) $ \nabla \times {E} = -\frac{\partial {B} }{\partial t}, {\qquad}\ \; \\ $

(1.1c) $ \nabla \cdot {E} = \frac{\rho}{\epsilon_0}, {\qquad}{\qquad}\ \\ $

(1.1d) $ \nabla \cdot {B} = 0, {\qquad}{\qquad}\ \ \\ $

这里的符号取它们通常表示的物理含义.描述一种无碰撞等离子的Vlasov方程为

(1.2) $ \begin{equation} \frac{\partial f({x}, {v}, t)}{\partial t} + {v}\cdot \nabla_{{x}} f({x}, {v}, t) + \frac{{F}}{m} \cdot \nabla_{{v}}f({x}, {v}, t) = 0, \end{equation} $

这里$ f $表示带电等离子体粒子在相空间的分布函数, $ {x} $, $ {v} $和$ m $分别是粒子的位置、速度和质量, $ {F} $是作用在带电粒子上的洛伦兹力. $ {v} = {\rm d}{x}/{\rm d}t \label{motionEq} $和$ {{{F}}\over{ m} } = {\rm d}{{v}}/{\rm d}t \label{LorentzForce} $是Vlasov方程(1.2)的特征方程.方程(1.1a)–(1.1d)和(1.2)的耦合是通过$ f $的前两阶矩,即电荷密度$ \rho({x}, t) $和电流密度$ {J}({x}, t) $实现的.具体地

(1.3) $ \begin{eqnarray} && \rho({x}, t) = q \int_{{{\Bbb R}} ^3} f({x}, {v}, t) {\rm d}v, \end{eqnarray} $

(1.4) $ \begin{eqnarray} &&{J}({x}, t) = q \int_{{{\Bbb R}} ^3}{v} f({x}, {v}, t) {\rm d}v, \end{eqnarray} $

这里$ q $是等离子体带电粒子的电荷.对于方程(1.1a)–(1.1d)来说,容易验证如果在初始时刻Maxwell方程组的后两个散度方程成立,并且$ \rho $和$ {J} $满足电流连续性方程

(1.5) $ \begin{equation} \frac{\partial \rho }{\partial t} + \nabla \cdot {J} = 0, \end{equation} $

那么这两个散度方程就一直成立,电磁场动力学问题仅需求解Maxwell方程组的前两个旋度方程即可.但是电流连续性方程(1.5)在离散情况下却不一定成立,电磁场量、电荷密度和电流密度定义在空间网格上,带电粒子可以取任意位置.采用某种分配方案,从带电粒子的运动信息得到的电荷密度和电流密度不一定满足离散电流连续性方程.本文关注满足电荷守恒的电流分配方案,即能够使得网格上的电荷和电流密度满足离散电流连续性方程的分配方案.

为了保证离散旋度算子和散度算子相容,离散电流连续性方程与Maxwell场求解方法是密切相关的^[15].当采用基于Yee网格的时域有限差分(FDTD)方法求解Maxwell场方程时,电磁场量、电荷密度及电流密度在二维Yee网格中一个元胞上的定义位置如图 1所示^[16].

对于指标为$ (i, \ j) $的元胞来说,相应的离散电流连续性方程为

(1.6) $ \begin{equation} \frac{\rho_{i, j}^{n+1} - \rho_{i, j}^{n}}{ \Delta t} + \frac{J^{n+1/2}_{x_{i+1/2, j}} - J^{n+1/2}_{x_{i-1/2, j}} }{\Delta x} + \frac{J^{n+1/2}_{y_{i, j+1/2}} - J^{n+1/2}_{y_{i, j-1/2}} }{\Delta y} = 0. \end{equation} $

当采用FDTD方法求解电磁场时,有三种常用的满足电荷守恒的电流分配方案^[16-18].具体地, Villasenor等^[16]假设一个电荷分布为常数的宏粒子在一个时间步内沿直线运动,任何被宏粒子所扫略过的元胞边界都会有一个电流通过.该方法在本文被称为“V方法”. Esirkepov等^[17]假设一个带有任意阶电荷分布的宏粒子在一个时间步内沿直线运动,同时给出了电流密度的一个特别的定义：电流密度的每一个分量都是几个“形因子(form-factor)”差的唯一线性组合.该方法在本文被称为“E方法”. Umeda等^[18]假定宏粒子在一个时间步内的运动轨迹是一条由两个直的子线段所组成的折线.网格上的电流密度由每一小段轨迹所对应的电荷运动通量相加得到.该方法在本文被称为“U方法”.基于这些方法, Barthelme^[19], Yu^[20]等人提出了基于更高阶形状函数的、满足电荷守恒的电流分配方案.

本文针对一个带有常数电荷分布的宏粒子在二维、三维Yee网格上移动所涉及的、电荷守恒的电流分配方案,用求解线性代数方程组的方式得出通用公式.对于Yee网格上的某一个元胞,当有电流从元胞边界流进或流出时,元胞内的电量会发生变化.因此,我们只关注那些在宏粒子运动前后电荷量发生变化的元胞及定义在这些元胞边界上的电流.当这些元胞及其内部电荷变化量被确定了之后,根据电荷守恒定律就能够建立起一个以所有元胞边界上电流密度为变量的线性代数方程组.求得线性代数方程组的一个解,即得到一种电流密度分配方案.求得线性代数方程组所有的解,即得到所有的电流密度分配方案.在二维Yee网格上,带有常数电荷分布的宏粒子在一个时间内的运动可能会涉及四个、六个或七个元胞共三种情形.针对每一种情形,本文列出并求解所对应的线性代数方程组,同时将V方法, E方法和U方法作为其特解详细分析.在三维Yee网格中,带有常数电荷分布的宏粒子在一个时间内的运动可能会涉及不同数目的元胞共八种情形,详见文献[18].本文将针对一种最常见、最简单涉及八个元胞的情形,作同样的分析.

在EMPIC模拟中,并不追踪所有的带电粒子,而是用相空间中$ N $个带有一定权重$ \omega_p $的宏粒子代替.每个宏粒子代表一些物理坐标和速度坐标都比较相近的带电粒子.带电粒子分布函数$ f({x}, {v}, t) $被近似为$ f_N({x}, {v}, t) $,即

$ f({x}, {v}, t) \approx f_N({x}, {v}, t) = \sum\limits_{p = 1}^{N} \omega_p g({x};\ {x}_p(t))\delta({v}-{v}_p(t)), $

其中, $ {x}_p(t) $是$ t $时刻第$ p $个宏粒子在物理空间中的坐标; $ g({x}; \ {x}_p(t)) $是定义在物理空间中的、用于描述宏粒子电荷分布的某个函数; $ \delta({v}, \ t) $是定义在速度空间的Dirac冲激函数; $ \omega_p $是第$ p $个宏粒子的权重.这种近似可以认为是将分布函数$ f $近似为相空间中$ N $个坐标为$ ({x}_p(t), {v}_p(t)) $、带有权重$ \omega_p $的宏粒子的总和.

在早期的模拟中,并没有赋予宏粒子形状和大小,函数$ g({x}; \ {x}_p(t)) $是一个Dirac冲激函数, $ f_N({x}, {v}, t) $近似为

$ f_{N, \delta}({x}, {v}, t) = \sum\limits_{p = 1}^{N} \omega_p \delta({x}-{x}_p(t))\delta({v}-{v}_p(t)). $

研究成果^[3]发现为每一个宏粒子(所带电荷)赋予一个有限大小和具体的形状,可以大大减小宏粒子的近距离相互作用、提高计算精度,这里用形状函数$ D({x}; {x_p}(t)) $来表示宏粒子(所带电荷)的大小和形状.假设所有的宏粒子有同样的大小和形状,则$ f_N({x}, {v}, t) $近似为

$ f_{N, D}({x}, {v}, t) = \sum\limits_{p = 1}^{N} \omega_p D({x};\ {x_p}(t))\delta({v}-{v_p}(t)). $

接下来继续介绍一个在EMPIC中经常用到的概念:形因子,用函数$ S_{({i})}({x}) $表示. $ qS_{({i})}({x}) $表示一个带电量为$ q $、中心位于$ {x} $的宏粒子对第$ {i} $个元胞的电量的贡献, $ S_{({i})}({x}) $表示宏粒子对元胞$ {i} $的电荷贡献量与宏粒子所带总电量的比值. 图 2给出了一个二维的Yee网格以及一个形状函数的支集为一个元胞大小的宏粒子.该Yee网格包含四个相邻的元胞, $ (i, \ j) $、$ (i+1, \ j) $、$ (i, \ j+1) $和$ (i+1, \ j+1) $分别是这些元胞的指标, $ \Omega_{ij} $是元胞$ (i, \ j) $的面积;阴影部分为宏粒子,它的中心坐标为$ (x_p, \ y_p) $,宏粒子形状函数的支集面积为$ \Omega_p $.

在任意一个元胞上对宏粒子的形状函数做积分,可以得到宏粒子相应于该元胞的形因子.假设图 2中宏粒子的带电量为$ q $,它对元胞$ (i, \ j) $的电量贡献为

$ qS_{(i, \ j)}(x_p, \ y_p) = q\int\!\!\!\!\int_{\Omega_{ij} \bigcap \Omega_p}D(x, \ y;\ x_p, \ y_p) {\rm d} x {\rm d} y, $

其中$ \Omega_{ij} $是元胞$ (i, \ j) $的区域, $ \Omega_p $是宏粒子的支集.宏粒子对该元胞的形因子为

$ S_{(i, \ j)}(x_p, \ y_p) = \int\!\!\!\!\int_{\Omega_{ij} \cap \Omega_p}D(x, \ y;\ x_p, \ y_p) {\rm d} x {\rm d} y. $

表 1列出了几个典型而常用的一维形状函数和形因子,通过张量积运算可以将它们自然推广至高维.

带有一定形状和大小的宏粒子在网格上移动,会引起网格上某些元胞内的电量发生变化.将宏粒子在一个时间内运动前后所涉及的元胞集合记为$ C = \{C_1, C_2, \cdots\} $,其中$ C_i $为涉及的第$ i $个元胞.由电荷守恒定律可知,单位时间内一个元胞的电荷增加量等于通过元胞边界流入的电流.对每一个元胞$ C_i $建立电荷增量与边界上电流的等式关系,将$ C $中所有元胞边界上的电流组成一个向量$ {J} $、所有元胞的电荷增量组成一个向量$ {b} $,则它们满足一个线性代数方程组: $ A {J} = {b} $,其中$ A $是由各个边界面积所构成的系数矩阵.

带有一阶形状函数的宏粒子在二维Yee网格上移动一个时间步,可能会引起四个、六个或七个元胞内的电荷量发生变化. 图 3是一个示意图,最左边的图是涉及四个元胞的情况, $ \xi_1, \cdots, \xi_4 $分别是需要计算的电流密度;中间的图是涉及六个元胞的情况, $ \xi_1, \cdots, \xi_7 $分别是需要计算的电流密度;最右边的是涉及七个元胞的情况, $ \xi_1, \cdots, \xi_{12} $分别是需要计算的电流密度.这里并没有将$ C $中所有元胞边界上的所有电流密度作为变量.对于元胞集合$ C = \{C_1, C_2, \cdots\} $的所有相邻外层元胞,若与$ C $中的所有元胞仅有一个公共边界,则可知穿过该边界的电流为$ 0 $.假设宏粒子在一个时间步内从$ (x^n, y^n) $移动到$ (x^{n+1}, y^{n+1}) $;坐标原点选在距离宏粒子起始位置$ (x^n, y^n) $最近的网格线的交叉点上;空间步长$ \delta x = \delta y = h = 1 $;宏粒子的形状函数为常数分布(一阶),形状函数的支集为一个元胞大小;宏粒子的电量为$ q $.下文将列出并详细分析与图 3相对应的线性代数方程组.

在涉及四个元胞时,相关形因子见附录表. $ rank(A) = 3 < 4 $, $ A{J} = {b} $的通解形式为$ {J} = k_1{j}^1 + {J}^0 $,其中$ k_1 \in {{\Bbb R}} $是参数, $ {j}^1 = \left(-1, \ 1, \ -1, \ 1\right)^T $是解空间的基向量, $ {J}^0 $是$ A{J} = {b} $的某一个特解.为了得到$ A{J} = {b} $的一个特解,可以为宏粒子指定一个特殊的路径:在一个时间步长内,先沿$ x $轴从$ (x^n, \ y^n) $运动到$ (x^{n+1}, \ y^n) $,再沿$ y $轴从$ (x^{n+1}, \ y^n) $运动到$ (x^{n+1}, \ y^{n+1}) $.由电荷守恒定律以及形因子的定义,可以得到特解$ {J}^0 $.具体写出通解的表达式为

(3.1) $ \begin{equation} {J} = k_1{j}^1 + {J}^0 = \left(\begin{array}{cc} { } -k_1 + \frac{q}{\Delta t}(x^{n+1}-x^n)(0.5-y^n) \\ { } k_1 + \frac{q}{\Delta t}(0.5-x^{n+1})(y^{n+1}-y^{n}) \\ { } -k_1 + \frac{q}{\Delta t}(0.5+x^{n+1})(y^{n+1}-y^{n}) \\ { } k_1 + \frac{q}{\Delta t}(x^{n+1}-x^n)(0.5+y^n) \end{array}\right). \end{equation} $

下面验证V方法、E方法和U方法给出的公式分别是通解公式(3.1)的三个特解.在四个元胞情形下, V方法、E方法和U方法给出的公式列在表 2中.

当参数$ k_1 = q (x^{n+1}-x^{n})(y^{n+1}-y^{n}) / (2\Delta t) $时,有$ {J} = {J}^V = {J}^E = {J}^U $成立.

在涉及六个元胞时,相关形因子见附录表 2. $ rank(A) = 5 < 7 $, $ A{J} = {b} $的通解为$ k_1{j}^1 + k_2{j}^2 + {J}^0 $,其中$ k_1 \in {{\Bbb R}} $和$ k_2 \in {{\Bbb R}} $是参数, $ {j}^1 $和$ {j}^2 $是解空间的两个基向量

$ \begin{eqnarray*} {j}^1 = \left(-1, 1, -1, 1, 0, 0, 0 \right)^T, \quad {j}^2 = \left( 0, 0, 1, 0, -1, -1, 1\right)^T. \end{eqnarray*} $

$ {J}^0 = \left(\xi^0_1, \xi^0_2, \cdots, \xi^0_{7}\right)^T $是特解,它是在为宏粒子指定了特殊路径后得到的(该特殊路径与在四个元胞的情形中所指定的路径一致).具体写出通解的表达式为

(3.2) $ \begin{equation} {J} = k_1{j}^1 + k_2{j}^2 + {J}^0 = \left(\begin{array}{cc} -k_1 + 0 + {q}(0.5-x^n)(0.5-y^n)/{\Delta t} \\ k_1 + 0 + 0 \\ -k_1 + k_2 + {q}(1.5-x^{n+1})(y^{n+1}-y^{n})/{\Delta t} \\ k_1 + 0 + {q}(0.5-x^n)(0.5+y^{n}) /{\Delta t} \\ -k_2 + {q}(x^{n+1}-0.5)(0.5-y^{n}) /{\Delta t} \\ -k_2 + {q}(x^{n+1}-0.5)(y^{n+1}-y^{n})/{\Delta t} \\ k_2 + {q}(x^{n+1}-0.5)(0.5+y^{n}) /{\Delta t} \end{array}\right). \end{equation} $

下面验证V方法、E方法和U方法给出的公式分别是通解公式(3.2)的三个特解.在六个元胞情形下, V方法、E方法和U方法给出的公式列在表 3中.

当参数$ k_1 = \frac{q}{(2\Delta t\Delta x)} (\Delta x_1)^2 \Delta y $, $ k_2 = \frac{q}{(2\Delta t)} \left(1+\frac{\Delta x_1}{\Delta x} \right)\Delta x_2 \Delta y $时,有$ {J}^V = {J} $成立.

当参数$ k_1 = \frac{q}{(2\Delta t)} \Delta x_1 \Delta y $, $ k_2 = \frac{q}{(2\Delta t)} \Delta x_2 \Delta y $时,有$ {J}^E = {J} $成立.

当参数$ k_1 = \frac{q}{(4\Delta t)} \Delta x_1 \Delta y $, $ k_2 = \frac{3q}{(4\Delta t)} \Delta x_2 \Delta y $时,有$ {J}^U = {J} $成立.

在涉及七个元胞时,相关形因子见附录表 3.注意这里,宏粒子在运动前后涉及的元胞集合为$ C = \{C_{i, j}, C_{i+1, j}, C_{i, j+1}, C_{i+1, j+1}, C_{i+2, j+1}, C_{i+1, j+2}, C_{i+2, j+2} \} $,分析所有与$ C_i $相邻的外层元胞可知,外层元胞$ C_{i+2, j} $、$ C_{i, j+2} $都与$ C $中的元胞有两条公共边界,故需要计算这两个元胞的电荷变化量. $ rank(A) = 8 < 12 $, $ A{J} = {b} $的通解为$ k_1{j}^1 + k_2{j}^2 + k_3{j}^3 + k_4{j}^4 + {J}^0 $,其中$ k_i \in {{\Bbb R}} \ (i = 1, 2, 3, 4) $是参数,有

$ {j}^1 = \left(1, -1, 1, -1, 0, 0, 0, 0, 0, 0, 0, 0 \right)^T, \quad {j}^2 = \left(0, 0, -1, 0, 1, 1, -1, 0, 0, 0, 0, 0 \right)^T, $

$ {j}^3 = \left(0, 0, 0, 0, 0, 0, 1, -1, 1, -1, 0, 0 \right)^T, \quad {j}^4 = \left(0, 0, 0, -1, 0, 0, 0, -1, 0, 0, 1, 1 \right)^T $

是解空间的四个基向量. $ {J}^0 = (\xi^0_1, \xi^0_2, \cdots, \xi^0_{12})^T $是特解,它是在为宏粒子指定了特殊路径后得到的(该特殊路径与在四个元胞的情形中所指定的路径一致).具体写出通解的表达式为

(3.3) $ \begin{equation} {J} = \sum\limits_{i = 1}^4 k_i{j}^i + {J}^0\\ = \left(\begin{array}{cccccc} { } k_1 + 0 + 0 + 0 + \frac{q}{\Delta t}(0.5-x^{n})(0.5-y^{n}) \\ { } -k_1 + 0 + 0 + 0 + 0 \\ { } k_1 -k_2 + 0 + 0 + \frac{q}{\Delta t}(1.5-x^{n+1})(0.5-y^{n}) \\ { } -k_1 + 0 + 0 -k_4 + \frac{q}{\Delta t}(0.5-x^{n})(0.5+y^{n}) \\ { } 0 + k_2 + 0 + 0 + \frac{q}{\Delta t}(x^{n+1}-0.5)(0.5-y^{n}) \\ { } 0 + k_2 + 0 + 0 + \frac{q}{\Delta t}(x^{n+1}-0.5)(0.5-y^{n}) \\ { } 0 -k_2 + k_3 + 0 + \frac{q}{\Delta t}(x^{n+1}-0.5)(0.5+y^{n}) \\ { } 0 + 0 -k_3 -k_4 + \frac{q}{\Delta t}(1.5-x^{n+1})(y^{n+1}-0.5) \\ { } 0 + 0 + k_3 + 0 + \frac{q}{\Delta t}(x^{n+1}-0.5)(y^{n+1}-0.5) \\ 0 + 0 -k_3 + 0 + 0 \\ 0 + 0 + 0 + k_4 + 0 \\ 0 + 0 + 0 + k_4 + 0 \end{array}\right). \end{equation} $

下面验证V方法、E方法和U方法给出的公式分别是通解公式(3.3)的三个特解.在七个元胞情形下, V方法、E方法和U方法给出的公式列在表 4中.

当参数

$ k_1 = -\frac{q}{2\Delta t} \Delta x_1 \Delta x_1 \frac{\Delta y}{\Delta x}, \quad k_2 = -\frac{q}{2\Delta t} \Delta x_2 \left( 0.5-y^{n}+\Delta y \frac{\Delta x_1}{\Delta x} \right), $

$ k_3 = \frac{q}{2\Delta t} (y^{n+1}-0.5) \left( \Delta x \frac{\Delta y_2}{\Delta y}-x^{n+1}+0.5 \right), {\quad} k_4 = 0 $

时,有$ {J} = {J}^V $成立.其中, $ \Delta x = x^{n+1}-x^{n} $, $ \Delta y = y^{n+1}-y^{n} $, $ \Delta x_1 = 0.5-x^{n} $, $ \Delta y_1 = \Delta y {\Delta x_1}/{\Delta x} $, $ y_{M_1} = y^{n}+\Delta y_1 $, $ \Delta y_2 = 0.5-y_{M_1} $, $ \Delta x_2 = \Delta x {\Delta y_2}/{\Delta y} $.

当参数

$ k_1 = -\frac{q}{2\Delta t} (0.5-x^{n})(0.5-y^{n}), \quad k_2 = -\frac{q}{2\Delta t} (x^{n+1}-0.5)(0.5-y^{n}), $

$ k_3 = -\frac{q}{2\Delta t} (x^{n+1}-0.5)(y^{n+1}-0.5), \quad k_4 = \frac{q}{2\Delta t} (0.5-x^{n})(y^{n+1}-0.5) $

时,有$ {J} = {J}^E $成立.

当参数

$ k_1 = -\frac{q}{2\Delta t} (0.5-x^{n})(0.5-y^{n}), \quad k_2 = -\frac{q}{\Delta t} (x^{n+1}-0.5)(0.5-y^{n}), $

$ k_3 = -\frac{q}{2\Delta t} (x^{n+1}-0.5)(y^{n+1}-0.5), {\quad} k_4 = 0 $

时,有$ {J} = {J}^U $成立.

带有一阶形状函数的宏粒子在满足CFL的条件下在三维Yee网格上移动一个时间步,会涉及更多的元胞.有八种不同的情况,具体情况列在了文献[18]的表 1中.本文以图 4所示的、涉及八个元胞的一种最常见且最简单的情形,作为例子进行阐述.列出并分析了相对应的线性代数方程组.假设条件在二维公式的假设条件下加上一个条件: $ \delta z = h = 1 $.

在该情形下,相关形因子见附录表 4, $ rank(A) = rank(A, b) = 7 < 12 $, $ A{J} = {b} $的通解

$ {J} = \sum\limits_{i = 1}^5 k_i{j}^i + {J}^0 , $

其中$ k_i \in {{\Bbb R}} \ (i = 1, \cdots, 5) $是参数,有

$ {j}^1 = \left(1, -1, 0, 0, -1, 1, 0, 0, 0, 0, 0, 0\right)^T, \quad {j}^2 = \left(0, 0, 1, -1, 0, 0, -1, 1, 0, 0, 0, 0\right)^T, $

$ {j}^3 = \left(1, 0, -1, 0, 0, 0, 0, 0, -1, 1, 0, 0\right)^T, \quad {j}^4 = \left(0, 0, 0, 0, 1, 0, -1, 0, -1, 0, 1, 0\right)^T, $

$ {j}^5 = \left(0, 1, 0, -1, 1, 0, -1, 0, -1, 0, 0, 1\right)^T $

是解空间的基向量, $ {J}^0 $是$ A{J} = {b} $的一个特解.为了得到$ A{J} = {b} $的一个特解,可以为宏粒子指定一个特殊的路径：在时间区间$ [t^n, t^{n+1}] $内,先沿$ x $轴从$ (x^n, y^n, z^n) $运动到$ (x^{n+1}, y^n, z^n) $,再沿$ y $轴从$ (x^{n+1}, y^n, z^n) $运动到$ (x^{n+1}, y^{n+1}, z^n) $,最后沿$ z $轴从$ (x^{n+1}, y^{n+1}, z^n) $运动到$ (x^{n+1}, y^{n+1}, z^{n+1}) $.由电荷守恒定律以及形因子的定义,得到特解

$ {J}^0 = (\xi^0_1, \xi^0_2, \ldots \xi^0_{11}, \xi^0_{12})^T. $

具体写出通解的表达式为

(4.1) $ \begin{equation} {J} = \sum\limits_{i = 1}^5 k_i{j}^i + {J}^0 = \left(\begin{array}{cccccc} { } k_1 + k_3 + \frac{q}{\Delta t}(x^{n+1}-x^n)(0.5-y^n)(0.5-z^n) \\ { } -k_1 + k_5 + \frac{q}{\Delta t}(x^{n+1}-x^n)(0.5+y^n)(0.5-z^n) \\ { } k_2 - k_3 + k_5 + \frac{q}{\Delta t}(x^{n+1}-x^n)(0.5-y^n)(0.5+z^n) \\ { } -k_2 - k_5 + \frac{q}{\Delta t}(x^{n+1}-x^n)(0.5+y^n)(0.5+z^n) \\ { } -k_1 + k_4 + k_5 + \frac{q}{\Delta t}(y^{n+1}-y^n)(0.5-x^{n+1})(0.5-z^n) \\ { } k_1 + \frac{q}{\Delta t}(y^{n+1}-y^n)(0.5+x^{n+1})(0.5-z^n) \\ { } -k_2 - k_4 -k_5 + \frac{q}{\Delta t}(y^{n+1}-y^n)(0.5-x^{n+1})(0.5+z^n) \\ { } k_2 + \frac{q}{\Delta t}(y^{n+1}-y^n)(0.5+x^{n+1})(0.5+z^n) \\ { } -k_3 - k_4 -k_5 + \frac{q}{\Delta t}(z^{n+1}-z^n)(0.5-x^{n+1})(0.5-y^{n+1}) \\ { } k_3 + \frac{q}{\Delta t}(z^{n+1}-z^n)(0.5+x^{n+1})(0.5-y^{n+1}) \\ { } k_4 + \frac{q}{\Delta t}(z^{n+1}-z^n)(0.5-x^{n+1})(0.5+y^{n+1}) \\ { } k_5 + \frac{q}{\Delta t}(z^{n+1}-z^n)(0.5+x^{n+1})(0.5+y^{n+1}) \end{array}\right) . \end{equation} $

上面已经给出了一种最常见、最简单的情形下,电流密度的通用公式.其他情形下电流密度的通用公式可以用同样的方法给出.下面给出V方法和E方法的公式,找出相应于通解(4.1)的参数.这里不再列出U方法,因为经过仔细推导发现, U方法给出的三维公式不再满足电荷守恒(具体推导见附录).

在该情形下, V方法给出的公式为

$ {J}^V = \frac{q}{\Delta t} \left(\begin{array}{cccccc} \Delta x [ (1-\bar{\eta})(1-\bar{\zeta}) + {\Delta y \Delta z}/{12} ] \\ \Delta x [ \bar{\eta} (1-\bar{\zeta}) - {\Delta y \Delta z}/{12} ] \\ \Delta x [ (1-\bar{\eta}) \bar{\zeta} - {\Delta y \Delta z}/{12} ] \\ \Delta x [ \bar{\eta} \bar{\zeta} + {\Delta y \Delta z}/{12} ] \\ \Delta y [ (1-\bar{\xi})(1-\bar{\zeta}) + {\Delta x \Delta z}/{12} ] \\ \Delta y [ \bar{\xi}(1-\bar{\zeta}) - {\Delta x \Delta z}/{12} ] \\ \Delta y [ (1-\bar{\xi})\bar{\zeta} - {\Delta x \Delta z}/{12} ] \\ \Delta y [ \bar{\xi}\bar{\zeta} + {\Delta x \Delta z}/{12} ] \\ \Delta z [ (1-\bar{\xi})(1-\bar{\eta}) + {\Delta x \Delta y}/{12} ] \\ \Delta z [ \bar{\xi}(1-\bar{\eta}) - {\Delta x \Delta y}/{12} ] \\ \Delta z [ (1-\bar{\xi})\bar{\eta} - {\Delta x \Delta y}/{12} ] \\ \Delta z [ \bar{\xi}\bar{\eta} + {\Delta x \Delta y}/{12} ] \\ \end{array}\right), $

其中, $ {\xi} = x-x_i, \ {\eta} = y-y_j, \ \zeta = z-z_k, \ (x_i, y_j, z_k) = (-0.5, -0.5, -0.5), $$ \Delta x = {\xi}^{n+1} - {\xi}^n, $$ \Delta y = {\eta}^{n+1} - {\eta}^n, $$ \Delta z = {\zeta}^{n+1} - {\zeta}^n, $$ \bar{\xi} = ({\xi}^n+{\xi}^{n+1})/2, \ \bar{\eta} = ({\eta}^n+{\eta}^{n+1})/2, \ \bar{\zeta} = ({\zeta}^n+{\zeta}^{n+1})/2 $.

在该情形下, E方法给出的公式为

$ {J}^E = \frac{q}{4\Delta t} \left(\begin{array}{cccccc} (x^{n+1}-x^n) \Big( (0.5-y^n)(1-\frac{4}{3}z^n-\frac{2}{3}z^{n+1}) + (0.5-y^{n+1})(1-\frac{2}{3}z^n-\frac{4}{3}z^{n+1}) \Big) \\ (x^{n+1}-x^n) \Big( (0.5+y^n)(1-\frac{4}{3}z^n-\frac{2}{3}z^{n+1}) + (0.5+y^{n+1})(1-\frac{2}{3}z^n-\frac{4}{3}z^{n+1}) \Big) \\ (x^{n+1}-x^n) \Big( (0.5-y^n)(1+\frac{4}{3}z^n+\frac{2}{3}z^{n+1}) + (0.5-y^{n+1})(1+\frac{2}{3}z^n+\frac{4}{3}z^{n+1}) \Big) \\ (x^{n+1}-x^n) \Big( (0.5+y^n)(1+\frac{4}{3}z^n+\frac{2}{3}z^{n+1}) + (0.5+y^{n+1})(1+\frac{2}{3}z^n+\frac{4}{3}z^{n+1}) \Big) \\ (y^{n+1}-y^n) \Big( (0.5-x^n)(1-\frac{4}{3}z^n-\frac{2}{3}z^{n+1}) + (0.5-x^{n+1})(1-\frac{2}{3}z^n-\frac{4}{3}z^{n+1}) \Big) \\ (y^{n+1}-y^n) \Big( (0.5+x^n)(1-\frac{4}{3}z^n-\frac{2}{3}z^{n+1}) + (0.5+x^{n+1})(1-\frac{2}{3}z^n-\frac{4}{3}z^{n+1}) \Big) \\ (y^{n+1}-y^n) \Big( (0.5-x^n)(1+\frac{4}{3}z^n+\frac{2}{3}z^{n+1}) + (0.5-x^{n+1})(1+\frac{2}{3}z^n+\frac{4}{3}z^{n+1}) \Big) \\ (y^{n+1}-y^n) \Big( (0.5+x^n)(1+\frac{4}{3}z^n+\frac{2}{3}z^{n+1}) + (0.5+x^{n+1})(1+\frac{2}{3}z^n+\frac{4}{3}z^{n+1}) \Big) \\ (z^{n+1}-z^n) \Big( (0.5-x^n)(1-\frac{4}{3}y^n-\frac{2}{3}y^{n+1}) + (0.5-x^{n+1})(1-\frac{2}{3}y^n-\frac{4}{3}y^{n+1}) \Big) \\ (z^{n+1}-z^n) \Big( (0.5+x^n)(1-\frac{4}{3}y^n-\frac{2}{3}y^{n+1}) + (0.5+x^{n+1})(1-\frac{2}{3}y^n-\frac{4}{3}y^{n+1}) \Big) \\ (z^{n+1}-z^n) \Big( (0.5-x^n)(1+\frac{4}{3}y^n+\frac{2}{3}y^{n+1}) + (0.5-x^{n+1})(1+\frac{2}{3}y^n+\frac{4}{3}y^{n+1}) \Big) \\ (z^{n+1}-z^n) \Big( (0.5+x^n)(1+\frac{4}{3}y^n+\frac{2}{3}y^{n+1}) + (0.5+x^{n+1})(1+\frac{2}{3}y^n+\frac{4}{3}y^{n+1}) \Big) \end{array}\right). $

$ {J} = {J}^V $,当参数

$ k_1 = q \Delta y \left[ -\Delta z \bar{\xi} - \Delta x (1-\bar{\zeta}) - 2\Delta x \Delta z/3 \right] / (2\Delta t), $

$ k_2 = q \Delta y \left[ \Delta z \bar{\xi} - \Delta x \bar{\zeta} + 2\Delta x \Delta z/3 \right] / (2\Delta t), $

$ k_3 = q \Delta z \left[ \Delta y \bar{\xi} - \Delta x (1-\bar{\eta}) + \Delta x \Delta y/3 \right] / (2\Delta t), $

$ k_4 = q \Delta z \left[ -\Delta y (1-\bar{\xi}) + \Delta x \eta + \Delta x \Delta y/3 \right] / (2\Delta t), $

$ k_5 = q \Delta z \left[ -\Delta y \bar{\xi} - \Delta x \eta - \Delta x \Delta y/3 \right] / (2\Delta t); $

$ {J} = {J}^E $,当参数

$ k_1 = \frac{q}{4\Delta t} (y^{n+1}-y^n) \Big( (0.5+x^n)( 1 - \frac{4}{3}z^n - \frac{2}{3}z^{n+1} ) + (0.5+x^{n+1})( -1 + \frac{10}{3}z^n - \frac{4}{3}z^{n+1} ) \Big), $

$ k_2 = \frac{q}{4\Delta t} (y^{n+1}-y^n) \Big( (0.5+x^n)( 1 + \frac{4}{3}z^n + \frac{2}{3}z^{n+1} ) + (0.5+x^{n+1})( -1 - \frac{10}{3}z^n + \frac{4}{3}z^{n+1} ) \Big), $

$ k_3 = \frac{q}{4\Delta t} (z^{n+1}-z^n) \Big( (0.5+x^n)( 1 - \frac{4}{3}y^n-\frac{2}{3}y^{n+1} ) + (0.5+x^{n+1})( -1 - \frac{2}{3}y^n + \frac{8}{3}y^{n+1} ) \Big), $

$ k_4 = \frac{q}{4\Delta t} (z^{n+1}-z^n) \Big( (0.5-x^n)( 1 + \frac{4}{3}y^n + \frac{2}{3}y^{n+1} ) + (0.5-x^{n+1})( -1 + \frac{2}{3}y^n - \frac{8}{3}y^{n+1} ) \Big), $

$ k_5 = \frac{q}{4\Delta t} (z^{n+1}-z^n) \Big( (0.5+x^n)( 1 + \frac{4}{3}y^n + \frac{2}{3}y^{n+1} ) + (0.5+x^{n+1})( -1 + \frac{2}{3}y^n - \frac{8}{3}y^{n+1} ) \Big). $

本文从代数角度提出了一个EMPIC模拟中、适用于Yee网格的、电荷守恒的电流分配方案所满足的统一公式.这个公式是通过求解建立在电荷守恒定律上的线性代数方程组而得到的.一个带有一阶形状函数的宏粒子单位时间内在二维Yee网格上移动,分为三种涉及不同元胞数目的情形.本文列出并分析了这三类情形所对应的线性代数方程组,同时研究了每个方程组所对应的三个特解：V方法、E方法和U方法所给的电流密度公式.一个带有一阶形状函数的宏粒子单位时间内在三维Yee网格上移动,分为八种涉及不同元胞数目的情形.针对其中最常见同时也最简单的一类情形,本文列出并分析了其所对应的线性代数方程组,同时研究了该方程组所对应的两个特解：V方法、E方法所给出的电流密度公式.值得一提的是U方法给出的三维公式,不满足电荷守恒.

从本文的分析可知,我们获得了一定条件下的通解,即：所有电荷守恒的电流分配方案.从这些通用公式的推导可以看出,公式对宏粒子形状函数的阶数没有限制并且可以推广到其它网格上去.同时从通解出发,我们可以想办法得到一个“最优解”,只要“最优”标准被确定.

在如图 4所示的、涉及八个元胞的最常见且最简单的三维情形下,根据U方法,有

$ x_r = (x^n+x^{n+1})/2, \ y_r = (y^n+y^{n+1})/2, \ z_r = (z^n+z^{n+1})/2, $

$ F_{x1} = \frac{q}{\Delta t} (x_r-x^n) = \frac{q}{\Delta t} \frac{x^{n+1}-x^n}{2}, \; \; F_{x2} = \frac{q}{\Delta t} (x^{n+1}-x_r) = \frac{q}{\Delta t} \frac{x^{n+1}-x^n}{2}, $

$ F_{y1} = \frac{q}{\Delta t} (y_r-x^n) = \frac{q}{\Delta t} \frac{y^{n+1}-y^n}{2}, \; \; F_{y2} = \frac{q}{\Delta t} (y^{n+1}-y_r) = \frac{q}{\Delta t} \frac{y^{n+1}-y^n}{2}, $

$ F_{z1} = \frac{q}{\Delta t} (z_r-z^n) = \frac{q}{\Delta t} \frac{z^{n+1}-z^n}{2}, \; \; F_{z2} = \frac{q}{\Delta t} (z^{n+1}-z_r) = \frac{q}{\Delta t} \frac{z^{n+1}-z^n}{2}, $

$ i_1 = i_2 = -0.5, \ j_1 = j_2 = -0.5\ k_1 = k_2 = -0.5, $

$ W_{x1} = \frac{x^n+x_r}{2\delta x}-i_1, \ W_{x2} = \frac{x^{n+1}+x_r}{2\delta x}-i_2, W_{y1} = \frac{y^n+y_r}{2\delta y}-j_1, $

$ W_{y2} = \frac{y^{n+1}+y_r}{2\delta y}-y_2, W_{z1} = \frac{z^n+z_r}{2\delta z}-k_1, \ W_{z2} = \frac{z^{n+1}+z_r}{2\delta z}-k_2, $

$ \begin{eqnarray*} {J}^U_1& = &J_x(i+\frac{1}{2}, j, k) = \frac{F_{x1}(1-W_{y1})(1-W_{z1}) + F_{x2}(1-W_{y2})(1-W_{z2})}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(x^{n+1}-x^n)(1-y^n-y^{n+1})(1-z^n-z^{n+1})\\ &&+ \frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \\ {J}^U_2& = &J_x(i+\frac{1}{2}, j+1, k) = \frac{F_{x1}W_{y1}(1-W_{z1}) + F_{x2}W_{y2}(1-W_{z2})}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(x^{n+1}-x^n)(1+y^n+y^{n+1})(1-z^n-z^{n+1})\\ &&-\frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \\ {J}^U_3& = &J_x(i+\frac{1}{2}, j, k+1) = \frac{F_{x1}(1-W_{y1})W_{z1} + F_{x2}(1-W_{y2})W_{z2}}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(x^{n+1}-x^n)(1-y^n-y^{n+1})(1+z^n+z^{n+1})\\ & &- \frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \\ {J}^U_4& = &J_x(i+\frac{1}{2}, j+1, k+1) = \frac{F_{x1}W_{y1}W_{z1} + F_{x2}W_{y2}W_{z2}}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(x^{n+1}-x^n)(1+y^n+y^{n+1})(1+z^n+z^{n+1})\\ & &+ \frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \\ {J}^U_5& = &J_y(i, j+\frac{1}{2}, k) = \frac{F_{y1}(1-W_{x1})(1-W_{z1}) + F_{y2}(1-W_{x2})(1-W_{z2})}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(y^{n+1}-y^n)(1-x^n-x^{n+1})(1-z^n-z^{n+1})\\ &&+ \frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \\ {J}^U_6& = &J_y(i+1, j+\frac{1}{2}, k) = \frac{F_{y1}W_{x1}(1-W_{z1}) + F_{y2}W_{x2}(1-W_{z2})}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(y^{n+1}-y^n)(1+x^n+x^{n+1})(1-z^n-z^{n+1})\\ & &- \frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \\ {J}^U_7& = &J_y(i, j+\frac{1}{2}, k+1) = \frac{F_{y1}(1-W_{x1})W_{z1} + F_{y2}(1-W_{x2})W_{z2}}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(y^{n+1}-y^n)(1-x^n-x^{n+1})(1+z^n+z^{n+1})\\ &&- \frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \\ {J}^U_8& = &J_y(i+1, j+\frac{1}{2}, k+1) = \frac{F_{y1}W_{x1}W_{z1} + F_{y2}W_{x2}W_{z2}}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(y^{n+1}-y^n)(1+x^n+x^{n+1})(1+z^n+z^{n+1})\\ &&+ \frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \\ {J}^U_9& = &J_z(i, j, k+\frac{1}{2}) = \frac{F_{z1}(1-W_{x1})(1-W_{y1}) + F_{z2}(1-W_{x2})(1-W_{y2})}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(z^{n+1}-z^n)(1-x^n-x^{n+1})(1-y^n-y^{n+1})\\ &&+ \frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \\ {J}^U_{10}& = &J_z(i+1, j, k+\frac{1}{2}) = \frac{F_{z1}W_{x1}(1-W_{y1}) + F_{z2}W_{x2}(1-W_{y2})}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(z^{n+1}-z^n)(1+x^n+x^{n+1})(1-y^n-y^{n+1})\\ &&- \frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \\ {J}^U_{11}& = &J_z(i, j+1, k+\frac{1}{2}) = \frac{F_{z1}(1-W_{x1})W_{y1} + F_{z2}(1-W_{x2})W_{y2}}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(z^{n+1}-z^n)(1-x^n-x^{n+1})(1+y^n+y^{n+1})\\ &&- \frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \\ {J}^U_{12}& = &J_z(i+1, j+1, k+\frac{1}{2}) = \frac{F_{z1}W_{x1}W_{y1} + F_{z2}W_{x2}W_{y2}}{\delta x \delta y \delta z } \\ & = &\frac{q}{4\Delta t}(z^{n+1}-z^n)(1+x^n+x^{n+1})(1+y^n+y^{n+1})\\ &&+ \frac{q}{16\Delta t}(x^{n+1}-x^n)(y^{n+1}-y^n)(z^{n+1}-z^n), \end{eqnarray*} $

对于图 4所示的某一元胞,如元胞$ (i, j, k) $来说,

$ {J}^U_{1} \delta y \delta z + {J}^U_{5} \delta z \delta z + {J}^U_{9} \delta x \delta y \neq -\frac{\Delta Q_{i, j, k} }{\Delta t} = -\frac{q}{\Delta t}(S_{i, j, k}(x^{n+1}, y^{n+1}, z^{n+1} - S_{i, j, k}(x^{n}, y^{n}, z^{n} ). $

即,不满足离散的电荷连续性方程.