# Analysis of the Calculation of a Plasma Sheath Using the Parallel SO-DGTD Method

Authors: Qian Yang; Bing Wei; Linqian Li; Debiao Ge
Journal: International Journal of Antennas and Propagation
 (2019)

---

## Abstract

The plasma sheath is known as a popular topic of computational electromagnetics, and the plasma case is more resource-intensive than the non-plasma case. In this paper, a parallel shift-operator discontinuous Galerkin time-domain method using the MPI (Message Passing Interface) library is proposed to solve the large-scale plasma problems. To demonstrate our algorithm, a plasma sheath model of the high-speed blunt cone was established based on the results of the multiphysics software, and our algorithm was used to extract the radar cross-section (RCS) versus different incident angles of the model.


---

## Body

## 1. Introduction

A high-speed aircraft in the atmosphere is surrounded by a plasma sheath particularly during hypersonic reentry into the earth’s atmosphere. As a result, aerodynamic heating has a great effect on the ionization of aircraft surface air [1–3]. The transmission signal between the aircraft and the ground base or the relay satellite could be interfered by the plasma sheath [4–7], resulting in the phenomenon of a “black barrier”. Thus, it is helpful to study the characteristics of the plasma sheath for tracking of and communication with the aircraft.The discontinuous Galerkin time-domain (DGTD) method offers superior versatility in geometry discrimination by using unstructured meshes; unlike the finite difference time-domain method (FDTD) [8, 9], the DGTD method has no numerical dispersion. DGTD is appropriate for solving complex electromagnetic problems [10–15]. The shift-operator DGTD (SO-DGTD) method was proposed for dispersive media [16–18]; to address the frequency domain constitutive relation and obtain the time-domain iteration formula, the shift-operator is used. Plasma sheath is a problem that requires a huge consumption of computer resources to solve, resulting in a huge resource consumption in time-domain simulation. In response, the parallel technology can improve the efficiency while allowing simulation of large-scale problems. References involving the study of the parallel DGTD algorithm based on MPI or CUDA can be found in the literature. In 2006, Marc Bernacki used MPI-DGTD to analyze the plane and the distribution of SAR of a human head [19]. The electromagnetic wave interaction with biological tissues has been simulated by GPU DGTD [20]. In [21], the MPI/GPU implementation of DGTD was proposed to simulate an SRR device, and the parallel efficiency of this algorithm was analyzed. Stylianos Dosopoulos proposed MPI-DGTD based on nonconformal meshes and the local time-stepping strategy for simulation of the IC package [22]. A parallel nonconforming DGTD based on MPI was presented for the simulation of metallic nanoparticles by R. Léger [23]. The above case studies do not have the detailed explanation in the MPI scheme, especially data exchange between processes. As far as the development of DGTD is concerned, an effective and extensible parallel scheme will provide certain help for later research work, such as complex media and the heterogeneous parallel system.According to the reference, the study of the parallel DGTD algorithm for the plasma sheath target has not been reported yet. In this paper, we focus on the study of DGTD in plasma medium. First, we introduce the SO-DGTD method; next, a specific parallel scheme for SO-DGTD is given, and numerical examples are given to illustrate parallel efficiency; finally, the proposed algorithm is applied to the analysis of the scattering characteristic of a hypersonic aircraft. The results show that the sheath of the ultrahigh-speed aircraft can significantly affect its RCS over a wide frequency band.


## 2. SO-DGTD Formulation For Plasma Medium

