1 Introduction

Upon to the proposal of the magnetohydrodynamic (MHD) power generator concept, it has attracted much attention and has been considered by many researchers as the most realistic and promising approach for thermoelectric direct conversion power generation. In contrast to conventional power generators, MHD generators may produce significant power output without mechanical moving parts and be integratable with coal-fired steam cycle power generation as the first stage of the steam cycle, which has been demonstrated that it has high efficiency and power density and low environmental pollution^[1-2].

The MHD generator is categorized as two types of structures^[3-4], disc and Faraday generation channels. Based on the shock wave tunnel equipment, the relationship among the total temperature, the generator performance, and the ionization characteristics of inert gas for the Faraday generator was investigated experimentally^{[1-2, 5-8]} with a high temperature inert gas as the working fluid. The experimental results showed that the plasma can be transferred from the non-uniform and unstable state to the uniform and steady state. The output power of the generator can be improved by increasing the inlet total temperature which varies for different inert gases. The inlet total temperature required for helium is higher than that of argon and that of xenon is the lowest. By carrying out numerical simulations on the experimental apparatus, xenon is proved to be the smallest ionization fluctuation, and helium has the highest enthalpy extraction rate at low temperature. Under the high temperature condition, the collision between ions and electrons becomes the main collision and the dependence of electrical conductivity on the electron number density declines. The MHD generator performance was numerically simulated for Ar/Cs mixture gas as the working fluid based on the experiment of shock tunnel equipment^[9]. The results showed that the current density near the electrode is larger than that of the mainstream. Plasma pattern with dark and bright alternating stripes is formed in the channel. The electron temperature ranges from 4 000 K to 6 000 K in the bright stripe and is about 2 500 K in the dark stripe. The electronic density is about 10¹⁹, which agreed with that of the experimentally observed phenomena. In cooperation with the American Marshall Space Flight Center^[10-11], Nagaoka University of Technology investigated the MHD generator system performance for a He/Xe mixture gas as a working medium under the total temperature of 1 800 K heated by nuclear reactors for the space energy system. They predicted that the mass power ratio can be reduced to 3 kg/kW or even below 2 kg/kW if the net power generation reached 1 MW or exceeded 3 MW, respectively. Harada et al.^[12-13] investigated the feasibility of replacing the alkali metal closed-loop with Xe seed for an MHD power generation system by two-dimensional (2D) numerical simulation. They considered that the performance of the generator with He/Xe mixture gas is far higher than that of He/Ce under the same boundary conditions. Moreover, the MHD generator performance with He/Xe mixture gas can be greatly improved and achieve 40% enthalpy extraction rate by improving the inlet electron temperature using pre-ionization and load factor. The effect of three-body recombination rate was studied on ionization stability for MHD generators^[14-17]. When the general Hivnov and Hirschberg curves of three-body recombination rate are applied, the discharge characteristics of the disc MHD generation channel are of a vortex form due to the instability of the ionization. At higher electron temperature, a steady enthalpy extraction rate can be obtained using the Biberman three-body recombination rate curve, where the discharge pattern and ionization process are fairly stable. The instability of the ionization can be suppressed by reducing the seed concentration or increasing the load factor and inlet rotation ratio. The dependence of different channel sizes on ionization degree of inlet and boundary layers, which influenced the disc power generator performance, was 2D numerically investigated by Veefkind et al.^[18] with pure He as the working fluid. The enthalpy extraction rate and entropy efficiency achieved 22.7% and 54.8%, respectively, by using a pre-ionization method to improve the inlet gas ionization degree for disc power generation. The pre-ionization degree of inlet is 4.79×10^-5 which costs 2% of the heat input power, while the load is 3 Ω under 4 T magnetic fields. When the ionization degree near the boundary layer is larger than that of the main flow, the full development of the boundary layer leads to the decrease of the generator performance due to the enhancement of the Lorenz force. While the ionization degree of the boundary layer is smaller than that of the main flow, the performance of the generator is improved. The relationship between plasma properties and MHD generator performance was performed experimentally by Veefkind^[19] and Litchford et al.^[20]. The instability of ionization leads to inhomogeneous distribution of conductivity, which makes the measured conductivity and Hall parameters deviate from the theoretical value. The plasma becomes quite unstable when the Hall parameter is below the critical value.

