Numerical Simulation of the Chemical Reaction on Faraday MHD Accelerator

Abstract

Magnetohydrodynamic (MHD) is one of the most promising novel propulsion technologies with the advantages of no pollution, high specific impulse, and high acceleration efficiency. As the carrier of this technology, the MHD accelerator has enormous potential for applications in hypersonic wind tunnels, supersonic ramjet engines, and deep space propulsion. In this study, a three-dimensional numerical simulation of an ideal Faraday magnetohydrodynamic (MHD) accelerator is conducted to assess the effect on performance with respect to applied potential and magnetic field intensity. The study is performed by employing a low magnetic Reynolds number MHD model coupled with a 7-component chemical reaction model to simplify the impact of real gas effects. The chemical reaction exhibits an increasing trend with rising applied potential and a decreasing trend with diminishing magnetic field strength. This variation influences the gas conductivity, subsequently affecting the velocity and thrust of the system. Specifically, at a magnetic field intensity of 2.0 T and an applied potential of 600 V, the accelerator exhibits maximum velocity and thrust growth rates of 18.6% and 59.8%, respectively.

1. Introduction

Magnetohydrodynamic (MHD) propulsion emerges as a highly promising and innovative technology, offering distinct advantages such as negligible environmental impact, elevated specific impulse, and exceptional acceleration efficiency. This novel approach holds significant potential across a range of applications, including hypersonic wind tunnels, trans-atmospheric vehicles, orbital transfer engines, and deep space exploration vehicles. However, the intricate coupling of electric fields, magnetic fields, and flow fields presents critical challenges, including an unclear understanding of the coupling mechanism, interference from induced magnetic fields, and prohibitive costs. Consequently, to establish itself as a viable carrier for this technology, it is imperative to undertake more extensive and systematic research on the MHD accelerator.

Research institutions in the United States, Russia, Japan, and China [1,2,3,4] have conducted comprehensive experimental investigations on the segmented Faraday MHD accelerator, establishing a rigorous experimental process to validate the practicality of MHD technology. However, the current devices are limited to measuring parameters at specific locations and generally have a short maintenance time. Consequently, relying solely on experiments proves insufficient to meet the demands of contemporary research. In tandem with scientific advancements, the increased power of computers has significantly enhanced the precision of numerical simulations, playing an increasingly crucial role in MHD accelerator research. Gaitonde et al. [5] performed a meticulous three-dimensional simulation of the flow characteristics of the magnetic fluid within the pipeline, developing a calculation program encompassing both induced magnetic fields and applied electric fields. Baranov et al. [6] analyzed the flow state of a water-cooled MHD accelerator, verifying the sustainability of hypersonic flow. Harada et al. [7,8] investigated the impact of basic parameters, such as electric field intensity, on the Faraday MHD accelerator through one-dimensional and three-dimensional numerical simulations. The performance of the linear MHD accelerator was also studied in their subsequent work [9], and it was identified that the maximum thrust is affected by the time interval of the current. Li et al. [10] concentrated on the influence of the magnetic field and analyzed the variation in electron density under magnets of different lengths. Sone et al. [11] studied the thrust performance of the diagonal-type MHD accelerator and identified that the relationship between magnetic flux density and maximum thrust increment is approximate to a steady-state model. Promson et al. [12] have conducted extensive research on the input current of the MHD accelerator and verified the importance of the discharge pulse width. Additionally, to explore the relevant properties of the MHD accelerator in a non-equilibrium ionization state, numerical simulations incorporating thermochemical models have also been undertaken [13].