The semidiscrete equations of DGTD are(1)Me∂∂tDe-SeHe-Fhe=0μMe∂∂tHe+SeEe+Fee=0where μ is the permeability and E and Η are the electric and magnetic field, respectively. D is the electric flux density. e is the index of element located in the computational region, M is the mass matrix, S is the stiffness matrix, and Fh and Fe are the flux matrices.The constitutive relation of plasma is defined by the Drude model(2)Dω=ε0εrωEωεrω=1+ωp2ωjνc-ωwhere ε0 is the permittivity of vacuum, εr(ω) is the relative permittivity (which is dependent on frequency), ωp is the plasma frequency, νc is the collision frequency, and ω is the angular frequency of electromagnetic waves. Equations (1) allow us to update D by using H and update H by using E. Obviously, we need a formula from D→E to complete the cycle of H→D→E→H. Bing Wei proved that the complex permittivity of the general dispersive media model can be described by a rational polynomial fraction in jω [24]:(3)Dω=ε0∑l=0Npljωl∑l=0MqljωlEωwhere M and N are the maximum powers of the fractional polynomial and pl and ql are the expansion coefficients of polynomials. According to (3), by using the relationship of frequency-time-domain (jω⇒∂/∂t), we can obtain(4)Dω=ε0∑l=0Npljωl∑l=0MqljωlEω→Dn=ε0∑l=0Npl∂/∂tl∑l=0Mql∂/∂tlEnwhere n denotes the time step. The relation of the time-domain differential operator and the shift-operator zt is as follows.(5)∂∂tl→2Δt·zt-1zt+1ll=1,2,…The final time-domain form corresponding to (3) is [16](6)∑l=0Mqlhlzt+1N-lzt-1lDn=ε0∑l=0Nplhlzt+1N-lzt-1lEnwhere zt represents the shift-operator, defined by ztfn=fn+1. The role of zt is to shift the time sampling from thenth step to the (n+1)-th step. h=2/Δt, where Δt is the time increment in the calculation.By using the characteristic ofzt, we obtain the iteration formula from (6)(7)Een=1b0∑l=0Malε0Den-l-∑l=1NblEen-lFor the Drude model (plasma), the expansion coefficients are described as follows.(8)M=N=2p0=ωp2,p1=νc,p2=1,q0=0,q1=νc,q2=1a0=q0+q1h+q2h2,a1=2q0-2q2h2,a2=q0-q1h+q2h2b0=p0+p1h+p2h2,b1=2p0-2p2h2,b2=p0-p1h+p2h2To match (7), the leap-frog scheme is used to solve (1).


## 3. MPI Parallelization of SO-DGTD

The proposed algorithm contains three types of variables, namely, electric fieldE, magnetic field H, and electric displacement field D. In our algorithm, E and H are the communications data between neighboring elements. {Fhe} and {Fee} in (1) can be written as below.(9)Fhe=kHeMfe·He+-He-vEeMge·Ee+-EeFee=kEeMfe·Ee+-Ee+vHeMge·He+-HeFor elemente, we denote the coefficients of flux by kHe, vEe, kEe, and vHe. Mfe and Mge denote the flux matrices. Ee and He denote the fields of element e; Ee+ and He+ denote the fields of adjacent elemente+. In parallel processing, each partition is mapped to an MPI process, and the main task is the communication of data. If e ande+ belong to the same process, then the data can be read from memory directly; otherwise, the communications must be handled by MPI function for the neighboring processes. The flow chart in Figure 1 shows the activities of the proposed algorithm.Figure 1
Flow chart.


### 3.1. Mesh Partition

As shown in Figure2, a sphere is partitioned into three subdomains corresponding to process 0, process 1 and process 2. Gmsh is used to provide the tetrahedral meshes and high quality partitions [25].Figure 2
Mesh partition.For each pair of neighboring submeshes, the boundary has no buffer elements. All communications between neighboring processes are handled with data packaging and transmitting technology, as shown in Figure3 and discussed in the next section.Figure 3
Data packaging and transmitting.


### 3.2. Communication and Calculation Setup

