Analysis of Radio Frequency Blackout for a Blunt-Body Capsule in Atmospheric Reentry Missions

Abstract

A numerical analysis of electromagnetic waves around the atmospheric reentry demonstrator (ARD) of the European Space Agency (ESA) in an atmospheric reentry mission was conducted. During the ARD mission, which involves a 70% scaled-down configuration capsule of the Apollo command module, radio frequency blackout and strong plasma attenuation of radio waves in communications with data relay satellites and air planes were observed. The electromagnetic interference was caused by highly dense plasma derived from a strong shock wave generated in front of the capsule because of orbital speed during reentry. In this study, the physical properties of the plasma flow in the shock layer and wake region of the ESA ARD were obtained using a computational fluid dynamics technique. Then, electromagnetic waves were expressed using a frequency-dependent finite-difference time-domain method using the plasma properties. The analysis model was validated based on experimental flight data. A comparison of the measured and predicted results showed good agreement. The distribution of charged particles around the ESA ARD and the complicated behavior of electromagnetic waves, with attenuation and reflection, are clarified in detail. It is suggested that the analysis model could be an effective tool for investigating radio frequency blackout and plasma attenuation in radio wave communication.

1. Introduction

One of the most important issues in the entry, descent and landing approach of spacecraft is radio frequency (RF) blackout; that is, the cutoff of communication between the reentry vehicle and ground stations and/or data-relay satellites. Figure 1 shows a schematic view of RF blackout between a reentry vehicle and a data relay satellite during atmospheric reentry. During reentry from a low Earth orbit, a vehicle flies at an orbital speed of more than 7 km/s, and a strong shock wave is formed at its front. Freestream air is rapidly dissociated by aerodynamic heating in the shock wave layer. Charged species, i.e., electrons and ions, are created, and these species then cover the vehicle’s surface. RF blackout is caused by a highly dense plasma layer on the vehicle surface, attenuating or reflecting electromagnetic waves used for telecommunication. This RF blackout has been observed in many reentry flights, including those of the RAM C, Apollo command module (NASA) and OREX (NASDA/NAL, currently JAXA).

Figure 1. Schematics of radio frequency blackout during atmospheric reentry.

Figure 1. Schematics of radio frequency blackout during atmospheric reentry.

To solve this problem, so far, several techniques for mitigating RF blackout and plasma attenuation have been proposed, including using lasers, magnetic fields and electrophilic substances [1,2,3]. For example, Kim et al. [4] reported that applying an electromagnetic $\mathbf{E}$×$\mathbf{B}$ layer to an antenna can control the plasma density and reduce RF blackouts. Belov et al. [5] experimentally investigated an aerodynamic shaping technique to mitigate RF blackout by mounting slender nose-cone probes ahead of a blunt-nosed vehicle. The thick plasma in front of the vehicle generated by the shock wave strongly prevents electromagnetic waves from propagating. On the other hand, because the plasma flow is rapidly cooled and rarefied by expansion when compressed gas in the shock layer flows toward the rear of the vehicle, propagation of electromagnetic waves is possible at the rear of the vehicle where the electron density is lower. In fact, it was suggested by Takahashi et al. [6] that an inflatable vehicle could passively reduce RF blackouts with data relay satellites during atmospheric reentry because of its low ballistic coefficient flight and distinctive vehicle configuration. Thus, it may be possible for a vehicle with such features and devices to communicate without RF blackouts during a reentry mission.

Because the RF blackout phenomenon is affected by the flow properties around the vehicle, specifically the number density of electrons, clarifying the distribution of charged particles is a key issue. In particular, determination of the plasma distribution around a blunt body reentry capsule is expected to be important from the point of view of their frequent use in reentry missions. However, it is difficult to completely reproduce the environment in which RF blackout is caused using ground testing facilities, such as a hypersonic wind tunnels, arc-heated wind tunnels, ballistic ranges, etc. Recently, remarkable progress has been made in the development of high-performance computers and algorithms for numerical simulation. Numerical prediction using high-performance computers has become a powerful tool for investigating the details of phenomena, such as flow fields or electromagnetic waves, and for the design and development of reentry vehicles. Computational fluid dynamics (CFD) techniques have been used widely to predict flow fields around reentry vehicles. So far, several analysis models to describe electromagnetic waves in reentry plasma have been proposed. Numerical simulations of communications in weakly-ionized plasma by using a ray-tracing solver and CFD were conducted by Vecchi et al. [7,8] and Delfino et al. [9] Moreover, Yucel et al. [10] proposed a volume-surface integral equation solver accelerated by a fast multiple method and fast Fourier transformation. Wave propagation in plasma around the reentry vehicle was numerically simulated using a higher-order FDTD method by White [11]. The frequency-dependent finite-difference time-domain (FD2TD) method [12] has also been an effective tool for simulating electromagnetic waves in a frequency-dependent medium [13]. An analysis approach to predict RF blackout and plasma attenuation for a reentry vehicle, based on a method combining CFD and FD2TD, has been suggested by Takahashi et al. [6,14,15]. However, the problem of uncertainties in using physical and numerical models for plasma properties in high temperature regions still remains. Validations of the combined model for plasma flow and electromagnetic waves have been conduced.