In summary, while significant progress has been achieved in the numerical simulation of the MHD accelerator, there remains a notable gap in systematic research and analysis of certain factors. The majority of researchers have concentrated on aspects including channel structure, ionization mode, properties of channel materials, distribution of magnetic fields, and discharge modes of electrodes. Numerous significant findings have been obtained, such as the observation that the conductivity achieved through electron beam ionization closely resembles a Gaussian distribution, the impact of current spacing on accelerator performance, and the adverse effects of excessive magnetic field intensity. Notably, limited attention has been devoted to the study of real gas effects, with existing research predominantly concentrating on the analysis of flow field parameters. Similarly, systematic discussions on species distribution and electric field parameters are lacking. The present work aims to address this gap by investigating the ideal Faraday MHD accelerator under the influence of chemical reactions, considering varying applied magnetic fields and external electric fields. This investigation employs a three-dimensional numerical simulation utilizing the low magnetic Reynolds number magnetohydrodynamic model and a 7-component chemical reaction model. The results obtained are comprehensively analyzed to discern the variation laws of component particles, electric field parameters, and flow field parameters within the channel. This analysis aims to establish a foundation for the optimal design of the MHD accelerator.

2. Numerical Model

2.1. Magneto-Hydrodynamic Equations

Typically, the conductivity achievable by the MHD accelerator is relatively low. The flow within the channel is commonly assumed to exhibit characteristics of low magnetic Reynolds number flow, where $R_{\mathrm{em}}$ ≪ 1. In addition to the fundamental Navier–Stokes equations, which include the electromagnetic source term, conductive fluids must also adhere to Maxwell’s equations. The N–S equations, when coupled with electromagnetic source terms, are expressed as follows [14]:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathit{u})=0$$

(1)

$$\frac{\partial (\rho \mathit{u})}{\partial t}+\nabla \cdot (\rho \mathit{u}\mathit{u})+\nabla P-\nabla \cdot \mathit{\tau }=\mathit{J}\times \mathit{B}$$

(2)

$$\frac{\partial }{\partial t}\left[\rho (c_{v}T+\frac{1}{2}{\mathit{u}}^2)\right]+\nabla \cdot \left[\rho \left(c_{v}T+\frac{1}{2}{\mathit{u}}^2\right)\mathit{u}\right]=\nabla \cdot (\tau \mathit{u})-\nabla \cdot \mathit{q}+\mathit{J}\cdot \mathit{E}$$

(3)

where $\rho$ is the gas density, $\mathit{u}$ is the fluid velocity, $P$ is the static pressure, $\tau$ is the viscous stress, $\mathit{J}$ is the current density, $\mathit{B}$ is the magnetic flux density, $c_v$ is the constant volume specific heat, $T$ is the static temperature, $\mathit{q}$ is the conductive heat flux, and $\mathit{E}$ is the electric field.

The state equation governing the behavior of the working gas is presented as

$$P=\rho RT$$

(4)

Given the constancy of the applied magnetic field in our study, the Maxwell equations are simplified as

$$\nabla \cdot \mathit{J}=0, \nabla \times \mathit{E}=0$$

(5)

Generalized Ohm’s law:

$$\mathit{J}=\sigma (\mathit{E}+\mathit{u}\times \mathit{B})-\frac{\beta }{\left|\mathit{B}\right|}(\mathit{J}\times \mathit{B})$$

(6)

where $\sigma$ is the conductivity, $\beta$ is the electron Hall parameter. Building upon Equations (5) and (6), we can obtain $\mathit{E}=-\nabla \phi$, where $\phi$ is the electric field potential and the potential Poisson equation can be expressed as

$$-\nabla \cdot (\sigma \nabla \phi )=-\sigma \left[\nabla \cdot \left(\mathit{u}\times \mathit{B}\right)+\beta B\nabla \cdot \mathit{u}\right]$$

(7)

2.2. 7-Component Chemical Reaction Model

In the study of chemical reaction phenomena in high-temperature ionized air, the current research commonly employs 5-component, 7-component, and 11-component models. The 5-component model, though relatively simple, is inadequately adapted to represent the complexities of the real atmospheric environment. The 7-component model, encompassing reactions between main particles, proves more effective in simulating real gas effects. On the other hand, the 11-component model, while the most comprehensive, demands substantial computing resources. Consequently, this paper opts for the 7-component reaction model, which includes N₂, O₂, NO, N, O, NO⁺, and $e$ as the components. The potential 21 chemical reactions among these components are as detailed below [15]:

$$N_2+M\rightleftarrows N+N+M$$

(8)

$$O_2+M\rightleftarrows O+O+M$$

(9)

$$NO+M\rightleftarrows N+O+M$$

(10)

$$N_2+O\rightleftarrows NO+N$$

(11)

$$NO+O\rightleftarrows N+O_2$$

(12)

$$N+O\rightleftarrows N{O}^{+}+e$$

(13)

where $M$ represents a generic particle.

The reaction rate R is controlled by the forward reaction rate $k_f$ and backward reaction rate $k_b$. The control law for this relationship is expressed as follows [16]:

$$R_1={\displaystyle \sum \left\{-k_f\left[N_2\right]\left[M\right]+k_b\left[N\right]\left[N\right]\left[M\right]\right\}}$$

(14)

$$R_2={\displaystyle \sum \left\{-k_f\left[O_2\right]\left[M\right]+k_b\left[O\right]\left[O\right]\left[M\right]\right\}}$$

(15)

$$R_3={\displaystyle \sum \left\{-k_f\left[NO\right]\left[M\right]+k_b\left[N\right]\left[O\right]\left[M\right]\right\}}$$

(16)

$$R_4=-k_f\left[NO\right]\left[O\right]+k_b\left[N\right]\left[O_2\right]$$

(17)

$$R_5=-k_f\left[O\right]\left[N_2\right]+k_b\left[N\right]\left[NO\right]$$

(18)

$$R_6=-k_f\left[N\right]\left[O\right]+k_b\left[N{O}^{+}\right]\left[e\right]$$

(19)

where $\left[X\right]$ is the molar concentration. The expressions of $k_f$ and $k_b$ are as follows:

$$k_f=C_f\cdot T^{n_f}\cdot e^{-E_r/(RT)}$$

(20)

$$k_b=k_f/k_{eq}$$

(21)

$$k_{eq}=\exp \left(A_1+A_{2}z+A_3{z}^2+A_4{z}^3+A_5{z}^4\right)$$

(22)

where $C_f$ is the pre-exponential factor, $n_f$ is the temperature exponent, $E_r$ is the activation energy, $R$ is the gas constant, z = 10,000/T, and A₁~A₅ are the fitting parameters [17].

In the context of a chemical reaction involving N components, additional characteristic parameters that require consideration include the isobaric-specific heat, viscosity coefficient, and heat transfer coefficient. The isobaric-specific heat of the i-th component $C_{p,i}$ can be determined through the fitting formula of T [18]:

$$C_{p,i}=\frac{R}{M_i}{\displaystyle \sum _{k=1}^5{A}_{k,i}{T}^{k-1}}$$

(23)

where $R=8314.4 \mathrm{J}/\left(\mathrm{kmol}\cdot \mathrm{K}\right)$, $M_i$ is the molar mass of the i-th component, $A_{k,i}$ is a constant.

The isobaric-specific heat of the mixed gas is determined according to the mixing law:

$$C_p={\displaystyle \sum _i{Y}_i{C}_{p,i}}$$

(24)

where $Y_i$ is the mass fraction of the i-th component.

The viscosity coefficient $\xi$ and heat transfer coefficient $\varsigma$ adhere to the Wilke semi-empirical formula [19], and the corresponding governing equations are as follows:

$$\xi ={\displaystyle \sum _{i=1}^N\left(\frac{{\xi }_i{X}_i}{{\displaystyle \sum _{j=1}^N{\Phi }_{ij}{X}_j}}\right)},\varsigma ={\displaystyle \sum _{i=1}^N\left(\frac{{\varsigma }_i{X}_i}{{\displaystyle \sum _{j=1}^N{\Phi }_{ij}{X}_j}}\right)}$$

(25)

where $X_i$ is the mole fraction, ${\Phi }_{ij}$ is the partition function between components, and ${\xi }_i$ and ${\varsigma }_i$ are the viscosity coefficient and heat transfer coefficient of the i-th component, respectively, which can be expressed as

