Modeling Continuous-Electrode Magnetohydrodynamic Channels with a Demonstration through the “Sakhalin” Generator

1. INTRODUCTION

A. BACKGROUND

Sakhalin is the name of a Russian island in Northeast Asia with a conspicuously elongated form, located near the Japanese island of Hokkaido ^1,2. Sakhalin is also the name of a powerful pulsed magnetohydrodynamic (MHD) generator that has been widely used and studied in the USSR (Union of Soviet Socialist Republics), with a demonstrated electric output of 510 MW during one of its performed tests. The Sakhalin pulsed MHD generator (PMHDG) was fueled by a solid propellant (pyrotechnic powder), and it contained a single linear continuous Faraday-type supersonic channel, with a length of 4.5 m ^3–5. The working medium for Sakhalin was hot-seeded weakly-ionized solid-propellant plasma (SPP). Continuous Faraday-type channels of MHG generators have two continuous electrodes placed opposite to each other, which is a simpler configuration compared to other electric connection designs ^6,7; namely, the segmented electrodes (segmented Faraday) design, the linear Hall design, and the diagonal design ⁸. The Sakhalin was the world's most powerful PMHDG, with a theoretical peak electric output power of 600 MW (this reduces to about 400 MW with a connected matched load) ⁹. It was a self-contained system, not requiring connections to external supply sources ¹⁰. The mass of the solid propellant charge was about 6000 kg (6 tonnes). The electrodes (negative anode and positive cathode) were graphite. A possible solid propellant is a mixture of powders having 64% Mg (magnesium), 35% KNO₃ (potassium nitrate), and 1% technological additives (percentages are mass-based), pressed to a density of about 1700 kg/m^{3 11}. The propellant can also be based on metal aluminum, rather than metal magnesium ¹². A small amount of alkali metal (cesium “Cs” or potassium “K”) is added to the solid propellant charge (the grain) to permit an appreciable level of ionization of the atoms of the seeded alkali metal (which changes phase to vapor at the high temperatures of the MHD channel) ^13,14. Potassium (K) vaporizes at about 762 °C (1035 K), while cesium (Cs) vaporizes at about 682 °C (955 K) ^15–17.

Pulsed magnetohydrodynamic generators (PMHDG) can be used for producing a powerful pulse of electricity, for a few seconds. PMHDG can be used in geoelectric prospecting, where the produced electric pulse allows for estimating the electric resistivity of underground formations at different depths ^18–20. PMHDG may also be used to power electromagnetic rail launchers (railguns), where projectiles can be accelerated to speeds exceeding 2.5 km/s ^21,22. PMHDG may also find applications in space systems ^23,24 and rocket propulsion ^25,26.

Besides Sakhalin, there were other solid-propellant plasma (SPP) PMHDG in the former USSR, but Sakhalin was by large the most powerful ^27–29. A comparison of different SPP-PMHDG that existed in the former USSR is presented in Table I ^30–32.

The table is not an exhaustive list. There were even additional SPP-PMHDG systems in the USSR, such as “Pamir”, and “Pamir-06”.

Aside from the pulsed mode of magnetohydrodynamic (MHD) generators, they can also operate in a continuous mode as a non-conventional topping power cycle in a dual-cycle thermal power plant ^33–37. MHD generators have a special feature of favoring elevated temperatures, which boost the electric conductivity in a rapid nonlinear fashion ^38–40. Other conventional power cycles (such as the Rankine-type steam cycle) can then be used as a bottoming cycle, where the enormous unexploited heat in the exhaust hot gases leaving the topping MHD power cycle can be used to generate steam that drives steam turbines, coupled with electromechanical generators, as happens in combined-cycle gas turbine (CCGT) power plants via a heat recovery steam generator (HRSG) to transfer heat from the higher-temperature topping gas cycle to the lower-temperature bottoming steam cycle ^41–44.

The use of MHD generators, either pulsed or continuous, has an environmental advantage over fossil-fuel-fired power plants, where reducing the emissions of greenhouse gases (GHG) becomes possible, and this reduction helps in attenuating global warming ^45–49. This possibility of GHG reduction when using MHD generators arises because the fuel and oxidizer combination in the case of combustion-plasma-based MHD generators can be selected such that the carbon dioxide (CO₂) dominates the exploited plasma gases leaving the MHD generator, which facilitates its capture through the oxy-fuel carbon capture (OFCC) technology ^50–53. The oxy-fuel technology is particularly attractive for MHD generators, because this specialized combustion technique largely raises the temperature of the combustion flame and the combustion products, thereby improving the electric conductivity of the thermal-equilibrium “hot” plasma gas, and thus increasing the electric output from the same volume flow rate of the plasma gas ^54,55. In thermal-equilibrium “hot” plasma gas, atoms are ionized due to the temperature alone (unlike non-thermal “cold” plasma), and thus higher temperatures increase the electric output from the same volume flow rate of the plasma gas ^56,57.

MHD generators also have the advantage of not having moving mechanical parts, which is unlike gas turbines, steam turbines, wind turbines, Stirling engines, and electromechanical generators that are commonly used in power generation ^58–62. This simplicity of the design and the lack of complex fluid-structure interaction (FSI) and solid-solid contact eliminate the need for lubrication, bearing, and sealing elements ^63–65. MHD generators are direct power extraction (DPE) systems, like solar photovoltaic (PV) modules and thermoelectric generators (TEG) modules, with no intermediate energy conversion processing needed ^66–70.

One more advantage of MHD generators operating with seeded combustion plasma is their flexibility with regard to the combusted fuel. Unlike conventional internal combustion engines (ICE), either reciprocating ones used in automobiles or turbomachinery ones used in jet airplanes, the MHD generator can be viewed as an external combustion engine (ECE), because the combustion of the fuel occurs outside the MHD channel, which is the primary energy conversion component (extracting thermal, pressure, and/or kinetic energy from the fluid and converting the fluid energy into direct-current “DC” electricity) ^71,72. Because the combustion products (flue gases) act mainly as a bulk high-speed medium that carries the ionized alkali metal species and the liberated charge-carrying electrons, a wide variety of fuels can be used in the plasma generator (PG) combustor unit ^73–75. For example, green hydrogen (hydrogen produced through water electrolysis powered by renewable energy sources) may be exploited for this type of power generation ^76–81. Gaseous fuels (such as natural gas) ^82–85. Also, liquid fuels (such as kerosene) ^86–88, and solid fuels (such as coal) are candidate sources of the plasma gases ^89,90. Despite this flexibility, gas-fired and liquid-fired MHD generators have the advantage of not forming a slag layer (a layer of molten ash material that is solidified on walls) over the electrodes of the MHD channel ^91–93.

B. OBJECTIVES

The first goal of this work is to present a proposed computational fluid dynamics (CFD) solver for continuous-electrode (continuous-Faraday) MHD channels, which is based on the OpenFOAM object-oriented programming (OOP) environment ⁹⁴. OpenFOAM is a popular collection of open-source CFD solvers for a wide variety of problems, which include compressible flows, incompressible flows, laminar flows, turbulent flows, multi-phase flows, and reacting flows ^95–99. OpenFOAM has been used in many applications in fluid dynamics ^100–102, combustion ¹⁰³, fluid-structure interaction (FSI) ¹⁰⁴, and heat transfer ¹⁰⁵. OpenFOAM is based on the finite volume method (FVM) for discretizing the governing equations of the flow, while attempting to enforce local conservation of mass, momentum, and energy at the level of a single computational cell (a single finite volume) ^106,107. The proposed MHD channel solver in this study is a valuable tool in the computer-aided design (CAD) of MHD channels with continuous electrodes, either operating in a pulsed mode or operating in a sustained mode. The solver may be further extended by the interested readers to add desired design alterations or explore different operating conditions such that one or more of the performance outputs can be controlled and optimized ^108,109. For example, the influence of the area cross-section profile of the supersonic channel (as a control variable) on the electric power output, the role of the externally applied magnetic field (applied magnetic-field flux density) on the plasma streamlines, and the use of seed alkali injection to boost the electric conductivity of the expanding plasma gas may be investigated in other studies by extending the model provided in this work ^110–112.