Flight experiment results for the atmospheric reentry demonstrator (ARD) [16] of the European Space Agency (ESA) have made a good target for validating the analysis model. In [14], a comparison with the analysis results and measured data from the ARD mission was made. However, that validation was insufficient because the actual flight environment was incompletely reproduced: for example, the angle of attack and a realistic configuration were not reproduced owing to the use of a two-dimensional axisymmetric flow solver for the plasma flow simulation. Takahashi [15] recently developed an RF blackout prediction tool combining RG-FaSTAR, which is a high-enthalpy flow solver applicable to unstructured grid systems and FD2TD simulation code. Using the combined model, it is possible to investigate the behavior of electromagnetic waves in reentry plasma. The present research objective is to reproduce the actual reentry environment encountered during the ARD flight and to clarify and acquire information about the plasma flow field and electromagnetic wave behavior. Because the ARD had a typical reentry capsule configuration and reentry trajectory, it is worth investigating further and providing the distributions of plasma and electromagnetic waves for the design and development of new reentry vehicles and for research and development concerning RF blackout mitigation techniques.

2. Reentry Vehicle

Detailed information about the reentry mission considered herein was reported in [16]. In this section, a brief explanation of the mission and vehicle configuration is presented. The atmospheric reentry demonstrator (ARD) of the European Space Agency (ESA) was launched by the Ariane-5 V503 rocket on 21 October 1998. The reentry capsule separated from the Ariane-5 rocket at an altitude of 218 km. ARD reached its apogee point at an altitude of 830 km and then reentered the Earth’s atmosphere. The reentry velocity at an altitude of 120 km was approximately 7450 m/s. A strong shock wave was generated in front of the vehicle owing to its high reentry velocity. The vehicle was expected to be surrounded by high-temperature plasma. After deceleration by aerodynamic forces, the ARD began its recovery sequence. Finally, the ARD splashed down at 133${}^{\circ }$ West and 3.7${}^{\circ }$ North in the Pacific Ocean and was successfully recovered by a recovery ship. The reentry trajectory, i.e., position, velocity and attitude, reconstructed by acceleration measurements was obtained using on-board inertial navigation system measurements and a GPS receiver. Figure 2a,b shows the altitude profiles of the reentry vehicle during the entry, decent and landing stage versus the longitude west of Greenwich and the latitude north, respectively.

The schematic configuration of the ESA ARD is shown in Figure 3. Its configuration was basically the same as the Apollo command module scaled down to 70%, which is typical for a blunt body reentry capsule. It comprised three key components: the front shield section, rear cone section and back cover section. The capsule was essentially a spherical cone with a nose radius of 3360 mm and a base radius of 2800 mm for the front shield section with a half-angle of 33${}^{\circ }$ for the rear section. To protect it from strong aerodynamic heating, an ablator material was coated onto the front heat shield. The ARD was equipped with six antennae and two reflectometers for a telemetry link and communication with a tracking and data relay satellite (TDRS), GPS satellites and airplanes. The antennae were positioned on the vehicle surface at the rear, 1.043 m from the top of the front heat shield. The link frequencies for the communication waves were approximately 2.267 GHz for TDRS, 1.575 GHz for GPS and 2.203 and 2.218 GHz for the two airplanes links. RF blackout for the TDRS link was not observed, while the blackout caused for the GPS and airplanes links was observed. However, for the TDRS link, strong plasma attenuation due to reentry plasma around the ARD was measured from an altitude of 86 km down to 44 km. These results on RF blackout and plasma attenuation of GPS and airplane links indicated that the blackout duration of the receivers in the forward position on the ARD is longer, whereas that of receivers on the rear position of the vehicle is shorter. This is because the shock layer in front of the vehicle strongly prevented electromagnetic waves used for communication from propagating in the forward direction relative to the ARD.

Figure 2. Reentry trajectory of the ESA atmospheric reentry demonstrator (ARD). (a) Altitude-longitude profile; (b) altitude-latitude profile.

Figure 2. Reentry trajectory of the ESA atmospheric reentry demonstrator (ARD). (a) Altitude-longitude profile; (b) altitude-latitude profile.

Figure 3. Analysis model of the ESA ARD. (a) Front side view; (b) rear side view.

Figure 3. Analysis model of the ESA ARD. (a) Front side view; (b) rear side view.

3. Flow Field Modeling

3.1. Governing Equation

In this study, the following assumptions for the flow field simulation are made: (1) The flow is laminar, steady and a continuum; (2) The chemical constituent of freestream gas is air; (3) The flow field is in thermochemical nonequilibrium, and the temperature is separated into translational, rotational, vibrational and electron temperatures. Note that the electronic excitation energy mode and the electron energy are assumed to be equilibrated.

The flow field around the ESA AED is expressed by the Navier–Stokes equations with a multi-temperature model. The governing equations include total mass, momentum, total energy, species mass, rotation energy, vibration energy and electron energy conservations. The equation system can be written in vector form as follows:

$$\frac{\partial \mathbf{Q}}{\partial t}+\frac{\partial \mathbf{F}}{\partial x_j}=\frac{\partial {\mathbf{F}}_v}{\partial x_j}+\mathbf{W}$$

(1)

The conservative vector is expressed as $\mathbf{Q}=\left[\rho ,\rho u_i,E,{\rho }_s,E_{\mathrm{rot}},E_{\mathrm{vib}},E_{\mathrm{e}}\right]$. Furthermore, $\mathbf{F}$ and ${\mathbf{F}}_v$ are the vectors of the inviscid and viscous terms, respectively, and $\mathbf{W}$ denotes the vector of the source term. The equation of state can be described by:

$$p=\sum _{s\ne \mathrm{e}}^{\mathrm{ns}-1}{\rho }_s{R}_s{T}_{\mathrm{trs}}+{\rho }_{\mathrm{e}}{R}_{\mathrm{e}}{T}_{\mathrm{ele}}=\rho \widehat{R}{T}_{\mathrm{tr}}+p_{\mathrm{e}}$$

(2)

The translational, rotational, vibrational and electron energies constituting the internal energy are given by:

$$\begin{array}{cc}& E_{\mathrm{trs}}=\sum _{s\ne \mathrm{e}}^{ns}\frac{3}{2}{\rho }_s{R}_s{T}_{\mathrm{trs}}\hfill \end{array}$$

(3)

$$\begin{array}{cc}& E_{\mathrm{rot}}=\sum _{s=\mathrm{M}}^{nm}{\rho }_s{R}_s{T}_{\mathrm{rot}}\hfill \end{array}$$

(4)

$$\begin{array}{cc}& E_{\mathrm{vib}}=\sum _{s=\mathrm{M}}^{nm}\frac{{\rho }_s{R}_s{\Theta }_{\mathrm{vib},s}}{\exp \left({\Theta }_{\mathrm{vib},s}/T_{\mathrm{vib}}\right)-1}\hfill \end{array}$$

(5)

$$\begin{array}{cc}& E_{\mathrm{ele}}=\frac{3}{2}{\rho }_{\mathrm{e}}{R}_{\mathrm{e}}{T}_{\mathrm{ele}}\hfill \end{array}$$

(6)

where $\Delta h_s^0$ and ${\Theta }_{\mathrm{vib},s}$ represent the enthalpy of formation and the vibrational characteristic temperatures, respectively. In addition, ${\delta }_{i,j}$, $S_{\mathrm{int}}$ and ${\dot{w}}_s$ present the Kronecker delta, the internal energy exchange rate and the mass production rate, respectively. The subscript “s” represents each chemical species, and “M” means the set of molecules ($M={\mathrm{N}}_2$, ${\mathrm{O}}_2$, $\mathrm{NO}$, ${\mathrm{N}}_2^{+}$, ${\mathrm{O}}_2^{+}$ and ${\mathrm{NO}}^{+}$). The total internal energy can be defined as:

$$\begin{array}{c}\hfill E=E_{\mathrm{trs}}+E_{\mathrm{rot}}+E_{\mathrm{vib}}+E_{\mathrm{ele}}+\sum _{s=1}^{ns}{\rho }_s\Delta h_s^0+\frac{1}{2}\rho u_i{u}_i.\end{array}$$

(7)

3.2. Transport Properties

Transport properties, such as the viscosity, thermal conductivity and binary diffusion coefficients, for a gas mixture are evaluated by Yos’s model [17], which is based on the first Chapman–Enskog approximation [18]. The collision cross-sections are given by Gupta [19] and Fertig [20,21]. The diffusion coefficients are expressed by Curtiss and Hirschfelder’s formula [22]. Ambipolar diffusion is assumed for charged species as $D_s^a=(1+T_{\mathrm{e}}/T_{\mathrm{tr}})D_s$, where $D_s$ is the effective diffusion coefficient of the ionic species.

3.3. Chemical Reactions

For chemical reactions in high-temperature air, the test gas is assumed to consist of 11 chemical species (${\mathrm{N}}_2$, ${\mathrm{O}}_2$, $\mathrm{NO}$, ${\mathrm{N}}_2^{+}$, ${\mathrm{O}}_2^{+}$, ${\mathrm{NO}}^{+}$, $\mathrm{N}$, $\mathrm{O}$, ${\mathrm{N}}^{+}$, ${\mathrm{O}}^{+}$ and ${\mathrm{e}}^{-}$), and 49 reactions are assumed to occur. The chemical reactions are composed of heavy-particle impact dissociation, electron impact dissociation, associative ionization, the nitric oxide (NO) exchange reaction, electron-impact ionization and charge exchange reactions in this analysis. The chemical reaction rate is determined using an Arrhenius-type form. The reaction rate coefficients are obtained from Park [23]. Moreover, the backward reaction rates are evaluated from the corresponding equilibrium constants. The equilibrium constants, being functions of temperature only, are calculated by the curve-fitting formula given in [19,24].

3.4. Internal Energy Exchange

We considered the energy transfer between each of the internal energy modes: translation-rotation (T-R) [25], translation-vibration (T-V) [26,27], translation-electron (T-e) [28,29,30], rotation-vibration (R-V) [31], rotation-electron (R-e) [32,33] and vibration-electron (V-e) [34]. The energy losses/releases for vibrations and rotations associated with chemical reactions reveal the dissociation energies for heavy-particle impact reactions. These can be obtained using a non-preferential dissociation model [30]. In addition, the electron energy losses/releases due to electron-impact dissociation and ionization are also considered.

3.5. Implementation