$${\Phi }_{ij}=\frac{{\left[1+{\left(\frac{{\xi }_i}{{\xi }_j}\right)}^{0.5}{\left(\frac{M_j}{M_i}\right)}^{0.25}\right]}^2}{{\left[8\left(1+\frac{M_i}{M_j}\right)\right]}^{0.5}}$$

(26)

$${\xi }_i=\left\{\begin{array}{ll}1.4584\times {10}^{-5}\times \frac{T^{1.5}}{T+110.3} ,& \frac{g}{cm\cdot s} (T<1000 K)\\ \left[\exp (C_{{\xi }_i})\right]T^{\left(A_{{\xi }_i}\ln T+B_{{\xi }_i}\right)} , & \frac{g}{cm\cdot s}\left(T\ge 1000 K\right)\end{array}\right.$$

(27)

$${\varsigma }_i=\left\{\begin{array}{ll}5.9776\times {10}^{-6}\times \frac{T^{1.5}}{T+194.4}& , \frac{cal}{cm\cdot s\cdot k}\left(T<1000 K\right)\\ \left[\exp (E_{{\varsigma }_i})\right]T^{\left[A_{{\varsigma }_i}{\left(\ln T\right)}^3+B_{{\varsigma }_i}{\left(\ln T\right)}^2+C_{{\varsigma }_i}(\ln T)+D_{{\varsigma }_i}\right]}& , \frac{cal}{cm\cdot s\cdot k}\left(T\ge 1000 K\right)\end{array}\right.$$

(28)

where $A_{{\xi }_i}$, $B_{{\xi }_i}$, $C_{{\xi }_i}$, $A_{{\varsigma }_i}$, $B_{{\varsigma }_i}$, $C_{{\varsigma }_i}$, $D_{{\varsigma }_i}$, and $E_{{\varsigma }_i}$ are constants [20].

In this paper, the chosen model for the calculation of conductivity is the multicomponent plasma conductivity model [21], and its expression is as follows:

$$\sigma =\frac{e^2{n}_e}{m_e{\displaystyle \sum _{s\ne e}{v}_{e,s}}}$$

(29)

where ${\nu }_{e,s}$ is the collision frequency of the electron and component s. When s is an ion, ${\nu }_{e,s}$ can be expressed as

$$v_{e,s}=6\pi {\left(\frac{e^2}{12\pi {\epsilon }_{0}k{T}_e}\right)}^2\ln \left[12\pi {(\frac{{\epsilon }_{0}k}{e^2})}^{1.5}\sqrt{\frac{T^3}{n_e}}\right]n_s\sqrt{\frac{8k{T}_e}{\pi m_e}}$$

(30)

where ${\epsilon }_0$ is the permittivity of the vacuum and k is the Boltzmann constant. When s is a neutral particle, ${\nu }_{e,s}$ can be written as

$${\nu }_{e,s}=\frac{4}{3}{\sigma }_{e,s}{n}_s\sqrt{\frac{8k{T}_e}{\pi m_e}}$$

(31)

where ${\sigma }_{e,s}$ is the collision cross-section of the neutral component and electron [18].

3. Simulation Procedure

3.1. Accelerator Configuration and Boundary Conditions

The channel structure utilized in this study is illustrated in Figure 1. The entrance size of the channel is $40\times 40$ mm, the exit size is $57\times 57$ mm, the total length is 500 mm, and the expansion angle of the channel is 1°. There are a total of 20 pairs of electrodes, with a width of 10 mm and a spacing of 5 mm between each pair. The initial pair of electrodes is positioned 25 mm away from the entrance.

To ensure a more accurate observation of changes in chemical composition closely resembling the actual situation, the inflow is assumed to be composed of N₂, O₂, NO⁺, and $e$. Furthermore, the particle number density of NO⁺ and $e$ is assumed to be equal. The inflow is considered electrically neutral on a macroscopic scale, as outlined in

Table 1. Initial conditions for species involved.

Species	Mass Fraction
O2	0.24
NO+	5.325 × 10−5
e	9.7625 × 10−10

. The boundary conditions are set as pressure inlet, pressure outlet, and non-slip isothermal walls, with the discretization method being a second-order upwind scheme. The simulation software is the commercial CFD package FLUENT, and the numerical calculation method is the finite volume method. The density-based solver, display solution format, and the Spalart–Allmaras turbulence model are used. In addition, user-defined functions are used to combine formulas and calculations. Specific calculation parameters are detailed in

Table 2. Calculation conditions.

Factor	Value
Static temperature, K	3240
Static pressure, Pa	3.1 × 105
Mach number	1.36
Velocity, m/s	1475
Wall temperature, K	1000

.

Table 3. Grid independence verification.

Mesh	Number of Grids	Vmax (m/s)	Relative Error
M1	$40\times 40\times 500$	1658.0	0.96‰
M2	$50\times 50\times 600$	1659.1	0.30‰
M3	$80\times 80\times 900$	1659.6	reference value

presents the maximum velocity and relative error at the channel exit using the parameters from previous work [22], considering both magnetic and electric fields at 0. The three types of grid numbers employed were refined at the boundary layer and channel entrance. After considering factors such as computational cost and accuracy, this paper exclusively utilizes the M2 grid for the calculations.

3.2. Computational Model Validation

Hartmann flow is a fundamental flow with an analytical solution in magnetohydrodynamics, providing a robust representation of the influence of the Lorentz force on flow. Consequently, the Hartmann flow was employed for validation purposes to assess the computational capability of the model for conductive fluid [14]. The numerical solution demonstrated a high degree of consistency with the analytical solution, thereby confirming the effectiveness of the model, as shown in Figure 2. Furthermore, program verification was also conducted with experiments in previous work [22]. The maximum error in speed is 7.87%, potentially attributed to the simplification of conductivity and neglecting the influence of the Hall effect. However, it effectively captures the velocity trend, affirming the feasibility of this algorithm.

Moreover, the chemical reaction model employed in this study incorporates 7 components and a total of 21 reaction formulas, constituting a complex mechanism that necessitates special verification. Figure 3 illustrates the mass fractions of N₂, O₂, NO, and O along the stagnation line. The accuracy of the chemical reaction model is verified with the strong agreement between the computed results and those reported in Ref. [23], indicating that the user-defined functions were implemented accurately and the parameters were appropriately configured.

4. Results and Discussion

4.1. Particles Distribution

The mole fraction of N₂ along the center axis under different magnetic field intensities and external potentials is illustrated in Figure 4. When the applied potential is between 200 V and 600 V, the mole fraction of N₂ gradually increases with the increase in magnetic field intensity but remains smaller than in the case of no electromagnetic effect (base). This is attributed to the influence of energy injection, where the rising temperature and intensification of chemical reactions lead to increased N₂ decomposition, subsequently reducing its mole fraction. When the applied potential is constant, the induced electromotive force $U=BLV$ grows with the increment of the magnetic field intensity. Correspondingly, the energy injection is inhibited, resulting in decreased N₂ decomposition and an increase in mole fraction. It is evident that the higher the applied potential, the higher the decomposition rate of N₂, and the applied potential plays a more significant role than the magnetic field in terms of the rate of the chemical reaction. Additionally, Figure 4e reveals that when U = 600 V and B = 0.5 T, the mole fraction of N₂ increases in the latter half of the channel. This is because the reaction is reversible; once the amount of N₂ is too low, the reverse reaction rate of N₂ becomes higher than the forward reaction rate, leading to an increase in the content of N₂.

The variation of the O₂ mole fraction is very similar to that of N₂, positively related to magnetic field intensity and negatively correlated with applied potential, as shown in Figure 5. However, since the dissociation temperature of O₂ is lower, the dissociation degree of O₂ is much higher than that of N₂ under the same conditions. With the magnetic field at 0.5 T and the potential at 200 V, the dissociation degrees of O₂ and N₂ at the exit of the channel are 4.25% and 0.830%, respectively. With the magnetic field at 0.5 T and the potential at 600 V, the dissociation degrees of O₂ and N₂ are 62.5% and 12.2%, respectively. Obviously, the difference in dissociation degree between the two also increases with the increase of the applied potential.

Besides N₂ and O₂, the other two particles are NO⁺ and $e$, which have a significant influence on the flow state inside the MHD accelerator. They share the same particle number density and variation trend. The mole fraction of $e$ along the center axis under different magnetic field intensities and external potentials is shown in Figure 6. It is indicated in Table 1 that the initial value of N and O is zero. Therefore, with the injection of energy, NO⁺ and e produce a large amount of N and O, resulting in a decrease in mole fraction similar to that of O₂ and N₂.

Figure 7 illustrates the distribution of the mole fraction of $e$ on the XOY plane at B = 1.5 T and U = 300 V. In the mainstream region, the mole fraction of $e$ decreases along the radial direction of the channel, which is consistent with the change along the center line. The Joule heating effect is not obvious near the wall when the potential is small, resulting in low temperature and weak reaction. On that account, the mole fraction of $e$ is larger than the mainstream region. As for N, O, and NO, since they play a secondary role in the flow, they will not be discussed in this chapter due to the limited space.

4.2. Electic Field Parameters

For the MHD accelerator with a chemical reaction model, the primary electric field parameters that require attention include conductivity, current, Lorentz force, and electric efficiency.

Figure 8 illustrates the distribution of conductivity along the central line of the channel. According to Equation (29), conductivity is directly proportional to the electron number density and inversely proportional to the particle collision frequency. In the absence of an external electromagnetic field, the supersonic flow accelerates in the dilated channel, causing a decrease in temperature and subsequently reducing particle collision frequency. Therefore, the channel’s conductivity increases despite the declining number of electrons.

Upon the application of an electromagnetic field, as observed in Figure 7 and Figure 8, the electron number density further decreases, and the particle collision frequency rises due to the injected energy. Consequently, the gas’s electrical conductivity remains smaller than that without electromagnetic action. When the applied potential is constant, higher magnetic field intensity leads to larger electrical conductivity, revealing the opposite tendency when the magnetic field intensity is constant. In the second half of the channel (X > 0.31 m), conductivity exhibits a significant upward trend. This is because, after the gas exits the electromagnetic action region, the influence of Joule heat continues to diminish with the drop in particle collision frequency, causing a slight increase in electrical conductivity.

Figure 9 displays the conductivity distribution in the XOY plane when B = 1.5 T and U = 300 V. In the core flow area, conductivity changes smoothly and symmetrically, consistent with the trend in Figure 8. At the same value of X, there is higher conductivity near the wall in the electromagnetic region due to low temperature, while the trend in the second half of the channel is the opposite.

Figure 10 depicts the distribution of current lines in the XOY plane under different magnetic field intensities and applied potentials. With the exception of Figure 10c, the other plots exhibit a similar trend. The current overflows at the front and back of the electromagnetic action region due to the conductivity of the working fluid throughout the channel, correctly flowing from the cathode to the anode. The overflow currents ultimately concentrate on the first and last pairs of electrodes, following the conductive wall condition of the cathode and anode and the insulation wall condition of the rest of the walls.

In Figure 10c, the current correctly flows from the cathode to the anode in the first half of the channel; however, it reverses direction after X = 0.25 m and forms a vortex current near X = 0.31 m. Several factors contribute to this situation, including the higher applied potential, continuous gas acceleration, and increased induced electromotive force along the flow direction in the first half of the channel. When the gas reaches X = 0.25 m, the induced electromotive force surpasses the applied potential, resulting in the generation of reverse current. Subsequently, the overflow current flows from the cathode to the anode, concentrating on the last pair of electrodes at the back end of the channel due to the disappearance of magnetic field intensity and induced electromotive force as the gas reaches X = 0.31 m. Consequently, currents with two directions occur at the last pair of electrodes, ultimately forming the vortex current.

At the same applied potential, a higher magnetic field intensity enhances the induced electromotive force, resulting in a decrease in current and even a change of direction, as depicted in Figure 11. The overflow current at the front end of the channel (X < 0.015 m) and the back end of the channel (X > 0.31 m) is independent of the magnetic field intensity; however, it increases with the rise in applied potential. The overflow current is more significant near the electrodes, with the maximum value occurring at the front end of the channel due to the smaller distance between the first pair of electrodes in this work compared to that of the last pair of electrodes.

The Lorentz force in the X direction is not only influenced by the current density $J_y$ but is also closely related to the magnetic field strength, as illustrated in Figure 12. The direction of $F_x$ is solely determined by $J_y$ since the magnetic fields used in this work are all uniform magnetic fields pointing in the positive direction of the Z-axis. When the current flows from the cathode to the anode ($J_{\mathrm{y}}>0$), the Lorentz force points in the positive direction of the X-axis, representing positive work. Examining Figure 12, it is evident that when the applied potential is small, $J_y$ decreases and may even assume negative values with the increase in the magnetic field strength, resulting in the intersection of $F_x$ in the diagram. Conversely, when the applied electric potential is large, the basic value of current density is strong and $F_x$ exhibits a certain regularity of increase with the rise in magnetic field strength.

Figure 13 illustrates the electrical efficiency under different applied potentials and magnetic field strengths. When the magnetic field strength is held constant, the channel’s electrical efficiency gradually decreases with the increase of the applied potential, as depicted in Figure 13a. The downward trend is nonlinear, and the decreasing rate diminishes as the potential rises. Notably, there is an observable phenomenon of energy extraction ($Jy<0$) at B = 2.0 T and U = 200 V. This explains why the electrical efficiency is lower than that at B = 2.0 T and U = 300 V. For applied potentials exceeding 300 V, the electrical efficiency exhibits nearly linear growth with the increase in the magnetic field intensity, as shown in Figure 13b. The variation trend of electrical efficiency is closely linked to conductivity. As the magnetic field intensity increases, the conductivity and the energy utilization rate of the channel also increase, leading to the growth of electrical efficiency. Additionally, the linear growth rate decreases with the increase of the applied potential.

4.3. Flow Field Parameters

The static temperature along the center line under different conditions is depicted in Figure 14. When the magnetic field intensity is constant, Joule dissipation becomes more intense with the increase of the applied potential, leading to an increment in channel static temperature. Conversely, when the applied potential is constant, the current density decreases with the increase in magnetic field strength, suppressing the Joule heating effect and causing a decrease in static temperature. In conjunction with Figure 8, it is evident that the Joule heating effect is weak when the applied potential is small, resulting in increased gas conductivity along the flow direction. Consequently, a significant portion of the gas’s internal energy is converted into kinetic energy, leading to a monotonous decrease in static temperature along the X direction after passing through the first pair of electrodes, as shown in Figure 14c. In contrast, the Joule heating effect is more pronounced at a larger applied potential, causing the static temperature to continuously increase in the electromagnetic action region, as shown in Figure 14d.

Figure 15 illustrates the distribution of static pressure along the center line of the channel. Following the application of the electromagnetic field, there is a significant increase in static pressure near the first pair of electrodes, indicating an MHD compression effect. When the magnetic field intensity is held constant, the static pressure rises with an increase in the applied potential. Conversely, when the applied potential is constant, an increase in magnetic field intensity helps alleviate the impact of the Joule heating effect, resulting in a decrease in static pressure. Additionally, in cases of low applied potential and high magnetic field strength, the static pressure curve closely resembles that without electromagnetic effect, as depicted in Figure 15c,d.

Figure 16 depicts the distribution of velocity along the center line of the channel. Under the influence of the MHD compression effect, the velocity decreases near the first pair of electrodes. This phenomenon becomes more pronounced with an increase in applied potential and a decrease in magnetic field intensity. Consequently, the gas must traverse a certain channel length, denoted as D, before achieving acceleration comparable to that without electromagnetic effect. When the magnetic field intensity is small, the value of D may exceed the entire electromagnetic region, as shown in Figure 16a. As the magnetic field intensity increases, the channel length gradually shortens, as demonstrated in Figure 16b,c, with maximum values of approximately 0.18 m and 0.1 m, respectively. When the magnetic field intensity is constant, an increase in applied potential also augments the value of D. However, the impact of external potential is considerably smaller than that of the magnetic field, particularly in cases of low magnetic field intensity. Combined with Figure 12, it becomes evident that when the applied potential is small, the Lorentz force performs negative work in the electromagnetic region with an increase in magnetic field intensity, leading to a decrease in velocity. Hence, the maximum exit velocity is observed at B = 1.0 T, while the velocity at B = 1.5 T falls between that at B = 1.0 T and B = 2.0 T in Figure 16d. With the magnetic field set at 2.0 T and the potential at 200 V, the exit velocity is 2081 m/s, representing only a 0.873% increase compared to that without electromagnetic effect. With the magnetic field at 2.0 T and the potential at 600 V, the exit velocity is 2448 m/s, which is 18.6% higher than that without electromagnetic effect.

The channel thrust under different applied potentials and magnetic field intensities is depicted in Figure 17. When the magnetic field intensity is constant, the velocity and pressure at the outlet increase with the rise of applied potential, resulting in enhanced thrust, as illustrated in Figure 17a. However, with constant applied potential, the channel thrust may decrease as the magnetic field strength increases, as shown in Figure 17b. This phenomenon is attributed to the strong magnetic field conditions, where the Lorentz force may perform negative work in the channel, leading to a reduction in electromagnetic power and velocity. Concurrently, the static pressure also decreases, contributing to a decline in thrust.

With the potential set at 200 V and magnetic fields at 0.5, 1.0, 1.5, and 2.0 T, the channel electromagnetic power is 0.0623, 0.0918, 0.0853, and 0.0445 MW/m³, respectively. The thrust values are 654, 659, 651, and 631 N, with thrust increasing rates of 7.74%, 8.57%, 7.25%, and 3.95% compared to that without electromagnetic action. When the potential is increased to 600 V and magnetic fields are at 0.5, 1.0, 1.5, and 2.0 T, the channel electromagnetic power becomes 0.185, 0.365, 0.525, and 0.651 MW/m³, respectively. The thrust values are 835, 890, 937, and 970 N, with thrust increasing rates of 37.6%, 46.6%, 54.4%, and 59.8%, respectively, compared to that without electromagnetic action.

5. Conclusions

In this work, a three-dimensional numerical simulation of an ideal Faraday MHD accelerator is conducted using the low magnetic Reynolds number MHD model and the 7-component chemical reaction model. The investigation analyzes the impact of magnetic field intensity and applied potential on channel particle distribution, electric field parameters, and flow field parameters, leading to the following conclusions:

(1)

With an increase in the applied potential, the mole fraction of N₂, O₂, NO⁺, and e decreases, and the collision frequency between particles increases, resulting in a reduction in electrical conductivity. Conversely, with an increment in magnetic field intensity, the magnetic field’s inhibitory effect on energy input strengthens, leading to a decrease in the degree of chemical reaction and an increase in electrical conductivity.

(2)

Since the working fluid exhibits conductivity throughout the channel, there is a current overflow at the front and back ends of the electromagnetic region, concentrating on the first pair of electrodes and the last pair of electrodes, respectively.

(3)

Under continuous acceleration, the electrical efficiency increases almost linearly with the rise of magnetic field intensity. However, the linear growth rate decreases with an increase in external potential.

(4)

When the magnetic field intensity remains constant, the velocity and pressure at the outlet increase with an increase in applied potential, resulting in thrust growth. However, with constant applied potential, the channel thrust decreases as the magnetic field strength increases, showing a monotonic increase only when U > 300 V.

In summary, the research findings offer valuable insights into the particle distribution patterns within MHD accelerators. However, the current research still encounters several challenges, including the impact of the Hall effect, temperature effects on wall surfaces, and voltage distribution between the anode and cathode. Building upon the existing foundation, these issues will be systematically addressed in future research endeavors, ultimately facilitating the advancement and application of MHD accelerators.