The second goal of this work is to present comprehensive results of three-dimensional computational fluid dynamics (CFD) modeling of the channel of the Sakhalin pulsed magnetohydrodynamic generator with solid-propellant plasma (PMHDG-SPP). Unlike various reported results in the literature, the current study shows a large number of three-dimensional and sampled one-dimensional variations of fluid fields and electric fields, as well as aggregate system-level scalar values. All these results help in understanding better the complex phenomena inside the channel. CFD modeling and post-processing of its results are very beneficial in this regard, due to the digital non-intrusive probing at any spatial location of the simulated domain. Although Sakhalin is an old system, it remains an exemplary MHD study case due to the details reported about it. Therefore, this goal is viewed as a beneficial contribution to the MHD generation field, as it can guide future attempts to construct a similar system, either pulsed for geological research or sustained for electric power generation.

The third goal of this work is to provide a brief review of former solid-propellant plasma (SPP) pulsed magnetohydrodynamic generators (PMHDG).

This study contributes to regaining interest in MHD power generation, aiming for it to be further investigated as a potential way for diversifying the energy mix locally and globally ^113–116, expanding the oxy-combustion carbon capture technology, benefiting from an exceptionally high volumetric power density ^117–119, utilizing concentrated direct power extraction, facilitating the process of electrification ^120–123, and pursuing unconventional energy solutions to help in achieving net-zero emissions ¹²⁴. Furthermore, the study can be helpful for educational purposes in physics, electric power engineering, numerical methods, or computer programming at the higher education level ^125–130.

2. RESEARCH METHOD

This section aims to describe the mathematical, computational, and geometric settings of the proposed simulation modeling for generic Faraday-type magnetohydrodynamic (MHD) channels having continuous electrodes ¹³¹.

By being continuous electrically-conducting electrodes (made of metal or graphite for example), either electrode is assigned a single voltage (a single electric potential). In this study, the anode is assigned a reference voltage of zero. Thus, the “cathode voltage” here refers to the voltage difference between the higher cathode voltage and the lower anode voltage.

A. FLUID FLOW EQUATIONS

We present here the governing partial differential equations for the plasma as a homogeneous compressible turbulent gas. These equations are extended versions of the classical Favre-averaged Navier-Stokes equations (FANS), by adding supplementary terms that mathematically account for the electromagnetic effects on the electrically-conductive plasma. The Favre-averaged Navier-Stokes equations are commonly referred to as RANS (Reynolds-averaged Navier-Stokes equations), although RANS strictly applies to incompressible (constant-density) fluids, whereas the Favre-averaged Navier-Stokes equations are extension of RANS to compressible (variable-density) fluids, which pertain to MHD generators operating on combustion plasma ^132–135. The RANS formulation utilizes simple (unweighted) time averaging, whereas the formulation FANS formulation utilizes density-weighted (mass-weighted) time averaging ^136–139.

The scalar mass conservation equation (the continuity equation) of the plasma gas is identical to the one governing ordinary gases (not electrically-conducting). The presence of a magnetic field and the existence of an electric conductivity of the plasma do not influence the mass conservation law of that plasma. Therefore, the mass conservation equation is

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla \cdot \left(\rho \stackrel{⃑}{u}\right)=0$$

(1)

where $\rho$ is the mass density of the plasma gas, $t$ is the time, $\stackrel{⃑}{u}$ is the velocity vector of the plasma (this is absolute velocity, relative to a fixed laboratory frame), and $\nabla \cdot$ is the divergence operator ^140,141.

The vector linear momentum conservation equation for the plasma flow within the MHD channel is an extension of the conventional linear momentum conservation equation for a non-reacting single-phase compressible turbulent flow, where a source term is added to account for the Lorentz force per unit volume, acting on the charged particles within the plasma gas as a body force exerted on a unit volume of the plasma flowing inside the MHD channel ^142,143. To explain the cause of this source term, it should be noted that any current-passing electrically-conducting medium that is in a translational motion while subject to a magnetic field experiences an induced force called the Lorentz force, which is perpendicular to both the translational motion and the applied magnetic field ^144,145.

Therefore, the extended plasma’s linear momentum equation within the MHD channel is

$${\displaystyle \frac{\partial \rho \stackrel{⃑}{u}}{\partial t}}+\nabla \cdot \left(\rho \stackrel{⃑}{u}⨂\stackrel{⃑}{u}\right)+\nabla p-\nabla \cdot \left(\stackrel{⃑}{\stackrel{⃑}{{\tau }_v}}+\stackrel{⃑}{\stackrel{⃑}{{\tau }_t}}\right)=\stackrel{⃑}{J}\times \stackrel{⃑}{B}$$

(2)

where the operator $⨂$ represents the outer product of two vectors, $p$ is the pressure (absolute pressure), $\nabla p$ is the pressure gradient vector, $\stackrel{⃑}{\stackrel{⃑}{{\tau }_v}}$ is the viscous shear stress tensor, $\stackrel{⃑}{\stackrel{⃑}{{\tau }_t}}$ is the turbulent shear stress tensor, $\stackrel{⃑}{J}$ is the electric-current density (the current density) vector, $\stackrel{⃑}{B}$ is the magnetic-field flux density (the magnetic field) vector, and the right-hand side term $\stackrel{⃑}{J}\times \stackrel{⃑}{B}$ is the added source term (the Lorentz force per unit volume of plasma) ^146–149.

The viscous shear stress tensor $\stackrel{⃑}{\stackrel{⃑}{{\tau }_v}}$ is a deviatoric (traceless) tensor, modeled using the following constitutive relation ^150,151:

$$\stackrel{⃑}{\stackrel{⃑}{{\tau }_v}}=\mu \left(\stackrel{⃑}{\stackrel{⃑}{\nabla u}}+{\left[\stackrel{⃑}{\stackrel{⃑}{\nabla u}}\right]}^T-{\displaystyle \frac{2}{3}}\left(\nabla \cdot \stackrel{⃑}{u}\right) \stackrel{⃑}{\stackrel{⃑}{I}}\right)$$

(3)

where $\mu$ is the absolute viscosity (the dynamic viscosity) of the plasma gas, $\stackrel{⃑}{\stackrel{⃑}{\nabla u}}$ is a tensor representing the gradient of the velocity vector, ${\left[\stackrel{⃑}{\stackrel{⃑}{\nabla u}}\right]}^T$is the transpose of $\stackrel{⃑}{\stackrel{⃑}{\nabla u}}$, and $\stackrel{⃑}{\stackrel{⃑}{I}}$ is the identity tensor ^152,153.

The turbulent (or Reynolds) stress tensor ${\stackrel{⃑}{\stackrel{⃑}{\tau }}}_t$ is modeled using the following gradient transport hypothesis (eddy-viscosity turbulence modeling) ¹⁵⁴:

$$\stackrel{⃑}{\stackrel{⃑}{{\tau }_t}}={\mu }_t\left(\stackrel{⃑}{\stackrel{⃑}{\nabla u}}+{\left[\stackrel{⃑}{\stackrel{⃑}{\nabla u}}\right]}^T-{\displaystyle \frac{2}{3}}\left(\nabla \cdot \stackrel{⃑}{u}\right) \stackrel{⃑}{\stackrel{⃑}{I}}\right)-{\displaystyle \frac{2}{3}}\rho k \stackrel{⃑}{\stackrel{⃑}{I}}$$

(4)

where ${\mu }_t$ is the turbulent viscosity, $k$ is the turbulent kinetic energy per unit mass, and the spherical tensor $-{\displaystyle \frac{2}{3}}\rho k\stackrel{⃑}{\stackrel{⃑}{I}}$ appears due to the eddy-viscosity turbulence modeling ^155,156.

