Research Article  Open Access
Analysis of the Calculation of a Plasma Sheath Using the Parallel SODGTD Method
Abstract
The plasma sheath is known as a popular topic of computational electromagnetics, and the plasma case is more resourceintensive than the nonplasma case. In this paper, a parallel shiftoperator discontinuous Galerkin timedomain method using the MPI (Message Passing Interface) library is proposed to solve the largescale plasma problems. To demonstrate our algorithm, a plasma sheath model of the highspeed blunt cone was established based on the results of the multiphysics software, and our algorithm was used to extract the radar crosssection (RCS) versus different incident angles of the model.
1. Introduction
A highspeed 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 timedomain (DGTD) method offers superior versatility in geometry discrimination by using unstructured meshes; unlike the finite difference timedomain method (FDTD) [8, 9], the DGTD method has no numerical dispersion. DGTD is appropriate for solving complex electromagnetic problems [10–15]. The shiftoperator DGTD (SODGTD) method was proposed for dispersive media [16–18]; to address the frequency domain constitutive relation and obtain the timedomain iteration formula, the shiftoperator is used. Plasma sheath is a problem that requires a huge consumption of computer resources to solve, resulting in a huge resource consumption in timedomain simulation. In response, the parallel technology can improve the efficiency while allowing simulation of largescale 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 MPIDGTD 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 MPIDGTD based on nonconformal meshes and the local timestepping 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 SODGTD method; next, a specific parallel scheme for SODGTD 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 ultrahighspeed aircraft can significantly affect its RCS over a wide frequency band.
2. SODGTD Formulation For Plasma Medium
The semidiscrete equations of DGTD arewhere is the permeability and and are the electric and magnetic field, respectively. is the electric flux density. is the index of element located in the computational region, is the mass matrix, is the stiffness matrix, and and are the flux matrices.
The constitutive relation of plasma is defined by the Drude modelwhere is the permittivity of vacuum, is the relative permittivity (which is dependent on frequency), is the plasma frequency, is the collision frequency, and is the angular frequency of electromagnetic waves. Equations (1) allow us to update by using and update by using . Obviously, we need a formula from to complete the cycle of . Bing Wei proved that the complex permittivity of the general dispersive media model can be described by a rational polynomial fraction in [24]:where and are the maximum powers of the fractional polynomial and and are the expansion coefficients of polynomials. According to (3), by using the relationship of frequencytimedomain (), we can obtainwhere denotes the time step. The relation of the timedomain differential operator and the shiftoperator is as follows.
The final timedomain form corresponding to (3) is [16]where represents the shiftoperator, defined by . The role of is to shift the time sampling from the nth step to the (n+1)th step. , where is the time increment in the calculation.
By using the characteristic of , we obtain the iteration formula from (6)
For the Drude model (plasma), the expansion coefficients are described as follows.
To match (7), the leapfrog scheme is used to solve (1).
3. MPI Parallelization of SODGTD
The proposed algorithm contains three types of variables, namely, electric field E, magnetic field , and electric displacement field . In our algorithm, E and are the communications data between neighboring elements. and in (1) can be written as below.
For element e, we denote the coefficients of flux by , , , and . and denote the flux matrices. and denote the fields of element ; and denote the fields of adjacent element e+. In parallel processing, each partition is mapped to an MPI process, and the main task is the communication of data. If and e+ 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.
3.1. Mesh Partition
As shown in Figure 2, 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].
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 Figure 3 and discussed in the next section.
3.2. Communication and Calculation Setup
The key problem of iterative calculation is addressing the issue of communications between neighboring processes. As shown in Figure 4, tetrahedron belongs to process M, and the neighboring tetrahedron also belongs to process M, but another neighboring tetrahedron belongs to process . Thus, the data among and can be read directly, whereas the data among and must be communicated between process and . Communication of parallel computing is implemented via the mapping of meshes.
Figure 4 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.

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_out n=N_out_NODE(i) CALL MPI_SENDRECV (E_cache_out(1:8,1:Num_NODE(n),n), 8Num_NODE(n), MPI_REAL, m1, 99+n, E_cache_in(1:8, 1:Num_NODE(n), n), 8Num_NODE(n), MPI_REAL, m1, 99+MYTIDS, MPI_COMM_WORLD, STATUS, IERR) Enddo
In 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 Figure 5. It should be noted that the parallel scheme is directed to dispersive media, and the scheme calculates the fields first, followed by the fields and fields lastly. This process repeats many times until the program ends. This scheme has the advantages of openness, convenience, and well expansibility.
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 (, ), 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 Figure 6. Compared with the Mie solution, the rootmean square error of our algorithm is 1.5370 dB.
In Figure 7, 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.
Example 2. The previous example shows that the proposed algorithm has good practicability and high efficiency. Next, we calculate the monostatic RCS of a plasmacoated PEC sphere. The sphere of radius 0.5 meters is coated with plasma, as shown in Figure 8. The plasma layer has a thickness of 0.2 meters. The plasma frequency and the collision frequency of the plasma are and , 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 Figure 9 (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 Figure 10, 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 SODGTD method can provide reliable results. The parallel SODGTD method can be used to analyze the plasma sheath target, as discussed in the next section.
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 Figure 12. Because of the very large span data, the plasma concentration is shown in terms of log_{10} in the plot. ne represents the plasma concentration (), 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 Figure 13, showing that the processes are wellbalanced. Three incident angles are used in the simulation (Figure 11). VV represent the polarized direction and receiving polarization vertical to the axis, and HH represent the polarized direction and receiving polarization horizontal to the axis. Moreover, we consider both the cone with a sheath and the cone without a sheath.
(a) Plasma concentration (log10(ne))
(b) Plasma temperature (K)
Table 2 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.

In Figure 14, 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 returned the same results, because the cone is axial symmetry.
(a) VV
(b) HH
(c) VV
(d) HH
(e) VV
(f) HH
For Figure 14(c), the results of the cone with and without a plasma sheath are almost the same. For , 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 SODGTD method by using MPI as a communication pattern, including mesh partition, communication between processes, and timedomain iteration. The algorithm was validated by the examples of a plasma sphere and a plasmacoated 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.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the National Natural Scientific Foundation of China (61571348, 61231003, 61401344).
References
 AIAA, “Radio frequency (RF) blackout during hypersonic reentry,” Annales Dhistoire Économique Et Sociale, vol. 7, no. 35, p. 514, 1935. View at: Publisher Site  Google Scholar
 M. Kim, M. Keidar, and I. D. Boyd, “Electrostatic manipulation of a hypersonic plasma layer: Images of the twodimensional sheath,” IEEE Transactions on Plasma Sciences, vol. 36, no. 4, pp. 11981199, 2008. View at: Publisher Site  Google Scholar
 O. I. Sukharevsky, G. S. Zalevsky, S. V. Nechitaylo, and I. O. Sukharevsky, “Simulation of scattering characteristics of aerial resonantsize objects in the VHF band,” Radioelectronics and Communications Systems, vol. 53, no. 4, pp. 213–218, 2010. View at: Publisher Site  Google Scholar
 J. P. Sheehan, I. D. Kaganovich, H. Wang, D. Sydorenko, Y. Raitses, and N. Hershkowitz, “Effects of emitted electron temperature on the plasma sheath,” Physics of Plasmas, vol. 21, no. 6, Article ID 063502, 2014. View at: Google Scholar
 G. He, Y. Zhan, N. Ge, Y. Pei, B. Wu, and Y. Zhao, “Channel characterization and finitestate Markov channel modeling for timevarying plasma sheath surrounding hypersonic vehicles,” Progress in Electromagnetics Research, vol. 145, no. 1933, pp. 299–308, 2014. View at: Publisher Site  Google Scholar
 B. Bai, X. Li, Y. Liu, J. Xu, L. Shi, and K. Xie, “Effects of reentry plasma sheath on the polarization properties of obliquely incident EM waves,” IEEE Transactions on Plasma Sciences, vol. 42, no. 10, pp. 3365–3372, 2014. View at: Publisher Site  Google Scholar
 M. Puttscher, A. Melzer, U. Konopka, S. LeBlanc, B. Lynch, and E. Thomas, “Vertical oscillations of dust particles in a strongly magnetized plasma sheath induced by horizontal laser manipulation,” Physics of Plasmas, vol. 24, no. 1, Article ID 013701, 2017. View at: Google Scholar
 K. Niu, Z. Huang, M. Li et al., “Optimization of the artificially anisotropic parameters in WCSFDTD method for reducing numerical dispersion,” IEEE Transactions on Antennas and Propagation, vol. 65, no. 12, pp. 7389–7394, 2017. View at: Publisher Site  Google Scholar
 K. Niu, Z. Huang, X. Ren, M. Li, B. Wu, and X. Wu, “An optimized 3D HIEFDTD method with reduced numerical dispersion,” IEEE Transactions on Antennas and Propagation, vol. 66, no. 11, pp. 6435–6440, 2018. View at: Publisher Site  Google Scholar
 I. Perugia and D. Schötzau, “The hplocal discontinuous Galerkin method for lowfrequency timeharmonic Maxwell equations,” Mathematics of Computation, vol. 72, no. 243, pp. 1179–1214, 2003. View at: Publisher Site  Google Scholar  MathSciNet
 P. Houston, I. Perugia, and D. Schötzau, “hpDGFEM for Maxwell's equations,” Numerical Mathematics and Advanced Applications, pp. 785–794, 2003. View at: Google Scholar
 B. Cockburn, F. Li, and C.W. Shu, “Locally divergencefree discontinuous Galerkin methods for the Maxwell equations,” Journal of Computational Physics, vol. 194, no. 2, pp. 588–610, 2004. View at: Publisher Site  Google Scholar  MathSciNet
 S. G. Garcia, M. F. Pantoja, C. M. D. J. V. Coevorden, A. R. Bretones, and R. G. Martín, “A new hybrid DGTD/FDTD method in 2D,” IEEE Microwave and Wireless Components Letters, vol. 18, no. 12, pp. 764–766, 2008. View at: Publisher Site  Google Scholar
 S. Lanteri, V. Dolean, and A. Catella, “An implicit discontinuous Galerkin timedomain method for twodimensional electromagnetic wave propagation,” COMPEL. The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, vol. 29, no. 3, pp. 602–625, 2010. View at: Publisher Site  Google Scholar  MathSciNet
 L. E. Tobón, Q. Ren, and Q. H. Liu, “A new efficient 3D discontinuous Galerkin time domain (DGTD) method for large and multiscale electromagnetic simulations,” Journal of Computational Physics, vol. 283, pp. 374–387, 2015. View at: Publisher Site  Google Scholar  MathSciNet
 Q. Yang, B. Wei, L. Li, R. Li, and D. Ge, “SODGTD method for dispersive media,” in Proceedings of the 2017 IEEE International Symposium on Antennas and Propagation and USNCURSI Radio Science Meeting, APSURSI 2017, pp. 24232424, USA, July 2017. View at: Google Scholar
 L. Li, B. Wei, Q. Yang, X. Yang, and D. Ge, “2D highorder DGTD method and its application in analysis of sheath propagation characteristics,” IEEE Transactions on Plasma Sciences, vol. 45, no. 9, pp. 2422–2430, 2017. View at: Publisher Site  Google Scholar
 B. Wei, L. Li, Q. Yang, and D. Ge, “Analysis of the transmission characteristics of radio waves in inhomogeneous weakly ionized dusty plasma sheath based on high order SODGTD,” Results in Physics, vol. 7, pp. 2582–2587, 2017. View at: Publisher Site  Google Scholar
 M. Bernacki, L. Fezoui, S. Lanteri, and S. Piperno, “Parallel discontinuous Galerkin unstructured mesh solvers for the calculation of threedimensional wave propagation problems,” Applied Mathematical Modelling, vol. 30, no. 8, pp. 744–763, 2006. View at: Publisher Site  Google Scholar
 T. Cabel, J. Charles, and S. Lanteri, MultiGPU Acceleration of a DGTD Method for Modeling Human Exposure to Electromagnetic Waves, INRIA, 2011.
 S. Dosopoulos, J. D. Gardiner, and J. F. Lee, “An MPI/GPU parallelization of an interior penalty discontinuous Galerkin time domain method for Maxwell's equations,” Radio Science, vol. 46, 2011. View at: Google Scholar
 S. Dosopoulos, B. Zhao, and J.F. Lee, “Nonconformal and parallel discontinuous Galerkin time domain method for Maxwell's equations: EM analysis of IC packages,” Journal of Computational Physics, vol. 238, pp. 48–70, 2013. View at: Publisher Site  Google Scholar  MathSciNet
 R. Léger, J. Viquerat, C. Durochat, C. Scheid, and S. Lanteri, “A parallel nonconforming multielement DGTD method for the simulation of electromagnetic wave interaction with metallic nanoparticles,” Journal of Computational and Applied Mathematics, vol. 270, pp. 330–342, 2014. View at: Publisher Site  Google Scholar  MathSciNet
 B. Wei, D.B. Ge, and F. Wang, “A general method for finite difference time domain modeling of wave propagation in frequencydispersive media,” Acta Physica Sinica, vol. 57, no. 10, pp. 6290–6297, 2008. View at: Google Scholar
 C. Geuzaine and J. F. Remacle, “Gmsh: A 3D finite element mesh generator with builtin pre and postprocessing facilities,” International Journal for Numerical Methods in Engineering, vol. 79, no. 11, pp. 1309–1331, 2009. View at: Publisher Site  Google Scholar  MathSciNet
Copyright
Copyright © 2019 Qian Yang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.