Abstract

The plasma sheath during reentry of hypersonic vehicle is an unmagnetized and weakly ionized nonuniform plasma flow, which causes radio frequency blackout and strong plasma attenuation of electromagnetic wave. The physical properties of the nonuniform plasma flow were obtained using computational fluid dynamics software with unstructured grids. In this study, a detailed computational model was reconstructed with the high-order Lagrange grids for the nonuniform plasma flow region and the high-order Serendipity grids for the homogeneous medium region. In order to calculate the numerical flux between the two types of grids in the discontinuous Galerkin time domain (DGTD) algorithm, an asymmetric high-order element is constructed as a transition unit. Finally, the simulation results in the plasma sphere show that the above method improves the computational accuracy and decreases calculation. The amplitude and scattering about electromagnetic wave in nonuniform plasma flow are clarified in detail. It is suggested that the presented method could be an effective tool for investigating interaction between electromagnetic waves and plasma flow.

1. Introduction

In recent years, hypersonic weapons near space have played an increasingly important strategic role in warfare. When a hypersonic vehicle reenters the atmosphere, ablative hypersonic viscous flows produces a plasma sheath [1, 2], which is a nonuniform plasma flow composed of free electrons, ions, and neutral particles around hypersonic vehicles [3, 4].

In early research, an approximate method was used to predict the plasma effect, assuming that the plasma flow had a certain electron density distribution. In [5], Rusch simulated the actual plasma sheath profile of the reentry vehicle using a parabolic electron density profile. In [6], a uniform plasma layer was set on the surface of the spherical model to calculate the interaction between the electromagnetic wave and the target. Later, with the development of computer technology, numerical calculation methods such as finite-difference time-domain (FDTD), integral equation method (IEM), and physical optics (PO) were widely used in the study of electromagnetic wave in plasma flow. Basically, the simulations of interactions between the plasma and electromagnetic wave are performed by using the particle-in-cell method, which is based on the finite difference time domain (FDTD) method [7, 8] or the FDTD method combined with the effective permittivity model.

In [9], the physical optics method was used to calculate the scattering characteristics of the electromagnetic wave in inhomogeneous plasma sheath from the S-band to the Ku-band frequencies. In [10], a well-conditioned internally combined volume surface integral equation (ICVSIE) for analyzing electromagnetic scattering from perfect electrically conducting (PEC) surfaces coated with negative permittivity plasmas is presented to be applied in various canonical scatterers and a plasma-engulfed reentry vehicle. Scarabosio et al. employ the Eikonal approximation in the large inhomogeneous plasma region [11] and compute radiation and scattering via the equivalence theorem, and the results show that significant radio link path losses can be associated with plasma spatial variations (gradients) and collisional losses. In [12], the Drude mode is proven to be adopted to characterize collisional plasma. In the calculation with FDTD, radar cross section (RCS) of nonuniform plasma spheres and conductor spheres coated by the plasma layer are calculated. A piecewise linear JE recursive convolution finite-difference time-domain (PLJERC-FDTD) was used to study the scattering characteristics of electromagnetic waves in time-varying and inhomogeneous plasma sheath [13]. In [14], the high-order FDTD method was used to simulate the RCS (radar cross section) from the plasma flow in a shock tube. To increase the computational efficiency, Chen and Wang et al. used the hybrid-implicit-explicit (HIE)-FDTD method [15, 16] to simulate the gyrotropic plasma in open regions [17]. Chen and Wang proposed to apply the weakly conditionally stable (WCS)-FDTD method [18] to simulate the plasma frequency selective surface [19]. Numerical experiments showed that both the HIE-FDTD and WCS-FDTD methods can enhance significantly the computational efficiency of plasma simulations. However, because of the used Yee cell, the FDTD methods suffer from serious accuracy deterioration when they come to irregular objects.

In [20], electromagnetic wave behavior in plasma with the combined CFD and FD2TD method was investigated. The distribution of charged particles around the ESA ARD and the complicated behavior of electromagnetic waves, with attenuation and reflection, are clarified in detail.