The sum of the laminar viscosity (the molecular viscosity) and the turbulent viscosity (the eddy viscosity) is the effective viscosity (${\mu }_{eff}$) ^157,158

$${\mu }_{eff}=\mu +{\mu }_t$$

(5)

Therefore, the sum of the viscous stress tensor and the turbulent stress tensor can be expressed as ^159,160

$$\stackrel{⃑}{\stackrel{⃑}{{\tau }_v}}+\stackrel{⃑}{\stackrel{⃑}{{\tau }_t}}={\mu }_{eff}\left(\stackrel{⃑}{\stackrel{⃑}{\nabla u}}+{\left[\stackrel{⃑}{\stackrel{⃑}{\nabla u}}\right]}^T-{\displaystyle \frac{2}{3}}\left(\nabla \cdot \stackrel{⃑}{u}\right) \stackrel{⃑}{\stackrel{⃑}{I}}\right)-{\displaystyle \frac{2}{3}}\rho k \stackrel{⃑}{\stackrel{⃑}{I}}$$

(6)

It is worth mentioning that the tensor

$\left(\stackrel{⃑}{\stackrel{⃑}{\nabla u}}+{\left[\stackrel{⃑}{\stackrel{⃑}{\nabla u}}\right]}^T-{\displaystyle \frac{2}{3}}\left(\nabla \cdot \stackrel{⃑}{u}\right)\stackrel{⃑}{\stackrel{⃑}{I}}\right)$ is the deviatoric (traceless) component of the velocity gradient tensor $\stackrel{⃑}{\stackrel{⃑}{\nabla u}}$, which has a trace of ${\displaystyle \frac{2}{3}}\left(\nabla \cdot \stackrel{⃑}{u}\right)$, and this trace quantity vanishes for incompressible fluids, but it is not zero here in the presented MHD plasma low model ^161–164.

The scalar equation of energy conservation is written here using the total (stagnation) specific internal energy, $\tilde{e}$, as the dependent variable, instead of the temperature $T$. This form is suitable for high-speed gases, as encountered in open-cycle MHD channels ^{165–168,168,169}. As with the case of the vector momentum equation, the scalar energy equation has a source term that is not present in the ordinary form (for non-electrically-conductive gases). This source is the power (per unit volume of plasma) extracted from the plasma and converted into electric power (per unit volume of plasma) to an external electric load connected to the electrodes of the MHD channel.

Thus, the extended energy equation for the plasma gas is

$${\displaystyle \frac{\partial \rho \tilde{e}}{\partial t}}+\nabla \cdot \left(\rho \stackrel{⃑}{u} \tilde{e}\right)+\nabla \cdot ({\kappa }_{eff }\nabla T)+\nabla \cdot \left(\stackrel{⃑}{u}\cdot \stackrel{⃑}{\stackrel{⃑}{\tau }}\right)=\stackrel{⃑}{J}\cdot \stackrel{⃑}{E_L}$$

(7)

where ${\kappa }_{eff}$ is the effective thermal conductivity (accounting for both laminar/molecular effect and turbulent/eddy effect in heat conduction), $-{\kappa }_{eff}\nabla T$ is the effective diffusion heat flux vector (expressed in terms of the temperature gradient $\nabla T$ according to Fourier’s law), $\nabla \cdot \left(\stackrel{⃑}{u}\cdot \stackrel{⃑}{\stackrel{⃑}{\tau }}\right)$ is the viscous-stresses power (rate of work by the plasma) per unit plasma volume due to viscous stresses, and $\stackrel{⃑}{J}\cdot \stackrel{⃑}{E_L}$ is the source term due to the electric properties of the plasma ^170,171. The vector $\stackrel{⃑}{E_L}$ is the electrostatic field (caused by the voltage difference between the MHD electrodes due to the externally applied load “if applicable” or due to the external electric circuit being open “the electrodes are not electrically connected electrically”). The vector field $\stackrel{⃑}{E_L}$ is an applied electric field on the plasma due to the voltage difference at the electrodes surrounding it, and these electrodes (grounded anode and positive cathode) are used to power an external electric load, thus we also refer to this electric field as the electrostatic field. The vector field $\stackrel{⃑}{E_L}$ is not directly dependent on the motion of the plasma, which explains the name “electrostatic”. The vector field $\stackrel{⃑}{E_L}$ can even be externally applied and ma manually controlled, by connecting the electrodes of the MHD channel to a DC power source (instead of an electric load). On the other hand; another electric field is induced due to the translational motion of the electric-conductive plasma gas with a velocity vector $\stackrel{⃑}{u}$ within an applied magnetic magnetic-field flux density $\stackrel{⃑}{B}$; this induced electric field is $\stackrel{⃑}{u}\times \stackrel{⃑}{B}$, and it can be referred to as an “electrodynamic” field. The load electric field $\stackrel{⃑}{E_L}$ is irrotational (as per Coulomb’s law) ^172–174. This means

$$\nabla \times \stackrel{⃑}{E_L}=0$$

(8)

The source term $\stackrel{⃑}{J}\cdot \stackrel{⃑}{E_L}$ actually has a negative value (the electrostatic field vector is directed against the electric-current density vector), making this source term act actually as an energy sink, which is the valid interpretation of it as an expression of the process of extracting energy from the plasma flow and supplying it as electric energy to an external electric load. In the idealized (most simplified) case of one-dimensional (unidirectional) constant plasma velocity through a constant-cross section MHD channel, one-dimensional (unidirectional) constant applied magnetic field and one-dimensional (unidirectional) constant electrostatic field being perpendicular to the parallel electrodes (from the grounded anode to the positive cathode), and with the three fields are mutually orthogonal to each other; this energy source term reduces to $-J E_L$. Thus, it is actually an energy sink term with a magnitude of $\left|J E_L\right|$, where $J$ and $E_L$ reduce to scalar values (but the current density $J$ is from the anode to the cathode, while the electrostatic field $E_L$ is from the cathode to the anode) ¹⁷⁵.

The total specific internal energy $\tilde{e}$ is the sum of the static specific internal energy $e$ and the specific kinetic energy $0.5\stackrel{⃑}{u}\cdot \stackrel{⃑}{u}$. Mathematically, this is expressed as ¹⁷⁶

$$\tilde{e}=e+{\displaystyle \frac{1}{2}}\stackrel{⃑}{u}\cdot \stackrel{⃑}{u}$$

(9)

The static specific internal energy $e$ is related to the static specific enthalpy $h$ as ^177,178

$$h=e+p/\rho$$

(10)

Likewise, the total specific internal energy $\tilde{e}$ is related to the total specific enthalpy $\tilde{h}$ as ^179,180

$$\tilde{h}=\tilde{e}+p/\rho$$

(11)

This means

$$\tilde{h}=h+{\displaystyle \frac{1}{2}}\stackrel{⃑}{u}\cdot \stackrel{⃑}{u}$$

(12)

The ideal gas law (ideal gas equation of state) is also needed in the MHD model, as an algebraic relation between the absolute pressure $p$, absolute temperature $T$, and density $\rho$. The ideal gas law is not affected by the electric characteristic of the weakly-ionized seeded combustion plasma, which is still treated as a gas thermodynamically. The ideal gas law has the form ^181,182

$$p=\rho R T$$

(13)

where $R$ is the specific gas constant for the plasma gas, which is related to the universal (molar) gas constant $\overline{R}$ and the plasma molecular weight $W$ as ¹⁸³

$$R=\overline{R}/W$$

(14)

The value of $\overline{R}$ is 8.314 J/(mol.K) or 8,314 J/(kmol.K) ^184,185. If the plasma is treated as a mixture of two of more gaseous species, then the molecular weight $W$ depends on the local chemical composition of the plasma (at a given computational cell) as ¹⁸⁶

$$W={\sum }_i^N{X}_i W_i$$

(15)