The governing equations for the flow field are transformed into the delta form and solved using a finite-volume approach. The flow properties are basically defined at the center of a control volume. The numerical fluxes of the advection term are calculated using the Harten, Lax, van Leer and Einfeldt (HLLE) scheme [35] using the monotonic upstream-centered scheme for interpolation of the conservation laws (MUSCL). The viscous terms of the equation are calculated using the second-order central differencing method. The spatial gradients of the flow variables are evaluated with the Green–Gauss method. A time integration is performed implicitly using the lower-upper symmetric Gauss–Seidel (LU-SGS) method [36] combined with the point-implicit method [37].

These physical and numerical models are implemented in the high-enthalpy flow solver RG-FaSTAR [15]. RG-FaSTAR is a version of the fast unstructured CFD code FaSTAR [38] originally developed by JAXA that incorporates real gas effects. These solvers can handle unstructured grid systems, such as tetrahedral, triangular prism, square pyramid and hexahedral meshes, and have high flexibility in generating the computational grid. For massive parallel computation on high-performance computers, the message-passing interface (MPI) technique with a domain partition approach is used.

3.6. Computational and Boundary Conditions

Numerical analyses of flow fields by RG-FaSTAR are performed for the cases as listed in Table 1. Note that AOA denotes the angle of attack of the ESA ARD. The freestream velocity is determined according to Tran’s report [16].

Table 1. Freestream parameters. AOA, angle of attack.

Altitude, km	Density, kg/m${}^{\mathbf{3}}$	Temperature, K	Velocity, m/s	AOA, Degree
85.0	8.18 × 10${}^{-6}$	191.0	7577	20.0
80.0	1.85 × 10${}^{-5}$	195.8	7609	20.0
75.0	4.07 × 10${}^{-5}$	201.7	7593	19.2
70.0	8.83 × 10${}^{-5}$	210.9	7542	19.2
60.0	3.40 × 10${}^{-4}$	242.0	6105	19.4
50.0	1.15 × 10${}^{-3}$	265.3	4567	20.0

Table 1. Freestream parameters. AOA, angle of attack.

Altitude, km	Density, kg/m${}^{\mathbf{3}}$	Temperature, K	Velocity, m/s	AOA, Degree
85.0	8.18 × 10${}^{-6}$	191.0	7577	20.0
80.0	1.85 × 10${}^{-5}$	195.8	7609	20.0
75.0	4.07 × 10${}^{-5}$	201.7	7593	19.2
70.0	8.83 × 10${}^{-5}$	210.9	7542	19.2
60.0	3.40 × 10${}^{-4}$	242.0	6105	19.4
50.0	1.15 × 10${}^{-3}$	265.3	4567	20.0

Computational grids for the ARD in RG-FaSTAR at an altitude of 85 km are composed of tetrahedral meshes, prismatic layer meshes near the wall, hexahedral meshes in the front region and squared pyramids, as shown in Figure 4. To finely resolve the shock layer near the ESA ARD, hexahedral meshes are inserted into its front region. In addition, prismatic meshes near walls are used to solve the boundary layer near the wall. Square pyramid meshes are used to connect the hexahedral meshes and the prism meshes. The numbers of cells and nodes of the grid system at an altitude of 85 km are 4,244,601 and 2,434,920, respectively. At other altitudes, the number of cells and nodes are different, but the composition of mesh types is basically similar.

Figure 4. Computational grids around the ESA ARD for plasma flow simulation. (a) The whole region; (b) the front region.

Figure 4. Computational grids around the ESA ARD for plasma flow simulation. (a) The whole region; (b) the front region.

For boundary conditions, inflow, outflow and wall conditions are considered as shown in Figure 5. The inflow and outlet boundaries in this figure are illustrated in a half-cut view. Freestream density, temperature and velocity, listed in Table 1, are given at the inflow boundary. Mass fractions of chemical species are always set to 0.765 for N${}_2$ and 0.235 for O${}_2$ at the inflow. An angle of attack of approximately 20 degrees is assumed for each case. At the outflow boundary, all flow parameters, e.g., density, velocity, temperature, etc., are determined by zeroth-order extrapolation. On the wall surface of the vehicle, there is no pressure gradient in the wall normal direction, and the no-slip wall condition applies. In addition, the wall temperature is fixed at 273 K as a cooling wall condition. The wall surface is assumed to be non-catalytic for chemical species.

Figure 5. Computational domains and boundary conditions around the ESA ARD for plasma flow simulation.

Figure 5. Computational domains and boundary conditions around the ESA ARD for plasma flow simulation.

As a study of grid dependency, plasma flow simulation at an altitude of 85 km was conducted using a fine mesh with cells of 7,943,085 and nodes of 4,190,013. Profiles of the number density of electron with the fine mesh in the rear region where it is important to predict electromagnetic wave behavior show good agreement with analytical results obtained using the present computational mesh.

4. Electromagnetic Wave Modeling

4.1. Maxwell’s Equations

To predict RF blackout and plasma attenuation during communication between the reentry vehicle and data relay satellite/ground station, it is important to clarify the behavior of electromagnetic waves in the plasma region. In general, Maxwell’s equations for describing the behavior of electromagnetic waves are expressed as follows:

$$\begin{array}{cc}& \nabla \times \mathbf{E}=-\frac{\partial \mathbf{B}}{\partial t}\hfill \end{array}$$

(8)

$$\begin{array}{cc}& \nabla \times \mathbf{H}=\frac{\partial \mathbf{D}}{\partial t}+\mathbf{J}\hfill \end{array}$$