For DGTD analysis, elements in tetrahedral form can be better fit to irregular objects [21, 22]. In [23], DGTD is used to solve electromagnetic scattering from hypersonic aircraft with plasma sheath; however, the flight status and the specific parameters of the plasma sheath were not given in the paper. In [24], an exponential-based DGTD method for the 3-D modeling of hypersonic vehicles with surrounded by a uniform plasma layer is presented, and this method enables the usage of larger time steps compared with the explicit time marching scheme. In the abovementioned DGTD methods, tetrahedral elements are used in the computational domain, which is inefficient when calculating complex geometric objects. More recently, the DGTD with hexahedron have attracted extensive attention for modeling with complex geometrical features [25]. In [26], DGTD with the high-order basis function is proposed to solve electromagnetic scattering from hypersonic aircraft with plasma sheath in which dielectric parameter changes with frequency, respectively. A hybrid mesh DGTD algorithm based on virtual elements for tetrahedra and hexahedra is proposed to improve the computational efficiency [27]. But tetrahedra and hexahedra are applied in different subdomains, respectively.

In this study, the Lagrange high-order element is used to read more plasma parameters in the plasma flow region, the efficient Serendicity high-order element is used in the uniform medium area, and the asymmetric high-order elements are constructed as the transition units between two types of elements to calculate the numerical flux. The high-order points of the Serendicity high-order element close to the shape of aircraft can be mapped to aircraft surface by isoparametric element. This article is organized as follows. In Section 2, time-domain Maxwell’s equations in the Drude dispersive medium are presented. Section 3 describes the physical and numerical models of asymmetric mixed higher-order mesh. As far as the authors’ knowledge, this is the first time that an asymmetric mixed higher-order mesh DGTD method is developed for the numerical simulation of nonuniform plasma flow target. Finally, numerical experiments of scattering problems are presented in Section 4, and the accuracy of the proposed method is verified through comparison with the results of Mie theory. It was confirmed that plasma flow under different flight altitudes has different effects on electromagnetic scattering.

2. Iterative Equation of DGTD

We recall the propagation of an electromagnetic wave in the dispersive medium. We can use the dielectric coefficient to express the characteristics of the nonmagnetized medium. The dielectric coefficient related to frequency can be given bywhere is the relative dielectric constant at infinity frequency, is the dielectric coefficient in vacuum, and is the polarizability function.

2.1. Drude Dispersive Medium

The nonuniform plasma flow of hypersonic vehicles is usually treated as cold plasma, which is expressed by the Drude dispersive medium. The thermochemical nonequilibrium model is used in computational fluid dynamics (CFD) to simulate the nonuniform plasma flow. The concentration of electrons, density of neutral particles, and temperature in the nonuniform plasma flow are important parameters for calculating the dielectric coefficient. In the Drude dispersive medium [28], is given bywhere is the plasma collision frequency and is the plasma frequency which are given bywhere is the electron mass, is the number density of electrons, is the number density of species s, e is the magnitude of the electronic charge equal to , and k is Boltzmann’s constant. The effective electron-neutral energy exchange cross section is defined by a curve fit of the form.

The constants for equation (4) are presented in Table 1.

2.2. Formulation in Computational Domain

The Maxwell equations are

Substituting equation (1) into equation (5), we obtain

Space is divided into nonoverlapping and continuous subregions. Each subregion is a separate hexahedral element, and its weight function (also called basis function) is (i is the serial number of the unit). The weighted integral of equation (6) in each element is taken and set equal to zero, such thatwhere is the weight function. The numerical fluxes are defined at the boundary between two elements aswhere is the unit outward normal vector in each surface of the element, and are the numerical fluxes, and are in the fields of the adjacent element, and and are in the fields of the element. , , , and are presented in Table 2.

Substituting (9) into (7), we obtain

3. Asymmetrically Mixed Higher-Order Mesh

The aircraft is a RAM-CII aircraft scaling model in this paper. The nonuniform flow is simulated by FASTRAN whose numerical algorithm is the finite volume method based on density [29].

