An improved immersed finte element particle-in-cell method for plasma simulation

TL;DR

This paper presents an enhanced immersed finite element particle-in-cell (IFE-PIC) method for plasma simulation, improving accuracy and charge conservation, suitable for complex interface problems with Cartesian meshes.

Contribution

The paper introduces a novel IFE-PIC method with a partially penalized IFE electric field solver and a new charge-conserving interpolation technique.

Findings

Enhanced accuracy of electric field computation.

Improved charge conservation in PIC implementation.

Effective handling of complex interface problems.

Abstract

The particle-in-cell (PIC) method has been widely used for plasma simulation, because of its noise-reduction capability and moderate computational cost. The immersed finite element (IFE) method is efficient for solving interface problems on Cartesian meshes, which is desirable for PIC method. The combination of these two methods provides an effective tool for plasma simulation with complex interface/boundary. This paper introduces an improved IFE-PIC method that enhances the performance in both IFE and PIC aspects. For the electric field solver, we adopt the newly developed partially penalized IFE method with enhanced accuracy. For PIC implementation, we introduce a new interpolation technique to ensure the conservation of the charge. Numerical examples are provided to demonstrate the features of the improved IFE-PIC method.

Tables

Table 1. Table 1: Comparison of the error of charge distribution.

$N$	Standard PIC		Improved PIC
	$\overline{\rho }$	$E_{\rho }$	$\overline{\rho }$	$E_{\rho }$
$1$	$-2.600000$	$35.00\%$	$-4.083332$	$2.08\%$
$4$	$-3.250000$	$18.75\%$	$-4.145832$	$3.65\%$
$16$	$-3.240625$	$18.98\%$	$-4.286457$	$7.16\%$
$64$	$-3.481812$	$12.95\%$	$-4.238281$	$5.96\%$
$256$	$-3.479996$	$13.00\%$	$-4.236654$	$5.92\%$
$1024$	$-3.482821$	$12.93\%$	$-4.240478$	$6.01\%$

Table 2. Table 2: Comparison of the L 2 superscript 𝐿 2 L^{2} norm error of electric potential.

$N$	IFE-PIC	Improved IFE-PIC
$1$	$1.868830\times {10}^{-3}$	$5.160630\times {10}^{-4}$
$4$	$1.190797\times {10}^{-3}$	$4.851096\times {10}^{-4}$
$16$	$1.208151\times {10}^{-3}$	$6.142589\times {10}^{-4}$
$64$	$1.281327\times {10}^{-3}$	$5.380597\times {10}^{-4}$
$256$	$1.284765\times {10}^{-3}$	$5.359054\times {10}^{-4}$
$1024$	$1.280203\times {10}^{-3}$	$5.403805\times {10}^{-4}$

Table 3. Table 3: The L 2 superscript 𝐿 2 L^{2} norm errors of electric potential for different mesh.

Mesh	Traditional IFE-PIC	Improved IFE-PIC
$10\times 10$	$1.112170\times {10}^{-2}$	$8.131783\times {10}^{-3}$
$20\times 20$	$3.953483\times {10}^{-3}$	$2.520218\times {10}^{-3}$
$40\times 40$	$1.280203\times {10}^{-3}$	$5.403805\times {10}^{-4}$
$80\times 80$	$5.275457\times {10}^{-4}$	$1.784664\times {10}^{-4}$
$160\times 160$	$3.553394\times {10}^{-4}$	$1.313621\times {10}^{-4}$
$320\times 320$	$1.372616\times {10}^{-4}$	$1.103237\times {10}^{-5}$
rate	$1.240231$	$1.771787$

Equations

$$m\frac{d\text{v}}{dt}=q(\mathbf{E}+\text{v}\times \mathbf{B}),$$

$$-\nabla \cdot (\beta \nabla \Phi )=\rho (\mathbf{x}).$$

$$\mathbf{E}(\mathbf{x})=-\nabla \Phi (\mathbf{x}).$$

$$-\nabla \cdot (\beta \nabla \Phi )$$

$$\Phi$$