(9)

$$\begin{array}{cc}& \nabla ·\mathbf{D}={\rho }_{\mathrm{em}}\hfill \end{array}$$

(10)

$$\begin{array}{cc}& \nabla ·\mathbf{B}=0\hfill \end{array}$$

(11)

In the present formulation of the FDTD method, Faraday’s law in Equation (8) and Ampere’s law in Equation (9) are used, while Gauss’s laws in Equations (10) and (11) are not used directly.

4.2. FD2TD Method

The FDTD method has been widely used for simulating electromagnetic wave propagation. However, in general, plasmas are frequency-dependent media whose dielectric constants and electrical conductivities vary over a frequency range. The FDTD formulation assumes that these medium properties are independent of frequency, and this leads to the risk of obtaining incorrect results when applying it to frequency-dependent media. On the other hand, the FD2TD method is effective for numerical simulation of Maxwell’s equations in such frequency-dependent media.

In the present FD2TD formulation, the complex relative permittivity, ${\epsilon }_{\mathrm{r}}^{*}$, assuming the first Drude dispersion, is given by:

$${\epsilon }_{\mathrm{r}}^{*}=1+\frac{{\omega }_{\mathrm{p}}^2}{\omega \left(i{\nu }_c-\omega \right)}=1+\chi \left(\omega \right)$$

(12)

where ${\omega }_{\mathrm{p}}$, ${\nu }_{\mathrm{c}}$ and $\chi \left(\omega \right)$ respectively represent the plasma angular frequency, electron collision frequency and electric susceptibility at a given angular frequency, ω. Note that i represents the imaginary unit ($i=\sqrt{-1}$). In the present analysis model, the plasma angular frequency ${\omega }_{\mathrm{p}}$ and electron collision frequency ${\nu }_{\mathrm{c}}$ are given by:

$$\begin{array}{cc}& {\omega }_{\mathrm{p}}=\sqrt{\frac{e^2{n}_{\mathrm{e}}}{m_{\mathrm{e}}{\epsilon }_0}}\hfill \end{array}$$

(13)

$$\begin{array}{cc}& {\nu }_{\mathrm{c}}=\sum _{s\ne \mathrm{e}}^{ns}{n}_s\pi {\Omega }_{s,\mathrm{e}}\sqrt{\frac{8k{T}_{\mathrm{ele}}}{\pi m_{\mathrm{e}}}}\hfill \end{array}$$

(14)

where e, $n_s$, $m_{\mathrm{e}}$ and k are the elementary charge, the number density of species s, the mass of electron and the Boltzmann constant, respectively. The collision cross-section Ω is obtained by Gupta’s model [19] and Fertig’s model [20,21], which are calculated by the same routine of the transport properties’ evaluation in RG-FaSTAR.

The complex relative permittivity ${\epsilon }_{\mathrm{r}}^{*}$ can be expressed with the relative permittivity ${\epsilon }_{\mathrm{r}}$ and dielectric conductivity σ as follows:

$${\epsilon }_{\mathrm{r}}^{*}={\epsilon }_{\mathrm{r}}+\frac{\sigma }{i\omega {\epsilon }_0}$$

(15)

where:

$$\begin{array}{c}\hfill {\epsilon }_{\mathrm{r}}=1-\frac{{\omega }_{\mathrm{p}}^2}{{\omega }^2+{\nu }_{\mathrm{c}}^2}\end{array}$$

(16)

$$\begin{array}{c}\hfill \sigma =\frac{{\epsilon }_0{\omega }_{\mathrm{p}}^2}{{\omega }^2+{\nu }_{\mathrm{c}}^2}{\nu }_{\mathrm{c}}\end{array}$$

(17)

Assuming that the electric field is constant during a time step ($\Delta t$), the electric flux density at the n-th step is expressed as follows:

$$\begin{array}{ccc}\hfill {\mathbf{D}}^n& =& {\epsilon }_0{\mathbf{E}}^n+{\epsilon }_0{\int }_0^{n\Delta t}\chi \left(\tau \right)\mathbf{E}(n\Delta t-\tau )d\tau \hfill \\ & =& {\epsilon }_0\left(1+{\chi }^0\right){\mathbf{E}}^n+{\epsilon }_0\sum _{m=0}^{n-2}{\mathbf{E}}^{n-1-m}{\chi }^{m+1}\hfill \end{array}$$

(18)

where ${\chi }^m$ is defined using the electric susceptibility $\chi \left(\tau \right)$ in the time region as follows:

$${\chi }^m={\int }_{m\Delta t}^{(m+1)\Delta t}\chi \left(\tau \right)d\tau$$

(19)

Applying the Fourier transformation to $\chi \left(\omega \right)$ yields:

$${\chi }^m=\frac{{\omega }_{\mathrm{p}}^2}{{\nu }_{\mathrm{c}}^2}\left\{{\nu }_{\mathrm{c}}\Delta t-\left(1-e^{-{\nu }_{\mathrm{c}}\Delta t}\right)e^{-{\nu }_{\mathrm{c}}m\Delta t}\right\}$$

(20)

Thus, the difference of the electric flux density in the time direction is given as follows:

$${\mathbf{D}}^n-{\mathbf{D}}^{n-1}={\epsilon }_0\left(1+{\chi }^0\right){\mathbf{E}}^n-{\epsilon }_0{\mathbf{E}}^{n-1}-{\epsilon }_0\sum _{m=0}^{n-2}{\mathbf{E}}^{n-1-m}\Delta {\chi }^m$$

(21)

where:

$$\Delta {\chi }^m={\chi }^m-{\chi }^{m+1}=-\frac{{\omega }_{\mathrm{p}}^2}{{\nu }_{\mathrm{c}}^2}{\left(1-e^{-{\nu }_{\mathrm{c}}\Delta t}\right)}^2{e}^{-{\nu }_{\mathrm{c}}m\Delta t}$$

(22)

Equation (9) can be expressed using Equation (21) as follows:

$${\mathbf{E}}^n=\frac{1}{1+{\chi }^0}{\mathbf{E}}^{n-1}+\frac{1}{1+{\chi }^0}{\Phi }^{n-1}+\frac{\Delta t/{\epsilon }_0}{1+{\chi }^0}\nabla \times {\mathbf{H}}^{n-\frac{1}{2}}$$

(23)

and:

$${\Phi }^{n-1}={\mathbf{E}}^{n-1}\Delta {\chi }^0+e^{-{\nu }_{\mathrm{c}}\Delta t}{\Phi }^{n-2},n\ge 2$$

(24)

As initial conditions, we set ${\mathbf{E}}^0=0$, ${\mathbf{H}}^{-\frac{1}{2}}=0$ and ${\Phi }^{-1}={\Phi }^0=0$ In the present study, the magnetic permeability (μ) is assumed not to be dispersive. Hence, Equation (8) is discretized as follows:

$${\mathbf{H}}^{n+\frac{1}{2}}={\mathbf{H}}^{n-\frac{1}{2}}-\frac{\Delta t}{\mu }\nabla \times {\mathbf{E}}^n.$$

(25)

The rotation terms in Equations (23) and (25) are discretized by the central difference method in the spatial direction. The components of the electric and magnetic field vectors are staggered on the Yee cell [39]. Therefore, the rotation of the electric fields in Equation (25) is calculated with coplanar magnetic fields, and that of the magnetic fields in Equation (23) is calculated with coplanar electric fields.

4.3. Computational Domain

Figure 6 shows the computational domain in the FD2TD approach. Note that the outer boundaries in the figure are illustrated by using a half-cut. The computational domain for electromagnetic waves is set to be a cuboid 7 m in length in the x direction, 6 m in width in the y direction and 6 m in height in the z direction. The computational grid system is fully composed of hexahedral meshes, and the grid nodes are set to 500 (x-direction), 300 (y-direction) and 300 (z-direction), i.e., 45,000,000 nodes. The ESA ARD is contained within the computational domain. The origin of the computational domain used in the CFD approach corresponds to that of the simulation of the electromagnetic waves. The physical properties (${\omega }_{\mathrm{p}}$ and ${\nu }_{\mathrm{c}}$) required in the present FD2TD method are obtained by reference to the computational results obtained by the CFD technique. On the other hand, the plasma parameters at a node in the FD2TD system that is positioned outside the CFD cell are obtained based on the parameters at the nearest boundary point in the CFD domain. The magnetic permeability (μ) is assumed to be the same as that of free space (${\mu }_0$) for simplicity.

Mur’s first-order absorption boundary [40] is imposed as a boundary condition at all six faces of the computational domain. On the other hand, a conduction boundary condition is imposed on the surface of the reentry vehicle.

Figure 6. Computational domains and boundary conditions around the ESA ARD for plasma flow simulation.

Figure 6. Computational domains and boundary conditions around the ESA ARD for plasma flow simulation.

On the antenna, transverse magnetic waves are emitted as a wave source by inputting a sinusoidal electric field.

The physical and numerical models are implemented in the electromagnetic wave simulation code “Arcflow/Arcwave” developed by Hokkaido University. For fast computation and efficient allocation of adequate memory, the FD2TD simulation code is parallelized using the MPI and OpenM techniques. A data transfer tool between the electromagnetic wave simulation code and RG-FaSTAR is also implemented. Both of the simulation codes using the CFD and FD2TD methods are packaged.

5. Results and Discussion

5.1. Stagnation Pressure

Comparison of the stagnation pressure coefficient in front of the ESA ARD between the flight experiment and numerical analysis is conducted to validate the present plasma flow analysis model. The pressure obtained in the flight experiment is reconstructed based on the flight information on velocity and density deduced by the on-board acceleration measurement technique. In this analysis, the pressure coefficient at stagnation $C_p$ is defined as the maximum pressure on the vehicle’s front surface, which is obtained by the following expression:

$$C_p=\frac{p_0}{\frac{1}{2}{\rho }_{\infty }{U}_{\infty }^2}$$

(26)

where $p_0$ is the maximum pressure on the vehicle’s front surface and ${\rho }_{\infty }$ and $U_{\infty }$ are the freestream density and velocity, respectively. Figure 7 shows a comparison of the stagnation pressure coefficient profiles with altitude. Error between the flight experiment results and the present numerical analysis is approximately 3% at most, except for an altitude of 85 km. This indicates that the analysis can correctly reproduce the aerodynamic tendencies in the shock layer.

Figure 7. Comparison of the pressure coefficient at stagnation at the vehicle’s front between the flight experiment and numerical analysis.

