3.1 Non-equilibrium magnetohydrodynamic equation

The governing equations in this paper are composed of thermochemical nonequilibrium coupling MHD equations. Consequently, the effect of magnetic interactions is described by incorporating a kinetic term into the momentum equation and a Joule heating term into the energy equation. The 2-dimensional axisymmetric thermochemical non-equilibrium MHD equation is derived based on the above premise. Due to the low magnetic Reynolds number (Rm << 1), the induced magnetic field can be ignored. The magnetic Reynolds number is defined by equation [Reference Nagata, Otsu and Yamada25]:

(1) \begin{equation}Rm = {\mu _0}\sigma {V_\infty }{L_{{\rm{ref}}}}\end{equation}

where μ ₀ is magnetic permeability in vacuum, σ is electrical conductivity, V _∞ is freestream velocity, and L _ref is the radius of the cone of the vehicle.

The mass, momentum, total energy, and vibrational electronic energy equations are given below as in Ref. (Reference Assonitis, Paciorri and Bonfiglioli30–Reference Luo, Wu and Chang33).

(2) \begin{equation}\frac{{\partial {\rho _s}}}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}\left( {{\rho _s}{u_j} + {J_{s,j}}} \right) = {\dot \omega _s}\end{equation}

(3) \begin{equation}\frac{\partial }{{\partial t}}\left( {\rho {u_i}} \right) + \frac{\partial }{{\partial {x_j}}}\left( {\rho {u_i}{u_j} + p{\delta _{ij}} - {\tau _{ij}}} \right) = {(J \times B)_i}\end{equation}

(4) \begin{equation}\frac{{\partial {E_t}}}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}\left( {({E_t} + p){u_j} - {\tau _{ij}}{u_i} + {q_j} + \sum\limits_{s = 1}^{ns} {{J_{s,j}}} {h_s}} \right) = J \cdot E\end{equation}

(5) \begin{equation}\frac{{\partial {E_{ve}}}}{{\partial t}} + \frac{\partial }{{\partial {x_j}}}\left( {{E_{ve}}{u_j} + {q_{ve,j}} + \sum\limits_{s = 1}^{ns} {{J_{s,j}}} {h_{ve,s}}} \right) + {p_e}\frac{{\partial {u_j}}}{{\partial {x_j}}} = {Q_{ve}} + {\gamma _{ev}}\frac{{J \cdot J}}{\sigma }\end{equation}

where $\rho$ _s represents the density of a species of s. u _i and u _j are velocity vectors in the i and j directions, respectively. p is the total pressure, and p _e is the partial pressure of the electron. The total and vibro-electronic energies per unit volume are denoted as E _t and E _ve , respectively. J _s,j is the mass diffusion flux of component s in the j direction. ${\dot \omega _s}$ is the source term of the chemical reactions. ${\tau _{ij}}$ is the viscous shear tensor. q _j and q _ve,j represent the total heat conduction vector and the vibration-electron heat conduction vector in the j direction, respectively. B is the magnetic field strength. Total enthalpy and vibration-electron enthalpy are expressed by h _s and h _ve,s , respectively. Q _ve is the vibration-electron energy source term. (J × B) _i represents the Lorentz force, J·E, and γ _ev (J·J)/σ represents Joule heat, and J is the current density. The parameter γ _ev represents the Joule thermal vibration energy ratio.

σ is the electrical conductivity of the plasma. The surface material is assumed to be entirely insulated, and the Hall effect is disregarded for simplicity. According to the Ohm’s law, the current can be calculated by the following Equation (6) [Reference Nagata, Otsu and Yamada25].

(6) \begin{equation}J = \sigma (E + u \times B)\end{equation}

The chemical reaction source term ${\dot \omega _s}$ is given by Equation (8) [Reference Yang and Liu34]:

(7) \begin{equation}{\dot \omega _{\rm{s}}} = {{\mathop{{\textit {M}}}\nolimits} _{\rm{s}}}\sum\limits_{r = 1}^{nr} {\left( {\nu _{s,r}^{\prime \prime } - \nu _{s,r}^\prime } \right)} \left[ {{{\mathop{\rm k}\nolimits} _{f,r}}\prod\limits_{{\rm{i}} = 1}^{n{\rm{i}}} {{{\left( {\frac{{{\rho _{\rm{s}}}}}{{{{\mathop{{\textit {M}}}\nolimits} _{\rm{s}}}}}} \right)}^{\nu _{s,r}^\prime }}} - {{\mathop{\rm k}\nolimits} _{b,r}}\prod\limits_{{\rm{i}} = 1}^{n{\rm{i}}} {{{\left( {\frac{{{\rho _s}}}{{{{\mathop{{\textit {M}}}\nolimits} _{\rm{s}}}}}} \right)}^{\nu _{s,r}^{\prime \prime }}}} } \right]\end{equation}