So far as we know, most of the studies on the MHD generator performance were accomplished by 2D numerical simulation. However, since the MHD effect is of temporal-spatial behaviors, some physical phenomena cannot be revealed by the 2D model. Therefore, for more comprehensive understanding, this paper investigates the performance of the power generation channel by three-dimensional (3D) numerical simulation for a Faraday channel with the helium/xenon mixture gas as the working fluid under different load factors and magnetic fields.

2 Physical model and numerical schemes 2.1 Physical and mathematical model

A Faraday channel model of the MHD power generator is shown in Fig. 1.

The channel length L is 0.12 m, the width a is 0.04 m, the height h is 0.03 m, the electrode width L₁ is 0.01 m, the width of the insulating wall L₂ between two electrodes is 0.005 m, and the channel involves 7 pairs of electrodes.

In the channel of MHD power generator, a fluid velocity which is much smaller than the speed of light results in a negligible displacement current term in the Maxwell equations. Due to the small magnetic Reynolds number for the weakly ionized gas, the effect of induced magnetic field can be neglected. Furthermore, for the applied static magnetic field, the Maxwell equations can be simplified as

(1)

(2)

(3)

(4)

The general Ohm's law for moving conducting fluids is

(5)

where B is the magnetic induction, H is the magnetic field, J is the current density vector, E is the electric field vector, U and σ are the velocity vector and electrical conductivity of ionized gas, respectively, and β is the Hall parameter.

As ∇ ×E=0, we have , and , in which φ is the electrical potential.

Derived from (1) and (5), the potential equation is written as follows:

(6)

The Hall effect can be neglected with β setting at 0 for the Faraday type channel with the segmented electrode which suppresses the generation of Hall current and Hall effect. Then, (6) can be rewritten as

(7)

After solving the current density vector J, the Lorentz force F, the electrical power density, and the other parameters can be obtained. Then, the source terms of the momentum equations and energy equation are obtained. For a supersonic viscous gas flow, the 3D compressible viscous MHD equations for electrical energy generation are as follows.

The continuity equation^[21] is

(8)

The momentum equations are

(9)

The energy equation is

(10)

where p is the static pressure, is the total energy of the mixed gas, e_s is the internal energy of gas, and (μ+μ_t)Q is the viscous dissipation term. μ is the fluid dynamic viscosity, and μ_t is the turbulent viscosity determined by which is calculated from the turbulent model of (11).

The Spalart-Allmaras (S-A) turbulent model, widely used in supersonic flows, is used to calculate the turbulence,

(11)

where is the turbulent kinematic viscosity, μ is the molecular viscosity, G_v is the product of the turbulent viscosity, Y_v is the dissipative term of the turbulent viscosity, and and C_b2 are constant.

2.2 Nonequilibrium ionization model

A dual temperature model is established to describe the nonequilibrium ionization process. The plasma is assumed to consist of two fluids, i.e., the electron flow and the heavy particle flow. The dual temperature model is implemented with the heavy particles at the gas temperature T_g and the electrons at another temperature T_e. The distribution of free electrons and binding electrons in the excited state is of thermodynamic equilibrium at the temperature T_e. Thus, the nonequilibrium ionization model of He/Xe mixture gas is written as follows:

Ion continuous equation^{[5-8, 11]}:

(12)

(13)

where n_i⁺ is the ionic concentration of He⁺ or Xe⁺, and n_i and n_e are number densities of atomic of He/Xe and electron, respectively. The ionization rate k_fi is

(14)