Figure 7. Comparison of the pressure coefficient at stagnation at the vehicle’s front between the flight experiment and numerical analysis.

5.2. Plasma Flow Field

A reentry vehicle is generally surrounded by high-temperature gas generated by strong shock waves encountered during atmospheric reentry. When freestream density is high and collisions between molecular particles are frequent, the flow field is closed in thermal and chemical equilibrium, whereas the flow field can be in a thermochemically nonequilibrium state when the freestream density is low. Because RF blackout tends to be caused in high altitude and high Mach number regions, the flow field is not necessarily in equilibrium. Thus, to accurately predict RF blackout and plasma attenuation, analysis considering nonequilibrium models may be required. In the present analysis methods implemented in the RG-FaSTAR code, thermochemically-nonequilibrium models are introduced as mentioned above. Figure 8a,b shows the distributions of translational and vibrational temperatures around the ESA ARD at an altitude of 85 km, respectively. In addition, distributions of pressure and velocity magnitude at an altitude of 85 km are presented in Figure 9 and Figure 10, respectively. A strong shock wave is generated in front of the vehicle, and a wide wake region appears in the rear. Temperatures rapidly increase across the shock wave at the vehicle’s front. An equilibrium region appears near the surface in the shock layer owing to highly dense gas. Then, the flow field passing through a region near the vehicle’s shoulder is rarefied and cooled due to expansion in the rear region. Translational and vibrational temperatures are completely separated in the whole region except for a region near the wall at the vehicle’s front, and strong nonequilibrium appears. In the wake region, the translational temperature decreases rapidly due to supersonic expansion, causing a transfer from thermal to kinetic energy while the vibrational temperature remains at high-temperature values. This is mainly because relaxation between the translational and vibrational energy is slowed by the low-density environment at the rear of the vehicle.

Figure 8. Temperature distributions around the ESA ARD at an altitude of 85 km (x-y plane). (a) Translational temperature; (b) vibrational temperature.

Figure 8. Temperature distributions around the ESA ARD at an altitude of 85 km (x-y plane). (a) Translational temperature; (b) vibrational temperature.

Figure 9. Pressure distributions around the ESA ARD at an altitude of 85 km (x-y plane).

Figure 9. Pressure distributions around the ESA ARD at an altitude of 85 km (x-y plane).

Figure 10. Velocity magnitude distributions around the ESA ARD at an altitude of 85 km (x-y plane).

Figure 10. Velocity magnitude distributions around the ESA ARD at an altitude of 85 km (x-y plane).

To investigate electromagnetic wave behavior in plasma, it is important to clarify the number density of electrons around the reentry vehicle, because the propagation of electromagnetic waves is strongly affected by the distribution of electrons. Validation of the analysis model of RG-FaSTAR in predicting the number density of electrons during atmospheric reentry was performed using measured results from electrostatic probes in the radio attenuation measurement (RAM C) program, as reported in [15]. The distributions of the electron number density around the ESA ARD at altitudes of 85, 70 and 60 km are shown in Figure 11. This figure shows the computed results on the x-y and x-z planes at various altitudes. Because the link frequency of the TDRS is 2.267 GHz, the critical number density of electrons for the link frequency is approximately 6.4 × 10${}^{16}$ 1/m${}^3$, i.e., approximately 16.8 in the common logarithmic scale. When electromagnetic waves pass through the regions where the electron density exceeds the critical value for the link frequency, attenuation, refraction or reflection occurs. In front of the vehicle, a region appears to exist in which the electron number density generated by the shock layer remains high at all altitudes. On the other hand, at the rear of the vehicle, the plasma density decreases to well below the critical value, except in the case of 70 km in altitude. This indicates that electromagnetic waves can propagate in a path through the wake. At an altitude of 70 km, the reentry vehicle is completely surrounded by highly dense plasma. There exists a high electron density region at the rear of the ESA ARD at an altitude of 60 km, as shown in Figure 11e. This is because of a recompression region in which the plasma and the wake flowing around the vehicle’s shoulder collide at the rear. It is indicated that the recompression gas at the rear plays an important role in RF blackout and plasma attenuation, increasing the electron density and blocking the propagation of electromagnetic waves.

Figure 11. Number density of electrons around the ESA ARD at each altitude. (a) Results at an altitude of 85 km (x–y plane); (b) results at an altitude of 85 km (x–z plane); (c) results at an altitude of 70 km (x–y plane); (d) results at an altitude of 70 km (x–z plane); (e) results at an altitude of 60 km (x–y plane); (f) Results at an altitude of 60 km (x–z plane).

Figure 11. Number density of electrons around the ESA ARD at each altitude. (a) Results at an altitude of 85 km (x–y plane); (b) results at an altitude of 85 km (x–z plane); (c) results at an altitude of 70 km (x–y plane); (d) results at an altitude of 70 km (x–z plane); (e) results at an altitude of 60 km (x–y plane); (f) Results at an altitude of 60 km (x–z plane).

5.3. Electromagnetic Waves