where $\nu _{s,r}^{\prime \prime }$ and $\nu _{s,r}^\prime $ are the stoichiometric coefficients of component s as products and reactants in reaction r, respectively. M _s is the molecular weight of species s. The Park 93 species-7 [Reference Park35] component reaction model is used to simulate the chemical processes in the aforementioned framework. This model primarily consists of the dissociation reaction of N₂ and O₂, the exchange reaction that produces NO, and the ionisation reaction that generates NO⁺ and e⁻.

The forward reaction rate coefficient k _f,r is determined by Arrhenius formulas as [Reference Casseau, Palharini and Scanlon36]:

(8) \begin{equation}{{{k}}_{{{f,r}}}}{{ = }}{{{A}}_{{f}}} \times {{T}}_{{{c,f}}}^{{b}}\exp \bigg ({{ - }}\frac{{{{{T}}_{{a}}}}}{{{{{T}}_{{{c,f}}}}}}\bigg)\end{equation}

where A _f is a pre-exponential factor, β is the temperature exponent, and T _a is the temperature of activation derived from the activation energy. T _c,f is the controlling temperature of the forward reaction.

Park et al. [Reference Zhang, W37] chose a least-squares curve fit of the equilibrium constant K _eq to calculate k _b,r [Reference Yang and Liu34]. The backward reaction rate coefficient k _b,r can be calculated by the following formula:

(9) \begin{equation}{{\mathop{ k}\nolimits} _{b,r}} = \frac{{{k_f}}}{{{K_{eq}}}}\end{equation}

where T _a is the activation temperature, T _c,f = T _tr ^0.7 T _ve ^0.3, and K _eq is the equilibrium constant (dimensionless) calculated by a function of temperature based on the curve-fitted method.

(10) \begin{equation}K_{\mathrm{eq}} = \exp \left[ a_1 \left( \frac{T}{10000} \right) + a_2 + a_3 \ln \left( \frac{10000}{T} \right) + a_4 \left( \frac{10000}{T} \right) + a_5 \left( \frac{10000}{T} \right)^2 \right]\end{equation}

where a ₁, a ₂, a ₃, a ₄ and a ₅ are the curve fit coefficients.

Bush [Reference Nagata, Otsu and Yamada25] provided a semi-empirical electrical conductivity formula in the form of a temperature fit based on the reference electrical conductivity and the reference temperature:

(11) \begin{equation}\sigma = {\sigma _0}{\bigg(\frac{{{T}}}{{{{{T}}_0}}}\bigg)^n}\end{equation}

where T ₀ is the reference temperature, T is the flow field temperature, and σ ₀ is the reference electrical conductivity.

The electrical conductivity reference values used in this paper are consistent with those in Ref. (Reference Nagata, Otsu and Yamada25). The reference parameters are: σ ₀ = 250 S·m⁻¹, T ₀ = 25000 K, and n = 2.

The computational schematic for the thermochemistry nonequilibrium MHD is shown in Fig. 2. The transport mixing rule is set by Wilke [Reference Candler38], and the multispecies transport model is Fick [Reference Kim, Kang and Park39]. The simulation starts with flow field initialisation, using initial conditions and applied transport and thermophysical properties. Subsequently, the initial conditions (such as ρ, u, T, etc.) are assigned to each grid cell. In the main solver, real gas effects are accounted for by introducing chemical reaction source terms into the continuity equation. Meanwhile, the momentum, total energy and vibrational energy equations are augmented with the Lorentz force, Joule heating, and vibrational energy source terms, respectively.

Figure 2.

Computational schematic for thermochemical nonequilibrium MHD.

3.2. Magnetic dipole field

The magnetic dipole field represents the magnetic field generated by a closed loop, which can be compared to the energised solenoid model. The magnetic dipole field distribution can be derived from the magnetic field strength at the feature points [Reference Ding, Jiang and Dong23, Reference Nagata, Otsu and Yamada26]:

(12) \begin{equation}\overrightarrow B = - \frac{{{B_0}}}{{2{{(r/{r_0})}^3}{r^2}}}\left( {3\overrightarrow r (\overrightarrow r \cdot \overrightarrow m ) - {r^2}\overrightarrow m } \right)\end{equation}

where B ₀ represents the magnitude of magnetic induction strength at the characteristic point. r ₀ denotes the distance of the feature point from the centre of the magnetic dipole and is assumed to be 175 mm here. $\overrightarrow r $ denotes the position vector from a point in the flow field to the centre of the magnetic dipole. r is the distance of the point from the centre of the magnetic dipole and $\overrightarrow m $ represents the vector denoting the angle of the magnetic dipole.