where g_i is the ground-state statistical weight of He⁺ or Xe⁺, g_h is the atomic ground-state statistical weight of He or Xe, m_e is the electron mass, h is the Planck constant, ε_i is the ionization potential of He or Xe atom, and k is the Boltzmann constant.

Using the Biberman model, the three-body recombination rate k_r is^[11]

(15)

where k_rh=1.09×10^-20T_e^-9/2, and .

The electron energy equations^{[5-8, 11]} are

(16)

(17)

where U_e is the electron energy, and v_en is the collision frequency among electrons, atoms, and ions.

(18)

where n_n is the quantity density of atoms and ions, Q_e-n is the cross-sectional area of the collision of electrons, atoms, and ions, and C_e is the average electronic velocity with the Maxwell distribution.

(19)

(20)

Here, e is the magnitude of charge by an electron, and ε₀ is the dielectric constant in vacuum. The electrical conductivity is expressed as .

2.3 Numerical schemes and boundary conditions

A coupled explicit scheme is used for the pressure-velocity solution. The finite difference method is adopted for the spatial discretization of the control equations. The diffusion term is discretized by the central difference scheme, and the convection term is discretized by the second-order upwind scheme and solved by SIMPLEC algorithm. The electrical potential equation is solved by the user-defined scalar (UDS) as transport equation solver. The inlet velocity of the channel is Ma=2.47. The total pressure is 0.6 MPa. The total inlet temperature is 2 000 K. The electron temperature is 8 000 K. The Xe seed concentration is 3×10^-5. The outlet employs the pressure outlet boundary condition. The bottom and top walls of the channel, no-slip and ideally insulated, have a constant temperature value at 600 K. The current density is suitable for the non-penetrating condition, i.e., J·n=0, where n is the unit normal vector of the wall. The electrode is a good conductor with zero voltage drop, and the amount of the electrical potential is given by the load factor.

3 Generator channel performance parameters

For measuring the performance of MHD generators, some parameters are characterized,

(21)

The power generation P₀ characterizes the amount of electric power generated per unit time of the generator channel (or the energy extracted from gas per unit time).

The electrical efficiency η, which indicates how much internal energy is converted into electrical energy in the channel, is an important parameter for energy conversion,

(22)

where P_F is the electromagnetic braking power, and P_g is the electric power generated by the channel.

The enthalpy effectiveness η_g is an important indictor to measure the conversion of fluid energy into electrical energy,

(23)

where W₀=m₀ (c_pT₁+u₀²/2) indicates the energy of the gas entering the MHD generator channel per unit time, in which m₀ is the inlet mass flow rate, c_p is the isothermal specific heat capacity, T₁ is the inlet static temperature, and u₀ is the inlet velocity.

4 Simulation results and discussion 4.1 Computational code verification

A non-uniform mesh of 120×40×80 is adopted in the computational domain, and all numerical results are checked for the grid independency. A finer mesh gradation is utilized in the regions near all walls.

In order to accurately simulate the MHD power generator performance with helium/xenon mixture gas as the working fluid, the ionization model of the gas and computational code must be verified. Therefore, the results of the ionization model are compared with the numerical results of Murakami and Okuno^[21] for a disc MHD generator.

From Fig. 2(a), the maximum deviation of the ionization rate occurs at 97 mm of radius where the result is 0.13% in the literature but 0.129% in our result. The maximum deviation of electron temperature appears at r=58 mm, where the literature and our simulation results are 5 500 K and 5 310 K, respectively. It is worth mentioning that the deviation is only 3.5%.

Figure 3 shows the distribution of the current lines at the y=0 cross section. The current line flows from the anode to the cathode as the power source. The current line does not pass through the insulating wall and does not form a loop in the channel cross section, which agrees with the analytical results for the insulating wall.

4.2 Effect of load factors on performance of MHD generator