where $X_i$ is the mole fraction of the i^th gaseous species in the plasma gas mixture (such as carbon dioxide CO₂ or water vapor H₂O), $W_i$ is the molecular weight of that i^th gaseous species, and N is the number of gaseous species in the plasma gas mixture.

B. ELECTRIC FIELDS EQUATIONS

In the previous subsection, it was shown that three electromagnetic vector fields are needed in order to integrate the system of partial differential equations governing the density, velocity, and the total specific internal energy of the plasma flow within the MHD channel. These three electromagnetic vector fields appear in the two source terms that arise in the extended vector momentum equation and the extended scalar energy equation. These three needed electromagnetic vector fields are (1) the electric-current density vector (or simply the current density) $\stackrel{⃑}{J}$, (2) the magnetic-field flux density vector (or simply the magnetic field) $\stackrel{⃑}{B}$, (3) and the electrostatic vector field $\stackrel{⃑}{E_L}$.

The magnetic-field flux density $\stackrel{⃑}{B}$ is a specified applied field. Thus, in the proposed MHD model here, only the current density field $\stackrel{⃑}{J}$ and the electrostatic field $\stackrel{⃑}{E_L}$ need to be estimated.

The approach proposed here for handling the electric aspect of the problem utilizes a scalar electric potential $\varphi$ field as a working variable that is directly solved for, and then its solution is used to indirectly compute the fields $\stackrel{⃑}{J}$ and $\stackrel{⃑}{E_L}$.

The electric potential scalar field $\varphi$ is related to the electrostatic vector field $\stackrel{⃑}{E_L}$ as follows:

$$\stackrel{⃑}{E_L}=-\nabla \varphi$$

(16)

Thus, the gradient of the electric potential is the electrostatic field, but in the opposite direction. The above mathematical relation is made possible given that the electrostatic field $\stackrel{⃑}{\stackrel{⃑}{E_L}}$ is irrotational ^187–190.

The electric potential $\varphi$ is governed by a Poisson-type parabolic partial different equation as ¹⁹¹

$$\nabla \cdot \left(\stackrel{⃑}{\stackrel{⃑}{{\sigma }_{eff}}} \cdot \nabla \varphi \right)=\nabla \cdot \left(\stackrel{⃑}{\stackrel{⃑}{{\sigma }_{eff}}}\cdot \left( \stackrel{⃑}{u}\times \stackrel{⃑}{B}\right)\right)$$

(17)

where $\stackrel{⃑}{\stackrel{⃑}{{\sigma }_{eff}}}$ is an effective electric conductivity tensor for the plasma gas. For the case of unidirectional magnetic-field flux density in the $Z$ direction only (thus, $\stackrel{⃑}{B}\equiv [0, 0, B]$), with $B$ being the magnitude of that magnetic-field flux density (which can be a function of space rather than being a uniform constant), the effective electric conductivity tensor has the following form:

$$\stackrel{⃑}{\stackrel{⃑}{{\sigma }_{eff}}}=\sigma \left[\begin{array}{ccc}1& -{\displaystyle \frac{\beta }{{\beta }^2+1}}& 0\\ {\displaystyle \frac{\beta }{{\beta }^2+1}}& 1& 0\\ 0& 0& 1\end{array}\right]$$

(18)

where $\sigma$ is the scalar electric conductivity of the plasma, and $\beta$ is the dimensionless Hall parameter (electron’s Hall parameter), which is the product of the electron mobility ${\mu }_e$ and the magnitude of the applied magnetic field $B$ ^192–194. Therefore,

$$\beta ={\mu }_{e}B$$

(19)

The effective electric conductivity tensor is used to express the generalized form (the vector form) of Ohm’s law for the electrically-conductive plasma, as ^195,175,196

$$\stackrel{⃑}{J}=\stackrel{⃑}{\stackrel{⃑}{{\sigma }_{eff}}}\cdot \left(\stackrel{⃑}{E_L}+\stackrel{⃑}{u}\times \stackrel{⃑}{B}\right)$$

(20)

It should be noted that the vector $\stackrel{⃑}{E_L}+\stackrel{⃑}{u}\times \stackrel{⃑}{B}$ is the net electric field acting internally on the plasma (as an electrically-conductive medium having its own internal electric resistivity), and it accounts for both the electrostatic field $\stackrel{⃑}{E_L}$ and the electrodynamic field $\stackrel{⃑}{u}\times \stackrel{⃑}{B}$. As an illustrative analogy with a simple battery-lamp circuit, the field $\stackrel{⃑}{u}\times \stackrel{⃑}{B}$ is likened to the battery voltage, the field $\stackrel{⃑}{E_L}$ is likened to the voltage drop across the lamp, and the field $\stackrel{⃑}{E_L}+\stackrel{⃑}{u}\times \stackrel{⃑}{B}$ is likened to the voltage drop across the battery itself due to its internal electric resistance. Furthermore, the absolute value of the energy equation source term $\left|\stackrel{⃑}{J}\cdot \stackrel{⃑}{E_L}\right|$ can be likened to the power dissipated in the lamp. In the plasma’s energy equation, this term has a negative value because, for the plasma, it is power extracted (removed).

It should also be noted that the proposed solution method for the electric fields assumes that the magnetic-field flux density $B$ is fixed in time (but it does not need to be fixed in space), and thus the magnetic-field flux density $B$ is not affected by the plasma flow. This corresponds to the low magnetic Reynolds number approximation ^197,198. The magnetic Reynolds number (Re_m) is a dimensionless quantity that helps in quantifying the interaction between the moving plasma and the magnetic field ^199,200. Under large magnetic Reynolds numbers (Re_m >> 1), the magnetic field lines can be likened to elastic bands that are frozen into the electrically-conductive plasma ²⁰¹. The velocity of the plasma induces a secondary induced magnetic-field flux density $\stackrel{⃑}{b}$, which is added to the applied magnetic-field flux density $\stackrel{⃑}{B}$, making the total magnetic-field flux density actually dependent on the plasma flow. This additional flow-dependent magnetic-field flux density then affects the plasma velocity through the source term in the plasma momentum equation, presented in (2). Thus, at high magnetic Reynolds numbers, a two-way coupling exists between the magnetic-field flux density seen by the traveling plasma and the velocity field of that plasma. On the other hand, under small magnetic Reynolds numbers (Re_m << 1), the movement of the plasma has a negligible effect on the magnetic-field flux density acting on it. Therefore, it is reasonable in that case to treat the magnetic-field flux density acting upon the plasma as the fixed externally-applied $\stackrel{⃑}{B}$ field only. This means that the low magnetic Reynolds number approximation treats the coupling between the magnetic-field flux density and the plasma motion as having a one-way direction, with the externally-applied magnetic-field flux density affecting the motion of the plasma while the induced current density within the moving plasma does not distort the magnetic field through producing an additional induced one ²⁰². The low magnetic Reynolds number (low-Re_m) approximation is also called an “inductionless” approximation ²⁰³. On the other hand, when the induced magnetic-field flux density is accounted for, this corresponds to a “self-excitation” or “self-excited” system ²⁰⁴.

Mathematically, the magnetic Reynolds number (Re_m) is defined as

$${Re}_m={\mu }_m^{*} {\sigma }^{*} u^{*} L^{*}$$

(21)

where ${\mu }_m^{\mathrm{*}}$ is a reference magnetic permeability. The magnetic permeability of vacuum (the permeability of free space or the magnetic constant) ${\mu }_0$ is a universal constant with the exact value of 4π × 10^–7 H/m ^205–207. In the above equation defining Re_m, ${\sigma }^{\mathrm{*}}$is a reference electric conductivity of the electrically-conductive medium, $u^{\mathrm{*}}$is a reference scalar speed for the electrically-conductive medium, and $L^{\mathrm{*}}$ is a reference length relevant to the geometric scales of the problem.

The definition of the magnetic Reynolds number (Re_m) in (21) can be reformulated as