Assuming the magnetic dipole is centred at (x ₀, y ₀, z ₀). The magnetic dipole angle ( $\theta $ ) is rotating counterclockwise around the Z-axis along the central position axis. The vectors $\overrightarrow m $ and $\overrightarrow r $ are expressed as:

(13) \begin{equation}\overrightarrow m = \left[ \begin{array}{l}\cos \theta \\{\rm{sin}}\,\theta \\0\end{array} \right]{\rm{ }}\overrightarrow r = \left[ \begin{array}{l}x - {x_0}\\y - {y_0}\\z - {z_0}\end{array} \right]\end{equation}

The magnetic field distribution of a magnetic dipole is centred at (x ₀, y ₀, z ₀) rotating around the Z-axis in the Cartesian coordinate system:

(14) \begin{equation}\left[ \begin{array}{l}{B_x}\\{B_y}\\{B_z}\end{array} \right] = - \frac{{{B_0}}}{{2{{(r/{r_0})}^3}{r^2}}}\left[ \begin{array}{c}2{(x - {x_0})^2}\cos \theta - {(y - {y_0})^2}\cos \theta - {(z - {z_0})^2}\cos \theta + 3(x - {x_0})(y - {y_0})\sin \theta \\[5pt]2{(y - {y_0})^2}\sin \theta - {(x - {x_0})^2}\sin \theta - {(z - {z_0})^2}\sin \theta + 3(x - {x_0})(y - {y_0})\cos \theta \\3(x - {x_0})(z - {z_0})\cos \theta + 3(y - {y_0})(z - {z_0})\sin \theta \end{array} \right]\end{equation}

The spatial magnetic field distribution of the model is shown in Fig. 3 [Reference Nagata, Otsu and Yamada26]. The magnetic field strength at point P is B ₀ and the distance from the magnetic dipole is r ₀. Point P is a feature point for the magnitude of the magnetic field strength.

Figure 3.

Schematic of dipole magnetic field distribution [Reference Nagata, Otsu and Yamada26].

Figure 4 illustrates the distribution of the magnetic field in the meridional plane for three angular configurations (θ = 0°, 45°, and 90°). Since other magnetic pole angles are similar to 45°, only three typical distributions of magnetic inductances are shown here. The electromagnetic coil is approximated as an idealised magnetic dipole aligned with the body axis for analytical clarity. The effects of MHD interactions are systematically analysed by varying θ (θ = 0°, 30°, 45°, 60° and 90°) between the vehicle axis and the applied dipole field orientation. The uniform magnetic field strength of B ₀ = 4 T is used for all simulations.

Figure 4.

Diagram of the magnetic field with different pole angles.

3.3 Boundary conditions

As practiced in Ref. (Reference Zhang, Wang and Zhang40), the boundary conditions for vehicles are established. At 293 seconds, the ELECTRE [Reference Muylaert, Walpot and Haeuser29] vehicle is flying at a velocity of 4230 m·s⁻¹ at an altitude of 53.3 km. The atmosphere temperature is 265 K, pressure is 52.8 Pa, viscosity is 1.68 × 10⁻⁵ (kg·(m/s)⁻¹) and density is 6.9 × 10⁻⁴ (kg/m³) [41]. The atmosphere composition consists of 22% O₂ and 78% N₂. The wall boundary conditions are applied without considering catalytic effects.

The computational domain is meshed based on the physical model presented in Section 2. A structured grid is employed to enhance convergence speed and improve computational accuracy. To precisely capture heat flux through the vehicle walls, the boundary layer grid is refined layer by layer. In accordance with the OpenFOAM solver framework, the axisymmetric computational domain is reformulated as a wedge-shaped structure. Figure 5 illustrates the grid structure of the fluid domain.

Figure 5.

ELECTRE vehicle boundary and grid model.

3.4 Code verification

The Hy2Foam solver, a density-based, compressible, non-equilibrium flow solver built on the OpenFOAM platform, was developed at Strathclyde University [Reference Casseau, Palharini and Scanlon36, Reference Zhang, Zhang and Yang42]. The magnetic dipole field is further integrated into the hy2Foam framework, forming the foundation for the subsequent research conducted in this study. The current study configures the solver with a maximum Courant number of 0.7.

Based on the RAM-C blunt cone model [Reference Ding, Jiang and Dong23], this section verifies the accuracy of the numerical model of thermochemical non-equilibrium MHD and the accuracy of the electrical conductivity model in this paper. The total length of the RAM-C blunt cone model is 1.295 m, the radius of the head is 0.1524 m and the angle between the cone and the central axis is 9°.