$$\beta(\mathbf{x})=\left\{\begin{array}[]{c}\beta^{-},~{}~{}\mathbf{x}\in\Omega^{-},\\ \beta^{+},~{}~{}\mathbf{x}\in\Omega^{+},\end{array}\right.$$

$${\left[\left[\Phi \right]\right]}_{\Gamma }=0,$$

$${\left[\left[\beta \frac{\partial \Phi }{\partial n}\right]\right]}_{\Gamma }=0.$$

$$\phi_{i,T}(\mathbf{x})=\left\{\begin{array}[]{lr}\phi_{i,T}^{+}(\mathbf{x})=a_{i}^{+}x+b_{i}^{+}y+c_{i}^{+},~{}~{}\mathbf{x}\in T^{+},&\\ \phi_{i,T}^{-}(\mathbf{x})=a_{i}^{-}x+b_{i}^{-}y+c_{i}^{-},~{}~{}\mathbf{x}\in T^{-},&\end{array}\right.\forall~{}i=1,2,3.$$

$$\phi_{i,T}(A_{j})=\delta_{ij}=\left\{\begin{array}[]{c}1,~{}~{}\text{if}~{}i=j,\\ 0,~{}~{}\text{if}~{}i\neq j,\\ \end{array}\right.~{}~{}~{}\forall~{}i,j=1,2,3.$$

$${\varphi }_{i,T}^{+}(D)={\varphi }_{i,T}^{-}(D),{\varphi }_{i,T}^{+}(E)={\varphi }_{i,T}^{-}(E),\forall i=1,2,3.$$

$${\beta }^{+}\frac{\partial {\varphi }_{i,T}^{+}}{\partial \text{n}}={\beta }^{-}\frac{\partial {\varphi }_{i,T}^{-}}{\partial \text{n}},\forall i=1,2,3,$$

$${\tilde{\mathbb{P}}}_1(T)=span\{{\varphi }_{1,T},{\varphi }_{2,T},{\varphi }_{3,T}\}.$$

$${\Phi }_i({\mathbf{x}}_j)={\delta }_{ij},\forall i,j=1,2,\cdots ,N,$$

$${\Phi }_i{|}_T\in {\tilde{\mathbb{P}}}_1(T),\forall T\in {\mathcal{T}}_h^i,\text{and}{\Phi }_i{|}_T\in {\mathbb{P}}_1(T)=span\{1,x,y\},\forall T\in {\mathcal{T}}_h^n.$$

$$a(u_h,v_h)=(f,v_h),\forall v_h\in S_h,$$

$$a(u,v)=\sum _{T\in {\mathcal{T}}_h}{\int }_T\beta \nabla u\cdot \nabla vd\mathbf{x},(f,v)={\int }_{\Omega }fvd\mathbf{x}.$$

$$q_{i,j}=q_P\frac{S(X_{i+1,j+1},X_P)}{S(X_{i,j},X_{i+1,j+1})},$$

$$\mathbf{E}(x,y)=-\nabla \Phi (x,y).$$

$$(E_x{)}_{i,j}=\frac{{\Phi }_{i-1,j}-{\Phi }_{i+1,j}}{2\Delta x},(E_y{)}_{i,j}=\frac{{\Phi }_{i,j-1}-{\Phi }_{i,j+1}}{2\Delta y},$$

$$\left\{\!\!\left\{u\right\}\!\!\right\}_{e}=\frac{1}{2}\Big{(}(u|_{T_{e,1}})_{e}+(u|_{T_{e,2}})_{e}\Big{)},~{}~{}~{}~{}\left[\!\left[u\right]\!\right]_{e}=(u|_{T_{e,1}})_{e}-(u|_{T_{e,2}})_{e}.$$

$$a_{ϵ}(u_h,v_h)=(f,v_h),\forall v_h\in S_h,$$

$$a_{ϵ}(u,v)$$

$$q_i^{(1)}=w_i{q}_{{}_P},i=1,2,3,4,$$

$$w_1=\frac{A_4}{A},w_2=\frac{A_3}{A},w_3=\frac{A_2}{A},w_4=\frac{A_1}{A}.$$

$$q_1^{(2)}$$

$$q_2^{(2)}$$

$$q_3^{(2)}$$

$$q_4^{(2)}$$

$$q_1^{(2)}$$

$$q_2^{(2)}$$

$$q_3^{(2)}$$

$$q_4^{(2)}$$

$$q_1^{(2)}$$

Full text

††journal: Computers $\&$ Mathematics with Applications

1 Introduction

There are two classes of methods for plasma simulation. The first one is the traditional dynamic simulation [33], which is mainly used to obtain the distribution function of particles change with time by solving the Vlasov equation. The second one is the particle simulation method [9, 10, 17], which is used to track a large number of individual particles and to obtain the trajectory parameters and characteristics of plasma by statistical methods. Due to the enormous number of particles to be tracked and the limited computational resources, the development of the particle simulation method was quite slow. The particle simulation method entered into a rapid developing period [2, 21, 32] since Birdsall and Langdon [22] introduced the particle-in-cell (PIC) method which utilizes the finite-sized particle (or cloud) instead of a huge number of real particles.

Immersed finite element (IFE) method is a finite element method for solving interface problems on uniform Cartesian meshes, which was first proposed by Li et al [23]. Different from classical finite element methods using body-fitted meshes, the mesh of IFE method is independent of the interface. However, the IFE basis functions around interfaces are modified to accommodate the interface jump conditions. The advantage of IFE method is that structured Cartesian meshes can be used to solve interface problems with arbitrary interface geometry. For problems with a moving interface, IFE methods are especially advantageous since there is no need to regenerated the solution meshes repeatedly [16, 26, 27]. The IFE methods have been developed for solving the second-order elliptic equations [14, 15, 24, 25, 30], elasticity equations [11, 29, 31], Stokes equations [1], to name only a few.

In the past decade, IFE method has been used together with PIC method for plasma simulations [5, 8, 19, 20]. IFE method used as an electric field solver is performed on well-structured Cartesian meshes. This is particularly desirable for PIC method because tracking a large number of plasma particles can be efficiently done in a uniform structured mesh. We refer to a few recent applications of IFE-PIC method for different types of particle simulations, such as ion thruster [3, 18], hall thruster [4], and lunar surface environment [12, 13]. Also, IFE-PIC method has been extended to handle unbounded interface problems with asymptotic boundary condition [7] and periodic boundary condition [6].

For the current IFE-PIC method, we noticed that there are two issues. First, the classical Galerkin IFE method is used as the field solver. As shown in [28], the classical IFE method is not accurate around the interface, because the IFE basis functions are discontinuous across the element boundaries, and the classical Galerkin formulation cannot control such discontinuity. Second, the particle interpolation method of PIC algorithm is imperfect. The conventional interpolation approach applied on the interface element often leads to the non conservation of charge, because it neglects the fact that some nodes of the interface cells are inside the conductors. Similar problems occur in applying the electric field force to the particles on the interface elements.

In this paper, we introduce an improved IFE-PIC method that focuses on overcoming the problems mentioned above. As a remedy of discontinuity of IFE field solver, we adopt the newly developed partially penalized immersed finite element (PPIFE) method [28] to improve the accuracy of IFE methods near interfaces. For PIC interpolation, we introduce a two-step approach for particle interpolation that preserves the charge conservation. Comparing with conventional charge distribution in PIC, we add a correction step that redistributes the quantity distributed to the nodes inside the conductor to the remaining nodes in order to maintain charge conservation. In addition, we use IFE basis functions to calculate the electric field and force on the interface elements. The new approach can calculate the motion of particles more accurately.

The rest of the article is organized as follows. In Section 2, we recall the classical IFE method and PIC interpolation. In Section 3, we present our improved IFE-PIC method. The improvement in the IFE solver part is the PPIFE method with additional penalty terms. The improvement in the PIC part includes the new particle interpolation scheme and the new method for the force deposit. In Section 4, we present some numerical experiments to compare the performance of traditional IFE-PIC method and improved IFE-PIC method. Brief conclusions will be given in Section 5.

2 Review of IFE-PIC Method

In this section, we first recall the main steps in a typical IFE-PIC computational cycle. Then we will recall the classical IFE method and PIC interpolation that are widely used in the literature.

2.1 Main Steps of IFE-PIC Method

Real plasma particles are modeled as many macro-particles in the PIC method, and they follow the evolution of the orbits of individual particles in the self-consistent electromagnetic field. The field is then updated by solving the governing elliptic equation with discontinuous dielectric coefficients. The IFE-PIC method is an iteration of solving for the electromagnetic field and particle motion until the steady state is achieved.

In general, an IFE-PIC computational cycle consists of the following five steps:

Step 1. Initialization

A series of initial settings of the simulation including domain, mesh, boundary condition, and initial position and velocity of particles must be set up.

Step 2. Particle Push

The motion of the particle is induced by the particles themselves and the applied external fields $\bm{E}$. The trajectory of an individual charged particle is obtained by integrating the Newton-Lorentz equation

[TABLE]

where $m$, $q$ and v are the mass, the charge, and the velocity of the particle, respectively. $\bm{B}$ denotes the static magnetic field.

Step 3. Charge Deposit

The change of the positions of particles leads to the change of the charge density $\rho$ on each node. Thus, we need to calculate the charge density at each node according to the new positions of particles. The process of interpolating the particle charges on the discrete mesh points is called weighting. In traditional PIC method, we use the ratio of the area of the rectangle formed by the opposite cell vertex and the particle to the area of element as the weighting.

Step 4. Solving for Potential

After obtaining the charging density $\rho$ on each node, the electric field should also be updated. The potential function $\Phi(\mathbf{x})$ can be described by the second-order Poisson’s equation with discontinuous dielectric coefficient $\beta(\mathbf{x})$, which represents different types of material:

[TABLE]

To solve this equation, we use IFE method as a field solver, because of its applicability of Cartesian mesh, which is desirable in the PIC simulation for fast tracking of particles’ locations.

Step 5. Force Deposit

After solving for the potential $\Phi(\mathbf{x})$, we calculate the electric field $\bm{E}$ at each node by

[TABLE]

Next, we need to deposit the electric field at nodes to the particles with arbitrary positions. The electric field of each node can be obtained using two-point difference method from the potential $\Phi(\mathbf{x})$ on each element node. The electric field of each particle can be obtained by interpolating.

We note that, in step 3, when the particle is located in an interface element, the current interpolation using the simple area-weighting will result in the non-conservation of charge. Also, in step 4, the classical IFE method may not be as accurate around interface as the rest of the domain due to the discontinuity of IFE basis functions. Finally, in step 5, the force deposit is not accurate in the interface element due to part of the element is the conductor. In this paper, our improvement of the IFE-PIC method will mainly focus on these three steps.

2.2 Classical IFE Method for Interface Problems

The electric field $\Phi(\mathbf{x})$ is assumed to be governed by the following second-order elliptic equation

[TABLE]

Here, we assume that $\Omega\subset\mathbb{R}^{2}$ is a rectangular domain separated by an interface curve $\Gamma$ into two sub-domains $\Omega^{-}$, $\Omega^{+}$ such that $\overline{\Omega}=\overline{\Omega^{+}\cup\Omega^{-}\cup\Gamma}$. See Figure 1 as an illustration.

The coefficient $\beta(\mathbf{x})$ is discontinuous across the interface. Without loss of generality, we assume $\beta$ is a piecewise constant function as follows

[TABLE]

where $\mathbf{x}=[x,y]^{t}$. Across the interface $\Gamma$, the following interface jump conditions are satisfied:

[TABLE]

[TABLE]

Let $\mathcal{T}_{h}=\{T\}$ be a uniform triangular or rectangular mesh of the domain $\Omega$ with size $h$. If an element $T$ is cut through by the interface, it is called an interface element; otherwise, it is said to be a noninterface element. The sets of interface elements and noninterface elements are denoted by $\mathcal{T}_{h}^{i}$ and $\mathcal{T}_{h}^{n}$, respectively. Standard linear or bilinear finite element functions are used on all noninterface elements. Special piecewise-polynomial basis functions are constructed on interface element to accommodate the interface conditions. To be more precise, we use the linear IFE method as an example. Assume that $T=\triangle A_{1}A_{2}A_{3}$ is an interface triangle, and the interface curve $\Gamma$ intersects $T$ at two points, denoted by $D$ and $E$. The element $T$ is divided into two sub-elements $T^{+}$ and $T^{-}$ by the line segment $\overline{DE}$. See Figure 2 as an illustration. The local IFE basis functions $\phi_{i,T}$, $i=1,2,3$ on $T$ are defined as follows:

[TABLE]

They satisfy the following condition:

1. nodal-value conditions

[TABLE]

2. function-value continuity

[TABLE]

3. flux continuity

[TABLE]

where n is the normal vector of $\overline{DE}$.

The local linear IFE space is defined as

[TABLE]

Denote the set of interior nodes on $\mathcal{T}_{h}$ by $\mathcal{N}_{h}$. On each node $\mathbf{x}_{i}$, $i=1,2,\cdots,N$, we define the global linear IFE basis function $\Phi_{i}$ such that

[TABLE]

and

[TABLE]

The global linear IFE space is formed as $S_{h}=span\{\Phi_{i}:i=1,2,\cdots,N\}$.

The classical (Galerkin) IFE method is to find $u\in S_{h}$ such that

[TABLE]

where

[TABLE]

The construction of the bilinear IFE spaces on rectangular meshes are similar, and we refer to [14, 15, 31] for more details.

2.3 Particle Interpolation in Traditional PIC Method

In the PIC method, particle interpolation is required at Step 3 and Step 5 described in Section 2.1. First, in charge deposit in Step 3, as shown in Figure 3(a), the portion of the total particle charge assigned to a certain cell vertex is proportional to the area of the rectangle formed by the opposite cell vertex and the particle. Thus, the interpolation of the particle $P$, located at position $X_{P}$, to the node $X_{i,j}$ can be calculated by

[TABLE]

where $q_{i,j}$ is the amount of charge of node $(i,j)$, $q_{P}$ is the amount of charge of particle $P$, and $S(A,B)$ is the area of the rectangle whose diagonal is $\overline{AB}$.

On non-interface elements, the interpolation (18) works well. However, on interface elements, it will face some difficulty, as shown in Figure 3(b), the charge should not be assigned to the node $X_{i+1,j}$, which is inside the conducting object. Directly applying (18) on interface elements will cause the non-conservation for the total charge and charge density.

Moreover, in Step 5 the electric field is required to be interpolated at the particles’ positions, so that the electric field force can be obtained. The potential $\Phi$ at grid points are solved by IFE method with the appropriate boundary conditions. Then, the electric field $\bm{E}=[E_{x},E_{y}]^{t}$ can be obtained from the potential $\Phi$ using the following equation

[TABLE]

The conventional approach for this interpolation uses the finite difference form for $E_{x}$ and $E_{y}$:

[TABLE]

where $\Delta x$ and $\Delta y$ are mesh sizes in the $x$- and $y$- directions, respectively. After the electric field is obtained, the forces caused by the field at mesh nodes can be deposited to the arbitrary particles positions.

When a particle is located in the interface element, the calculation (20) is apparently inaccurate, because of the electric field of the four points that located at the different locations of the object are discontinuous.

3 Improved IFE-PIC Method

In this section, we present an improved IFE-PIC method. The improvement involves both the IFE solver and the PIC particle interpolation.

3.1 Partially Penalized IFE Method

Due to the discontinuity of IFE basis functions, classical IFE method using Galerkin formulation may generate large errors around the interface. Partial penalized immersed finite element (PPIFE) method is introduced in [28] that is known to greatly improve the accuracy of IFE solution around interface. The main idea of this method is to add penalty terms on interface edges to reduce the negative impact of the discontinuity introduced by IFE functions. In our improved IFE-PIC method, we adopt this new PPIFE scheme as the new field solver.

To present the PPIFE method, we need to introduce a few notations. Let $\mathcal{E}_{h}$ be the set of all interior edges of the mesh $\mathcal{T}_{h}$. If an edge $e\in\mathcal{E}_{h}$ intersects with the interface curve $\Gamma$, we call it an interface edge; otherwise a noninterface edge. The sets of interface edges and noninterface edges are denoted by $\mathcal{E}_{h}^{i}$ and $\mathcal{E}_{h}^{n}$, respectively. For each interior edge $e\in\mathcal{E}_{h}$, it must be shared by two adjacent elements, denoted by $T_{e,1}$ and $T_{e,2}$. For a function $u$ defined on $T_{e,1}\cup T_{e,2}$, the average and jump of $u$ on $e$ are defined as follows

[TABLE]

The PPIFE method for solving (4) - (8) is to find $u_{h}\in S_{h}$ such that

[TABLE]

where

[TABLE]

Here, $\epsilon$ can be chosen as $-1$, [math], or $1$, which corresponds to symmetric, incomplete, or non-symmetric PPIFE methods, respectively. $\sigma_{e}$ is a positive penalty parameter. As shown in [28], for non-symmetric PPIFE method, the scheme (22) is stable as long as $\sigma_{e}>0$ on every edge $e$. For symmetric and incomplete PPIFE methods, (22) is stable when $\sigma_{e}$ is large enough.

3.2 Charge-Conservative Particle Interpolation

As shown in Section 2, the particle interpolation technique (18) is not accurate for particles in the interface elements. Also, the total charge is not conservative due to part of the nodes are inside the objects such as metal or ceramic materials. In this subsection, we introduce a two-step interpolation that can maintain the charge conservation.

3.2.1 Improved Algorithm for Charge Deposit

We note that if a particle is in a noninterface element, it is suitable to use the standard charge distribution (18) with area weight method. If a particle $P$ is in an interface element, we introduce a correction step to redistribute the charge assigned to nodes inside the conductor to other nodes. This step ensures the conservation of the charge in the simulation domain. The complete two-step interpolation is as follows.

Step 1: Initial distribution

Interpolate the charge $q_{{}_{P}}$ of the particle to four nodes of the element using the standard area-weight approach:

[TABLE]

where

[TABLE]

Here $A_{i}$, $i=1,2,3,4$ are the areas of small rectangles, illustrated in Figure 4.

Step 2: Re-distribution

Correct the charge distribution to nodes inside the object.

• 1.

Case 1: There is only one node (e.g. node $\#$2) in the object, as shown in Figure 4(a). The charge initially distributed to node $\#$2 should be redistributed to the other three nodes. The correction procedure is the following:

[TABLE]

• 2.

Case 2: If two nodes are inside the object, as shown in Figure 4(b), then the correction procedure becomes

[TABLE]

• 3.

Case 3: If three nodes are inside the object, as shown in Figure 4(c), then

[TABLE]

We note that if we use the traditional PIC charge interpolation (without correction) on interface elements, the total charge of nodes is less than the charge of particle, that is

[TABLE]

where $\Lambda$ denotes the indices of nodes outside the object. However, adding the correction, the interpolation satisfies that

[TABLE]

This clearly shows that the new interpolation scheme preserves the charge-conservation.

Moreover, this new interpolation scheme is robust regarding the location of the particle. That is, if the particle touches the boundary of the interface element, the algorithm is still valid.

3.2.2 Improved Algorithm for Force Deposit

The potential distribution at grid points is obtained by solving linear algebraic equations contained in the immersed finite element method, and then for the electric potential at point $(x,y)$ in an element $T$ can be approximated using the following equation:

[TABLE]

Here, $u_{i}$ are the numerical solutions of the potential at the vertices of $T$, and $\phi_{i,T}$ are local IFE basis functions. The degree of freedom $d=3$ for linear IFE method, and $d=4$ for bilinear IFE method.

Suppose the electrical field is in $\Omega^{+}$, and the conductor is in $\Omega^{-}$. Then the electric field E at a particle $P(x,y)$ can be obtained by

[TABLE]

For points located in $T^{-}$, the electric field can be calculated in the same way. This method is more accurate than linear interpolation when the internal and external potential of objects are discontinuous. With the new PPIFE field solver, and the improved PIC interpolation technique on particles and force, the workflow of our improved IFE-PIC algorithm can be summarized in Figure 5.

4 Numerical Examples

In this section, we present some numerical examples to demonstrate the features of the improved IFE-PIC method.

We set up a test problem of conducting cylinder with background plasma particles. Let the two-dimensional simulation domain be $\Omega_{s}=[-1,1]\times[-1,1]$. The center and the radius of the cylinder are set to be $O=[0,0]$ and $r_{0}=\pi/12$, respectively, as shown in Figure 6. For plasma particles, we load $1635841~{}(=1279^{2})$ uniformly distributed particles $P_{i,j}$ into the simulation domain and the positions of particles are given by

[TABLE]

Assume that the loaded particles are electrons, and the charge density on the mesh point is $\rho=-4$, so we can get the charge of each particle. Then we remove all the particles inside the cylinder. We choose appropriate Dirichlet boundary conditions and source function so that the analytical solution of this simulation problem is

[TABLE]

where $r=\sqrt{x^{2}+y^{2}}$, and $[\beta^{-},\beta^{+}]=[1,10]$. The analytical solution (30) of the potential in the simulation domain is shown in Figure 6.

4.1 Comparison of PIC Interpolations

We first compare the performance of traditional PIC charge interpolation with the new interpolation method. The charge densities on mesh points can be obtained by depositing the physical quantities (charge of electrons) from the particle locations onto mesh nodes. Figure 7 shows the charge density distribution with standard charge deposit algorithm on a $20\times 20$ Cartesian mesh. It can be shown that all the charge of interface element is less than the prescribed value $\rho=-4$, because the charge of particles are deposited to several mesh points inside the object.

Using the improved PIC algorithm introduced in Section 3.2, the charge density distribution in simulation domain $\Omega_{s}$ is shown in Figure 8. On the right side of Figure 8, we show a zoom-in plot around interface. It can be shown that the charge on the nodes of interface elements are not unanimously less than $-4$. On some nodes (e.g. [-0.3,0.1], and [-0.3,0.2]), the change quantities are actually greater than $-4$. This is because of the redistribution step we added, so that the charge of particles will not be deposited to any mesh points inside the conductor. Moreover, this new interpolation algorithm ensures the conservation of the charge in the calculation domain.

4.2 Comparison of IFE-PIC Performance

Next, we combine the new PIC algorithm with the new electrical field solver, PPIFE scheme. We compare the performance of traditional and improved IFE-PIC schemes. In Figure 9 and Figure 10, we plot the numerical errors of electrical field by traditional and improved IFE-PIC methods, respectively. It can be easy observed that the new IFE-PIC method has much smaller error than traditional IFE-PIC method especially around the interface.

4.3 Sensitivity of the Number of Particles

In this test, we focus on the sensitivity of our IFE-PIC method to the number of particles. Note that on non-interface elements, regardless of the number of particles within each element, the final charge of interpolation to nodes are same, i.e., $\rho=-4$. However, on interface elements, the element is split into two parts. The particles located in the conductor will not contribute to the charge distribution, as a result, the number of particles in an interface element should have greater influence on the charge distribution results.

Figure 11 demonstrate an interface element with $1$, or $16$, or $256$ equally distributed particles, respectively. As shown in Figure 11, if the number of particles is small, the charge distribution will rely more on the random fall of individual particles. On the other hand, if the number of particles is large, the effect of a single particle could be neglected. This is consistent with the actual situation.

Next, we compute the average density for different number of particles with traditional and improved PIC interpolations. We use $N$ to denote the number of particles in each interface element. We let $\mathbf{x}_{i}$, $i=1,\cdots N_{I}$ be the nodes that belongs to at least one interface element. The mean value of density at interface nodes is computed as follows

[TABLE]

where the $\rho_{h}(\mathbf{x}_{i})$ is the interpolated charging density at the point $\mathbf{x}_{i}$. The percentage error of the density is defined by

[TABLE]

The $L^{2}$ norm error of the potential function is defined by

[TABLE]

Table 1 and Table 2 report the error of charging density and error of electric potential using two IFE-PIC methods. These tests are conducted on a uniform $40\times 40$ Cartesian mesh. Table 1 clearly indicates that as we increase the number of particles belong each element, the traditional PIC interpolation become more accurate, this indicates that the traditional approach is very sensitive to the number of particles. On the other hand, the new interpolation technique is much more accurate than the traditional approach, and are more robust with respect to the number of particles. Table 2 shows that the improved IFE-PIC scheme is about three times more accurate than the classical scheme in the $L^{2}$ norm of the electric potential.

Remark 4.1

Based on our current numerical experiments, simply improving only one part of IFE solver or PIC interpolation is less significant, and the behavior on accuracy varies in different choices of interface geometry, conductivity coefficient, and the number of particles.

4.4 Convergence of Potential on a Sequence of Meshes

In this example, we test convergence of potential on a family of uniform meshes with the same number of particles in each element. We consider the number of particles is the same as in Section 4.1. Table 3 reports the $L^{2}$ norm error of the potential. From the comparison, we can see that our improved IFE-PIC scheme is much more accurate on every mesh level than the widely used scheme.

5 Conclusions

In this paper, we proposed a new IFE-PIC method for plasma simulation. The new method has improvement in both IFE solver and PIC interpolation. One prominent feature is that our new PIC interpolation has the charge-conservation property. Moreover, the improved IFE solver with partial penalty terms produce more accurate approximation around interface.