International Journal of Antennas and Propagation
Volume 2019 (2019), Article ID 7160913, 9 pages
https://doi.org/10.1155/2019/7160913

Research Article

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

Qian Yang,1,2 Bing Wei,1,2 Linqian Li,1,2 and Debiao Ge1,2

¹School of Physics and Optoelectronic Engineering, Xidian University, Xian 710071, China
²Collaborative Innovation Center of Information Sensing and Understanding, Xidian University, Xian 710071, China

Correspondence should be addressed to Qian Yang; qyang@xidian.edu.cn

Received 16 January 2019; Accepted 18 March 2019; Published 21 April 2019

Academic Editor: Ping Li

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.

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.

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

where 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 model

where 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 frequency-time-domain (), we can obtain

where denotes the time step. The relation of the time-domain differential operator and the shift-operator is as follows.

The final time-domain form corresponding to (3) is [16]

where represents the shift-operator, 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 leap-frog scheme is used to solve (1).

3. MPI Parallelization of SO-DGTD

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.

Figure 1: Flow chart.

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].

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 Figure 3 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 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: 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.

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_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.

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 (, ), 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 root-mean 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.

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 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 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 Figure 12. Because of the very large span data, the plasma concentration is shown in terms of log₁₀ 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 well-balanced. 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.

Figure 11: Cone model.

(a) Plasma concentration (log₁₀(ne))

(b) Plasma temperature (K)

Figure 12: Plasma parameters.

Figure 13: Number of tetrahedrons per process.

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.

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 GB

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

Figure 14: Monostatic RCS.

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 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.

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.

M. Kim, M. Keidar, and I. D. Boyd, “Electrostatic manipulation of a hypersonic plasma layer: Images of the two-dimensional sheath,” IEEE Transactions on Plasma Sciences, vol. 36, no. 4, pp. 1198-1199, 2008.

O. I. Sukharevsky, G. S. Zalevsky, S. V. Nechitaylo, and I. O. Sukharevsky, “Simulation of scattering characteristics of aerial resonant-size objects in the VHF band,” Radioelectronics and Communications Systems, vol. 53, no. 4, pp. 213–218, 2010.

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.

G. He, Y. Zhan, N. Ge, Y. Pei, B. Wu, and Y. Zhao, “Channel characterization and finite-state Markov channel modeling for time-varying plasma sheath surrounding hypersonic vehicles,” Progress in Electromagnetics Research, vol. 145, no. 1933, pp. 299–308, 2014.

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.

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.

K. Niu, Z. Huang, M. Li et al., “Optimization of the artificially anisotropic parameters in WCS-FDTD method for reducing numerical dispersion,” IEEE Transactions on Antennas and Propagation, vol. 65, no. 12, pp. 7389–7394, 2017.

K. Niu, Z. Huang, X. Ren, M. Li, B. Wu, and X. Wu, “An optimized 3-D HIE-FDTD method with reduced numerical dispersion,” IEEE Transactions on Antennas and Propagation, vol. 66, no. 11, pp. 6435–6440, 2018.

I. Perugia and D. Schötzau, “The hp-local discontinuous Galerkin method for low-frequency time-harmonic Maxwell equations,” Mathematics of Computation, vol. 72, no. 243, pp. 1179–1214, 2003.

P. Houston, I. Perugia, and D. Schötzau, “hp-DGFEM for Maxwell's equations,” Numerical Mathematics and Advanced Applications, pp. 785–794, 2003.

B. Cockburn, F. Li, and C.-W. Shu, “Locally divergence-free discontinuous Galerkin methods for the Maxwell equations,” Journal of Computational Physics, vol. 194, no. 2, pp. 588–610, 2004.

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 2-D,” IEEE Microwave and Wireless Components Letters, vol. 18, no. 12, pp. 764–766, 2008.

S. Lanteri, V. Dolean, and A. Catella, “An implicit discontinuous Galerkin time-domain method for two-dimensional electromagnetic wave propagation,” COMPEL. The International Journal for Computation and Mathematics in Electrical and Electronic Engineering, vol. 29, no. 3, pp. 602–625, 2010.

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.

Q. Yang, B. Wei, L. Li, R. Li, and D. Ge, “SO-DGTD method for dispersive media,” in Proceedings of the 2017 IEEE International Symposium on Antennas and Propagation and USNC-URSI Radio Science Meeting, APSURSI 2017, pp. 2423-2424, USA, July 2017.

L. Li, B. Wei, Q. Yang, X. Yang, and D. Ge, “2-D high-order DGTD method and its application in analysis of sheath propagation characteristics,” IEEE Transactions on Plasma Sciences, vol. 45, no. 9, pp. 2422–2430, 2017.

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 SO-DGTD,” Results in Physics, vol. 7, pp. 2582–2587, 2017.

M. Bernacki, L. Fezoui, S. Lanteri, and S. Piperno, “Parallel discontinuous Galerkin unstructured mesh solvers for the calculation of three-dimensional wave propagation problems,” Applied Mathematical Modelling, vol. 30, no. 8, pp. 744–763, 2006.

T. Cabel, J. Charles, and S. Lanteri, Multi-GPU 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.

S. Dosopoulos, B. Zhao, and J.-F. Lee, “Non-conformal 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.

R. Léger, J. Viquerat, C. Durochat, C. Scheid, and S. Lanteri, “A parallel non-conforming multi-element DGTD method for the simulation of electromagnetic wave interaction with metallic nanoparticles,” Journal of Computational and Applied Mathematics, vol. 270, pp. 330–342, 2014.

B. Wei, D.-B. Ge, and F. Wang, “A general method for finite difference time domain modeling of wave propagation in frequency-dispersive media,” Acta Physica Sinica, vol. 57, no. 10, pp. 6290–6297, 2008.

C. Geuzaine and J. F. Remacle, “Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities,” International Journal for Numerical Methods in Engineering, vol. 79, no. 11, pp. 1309–1331, 2009.