$${Re}_m={\displaystyle \frac{u^{*} L^{*}}{{\mu }_m^{*} {\sigma }^{*}}}$$

(22)

where the product ${\mu }_m^{\mathrm{*}}{\sigma }^{\mathrm{*}}$ has the dimension of (area/time), which is the same as the momentum diffusivity (the kinematic viscosity) for fluids. Therefore, the quantity ${\mu }_m^{\mathrm{*}}{\sigma }^{\mathrm{*}}$ may be called a reference “magnetic diffusivity”, and assigned the symbol ${\nu }_m^{\mathrm{*}}$ that resembles the symbol ${\nu }^{*}$ used for designating a reference momentum diffusivity (a reference kinematic viscosity) in fluid mechanics ²⁰⁸.

Equation (22) thus can be written as

$${Re}_m={\displaystyle \frac{u^{*} L^{*}}{{\nu }_m^{*}}}$$

(23)

The above form is very similar to the one used for defining the conventional Reynolds number (Re) in fluid mechanics applications, which is ^209,210

$$Re={\displaystyle \frac{u^{*} L^{*}}{{\nu }^{*}}}$$

(24)

This conventional Reynolds number (Re) in fluid mechanics is used as a nondimensional criterion to infer whether the fluid flow is expected to be in a laminar (regular and smooth) regime or a turbulent (irregular with a lot of mixing) regime. However, unlike the magnetic Reynolds number (Re_m); the threshold of regime transitioning for the conventional Reynolds number (Re) tends to be much larger than unity, and this threshold is dependent on the type of flow (for example, internal flows inside a cylindrical pipe have a lower threshold than external flows over the same cylindrical pipe).

The low magnetic Reynolds number (low-Re_m) approximation is relevant in MHD generators, and it is adopted in the proposed MHD model here ^211,212. To demonstrate this relevance, we consider a representative situation of a weakly-ionized plasma gas for MHD channels, formed by seeding hot air at 3,000 °C (3,273 K) with a small amount of potassium vapor (2% fraction). Under a reasonable reference length of 1 m and a reasonable reference speed of 1,000 m/s, the magnetic Reynolds number (Re_m) is near 0.13 only, while the conventional Reynolds number (Re) is about 3.5 × 10^{6 213}. For the Sakhalin pulsed MHD generator, the induced magnetic-field flux density $\stackrel{⃑}{b}$ is small compared to the external magnetic-field flux density $\stackrel{⃑}{B}$, with the peak of $\stackrel{⃑}{b}$ (about 0.19 T) being less than 10% of the peak of $\stackrel{⃑}{B}$ (about 2.1 T). This induced magnetic-field flux density causes the overall magnetic-field flux density (the sum of $\stackrel{⃑}{B}$ and $\stackrel{⃑}{b}$) to be less than the externally applied $\stackrel{⃑}{B}$ at the first region of the channel (nearly in its first half), but more than the externally applied $\stackrel{⃑}{B}$ in the second region of the channel (nearly in its second half), and the alteration due to $\stackrel{⃑}{b}$ are smallest near the middle of the longitudinal dimension of MHD channel (around $X$ = 2.25 m). The influence of the induced $\stackrel{⃑}{b}$ on the estimated electric power generation is only 2.5% ²¹⁴. So, our low magnetic Reynolds number (low-Re_m) approximation is considered suitable, despite that the magnetic Reynolds number (Re_m) for the Sakhalin PMHD generator is not very small compared to unity, but it is still less than one (about 0.5).

For a solid conductor, the generalized form (vector form) of Ohm’s law, as presented in (20), reduces to the classical scalar form of Ohm’s law, which is

$$I=V/R$$

(25)

where $I$ is the electric current, $V$ is the voltage difference (voltage drop across the electric load), and $R$ is the electric resistance of the electric load.

Using (18) in (20), and expanding the resulting expression from a vector equation into three scalar equations gives the following three Cartesian components of the electric density vector, $\stackrel{⃑}{J}\equiv [J_x, J_y, J_z$]:

$$J_x={\displaystyle \frac{\sigma }{{\beta }^2+1}}\left[\left(E_{L,x}+v B\right)- \beta \left(E_{L,y}-u B\right)\right]$$

(26)

$$J_y={\displaystyle \frac{\sigma }{{\beta }^2+1}}\left[\left(E_{L,y}-u B\right)+\beta \left(E_{L,x}+v B\right)\right]$$

(27)

$$J_z=\sigma E_{L,z}$$

(28)

where $u,v,\mathrm{a}\mathrm{n}\mathrm{d}w$ are the three Cartesian components of the plasma velocity vector $\stackrel{⃑}{u}$, and $E_{L,x},E_{L,y},\mathrm{a}\mathrm{n}\mathrm{d}{E}_{L,z}$ are the three Cartesian components of the electrostatic vector $\stackrel{⃑}{E_L}$. Thus, $\stackrel{⃑}{u}\equiv [u,v,w]$. and $\stackrel{⃑}{E_L}\equiv [E_{L,x},E_{L,y},E_{L,z}]$. It can be noticed that while there is a coupling between $J_x$ and $J_y$, they are decoupled from $J_z$. This decoupling is a result of assuming the magnetic-field flux density $\stackrel{⃑}{B}$ to be unidirectional, thus $\stackrel{⃑}{B}\equiv [0, 0, B]$ in the proposed MHD model.

The Poisson equation (17) for the electric potential $\varphi$ has the electric potential in the left-hand side only, which is treated with implicit discretization, while the right-hand side has the flow velocity, which is discretized explicitly.

When the electric potential $\varphi$ is obtained by integrating the Poisson equation (17) numerically, a Dirichlet boundary condition (a first-type boundary condition) is applied at the anode (the top MHD boundary), by enforcing the value of $\varphi$ to be zero, which is an arbitrary value that represents electric grounding ^215–217. Another Dirichlet boundary condition is applied at the positive cathode electrode (the bottom MHD boundary), by enforcing the value of $\varphi$ there to be equal to the aimed value for the MHD operation (the target voltage drop needed for powering the external electric load, or the target open-circuit voltage if no load is connected). Therefore, the two boundary conditions for the electric potential $\varphi$ at the two electrodes are

$${\varphi }_{anode}=0$$

(29)

$${\varphi }_{cathode}={\Delta V}_L$$

(30)

where ${\Delta V}_L$ is the target voltage difference across the electrodes and thus across the external connected load (or the open-circuit voltage).

For the two other boundaries of the MHD channel (the non-electrically-conductive walls), a Neumann boundary condition (a second-type boundary condition) is applied, which has the form ^218–220

$${\displaystyle \frac{\partial \varphi }{\partial n}}\left(or \nabla \varphi \cdot \widehat{n}\right)={\left(\stackrel{⃑}{u}\times \stackrel{⃑}{B}\right)}_{boundary}\cdot \widehat{n}$$

(31)

where $\widehat{n}$ is a unit normal vector pointing perpendicularly away from the MHD boundary wall (as if leaving the channel laterally).

After obtaining the electric potential field $\mathit{\varphi }$, the electrostatic field $\stackrel{⃑}{{\mathit{E}}_{\mathit{L}}}$ is computed as the negative gradient of $\mathit{\varphi }$ using (16), which is repeated below for a convenient flow of the analysis.

$$\stackrel{⃑}{E_L}=-\nabla \varphi$$

(16)

Finally, the generalized Ohm’s law is used to compute the electric-current density $\stackrel{⃑}{J}$ according to (20), which is repeated below for a convenient flow of the analysis.

$$\stackrel{⃑}{J}=\stackrel{⃑}{\stackrel{⃑}{{\sigma }_{eff}}}\cdot \left(\stackrel{⃑}{E_L}+\stackrel{⃑}{u}\times \stackrel{⃑}{B}\right)$$

(20)

C. GEOMETRY AND MESH OF THE SAKHALIN MHD CHANNEL