The key problem of iterative calculation is addressing the issue of communications between neighboring processes. As shown in Figure4, tetrahedron P belongs to process M, and the neighboring tetrahedron R also belongs to process M, but another neighboring tetrahedron Q belongs to process N. Thus, the data among P and R can be read directly, whereas the data among P and Q must be communicated between process M and N. Communication of parallel computing is implemented via the mapping of meshes.Figure 4
Mapping of meshes.Figure4 illustrates the communication between neighboring processes. In addition, the adjacency relations among all processes and elements (tetrahedrons) must be established, as shown in Table 1.Table 1
Topology table.
Current process Neighboring process Neighboring elements M N (P, Q)  … … …In the calculation of DGTD, we need to establish a topology table for each processes, and the topology tables are extracted from meshes. The send buffers and receive buffers are established based on the topology tables. This part is particularly important for parallel computation. MPI_SENDRECV() function is used to implement the communication. Part of the code is as follows:Do i=1, N_outn=N_out_NODE(i)CALL MPI_SENDRECV (E_cache_out(1:8,1:Num_NODE(n),n), 8∗Num_NODE(n),MPI_REAL, m1, 99+n, E_cache_in(1:8, 1:Num_NODE(n), n), 8∗Num_NODE(n),MPI_REAL, m1, 99+MYTIDS, MPI_COMM_WORLD, STATUS, IERR)EnddoIn the code above, MPI_SENDRECV() function is only running on the neighboring processes. This code can achieve the desired result and efficiency.Subsequently, these procedures are packaged into a solution wrapper, as shown in Figure5. It should be noted that the parallel scheme is directed to dispersive media, and the scheme calculates the D fields first, followed by the E fields and H fields lastly. This process repeats many times until the program ends. This scheme has the advantages of openness, convenience, and well expansibility.Figure 5
Calculation strategy.


### 3.3. Parallel Efficiency

Example 1.
The monostatic RCS of a plasma sphere is simulated to show the parallel efficiency. The plasma sphere has a radius of 1 m (ωp=6.75×108rad/s, νc=7.5×107Hz), and the mesh size is 0.05 m. Three sets of platform with different numbers of processes are used to test the program (Figure 6). The computing platform is the Xidian high performance computing center.
Good agreements are observed between the numerical results and the Mie analytical solution, as shown in Figure6. Compared with the Mie solution, the root-mean square error of our algorithm is 1.5370 dB.
In Figure7, as a benchmark we chose 24 processes, and the parallel efficiency of 72 processes is 90.1%. The proposed algorithm has good performance in this example.Figure 6
Monostatic RCS of plasma sphere.Figure 7
Parallel efficiency.Example 2.
The previous example shows that the proposed algorithm has good practicability and high efficiency. Next, we calculate the monostatic RCS of a plasma-coated PEC sphere. The sphere of radius 0.5 meters is coated with plasma, as shown in Figure8. The plasma layer has a thickness of 0.2 meters. The plasma frequency and the collision frequency of the plasma are ωp=2.0×109rad/s and νc=0.1×109Hz, respectively.
The model is discretized into 3276319 tetrahedrons. Three sets of platform with different numbers of processes are used, and the numbers are 24, 48, and 72, respectively. The computing platform is the Xidian high performance computing center.
The Mie theory is an analytical method for the sphere target, and the Mie solution is used as a reference value in this example, as shown in Figure9 (cyan line). The monostatic RCS results of different sets are also shown in Figure 9, where the black line, red dash, and blue dot represent the results of 24 processes, 48 processes, and 72 processes, respectively. The Mie solution is consistent with the numerical results.
As seen in Figure10, the parallel efficiency of 72 processes is 91.50%. The proposed algorithm still has good performance in this example.
The calculation results illustrate the proposed parallel scheme is effective in parallel processing, and the results suggest the parallel SO-DGTD method can provide reliable results. The parallel SO-DGTD method can be used to analyze the plasma sheath target, as discussed in the next section.Figure 8
Plasma-coated PEC sphere.Figure 9
Monostatic RCS of the plasma-coated PEC sphere.Figure 10
Parallel efficiency.


## 3.1. Mesh Partition

As shown in Figure2, a sphere is partitioned into three subdomains corresponding to process 0, process 1 and process 2. Gmsh is used to provide the tetrahedral meshes and high quality partitions [25].Figure 2
Mesh partition.For each pair of neighboring submeshes, the boundary has no buffer elements. All communications between neighboring processes are handled with data packaging and transmitting technology, as shown in Figure3 and discussed in the next section.Figure 3
Data packaging and transmitting.