Figure 4(a) presents the center line (x, 0, 0) velocity distribution along the channel for different load factors. The tendency indicates a slower and slower declination of the channel velocity with the increasing load factor. Figure 4(b) displays the comparison of Lorenz force F_x distributions along the channel center line with that near the channel wall while k=0.5. The Lorentz force F_x close to the channel wall behaves a periodic distribution along the flow direction because current concentrates at the edge of electrodes and flows primarily from the cathode toward the anodes, although a relatively small amount of current flows into insulators due to so-called numerical leakage. The current density at the electrode region is high, and there is no current passing through the insulating wall. Therefore, the Lorentz force F_x which is proportional to the current density is high at the electrode region and low at the insulating wall. Far from the walls, the current line becomes uniform (see Fig. 3). Then, the F_x distribution along the center line has hardly any fluctuation. The average of F_x near the wall decreases slightly first and then gradually increases along the channel wall to be the same as the distribution along the center line. The magnitude of Lorenz force F_x near-wall surfaces at the electrode region, for example at x=0.04 m, is larger than that in the main flow. The behavior is just on the contrary in the immediate vicinity insulating walls. Therefore, the contours of the Mach number in the channel show a wave state (see Fig. 5).

Figure 5 shows that the Mach number in the y=0 cross section decreases as the load factor increases. Because of the addition of load factor, the decrement of ionization rate in the channel makes the electromagnetic braking power decrease and thus the ratio of kinetic energy conversion to electrical power gradually decreases. From Fig. 5, the Mach number distribution at the y=0 plane obviously has a striking wave phenomenon due to the distorted current line shown at Fig. 3. At the core of the channel, the Lorentz force vectors point opposite to the flow direction rather regularly similar to that of the 2D case, altering their directions near the boundary layers, especially on insulators. The vertical component of Lorentz force accelerates the secondary flow which is responsible for the distortion of the current lines. As a result, the vertical component of the Lorentz forces will be strengthened locally due to the interaction between the flow and electromagnetic fields. The phenomenon becomes more notable by increasing the magnetic fields intensity. It may even show that the current lines loop at the center of the channel.

The variety of ionization degrees corresponds with the change of electron temperature. When k > 0.6, with the increase in the load factor, the ionization rate along the center line falls down faster and faster (see Fig. 6). As the load factor increases, the current density decreases, and the channel temperature rises a little bit, as shown in Fig. 7(a). Then, the three-body recombination rate is greater than the ionization rate, which makes the plasma keep combination, and the ionization rate decreases. When the load factor k≤0.6, the ionization ratio along the center line remains unchanged, i.e., the ionization rate equals that of the three-body recombination, and the plasma structure in the main flow is relatively stable.

Figure 7 shows the distributions of static temperature and static pressure along the center line for different load factors. As a result of the power generated by the Lorenz force, part of the energy is converted into the gas internal energy by Joule heating, which increases the static temperature and the static pressure. The increment degree of static temperature and static pressure is higher for smaller load factors with a lower resistance of the closed circuit. Consequently, the current density through the MHD generator channel decreases and induces the reduction in the electromagnetic braking power which reduces the kinetic energy conversion rate, which results in the increment degree of static temperature and static pressure reducing with the load factor increasing.

From Fig. 8, the temperature of the electrons near the electrode wall is higher than that of the main flow and the temperature near the insulating wall is the lowest one. According to the numerical simulation results of Tanaka and Okuno^[6], a higher current density induces a higher electron temperature and a greater collision frequency between electrons and neutral particles. The electrical conductivity is inversely proportional to the collision frequency, which results in the opposite change for the conductivity and electron temperatures.