The geometry and operational settings for the modeled Sakhalin MHD channel are illustrated in Figure 1. The bulk flow direction of the supersonic plasma is along the $X$-axis. The unidirectional magnetic field depends only on the axial distance ($X$) from the channel inlet, and the range of $X$ is from 0 to 4.5 m. The anode is the upper surface of the channel (having positive linearly-increasing $Y$ coordinates), while the cathode is the lower surface of the channel (having negative linearly-decreasing $Y$ coordinates). The width of the channel (along the Z direction) is constant, with a value of 1 m. The height (along the Y direction) of the channel increases from 0.9 m at the inlet to 1.6 m at the outlet.

Projected views (front view and top view) of the MHD channel are provided in Figure 2. The front view is a trapezium (a trapezoid), while the top view is a rectangle.

The spatial discretization of the Sakhalin MHD channel is illustrated in Figure 3. Due to its geometric simplicity, the channel is discretized spatially as a single block, with a total of 376,380 hex cells (hexahedral computational cells). We used 180 cells along the axial $X$ dimension (the length), 51 cells along the vertical $Y$ dimension (the height), and 41 cells along the lateral $Z$ dimension (the width). Hexahedral (structured) meshes tend to be favored over tetrahedral (unstructured) meshes due to the higher accuracy, and the ability to have a larger aspect ratio ^221–223. Similarly, in two-dimensional problems, quadrilateral cells are computationally preferred over triangular cells ^224,225. Therefore, when possible, a tetrahedral mesh may be limited to complex geometries that cannot be represented well through hexahedral mesh elements. In our CFD simulation of the Sakhalin MHD channel, tetrahedral cells were not necessary.

Our selected spatial discretization makes each cell nearly a cube (cells have a comparable size of their edges, near 0.025 m), and this rough regularity in the mesh is preferred in CFD simulations compared to highly skewed or stretched cells ^226,227. Furthermore, our selection of an odd number of cells for the height (51 cells) and the width (41 cells) ensures that the longitudinal centerline of the channel coincides with a sequence of cell centers, which are also equally spaced. This is helpful because we sample the simulation results along that longitudinal centerline, which starts from the $[X, Y, Z]$ coordinates [0,0,0] and ends at the point [4.5,0,0]. Thus, exact cell-center values are readily available for sampling, without approximating interpolation.

D. EMPIRICAL EXPRESSIONS FOR SAKHALIN

The previous mathematical description of the governing equations and computational modeling for the MHD plasma is not restricted to the Sakhalin pulsed magnetohydrodynamic generator (PMHDG). However, a number of thermofluidic and electromagnetic properties appear in these governing equations, which should be modeled in the computational fluid dynamics (CFD) simulations using customized submodels ^228–231. Such properties include the specific heat capacity at constant pressure $C_p$, the magnitude of the unidirectional magnetic-field flux density $B$ the Hall parameter $\beta$, and the scalar electric conductivity $\sigma$.

The Sakhalin PMHDG has empirical expressions published in the literature that describe the aforementioned four quantities ²³². Such expressions are listed in Table II. In addition, the same table lists other operational conditions used in our simulation based on published data in the literature.

For the magnitude of the magnetic-field flux density $B$, the expression listed in the table is our seventh-degree polynomial fitting, which we obtained using the (curve_fit) function of the (optimize) package for optimization and root finding, provided with the (SciPy) open-source Python library ^233–235. The function (curve_fit) utilizes a non-linear least squares method. Our proposed fitting function is visualized in Figure 4.

E. COMPUTATIONAL SETTINGS

As with other computational fluid dynamics (CFD) simulations, a number of choices need to be made with regard to discretization schemes and/or numerical parameters. In this subsection, we cover this aspect of the OpenFOAM model presented here for magnetohydrodynamics (MHD) channels by listing some configuration settings adopted in the simulation in Table III.

The simulated time of 0.04 s was selected because the solution reached a quasi-steady state, with negligible transience. Although the governing equations are unsteady, the simulation solution reached a stable stage.

OpenFOAM Code Lines

We provide below code lines from the proposed OpenFOAM magnetohydrodynamic (MHD) solver, based on our mathematical description provided earlier.

For example, the first part of these code lines shows how the Poisson equation (17) for the electric potential $\varphi$ is solved using the high-level OpenFOAM programming language, through an iterative loop with (nCorr) iterations. In the code, the electric potential $\varphi$ is assigned the variable name (phi_e).

The remaining parts of the code lines show how the scalar mass conservation equation is solved, followed by the vector momentum equation, and finally the scalar energy equation ^255,256.

for (int i = 0; i < nCorr; i++)

  {

   fvScalarMatrix phiEqn

   (

    fvc::div( sigmaEff & (U ^ B) ) - fvm::laplacian(sigmaEff,phi_e)

   );

   phiEqn.relax();

   phiEqn.solve();

   E = -fvc::grad(phi_e);

   J = sigmaEff & ( E + (U^B) );

   lorentzForceDensity = J ^ B;

   powerDensity = E & J;

  }

solve(fvm::ddt(rho) + fvc::div(phi));

solve

   (

   fvm::ddt(rho, U) - fvc::ddt(rho, U)

   - fvm::laplacian(muEff, U)

   - fvc::div(tauEffExplicit)

   );

surfaceScalarField TauEffDotUDotSf

   (

   (

    fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)

   + (mesh.Sf() & fvc::interpolate(tauEffExplicit))

   )

   & (frac_ap*U_ifAllFluxesOut + frac_am*U_ifAllFluxesIn)

  );

solve

   (

   fvm::ddt(rhoE)

   + fvc::div(phiE_pU)

   - fvc::div(TauEffDotUDotSf)

   - powerDensity

   );

e = rhoE/rho - 0.5*magSqr(U);

e.correctBoundaryConditions();

volScalarField kEff("kEff", thermo.Cp()*turbulence->alphaEff());

solve

   (

    fvm::ddt(rho, e) - fvc::ddt(rho, e)

   - fvc::laplacian(kEff, T)

   );

rhoE = rho*(e + 0.5*magSqr(U));

3. VALIDATION AND VERIFICATION

After presenting the mathematical/computational formulation of the OpenFOAM solver developed for handing continuous Faraday channels in MHD generators, this section is provided to compare some results we obtained using the solver with values previously published in the literature, and this is a validation step for our CFD solution of the Sakhalin pulsed magnetohydrodynamics generator (PMHG), through external benchmarking. In the verification step, we perform internal benchmarking by comparing the CFD solution at the nominal spatial and temporal resolutions with another CFD solution using the same solver but at finer spatial and temporal resolutions; this verifies the ability of the CFD model to converge at the nominal resolutions ^257–259.

Before presenting the validation step; we point out that due to the uncertainty in reported data of the Sakhalin PMHDG where a range for a variable is reported, rather than a single value for it (such as the mass flow rate), the discrepancy of some reported operational characteristics (such as the MHD channel inlet velocity), the time-unsteadiness nature of the actual problem, our simplification introduced in the inlet conditions, and the empirical fitting functions used in modeling some parameters; we do not aim at achieving a close quantitative agreement with the literature information ²⁶⁰. Instead, we seek to obtain values of aggregate quantities that are comparable to those reported by others. We consider three key benchmarking quantities for the validation step, as summarized in Table IV. Based on the acceptable small deviation in this table, our simulation is considered valid.

Table V summarizes the verification step for our CFD simulation of the Sakhalin PMHDG channel using the proposed OpenFOAM solver. In this step, we generated an internal benchmarking solution using the same proposed OpenFOAM solver but at a finer spatial resolution and a finer temporal resolution. This auxiliary fine-resolution solution acts as a higher-fidelity reference to assess the nominal (the solution obtained with the nominal resolutions). Based on the observed small changes in the compared characteristic of the MHD channel, as shown in the table, our verification step is considered successful, and our nominal OpenFOAM CFD solution is considered satisfactorily resolution-independent. All the remaining results in this study correspond to the nominal resolutions (which are 376,380 cells and 5 × 10^–6 s time-step).