The behaviors of electromagnetic waves determined using the FD2TD analysis code are presented in this section. As mentioned above, antennae emitting communication waves are positioned on the surface of the rear cone section of the ARD. Propagation of communication waves is strongly prevented in the vehicle’s front region by highly dense plasma. Thus, the vehicle basically conducts communication with satellites and airplanes through the rear region, where the plasma density becomes relatively low. Figure 12 shows the time sequence of isosurfaces of the electric field magnitude from wave emission around the ESA ARD at an altitude of 85 km. At an elapsed time of 1.875 × 10${}^{-9}$ s, the electromagnetic waves can propagate in all directions. However, reflection and attenuation by the highly dense plasma start to appear at an elapsed time of 3.75 × 10${}^{-9}$ s. From an elapsed time of 7.5 × 10${}^{-9}$ s, the waves partly reach the boundary surface of the computational domain; then, the isosurface partly vanishes. It is confirmed that the electromagnetic waves diffract on the edge of the front shock wave at 9.375 × 10${}^{-9}$ s and partly propagate along the frontal plasma. On the other hand, electromagnetic waves propagate in the rear region without strong attenuation.

Figure 12. Time sequence of isosurfaces of the electric field around the ESA ARD at an altitude of 85 km.

Figure 12. Time sequence of isosurfaces of the electric field around the ESA ARD at an altitude of 85 km.

Figure 13a–d shows the electric field magnitudes of the electromagnetic waves for the cases of free space, an altitude of 85 km, an altitude of 60 km and an altitude of 50 km, respectively. Note that the magnitude of the electric field is normalized to 1.0 on the antenna here. All values of the electric field in the figure are expressed on a logarithmic scale. The electric field distribution in free space where there is no plasma is presented to compare with the results for the cases of other altitudes. For the free-space condition, the electromagnetic waves can propagate in all directions, diffracting and reflecting on the ARD surface. On the other hand, as confirmed in the time sequence of electromagnetic wave propagation, electromagnetic waves are completely reflected below an altitude of 85 km and cannot propagate on the front side of the ARD along the shock layer surface. For an altitude of 70 km in the computational results, electromagnetic waves cannot propagate in any direction, and it is observed that a complete RF blackout occurs, although this is not shown here. However, according to Tran’s work in [16], this RF blackout was not observed in communication between the ARD and the TDRS using a link frequency of 2.267 GHz during atmospheric reentry, although strong plasma attenuation appeared at an altitude of 70 km. Uncertainties in chemical reaction models and internal energy exchange models used in this simulation code may cause this discrepancy, as pointed out by Takahashi et al. [14]. On the other hand, because uncertainty analysis incurs a large computational cost when using a three-dimensional high-enthalpy flow solver, this issue is left for future works. At an altitude of 60 km, wave propagation recovers from the RF blackout, while strong attenuation and reflection in the plasma wake are caused by the recompression gas region. The number density of electrons becomes low in the ARD’s rear region at an altitude of 50 km, and electromagnetic waves propagate in a similar fashion to the computed results for the free-space condition.

Figure 13. Electromagnetic waves (magnitude of electric field) around the ESA ARD at an elapsed time of 3.75 × 10${}^{-8}$ s. (a) Results for free space (x–y plane); (b) results at an altitude of 85 km (x–y plane); (c) results at an altitude of 60 km (x–y plane; (d) results at an altitude of 50 km (x–y plane).

Figure 13. Electromagnetic waves (magnitude of electric field) around the ESA ARD at an elapsed time of 3.75 × 10${}^{-8}$ s. (a) Results for free space (x–y plane); (b) results at an altitude of 85 km (x–y plane); (c) results at an altitude of 60 km (x–y plane; (d) results at an altitude of 50 km (x–y plane).

6. Conclusions

Numerical analyses of plasma flow fields and communication waves around the Atmospheric Reentry Demonstrator(ARD) of the European Space Agency (ESA) were conducted for several cases along its reentry trajectory for the prediction of radio frequency blackout and plasma attenuation. The capsule configuration and reentry trajectory used were typical for reentry missions with blunt body reentry capsules. Formulations of nonequilibrium high-enthalpy flow based on the Navier–Stokes equations and electromagnetic waves using Maxwell’s equations were presented. These physical models were implemented with the high-enthalpy flow solver “RG-FaSTAR” based on the computational fluid dynamics (CFD) and electromagnetic wave simulation code “Arcflow/Arcwave”, which uses the frequency-dependent finite-difference time-domain (FD2TD) method. A detailed computational model to reproduce the reentry environment of the ESA ARD was reconstructed with unstructured grids for the CFD code and rectangular grids for the FD2TD code. Plasma flow properties obtained using the CFD code were used in the FD2TD simulation to express the features of the dispersion medium, i.e., plasma, around the ESA ARD. Computed results were validated with experimental data for stagnation pressure profiles during atmospheric reentry, and good agreement was observed. Plasma flow properties, i.e., temperature, pressure and velocity around the reentry vehicle, were investigated. Then, distributions of the electron number density were presented, and radio frequency blackout and plasma attenuation were discussed using information from the plasma flow field. The time sequences of electromagnetic wave propagation behavior at various altitudes were clarified. It was confirmed that reflection, diffraction and attenuation on the strong shock wave surface in front of the vehicle occur.

In the present study, electromagnetic wave behavior in plasma with the combined CFD and FD2TD method were investigated in detail. The computed results and analysis model may provide effective insight into the design of antennae and for research on RF blackout mitigation techniques for reentry vehicles. It was indicated that the radio frequency blackout and plasma attenuation situation for the blunt body capsules widely used in reentry missions could be investigated with the present model.