As shown in Fig. 9, the volume average electrical conductivity in the channel continuously increases and the enthalpy extraction rate increases first and then decreases with the addition of the load factor. When the load factor increases, the static pressure rises at a decreasing speed and results in a decrement in the increasing amplitude of the neutral particles. Therefore, the larger the load factor is, the lower the collision frequency between the electrons and the neutral particles is. The opposite of electrical conductivity and the collision frequency results in the linear increase in the electrical conductivity with the load factor, as shown in Fig. 9. When the load factor k=0.625, the maximum of enthalpy extraction rate is about 4.37%. This optimum load factor k=0.625 is not the same as the best load factor k=0.5 for the liquid metal power generation channel. The liquid metal has a constant electrical conductivity, while the electrical conductivity of the plasma increases with the increase of the load factor. Therefore, it is helpful to improve the conductivity and enthalpy extraction rate. The power generator performance with the load factor is increased properly.

Table 1 lists the performance parameters of the MHD generator under different load factors. It indicates that the power output and enthalpy extraction increase first and then decrease with the increase in the load factors. The Joule heating and the electromagnetic braking power decrease with the increase in the load factor, which makes the electrical conductivity distribution gradually uniform. Therefore, the Joule heating dissipation is reduced which gives rise to improvement on the performance of the MHD generator.

4.3 Effect of magnetic fields on performance of MHD generator

Figure 10 shows the distribution of electrical conductivity along the center line under different magnetic fields. Obviously, when x < 0.02 m, the electrical conductivity decreases along the flow direction. When x > 0.02 m, the electrical conductivity increases with the magnetic field intensity while B≥2 T, but still keeps a downward tendency at B=1 T. When x < 0.02 m, the channel static temperature is low, and the ionization rate is less than the three-body recombination rate which induces the conductivity decreasing. With the channel static temperature increasing, the ionization rate increases, the three-body recombination rate reduces so that the electrical conductivity gradually increases along the flow direction. At B=1 T, the static temperature increment magnitude is small, and the electrical conductivity keeps declining as the ionization rate is still less than that of the three-body recombination.

The distribution in Fig. 11 indicates that the ionization rate decreases more and more slowly along the center line as the magnetic field B≥2 T but remains essentially unchanged along the center line. When B < 2 T, the ionization rate decreases faster because the magnetic field decreases. The current density decreases, and the increase degree of the static temperature in the channel is low. At this time, the three-body recombination rate is greater than the ionization rate, and the plasma behaves combination.

With the increase in the magnetic field, the average conductivity in the channel increases first and then decreases. When B=1.75 T, the conductivity reaches the maximum, about 11.05 S/m. The enthalpy extraction rate linearly increases with the magnetic field. When the magnetic field increases, the enthalpy extraction rate is remarkably improved although the conductivity is reduced. These interesting phenomena can be explained by the induced current distribution in the channel. The current directs itself from cathode to anode due to the fact that the electric field induced by the flow, proportional to U×B, is larger than the load application field E, which is related with the load factor, for the applied high magnetic field. On the contrary, the induced electric field will be less than the load application field E.

From Fig. 13, with the magnetic field strength increasing, the electromagnetic braking power, the power generation, and the Joule heat increase almost linearly. The electromagnetic brake power increases the fastest.

5 Conclusions

In this paper, the effect of the load factor and magnetic field strength on the power generation performance is 3D numerically studied for the Faraday type MHD generator with the He/Xe mixture gas as the working fluid. The conclusions can be summerized as follows.

(ⅰ) The volume average electrical conductivity in the channel increases with the load factor. When the average electrical conductivity increases, the Joule heating dissipation decreases. Therefore, the electrical efficiency increases.

(ⅱ) With the increase in the load factor, the enthalpy extraction increases first and then decreases. At the optimum load factor k=0.625, the enthalpy extraction rate reaches the maximum, and the power generation is also the largest.

(ⅲ) With the increase in the magnetic field, the volume average conductivity in the MHD power generation channel increases first and then decreases. At B=1.75 T, the conductivity reaches the maximum, about 11.05 S/m. The enthalpy extraction rate of the MHD power generation channel increases with the increase in the magnetic field strength. Therefore, the power generation increases.