## 3.2. Communication and Calculation Setup

The key problem of iterative calculation is addressing the issue of communications between neighboring processes. As shown in Figure4, tetrahedron P belongs to process M, and the neighboring tetrahedron R also belongs to process M, but another neighboring tetrahedron Q belongs to process N. Thus, the data among P and R can be read directly, whereas the data among P and Q must be communicated between process M and N. Communication of parallel computing is implemented via the mapping of meshes.Figure 4
Mapping of meshes.Figure4 illustrates the communication between neighboring processes. In addition, the adjacency relations among all processes and elements (tetrahedrons) must be established, as shown in Table 1.Table 1
Topology table.
Current process Neighboring process Neighboring elements M N (P, Q)  … … …In the calculation of DGTD, we need to establish a topology table for each processes, and the topology tables are extracted from meshes. The send buffers and receive buffers are established based on the topology tables. This part is particularly important for parallel computation. MPI_SENDRECV() function is used to implement the communication. Part of the code is as follows:Do i=1, N_outn=N_out_NODE(i)CALL MPI_SENDRECV (E_cache_out(1:8,1:Num_NODE(n),n), 8∗Num_NODE(n),MPI_REAL, m1, 99+n, E_cache_in(1:8, 1:Num_NODE(n), n), 8∗Num_NODE(n),MPI_REAL, m1, 99+MYTIDS, MPI_COMM_WORLD, STATUS, IERR)EnddoIn the code above, MPI_SENDRECV() function is only running on the neighboring processes. This code can achieve the desired result and efficiency.Subsequently, these procedures are packaged into a solution wrapper, as shown in Figure5. It should be noted that the parallel scheme is directed to dispersive media, and the scheme calculates the D fields first, followed by the E fields and H fields lastly. This process repeats many times until the program ends. This scheme has the advantages of openness, convenience, and well expansibility.Figure 5
Calculation strategy.


## 3.3. Parallel Efficiency

Example 1.
The monostatic RCS of a plasma sphere is simulated to show the parallel efficiency. The plasma sphere has a radius of 1 m (ωp=6.75×108rad/s, νc=7.5×107Hz), and the mesh size is 0.05 m. Three sets of platform with different numbers of processes are used to test the program (Figure 6). The computing platform is the Xidian high performance computing center.
Good agreements are observed between the numerical results and the Mie analytical solution, as shown in Figure6. Compared with the Mie solution, the root-mean square error of our algorithm is 1.5370 dB.
In Figure7, as a benchmark we chose 24 processes, and the parallel efficiency of 72 processes is 90.1%. The proposed algorithm has good performance in this example.Figure 6
Monostatic RCS of plasma sphere.Figure 7
Parallel efficiency.Example 2.
The previous example shows that the proposed algorithm has good practicability and high efficiency. Next, we calculate the monostatic RCS of a plasma-coated PEC sphere. The sphere of radius 0.5 meters is coated with plasma, as shown in Figure8. The plasma layer has a thickness of 0.2 meters. The plasma frequency and the collision frequency of the plasma are ωp=2.0×109rad/s and νc=0.1×109Hz, respectively.
The model is discretized into 3276319 tetrahedrons. Three sets of platform with different numbers of processes are used, and the numbers are 24, 48, and 72, respectively. The computing platform is the Xidian high performance computing center.
The Mie theory is an analytical method for the sphere target, and the Mie solution is used as a reference value in this example, as shown in Figure9 (cyan line). The monostatic RCS results of different sets are also shown in Figure 9, where the black line, red dash, and blue dot represent the results of 24 processes, 48 processes, and 72 processes, respectively. The Mie solution is consistent with the numerical results.
As seen in Figure10, the parallel efficiency of 72 processes is 91.50%. The proposed algorithm still has good performance in this example.
The calculation results illustrate the proposed parallel scheme is effective in parallel processing, and the results suggest the parallel SO-DGTD method can provide reliable results. The parallel SO-DGTD method can be used to analyze the plasma sheath target, as discussed in the next section.Figure 8
Plasma-coated PEC sphere.Figure 9
Monostatic RCS of the plasma-coated PEC sphere.Figure 10
Parallel efficiency.