The above verification table also has indications about the validity of our OpenFOAM model, through the correct sign of the $Y$-component of the volumetric Lorentz force density (force per unit volume), which is pointing from the lower cathode to the upper anode for our case where the plasma’s bulk motion is in the positive $X$ direction, and the magnetic-field flux density is unidirectional in the positive $Z$ direction. The Lorentz force vector $\stackrel{⃑}{F}$ on a general electrically charged particle with a positive charge $q$, when this particle is carried with the plasma gas (moving at the same bulk velocity $\stackrel{⃑}{u}$) is ²⁶⁵

$$\stackrel{⃑}{F}=q \left({\stackrel{⃑}{E}}_L+\stackrel{⃑}{u}\times \stackrel{⃑}{B}\right)$$

(32)

Although electrostatic field ${\stackrel{⃑}{E}}_L$ reduces or alters the electrodynamic fields $\stackrel{⃑}{u}\times \stackrel{⃑}{B}$, the resultant overall electric field ${\stackrel{⃑}{E}}_L+\stackrel{⃑}{u}\times \stackrel{⃑}{B}$ is like electrodynamic fields $\stackrel{⃑}{u}\times \stackrel{⃑}{B}$ in terms of pointing vertically down (in the negative $Y$ direction, from the upper anode to the lower cathode).

Because the effective mobile charge carriers in the plasma are the negatively-charged electrons (not the positively-charged seed ions), the charge $q$ is actually negative, and thus the Lorentz force actually points upward, and this is in agreement with our CFD predictions ^266,267.

The finding that the averaged electric-current density vector is pointing opposite to the averaged Lorentz force vector is also a favorable indication that the CFD model is able to predict correct vector fields. While the collected conventional electric current flows externally through the electric load from the positive cathode to the grounded anode, the electron electric current flows externally in the opposite direction from the grounded anode to the positive cathode. In the presented OpenFOAM model, the electric-current density refers to the conventional electric current. Therefore, its averaged direction within the MHD channel should be pointing down (thus, having a negative $Y$-component), which is correctly predicted in out simulation.

Finally, the obtained zero $Z$-component of either the averaged Lorentz force density vector or the averaged electric-current density vector is a logical outcome, given the unidirectional $Z$-oriented nature of the applied magnetic-field flux density, and this does not cause any electric-current density or Lorentz force along the $Z$ direction. There is no reason for having an averaged $Z$-component of the electrostatic field (thus, $E_{L,z}$ = 0).

4. RESULTS

We present the results of the proposed OpenFOAM solver when applied to the Sakhalin magnetohydrodynamic (MHD) channel in the next three subsections.

In the first subsection, we present scalar aggregate quantitative values. In the second subsection, we present one-dimensional graphical profiles. In the third subsection, we present three-dimensional visualizations of selected fields in the entire MHD channel.

Scalar Values at the MHD Channel Outlet

By post-processing our simulation results, we obtain the mass-averaged values of six key selected properties of the plasma as it exits the MHD channel. These values are summarized in Table VI.

A generic mass-averaged quantity (say $\psi$) at the outlet is computed through density-weighted integration of that quantity ${\psi }_{\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}-\mathrm{a}\mathrm{v}\mathrm{g}}$ at the outlet area, and then dividing the mass flow $\dot{m}$ leaving the channel form that outlet area. Thus,

$${\psi }_{mass-avg}={\displaystyle \frac{1}{\dot{m}}}{\int }_{outlet area}\rho u \psi dA$$

(33)

where $u$ is the axial component of the exiting plasma (this component is perpendicular to the outlet surface), $dA$ is an infinitesimal area element at the outlet (which has a total finite area of 1.6 m²).

To better demonstrate the change in the plasma flow due to the direct power extraction (DPE) process within the MHD channel, the inlet values in our OpenFOAM simulation for the same six plasma properties are also listed in the table, which helps in contrasting them with the computed outlet values.

It can be seen that one major influence of the power extraction process from the plasma is an occurring deceleration process, although supersonic flows accelerate in divergent channels (as in the case of the Sakhalin channel) ^268,269.

On the other hand, the temperature and the pressure nearly did not change from their inlet values. This means that the main source of the electricity extracted from the flowing plasma is its kinetic energy, rather than its internal energy (thermal energy) or its pressure energy (flow energy).

Due to the weak change in the pressure and temperature, and because the electric conductivity and the Hall parameter are fitted as dependent functions of these two thermodynamic properties (the Hall parameter depends additionally on the magnitude of the applied magnetic-field flux density), the electric conductivity and the Hall parameter of the plasma do not change largely when the outlet location is compared with the inlet location.

Because the Mach number depends largely on the axial velocity of the plasma, and because the speed of sound is dependent on the temperature which does not show a large change between the inlet and the outlet locations of the MHD channel; the Mach number decreases strongly between the inlet and the outlet locations, similar to the change in the axial plasma velocity. This decrease of a supersonic Mach number for a flow within a divergent channel or nozzle also emphasizes the loss of kinetic energy from the plasma, which is actually decelerated within the channel because of the electric power extraction.

Centerline Profiles

Through one-dimensional sampling along the centerline of the MHD channel (thus, along all $X$ values from 0 to 4.5 m, but with $Y$ = 0 and $Z$ = 0), the variations of ten thermodynamic, fluidic, or electric properties of the simulated Sakhalin channel are illustrated in the current subsection. None of the investigated profiles has a discontinuity (abrupt change), which is an implicit indication of the lack of internal instabilities in the performed simulation, which agrees with the expectation for this problem and supports the correctness of the presented model.

We start with the electric conductivity $\sigma$, and display its centerline variation in Figure 5. We point out here that this plasma electric variable is described through an empirical nonlinear regression model as a function of the local temperature and the local pressure. Therefore, its variation reflects variation in either of these two independent variables. Despite the apparent non-uniformity of the electric conductivity, its variation is smooth and bounded within a narrow range of about 4 S/m only (from 48.6 S/m to 52.8 S/m).

Next, we display the variations of the Hall parameter $\beta$ along the centerline of the MHD channel in Figure 6. Like the electric conductivity, the Hall parameter is described through an empirical nonlinear regression model as a function of the local temperature and the local pressure. In addition, the Hall parameter depends linearly on the magnitude of the applied magnetic-field flux density $B$. The magnitude of the applied magnetic-field flux density declines rapidly near the inlet and outlet of the MHD channel, and this is clearly shown in the profile of the Hall parameter. In the internal region of the channel, the Hall parameter qualitatively follows a similar profile to that of the unidirectional applied magnetic-field flux density, which can be explained by the weak influence of the temperature and the pressure because they do not change significantly within the MHD channel. The maximum encountered value of the Hall parameter is 0.39, while the lowest value is at the MHD outlet, and it is 0.26.

Figure 7 shows the centerline profile of the Mach number, which declines in a nearly monotonic pattern, leaving the MHD channel at a noticeably smaller value of 1.26 than its inlet value of 2.05. However, the plasma remains supersonic even at the outlet (the exit Mach number is still above unity). Therefore, the exiting plasma still has a lot of kinetic energy that can be exploited in a subsequent energy conversion system.

Figure 8 shows the centerline profile of the absolute pressure, which manifests slow variations within 0.6 bar. The leaving plasma still has a high pressure close to the inlet pressure (in addition to the still-high speed), thus it is subject to further energy extraction through another energy system (such as gas turbines) after deceleration to proper subsonic speeds using supersonic diffusers ^244,270.