The flight altitude of the RAM-C blunt cone model is 71 km, and the speed is 7650 m/s. The specific boundary conditions are set in Table 1. The static pressure, static temperature and dynamic viscosity of the freestream are 4.737 Pa, 217.35 K and 1.4377 × 10⁻⁵ (kg·(m/s)⁻¹), respectively. The Reynolds number of the freestream is 12316. The wall is set slip with complete diffuse reflection. The momentum adjustment coefficient (α_M) is 1.0, the energy adjustment coefficient (α_T) is 1.0, and the wall temperature is 1500 K. The wall catalytic reaction is not considered. The magnetic field distribution is consistent with Ref. (Reference Nagata, Otsu and Yamada26), where the magnetic field is a dipole field. The centre of the magnetic dipole overlaps with the center of the head of the RAM-C blunt cone. The magnetic field strength at the feature point is B ₀ = 0.5 T, and the distance between the feature point and the magnetic dipole center is r ₀ = 0.1524 m.

Table 1.

RAM-C blunt cone boundary conditions

Figure 6 shows the validation of the trans-rotational temperature curve for the RAM-C vehicle along the stagnation point line. As shown in Fig. 6, the trans-rotational temperature and shock standoff distance calculated in this study agree well with the results reported in Ref. (Reference Ding, Jiang and Dong23).

Figure 6.

Case validation of MHD.

3.5 Grid independence

Following the grid generation methodology established in Ref. (Reference Zhang, Zhang and Yang42), the grid convergence study is conducted using three different grid numbers, as shown in Table 2. The terms n _x , n _y and y _wall refer to the number of grid nodes in the direction normal to the wall, the number of grid nodes along the wall, and the thickness of the first layer of the grid closest to the wall, respectively [Reference Zhang, Zhang and Yang42]. The thickness of the grid at the location near the wall decreases from Coarse to Fine, and the number of grids rises.

Table 2.

Parameters for different grid sizes

Wall heat flux predictions from three grid resolutions are compared in Fig. 7. The spatial distribution of the wall heat flux is characterised using polar coordinates, where the origin coincides with the geometric centre of the vehicle (point O in Fig. 2) and the polar axis aligns with the stagnation point line. The position of the wall with respect to the polar axis is defined by the angular coordinate α. Grid resolution demonstrates a significant influence on heat flux distribution characteristics. Sensitivity to grid parameters diminishes with grid refinement, both through increased cell count and reduced first-layer thickness. The Medium grid produces wall heat flux distributions that closely match those of the Fine. Despite a 70% reduction in grid cells, the spacing of the first layer increased by 50% (as shown in Table 2). Considering the balance between computational efficiency and numerical accuracy, the medium-resolution grid is adopted for all subsequent simulations.

Figure 7.

Verification of grid independence. (a) Wall heat flux, and (b) wall heat flux at the blunt position.

The experimental study of ELECTRE [Reference Muylaert, Walpot and Haeuser29] has been selected to validate the solver. A comparison between the calculated wall heat flux and the reference data is shown in Fig. 7. Overall, the computed wall heat flux agrees well with the experimentally measured results. Notably, at the stagnation point position, the predicted values show good agreement with the experimental values, as shown in Fig. 7(b). However, a discrepancy exists where the experimental values of wall heat flux surpass those of the prediction at the shoulder location. Such deviation may be because the wall catalysed is ignored [Reference Zhang, Zhang and Yang42]. Therefore, the validated code forms a reliable foundation for subsequent hypersonic MHD numerical simulations.

A domain independence study is conducted to ensure that the computational domain is sufficiently large to accurately capture the shock wave. The optimal computational domain size is determined by evaluating different domain sizes and comparing the resulting shock wave positions. Figure 8 presents the shock wave positions for various domain configurations.

Figure 8.

Mach number distribution under different fluid domains.

In Case 1, the grid configuration from Section 3.5 is used, with a magnetic field strength of 4 T and a magnetic pole angle of 0°. Case 2 expands the fluid domain by a factor of 1.6 in the axial direction and 1.35 in the radial direction compared to Case 1. Similarly, Case 3 expands the fluid domain by a factor of 2 in the axial direction and 1.5 in the radial direction. Apart from mesh resolution, Cases 2 and 3 maintain the same conditions as Case 1.

As shown in Fig. 8, variations in the fluid domain have an insignificant impact on the computed results. In particular, the predicted shock standoff distance remains nearly consistent. This confirms that the computational domain of Case 1 is sufficient for accurately capturing the shock wave.