## 4. Simulation Results of Hypersonic Target

The target is simplified as a PEC blunt cone. The cone is a rotationally symmetric target with a bottom radius of 0.05 meters and height of 0.3 meters. Let the height of the cone be 50 kilometers and the velocity be Mach 20. The parameters of the plasma sheath around the target at the certain altitude and velocity were calculated by COMSOL.The section view is shown in Figure12. Because of the very large span data, the plasma concentration is shown in terms of log10 in the plot.ne represents the plasma concentration (m-3), and the unit of plasma temperature is K. The Lagrange interpolation technique was used to obtain the parameters at the sampling point for the DGTD calculation.The meshes are assigned to a cluster containing 240 processes. The number of tetrahedrons per process is given in Figure13, showing that the processes are well-balanced. Three incident angles are used in the simulation (Figure 11). VV represent the polarized direction and receiving polarization vertical to the X axis, and HH represent the polarized direction and receiving polarization horizontal to the X axis. Moreover, we consider both the cone with a sheath and the cone without a sheath.Figure 11
Cone model.Figure 12
Plasma parameters.
(a)
Plasma concentration (log10(ne)) (b)
Plasma temperature (K)Figure 13
Number of tetrahedrons per process.Table2 provides the resource consumption of this example. There were 12 sets of parameters in this example, and the average CPU time is mean of 12 tests.Table 2
Calculation statistics (the computing platform is the Xidian high performance computing center).
Tetrahedron Processes Average CPU time Overall memory requirement  12805124  (Mesh size: 0.003 meters) 240 210 minutes 108.61 GBIn Figure14, the black line represents the result of the cone with a sheath; on the contrary, the red dash represents the result of the cone without a sheath. For Figure 14(a), because the plasma sheath is mainly distributed in the head part of the cone, the wave is strongly attenuated in this case. Thus, the RCS of the cone with a sheath is less than that without a sheath. It can be noticed that VV and HH with the incident direction k^:(θ=90∘,φ=90∘) returned the same results, because the cone is axial symmetry.Figure 14
Monostatic RCS.
(a)
k ^ : ( θ = 9 0 ∘ , φ = 9 0 ∘ ) VV (b)
k ^ : ( θ = 9 0 ∘ , φ = 9 0 ∘ ) HH (c)
k ^ : ( θ = 4 5 ∘ , φ = 9 0 ∘ ) VV (d)
k ^ : ( θ = 4 5 ∘ , φ = 9 0 ∘ ) HH (e)
k ^ : ( θ = 0 ∘ , φ = 9 0 ∘ ) VV (f)
k ^ : ( θ = 0 ∘ , φ = 9 0 ∘ ) HHFor Figure14(c), the results of the cone with and without a plasma sheath are almost the same. For k^:(θ=45∘,φ=90∘), in the HH case (Figure 14(d)), the RCS of the cone with a sheath is less than that without sheath below 8 GHz, above which they are similar. For Figures 14(e) and 14(f), among the two types of results (with and without a sheath) of VV or HH case, there are few differences. Considering that the thin plasma sheath lies at the tail of the cone, for the latter two cases, the sheath is not effective.


## 5. Conclusions

In this paper, our study provided a parallel framework for the SO-DGTD method by using MPI as a communication pattern, including mesh partition, communication between processes, and time-domain iteration. The algorithm was validated by the examples of a plasma sphere and a plasma-coated sphere. For a hypersonic target, the plasma is mainly distributed in the head part of the target; the results showed that when the electromagnetic wave vector is along the head of the cone, namely, the main distribution area of the plasma sheath, the reflection decline is influenced by the sheath. When the change in incident angle has distanced wave vector from the main distribution area of the plasma sheath, then the difference in the results of the cone with and without the sheath is declined.



---

*Source: 7160913-2019-04-21.xml*