Figure 9 and Figure 10 show the centerline profiles of the plasma density, and the plasma absolute temperature, respectively. It is interesting that the temperature is practically constant near its inlet value of 2,750 K. Following the ideal gas law, the density in such a case (a nearly isothermal process) becomes essentially a scaled version of the absolute pressure, which is a satisfied feature that can be confirmed by comparing the previous profile of the absolute pressure with the profile of the density. The high temperature of the leaving plasma indicates an enormous amount of thermal energy that can be transferred through a heat exchanger to another fluid, for example, to generate steam that can be used in a steam power cycle. Using the empirical formula for the specific heat capacity at constant pressure for the plasma of the Sakhalin PMHDG, namely $C_p\left(T\left[\mathrm{K}\right]\right)\left[\mathrm{J}/(\mathrm{k}\mathrm{g}.\mathrm{K})\right]=649.863+0.620932T$, the change in the specific enthalpy, designated by $\mathsf{\Delta }h\equiv h_2-h_1$, between two states (a lower-temperature state 1 and a higher-temperature state 2) can be obtained by integrating this $C_p$ expression as

$$\Delta h\equiv h_2-h_1={\int }_{T_1}^{T_2}{C}_p\left(T\right)dT$$

(34)

This gives

$$\Delta h=649.863 \left(T_2-T_1\right)+0.310466 (T_2^2-T_1^2)$$

(35)

If the higher-temperature state has $T_2$ = 2,750 K (nearly the plasma exit temperature), and the lower-temperature state has $T_1$ = 300 K (nearly the ambient temperature), then the enthalpic energy available per kg of plasma is 3.912 × 10⁶ J/kg (3.912 MJ/kg). This is more than enough to generate superheated steam at extremely high temperatures (even exceeding 700 °C), which is suitable to operate steam power plants ^271,272.

Figure 11 and Figure 12 show the centerline profiles of the axial and upward velocity components of the plasma, respectively. As was the case for the Mach number, the centerline axial velocity declines as the plasma progresses through the MHD channel, reflecting a loss of kinetic energy that is converted into electric energy. The upward velocity component (the plasma velocity in the $Y$ direction, being positive if the plasma is moving from the downward cathode to the upward anode) varies mildly around the zero value, bounded between –50 m/s and 30 m/s.

Figure 13 shows the centerline profile of the third velocity component of the plasma (the $Z$-component), which is along the unidirectional applied magnetic-field flux density. This component is practically zero at the centerline, reflecting perfect symmetry and balance of lateral forces. This is reasonable given the lack of gravitational effect along the $Z$ direction, as well as the constant channel width along the same direction.

Figure 14 shows the centerline profile of the downward component of the electric-current density. This component has negative values; thus the conventional electric-current density is directed toward the bottom cathode (although the physical charge-carrying electrons are forced upward toward the upper anode). In the figure, the negative value of that already-negative conventional component is displayed, which is equivalent to displaying the absolute value of that negative conventional component. The displayed absolute conventional current density component increases nonlinearly until reaching a maximum value of 6.1 A/cm² at a distance of 1 m from the MHD channel inlet; it then declines nonlinearly until reaching a zero value at the channel outlet.

B. THREE-DIMENSIONAL DISTRIBUTIONS

In this subsection, we present three-dimensional views for eight selected fields as per the solution obtained in our simulation of the Sakhalin MHD channel. These views are instrumental in demonstrating the variation of each field in the entire channel.

Figure 15 displays the spatial distribution of the electric conductivity, which changes mildly between 43.9 S/m and 53.6 S/m. Again, the variations in this electric conductivity are determined by the variations in the temperature or pressure.

Figure 16 displays the spatial distribution of the Hall parameter, which changes between 0.253 and 0.424. The variations in the Hall parameter are determined by variations in the temperature, pressure, and the applied magnetic-field flux density.

Figure 17 displays the spatial distribution of the absolute pressure, which changes between 2.56 bar and 3.62 bar. A pressure recovery phenomenon is observed near the outlet, where the pressure increases back after its initial drop after entering the channel.

Figure 18 displays the spatial distribution of the axial velocity component of the plasma, which changes between 832.9 m/s and 2,089 m/s. The deceleration process is efficiently visible through this figure, with the axial velocity generally dropping as the plasma moves from the inlet toward the outlet.

Figure 19 displays the spatial distribution of the specific turbulent kinetic energy (the turbulent kinetic energy per unit mass, designated by the symbol $k$), which changes between 1,803 m²/s² and 16,262 m²/s². This variable represents the intensity of the modeled turbulence in the plasma ²⁷³. This turbulence quantity is higher in the intermediate part of the channel compared to the inlet and the outlet.

Figure 20 displays the spatial distribution of the turbulence dissipation or the turbulence dissipation rate (designated by the symbol $ϵ$ or the name “epsilon”), which changes between 2.05 × 10⁵ m²/s³ and 7.2 × 10⁷ m²/s³. This variable describes the loss of turbulent motion that is converted into thermal internal heat due to viscous effects ²⁷⁴. This turbulence quantity is higher in the intermediate part of the channel than in the inlet and outlet. This was the case for the turbulent kinetic energy per unit mass, and it is as reasonable behavior that as the turbulence intensity increases, the turbulence dissipation also increases ²⁷⁵.

Figure 21 and Figure 22 display the spatial distribution of the axial and downward components of the electric-current density, respectively. Thus, the variables displayed in these figures are $J_x$ and $-J_y$, respectively. The axial component changes between –0.008 A/cm² and 10.5 A/cm². Thus, the axial electric-current density is mostly from the channel inlet to the channel outlet. The downward component changes between –0.310 A/cm² and 21.9 A/cm². Thus, the vertical electric-current density is mostly from the top anode to the bottom cathode.

5. CONCLUSIONS

In the current study, we presented a mathematical and computational model for simulating linear magnetohydrodynamic (MHD) Faraday channels with continuous electrodes.

The computational fluid dynamics (CFD) computer language OpenFOAM is used to develop an MHD solver which we applied to the Sakhalin pulsed solid-propellant plasma (SPP) power system, with a supersonic divergent channel. The proposed solver corresponds to the low magnetic Reynolds number (Re_m), and it utilizes the scalar electric potential as the working variable for handling the electric aspect of the problem.

After validation and verification, we used the solver to explore various aspects of the thermal, fluidic, and electric features of the problem. These include the pressure, temperature, density, velocity, Mach number, electric-current density, specific turbulent kinetic energy, and turbulence dissipation.

Despite earlier studies about this Sakhalin pulsed magnetohydrodynamic generator (PMHDG), we do not know of any previous that that provided the same comprehensive aspects as our study, and we do not know of any previous study that presented the OpenFOAM MHD solver we propose.

The following items can be considered as contributions made in this study:

• A review of real examples of pulsed MHG generators (PMHD) was made, with emphasis on the Sakhalin model that was described with quantitative values and empirical regression models for its operational settings that help in modeling it by others. A large number of literature research articles were used in preparing this review.
• A mathematical/computational model is proposed for simulating the three-dimensional fluid dynamic and electric fields for linear continuous-electrode MHD channels. The proposed solver was found to yield satisfactory outputs when compared to other independent studies. The ability of the solver to give resolution-independent results was verified.
• Selected results from our CFD simulation for the Sakhalin MHD channel were presented and discussed. Using the power of CFD and data post-processing, it was possible to reveal various attributes of the channel of the Sakhalin PMHDG, such as the nearly-isothermal operation, the extraction of energy mainly from the kinetic energy of the plasma (rather than the thermal or the pressure energies), the unusual deceleration (opposite to the conventional situation) of supersonic flow despite flowing in a divergent channel, the weak change in the electric conductivity, and the pressure recovery phenomenon.

The MHD solver presented here can be expanded in subsequent studies to make it even more useful as a generic tool for designing MHD channels for open-cycle MHD generators, without being customized to the Sakhalin channel. One way to improve the applicability of the presented model is to replace the prescribed regression function for the local plasma electric conductivity with a submodel that estimates the local plasma electric conductivity for an arbitrary thermal-equilibrium plasma medium, given its (1) chemical composition, (2) temperature, and (3) pressure.