In Figure 1, the plasma flow forms a viscous boundary layer on the aircraft surface when flight altitude is 55 km and velocity altitude is 20 Mach, and the plasma parameters have a large gradient. The maximum value of plasma frequency appears in the area of the blunt cone head (with an order of magnitude of ), and the closer to the blunt cone head, the greater the value of plasma frequency (as shown in Figure 1(a). The plasma frequency in the other regions decreased approximately one order of magnitude, about . The maximum value of the collision frequency of appears in the area of the blunt cone head, with an order of magnitude of (as shown in Figure 1(b)). The collision frequency in the other regions decreased approximately one order of magnitude around .

A detailed computational model was reconstructed with mixed high-order hexahedral grids for the DGTD code. In order to read more parameters of the nonuniform plasma flow, the electromagnetic model is divided by mixed high-order hexahedral grids.

A hexahedral element constructs the basis function through isoparametric elements [30, 31]. The nodes of the actual mesh and isoparametric cell are in one-to-one correspondence, as shown in Figure 2. Construct three isoparametric unctions (as shown in Figures 2(a)2(c)) for three actual unctions (as shown in Figures 2(d)2(f)) defined over the square .

The isoparametric element is the standard element in coordinates. If there are n basis functions, they can be given by

The basis function satisfies the abovementioned conditions, and the mapping formula from local coordinates to actual coordinates is as follows:

The isoparametric element can be mapped to element with arbitrary shape using the abovementioned formula. In the higher-order element, the right-angle edge of the standard element can be mapped to the parabolic boundary of the actual element.

In Figure 3, the 27 nodes high-order Lagrange element is used in the nonuniform plasma flow, which can read more plasma parameters by using more internal nodes. The efficient 20 nodes Serendipity element is used in the homogeneous medium region, which has no internal nodes, thus ensuring the accuracy of calculation and reducing the amount of computation. In order to calculate the numerical flux between the two types of elements in the DGTD algorithm, an asymmetric high-order element is constructed as a transition unit.

3.1. NonUniform Plasma Flow Region Mesh

Lagrange high-order element is used in the nonuniform plasma flow region. It is assumed that the isoparametric element is divided into groups, groups, and groups of nodes in the x, y, and z directions, respectively, and the nodes are located at the intersection lines of the element. Node i is located in groups m, n, and k in the x, y, and z directions, respectively. Its interpolation functions are

The basis function of node i can be obtained by multiplying the Lagrange polynomials in the above:

Grid-generation software usually provides only the first order nodes of an element. In order to effectively fit the higher order element nodes near the boundary of the plasma flow to the contour of the plasma flow, we set more surface points of the plasma flow region as mapping points of the higher-order element by constructing a cylindrical coordinate system.

Firstly, the surface mesh of the flow region is further subdivided separately. The further subdivided mesh is denser than the working mesh as shown in Figure 4(b).

Then, the cylindrical coordinate system rθZ in Figure 4 is established in the further mesh. The actual coordinate points and the higher-order grid nodes close to the flow field have the same values on the Z and θ axes (as shown in Figure 5). The transformation relationship between the cylindrical coordinate system and the rectangular coordinate system is as follows:

Finally, a suitable cylindrical surface point is found according to the value of Z axis and θ axis corresponding to the higher-order node.

3.2. Uniform Medium Region Mesh

The uniform medium region is divided by the Serendipity high-order element. The Serendipity higher-order element is shown in Figure 3(c), whose basis functions are

3.3. Transition Unit (Asymmetric High Order Element)

The Lagrange high-order element with 27 nodes and the Serendipity high-order element with 20 nodes are used in this paper. There are 9 nodes at each face of the Lagrange high-order element and 8 nodes at each face of the Serendipity high-order element. An asymmetric high-order element with 21 nodes is constructed as a transition element (as shown in Figure 6).

To determine the arbitrary basis function as follows, we divide n planes that do not pass through node i but pass through all other nodes of the element.

Next, the basis function of the asymmetric higher-order element is constructed. First, high-order node 21 is constructed on the basis of 20 nodes Serendipity high-order element. is given by

Firstly, we assume that the basic functions of a 20-nodes Serendipity higher-order grid are . Substitute and into the displacement function [32]:

Substitute into equation (19)

Then, we can get as

Finally, substitute instead of into the displacement function:

All basis functions can be deduced in terms of the abovementioned way.

4. Examples

In the following, two numerical experiments are used to verify the effectiveness of the asymmetrically mixed high-order mesh DGTD method, and the influence of nonuniform plasma flow on electromagnetic scattering is discussed.

4.1. RCS of a Plasma Sphere

We calculate the monostatic radar cross section (RCS) of a plasma ball . In the case of mixed higher order element, the plasma sphere mesh is mapped to 27 nodes second-order Lagrange element, and the other mesh is mapped to 20 nodes second-order Serendipity element, an asymmetric high order element with 21 nodes is constructed as a transition unit between two types of elements.

As shown in Figure 7, a very good agreement is achieved among the results obtained using the mixed higher-order element method, second-order Lagrange element method, and the Mie theory. Compared with DGTD with the first order element, this method has less RCS error, especially when the frequency is 76.2 GHz, 93.5 GHz (as shown in Table 3).

4.2. RCS of Inhomogeneous Plasma Flow around Blunt Cone

In this example, monostatic and bistatic radar cross section (RCS) of the blunted cone surrounded by nonuniform plasma flow is illuminated. We consider a Gaussian pulse incidents from the front of the blunted cone . The blunted cone has a blunted nose of radius 0.01 m, a conical cross section with half-angle of , and a whole length 0.1 m.

As shown in Figure 8, the 3D geometric model is divided by the hexahedral grids consists of 1067631 hexahedrons. Table 4 shows the memory required and calculation times of 27 nodes Lagrange mesh and mix-order mesh with 8 parallel threads algorithm. The results demonstrate that mix-order mesh is more efficient, with the required memory being reduced by . To truncate the open boundary of the simulation, the uniaxial perfectly matched layer (UPML) with a thickness of 0.0032 m is placed at least 0.006 m away from the plasma sheath.

In Figure 9, the influence of flight altitude and Mach number on the backward RCS of blunt cones at altitudes of 50 and 55 km and speeds of 20 Mach are analyzed. When the frequency is below 1.8 GHz, the plasma flow enhances the backward RCS of the blunt cone greatly. When the collision frequency is greater than the incident wave frequency, reducing the incident wave frequency will increase reflection [33]. When the frequency is 1.8 GHz13 GHz, the plasma flow will reduce the backward RCS of blunt cone. When the frequency is above 13 GHz, the backward RCS of blunt cone with plasma is consistent with that without plasma, indicating that the plasma flow has little influence on the electromagnetic scattering of the reentry target in the high frequency band.

Aiming at abovementioned characteristics, we discuss the bistatic RCS and the electric field amplitude of the nonuniform flow field target, which is at the altitude of 55 km and the velocity of 20 Mach with three frequency bands of 1.0 GHz, 4.0 GHz, and 15.0 GHz.

The amplitude at 1.0 GHz without and with the plasma flow of XOZ section is given in Figures 10(a) and 10(b). When the frequency is 1.0 GHz, the incident wave cannot penetrate the plasma flow near the nose as show in the partial enlarged view. In the plasma layer near the nose, the blue shadow in Figure 10(b) increases compared to Figure 10(a). Therefore, the backward and lateral RCS of the blunt cone with plasma flow, particularly in 50 km 20 Mach condition, are greater than that of the blunt cone without plasma flow in Figure 10(c). The interaction between the plasma flow and the incident wave is more obvious, as the altitude decreases.

When the frequency of the incident wave increases to 4.0 GHz, the lateral RCS of the blunt cone with plasma flow are smaller than that of the blunt cone without plasma flow, particularly in 50 km 20 Mach condition, the backward RCS is the reverse of the lateral RCS (as shown Figure 11(c)). Compared with Figure 11(a), the effect of plasma flow on incident waves is relatively strong in Figure 11(b).

When the frequency of the incident wave is 15.0 GHz, the incident wave can penetrate the plasma flow. The influence of plasma flow on amplitude decreases in Figure 12(b). The RCS of the blunt cone with plasma flow is close to that without plasma flow, particularly in 55 km 20 Mach condition (as shown in Figure 12(c)).

5. Conclusion

Most current calculation schemes for DGTD analysis are based on a single type of mesh scheme, which is inefficient when calculating no-uniform flow field. Here, we propose an asymmetrically mixed higher-order mesh DGTD method to address this problem. This method took into account the higher-order computational accuracy while reducing computation. The reliability of the asymmetrically mixed high-order mesh DGTD method is verified by plasma ball experiments. Compared with the single-type high-order mesh, the proposed method reduces the computational effort and has the advantage over the first-order mesh to further profile the target. We study the influence between plasma flow of hypersonic vehicle and electromagnetic wave from the perspective of amplitude and scattering for the first time. More research ideas are proposed for studying the relationship between flight height and electromagnetic wave frequency.

This article has considered an asymmetrically mixed higher-order mesh DGTD method for only two types of mesh, the Lagrange high-order mesh with 27 nodes and the Serendipity high-order mesh with 20 nodes, for a type of transition unit. This method currently applies only to these two kinds of mesh. Nevertheless, it can also be adapted for other types of Lagrange grids and Serendipity grids and constructs corresponding transition unit. Moreover, when the order of the magnitude gradient of nonuniform medium parameters is relatively large, different types of mesh need to be constructed according to the gradient values. This will be investigated in future work.

Data Availability

The processed data are available from the corresponding author upon request. The data for other figures in this paper used to support the findings of this study are also available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was supported by the National Key Research and Development Program of China under grant no. 6142411203109.