Multi-Physics Mathematical Model of Weakly-Ionized Plasma Flows

This work presents a multidisciplinary mathematical model, as a set of coupled governing equations and auxiliary relations describing the fluid-flow, thermal, and electric fields of partially-ionized plasma with low magnetic Reynolds numbers. The model is generic enough to handle three-dimensionality, Hall effect, compressibility, and variability of fluid, thermal, and electric properties of the plasma. The model can be of interest to computational modelers aiming to build a solver that quantitatively assesses direct extraction of electric energy from a plasma flow. Three different approaches are proposed to solve numerically for the electric fields with different levels of tolerance toward possible numerical instability encountered at a large Hall parameter, where the effective conductivity tensor loses diagonal dominance and becomes close to singular. A submodel for calculating the local electric properties of the plasma is presented in detail and is applied to demonstrate the effect of different factors on the electric conductivity, including the fuel’s carbon/hydrogen ratio and the alkaline seed element that acts as the ionizing species. An analytical expression for the collision cross-section for argon is developed, such that this noble gas can be included as one of the gaseous species comprising the plasma.


Introduction
Ionization is an endothermic process in which a sufficient amount of energy exceeding the ionization energy of an atom or molecule is supplied to it, leading to the liberation of an electron. First ionization refers to the removal of the first electron. The removal of the second electron is called second ionization. We here consider only first ionization because it requires the least amount of energy and it is much more encountered in practical and engineering applications.
Ionization is classified into 3 main types based on the source of the supplied energy to overcome the ionization energy. The first method is electron ionization, where the gas is bombarded with an energetic electron beam, stripping off other electrons from the atoms or molecules upon collision. This type is common in mass spectrometry [1]. The second type is ionizing radiation, where sufficiently-powerful phonons (e.g., in the ultraviolet or X-ray radiation spectrum) overcome the binding energy of electrons in atoms or molecules and free them. This type finds uses in the medical field through diagnostic and therapeutic applications [2]. The third type is thermal ionization, which is due to the collision of thermally-agitated atoms or molecules, constituting a hot gas. When these atoms or molecules are very hot, their random-motion velocity becomes high enough so that the kinetic energy transferred in a collision between two atoms or molecules is sufficient to ionize one of them. This type occurs in electric arcs, high-temperature flames, and behind strong shock waves in hypersonic flow fields [3]. We here consider thermal ionization, which is relevant to the application we are concerned with, namely extracting electric energy from weakly-ionized hot-gas plasma in a magnetohydrodynamic (MHD) generator channel. This plasma is formed when a small amount (e.g., 1% by mole) of an alkali metal vapor, especially cesium or potassium [4], is present in the combustion products. This alkali metal can be introduced into the combustor as an alkaline salt compound added to the reactants. For example, potassium carbonate (K 2 CO 3 ) can be used for potassium [5]. Alkali metals have low ionization energies, making them much more ionizable than conventional combustion gases (namely CO 2 and H 2 O), which in turn makes them practically the only source of electrons in the hot-gas plasma.
MHD generators are thermoelectric devices that convert chemical energy in the form of a fuel into electricity. The fuel is burned in a combustor using a suitable oxidizer and seeded alkali metal compound. The electric conductivity of alkaliseeded slightly-ionized plasma is strongly dependent on the absolute temperature [4], taking roughly the following formula [6], Eq. (1): σ (1) where and are coefficients that depend on the plasma gas composition. Thus, a high temperature is critical to the success of the MHD concept. This favors the use of oxygen as the oxidizer rather than air to avoid the dilution effect of nitrogen or other recirculated flue gases (RFG) [7]. The plasma then flows downstream the channel which is subject to an externally applied lateral magnetic field of magnitude (B) and has a top electrode (anode) with low voltage and a bottom electrode (cathode) with a high voltage. These voltages can be controlled. For example, the low voltage is zero if the anode is grounded; then the high voltage will be at its largest value if the external electric circuitry is open, but it will be at its smallest value if the external electric circuitry is shorted. Optionally, the combustor exit is fitted with a convergent nozzle ending with a throat upstream of the channel. This makes the flow in the channel supersonic, which critically intensifies the electric output from the channel. In that case, a diffuser (convergent passage) is fitted downstream of the channel to reduce the high velocity of plasma. Figure 1. elucidates the aforementioned components of an MHD generator. The interest in MHD electric power is not new, and it goes back to the early 1970s, where the former USSR has built different MHD generator units as portable, short-term (pulse mode), high-output electric sources for area and deep electrical geophysical surveying [8,9]. Later, the USA has launched an experimental proof-of-concept (POC) [10] study to explore the technical and economical aspects of coal-fired MHD power generation for retrofitting a plant of 300 MWth (thermal input power). The program ran during 1987 to 1993. Aside from the experimental studies, computational simulations have been performed to predict the performance of potential designs for MHD channels. These studies were contemporary with that phase of interest in MHD generators, and thus were based on simplified mathematical modeling (e.g., steady problem with no time variations [11], axisymmetric [12], one dimensional [13], two-dimensional [14], or three-dimensional but with two-dimensional electric fields and parabolized form of Navier-Stokes equations [15] which eliminates the second-order derivatives in order to allow marching in the axial direction in an technique similar to the boundary-layer techniques [16]) of the plasma motion and electric fields. Recent interest in MHD generation has aroused due to different factors, including stronger magnets which boost the output power from MHD channels, enhanced combustor designs that help lengthening the electrode life, and the advance in computational fluid dynamics (CFD) enabling more reliable prediction of the performance of MHD units [17]. It is the aim of this work to provide a complete mathematical model that describes the temporal and spatial evolution of plasma-gas fields and electric fields within the channel, such that computational researchers can pursue further and use this model in numerical simulations (e.g., based on finite volume discretization) to predict the performance of arbitrary MHD channel designs. The model consists of 3 conservation equations plus the equation of state for the 4 flow fields (density, pressure, stagnation enthalpy, and the velocity vector), and 3 equations for the electric fields (electric potential, current density vector, and electrostatic field vector). Thus, there is a total of 13 unknown scalar fields described by 13 equations. The model is augmented with a submodel to calculate the plasma electric properties, such as the electric conductivity. We used this submodel to analyze multiple factors affecting the electric conductivity. Three different computational approaches are discussed to resolve the electric fields, coping with poor numerical behavior at plasma conditions characterized by a large Hall parameter.

Elementary Electromagnetic Equations
Before 1864, the laws of electromagnetism were stated in integral forms until James Maxwell reformulated them in differential forms. Maxwell equations constitute of Gauss' law of electricity (Gauss' law), Gauss' law of magnetism (solenoidal law of the magnetic field), Faraday's law of induction, and Ampère-Maxwell equation. Gauss' law, the solenoidal law of a magnetic field, and Ampère-Maxwell equation are not very instrumental in our work. In the following, we present the remaining law among Maxwell equations, namely Faraday's law, plus the generalized Ohm's law, and combine both with the charge conservation law to devise the equations governing the electric fields in an MHD generator channel. This will be manipulated further in section 0.

Faraday's Law
Verbally, Faraday's law mentions that if a conductor (plasma in our case) moves within a magnetic field, where a closed loop is moving with the plasma, an electromotive force will be generated that is equal to minus the rate of change of the magnetic flux through the loop. In an integral form, this law has the form in Eq. (2). The differential form follows from applying Stokes theorem and is given in Eq.
(3). Equation (2) implies that an emf is induced if the applied magnetic field is time-dependent. It should be noted that if the magnetic field is not time-dependent but the conductor is moving, then an emf is still induced in the closed path because there is still a change in the magnetic flux crossing perpendicularly the area spanned by the loop.
The electric field in Eq. (3) can be decomposed into two fields, namely the electrostatic field & and the unsteady electric field ' induced due to any unsteadiness of the magnetic field [18]. The first electric field is irrotational, by Coulomb's law. Only the second electric is tied to the time rate of change in the magnetic field. Because the magnetic field is not changing in the MHD channel, we are only concerned about the electrostatic field, and we have: The above equation means that an electric potential function ) can be introduced such that There is another electric field induced due to the motion of the plasma (the conductor) with velocity * within the applied magnetic field ! . This motional electric field is equal to + * × ! . This can be shown from Eq. (2), and it is the working principle for the turbo-generators. When applying Eq. (3) for that case, one should note that the change in the magnetic field is with respect to a moving observer traveling with the plasma, thus the closed loop in Eq. (2) is moving rather than being stationary. Moreover, the drift of the electrons (as a conductor) within the plasma also induces another electric field, manifesting the Hall effect. So, in general we have The second term is zero for our applications because we consider only the case a time-fixed (steady) magnetic field.
Theoretically, a magnetic field (/ ) due to the currents within the plasma is induced but it is safely neglected here. This is a very reasonable assumption for an MHD channel, called low magnetic Reynolds number, and is described in Appendix A.

Generalized Ohm's Law
The (simple) Ohm's law for a conducting medium subject to an electrostatic field simply establishes a linear relation between the current density within the medium and that electrostatic field, thus It is to be noted that the current density vector is collinear with the electrostatic field and has the same direction. However, in the presence of other induced electric fields as a result of any unsteadiness of the magnetic field, the motion of the plasma under the effect of an applied magnetic field, or the drift of the electrons within the plasma under the effect of an applied magnetic field, the law is generalized as: The three terms in the RHS represent electrostatic current density, Faraday current density, and Hall current density, respectively. The derivation of Eq. (8) is given in Appendix B.

Charge Conservation Law
Because the charge is a conserved quantity, the net flux of 0 leaving out of an arbitrary closed surface " enclosing a volume V is equal to the time rate of decrease of the charge contained in that volume. In its integral form, the charge conservation law is expressed as Using Gauss' theorem, the differential form of this law is For our purposes, the RHS of Eq. (10) is zero since the plasma is macroscopically neutral (in general, this term is negligible in conductors [18]), and we have

Derived Electromagnetic Equations
In this section, we manipulate the elementary electromagnetic equations presented in the previous section to derive partial differential equations that are more suitable for implementation computationally, by augmenting them to an arbitrary fluid-dynamic solver.

Generalized Ohm's Law with Tonsorial Conductivity
The generalized Ohm's law, Eq. (8), can be expressed in a more compact form that is explicit in the current density (i.e., 0 appears only in the LHS). To achieve this, we recall that the vector product 10 × ! 2 can be expressed as the inner product of a skew-symmetric tensor ! and 0 (as a first-rank tensor).
The tensor ! has the following form: where B E , B D , B C are the Cartesian components of the magnetic field vector ! . Using this vector rule into the generalized Ohm's law, Eq. (8), and slightly manipulating the terms, we obtain where J is the identity tensor. Denoting the coefficient tensor in the LHS as L , and multiplying the above equation by the inverse of this tensor, we obtain or where the tensor 4NN represents an effective electric conductivity, which is calculated as Using the symbolic manipulation package SymPy [19] (a Python [20] library for symbolic mathematics), we determined the analytic expression of 4NN , and this is given in Appendix C.

Poisson Equation for Electric Potential O
Using the definition of the electric potential, Eq. (5), into the modified generalized Ohm's law, Eq. (15), and then applying the divergence operator to both sides, and recalling that the divergence of the LHS vanishes by the charge conservation law, Eq. (11), we obtain an elliptic Poisson-type partial differential equation for the electric potential )

Calculating the Electric Properties
In order to handle the equations described in the previous section, it is apparent that we need a submodel to calculate the local electric conductivity and the electron mobility 3 4 as a function of other local plasma parameters, such as the temperature and chemical composition. We propose a submodel based on the equilibrium ionization model of Saha [21,22], which predicts the number density of free electrons in the plasma at a specified equilibrium temperature. Saha equation for alkali metals gives that, under equilibrium ionization, we have the following expression for the equilibrium constant .X (18) where S 4 is the number density of the free electrons after ionization (per m 3 ), S &,NTU is the number density of seed alkaline atoms after ionization (per m 3 ), and V #W;W is a lumped constant which has the value of 2.4247×10 21 (in m -3 ·K -1.5 ). The ionization energy for cesium (Cs) and potassium (K) are 3.8939 eV and 4.3407 eV, respectively [23]. From the conservation of alkaline nuclei, we have where S &,TUT is the number density of seed alkaline atoms before ionization (per m 3 ), which is calculated as [24] where i &,TUT is the mole fraction of the alkaline atoms before ionization (dimensionless), TUT is the total (not partial) pressure of the gases before ionization (in Pa = N/m 2 ). In the Eq. (19), the number density of electrons S 4 is used instead of the number density if alkaline ions in the RHS, because both have the same values under single-ionization. Combining Eq. (18) and Eq. (19), we obtain a quadratic equation for the electron number density after ionization S 4 , whose only acceptable (positive) solution is Once this electron number density is known, the composition of the plasma is known. One can then use the collision cross sections of the neutral species in the plasma along with the number density of each of them and obtain the mean collision frequency m n due to these gaseous species according to Eq. (22) [25] m n = 4,+4WU ∑ S p q 4r where 4 is the mean of the magnitude of electron velocity based on Maxwell speed distribution (in m/s), S r is the number density (in 1/m 3 ) of the species u, q 4r is the average momentum-transfer cross section for electron collisions with the molecules of the species u (in m 2 ), ]^is Boltzmann constant, is the absolute temperature (in kelvins). The lumped constant V 4 has the value 6212 It should be noted that Eq. (22) is suitable for slightly-ionized plasma rather than fully-ionized plasma, as it does not account for electron-ion and electron-electron collisions.
We used analytical fits (based on experimental data) for the collision cross sections, as a function of the electron energy in eV (electron volt) made by Frost [26]. For Argon, we deduced a fitting function from the data in Ref. [27], because argon is not covered in the work of Frost. We present the details of these analytical fits in Appendix D.
The mean collision frequency due to neutrals only is augmented by the mean collision frequency due to ions (and electrons) scattering m (in Hz), to yield the total mean collision frequency m x . where In the above set of equations, † n is the electric permittivity of vacuum (8.854 × 10 -12 F/m), the lumped constant V~ has the value 7.727×10 -12 (when using SI units), and the lumped constant V Š has the value 1.549×10 13 (when using SI units). Once the m x frequency is calculated, we proceed with calculating the electron mobility as described later in the Nomenclature section, in the part devoted there to the electron mobility, but with m x replacing m 4 .
Comparisons with reference data [6] support the validity of our implementation of this submodel, which we carried out using the Python platform. Another comparison is made with a case of air at 3273 K with 2% potassium [28]. The cited value of the magnetic Reynolds number (Re m ) for this case is 1.3×10 -5 , under a reference velocity of 1 m/s and a reference length of 1 m. Taking the magnetic permeability of air to be the same as the one of vacuum, i.e., 4… × 10 M• H/m, (because there are very close to each other), and from the definition of Re m in Eq. (A-1), we infer that the corresponding electric conductivity for this case is 103 S/m. Our calculation estimates the electric conductivity for this case to be 192 S/m. Although the difference is not small, the agreement is thought to be favorable in the light of large uncertainty in the analytical fits, and the lack of details about the limitations and calculation method in the source reference. Relatively large deviations in across different sources are acceptably waived [4]. In our calculations, we approximated the air as a mixture of 21% O 2 and 79% N 2 , by mole. We also assumed that the pressure is atmospheric; the pressure value in the source reference is not disclosed. In fact, if the pressure in the reference source is high, then our calculation will approach the cited value. For example, if the pressure is twice the atmospheric value, our predicted electric conductivity drops to 158 S/m, and drops further to 139 S/m if the pressure is increased by another atm. We performed another validation against an experimental value of [29] for combustion-plasma of a paraffin fuel (kerosene) having a C/H ratio of 0.5 (1 carbon atom per 2 hydrogen atoms [30]) which is typical for transportation fuels [31], in oxygen at two temperatures with 1% potassium by mole. We also compared our to a computational value [4]. The comparison is summarized in Table 1, and our results fairly agree with the measurements. The submodel should be applied at each spatial point where the electric properties of the plasma are sought. The inputs to the submodel include the local volume fraction of the alkali metal vapor (either cesium or potassium) and the mole fractions of potential gaseous molecules which may appear in the combustion products. Due to the relatively high ionization energies of all species except the alkali metal, the free electrons are taken to be due to the ionization of the alkali metal only. However, the other species affect the electron mobility because they form inhibit a collisionless electron motion to different degrees, depending on their momentum cross sections. This submodel fits well the plasma conditions in an MHD channel, where the ionization is limited to a small fraction (≈1%) of the alkaline vapor, which is a small fraction (≈1%) of the plasma gas.

Results and Discussion
This section is divided into three main parts; the first and second parts are devoted to presenting the equations governing the thermo-fluid and electric field variables of the plasma, with proposed algorithms to solve for the electric fields. In the last part, we shed light on important factors that affect the electric conductivity of the plasma, through applying the electric submodel described in the previous section.

Governing Equations -Fluid-Flow Part
In this section, we put the partial differential equations describing the evolution of the velocity vector, temperature, pressure, and density of the plasma-gas [32], with suitable source terms to couple the electric fields to the fluid-flow fields.
There is no major restricting assumption about the fluid regime or dimensionality, and a generic compressible threedimensional flow with non-constant momentum and thermal properties is treated. Similarly, the governing equations for the electric fields do not stipulate a restriction on the dimensionality of these fields or a predefined spatial distribution for the electric conductivity or electron mobility.

Mass Conservation
The mass conservation equation of the plasma gas is not affected by the presence of electric fields. The equation preserves its differential form for a non-conducting fluid, being where : is the mass density of the plasma gas.

Momentum Conservation
The momentum conservation equation of the plasma gas derives from the original vectorial Navier-Stokes equation, but a source term is added representing the force (per unit volume) acting on the charged particles. Any translating conductor carrying a current while subject to a magnetic field will experience a force called the Lorentz force [18], perpendicular to both the motion and the magnetic field. Its origin comes from the force acting on an individual charged particle with a charge ' (in coulombs), where we have where ' = − for an electron, but ' = for an ion.
Applying this equation to the electrons and ions within a unit volume of the plasma, we obtain Given the large ion-to-electron mass ratio, the ion is considered to move with the same bulk velocity of the plasma-gas, thus * -xU = * . This allows us to re-write the above equation as But (from the Nomenclature section provided later, under the part devoted there to the electric-current density, 0 ) this means This is the Lorentz force per unit volume, which is a body force exerted on a unit volume of the plasma in the MHD channel due to the presence of electromagnetic fields. It is the means of coupling the electric field into the motion of the plasma. The modified momentum equation with the additional RHS source term is thus where is the pressure of the plasma gas, and š is the viscous stress tensor. We point out that in the limiting case of one-dimensional plasma velocity and magnetic field and electrostatic field (mutually orthogonal to each other), the Lorentz source term reduces to −|0|! . The minus sign of this term reveals the dissipative nature of this source, which is in fact a sink of plasma momentum, trying to decelerate the plasma and/or decrease its pressure as it progresses downstream the MHD channel.
As for any Newtonian fluid, the viscous stress tensor š is modeled using the following constitutive relation [33] where 3 is the dynamic viscosity of the plasma-gas, ˜* is the velocity gradient tensor, I˜* K is its transpose, and J is the identity tensor.

Energy Conservation
There are several ways to express the energy conservation of a fluid. We here use one form convenient to high-speed gases, which is demanded for MHD channels. This form utilizes the stagnation specific enthalpy oe as the dependent variable rather than the temperature. A nonlinear solver is needed to extract the temperature from oe , where we first compute the static specific enthalpy ℎ as ℎ = oe − * ⋅ * (33) and then solve the following nonlinear equation for the temperature The specific heat at constant pressure V h ( ) is customarily expressed as a piecewise polynomial of temperature [34], with a reference specific enthalpy ℎ ž4N at a reference temperature ž4N for each species.
The partial differential equation for the energy conservation is where ' is the diffusion heat flux vector, which is typically expressed using the following constitutive relation (Fourier's law in this case) [35]: where | is the thermal conductivity (in W/m·K). The MHD scalar source term in Eq. (35) expresses the power (per unit volume) extracted from the plasma and converted into electric power to an external load connected to the channel electrodes. This term is negative, acting actually as an energy sink (not source), as expected. In the limiting case of one-dimensional plasma velocity and magnetic field and electrostatic field (mutually orthogonal to each other), this source term reduces to −|0| .

Ideal Gas Law
With ¢ being the specific gas constant (in J/kg·K), the relation between the pressure, density, and temperature follows the following form of the ideal gas law Because combustion plasma is composed of more than one gaseous species, the ¢ value in the above equation is an average value, calculated by where ¦ ¥ is the molecular weight of the mixture, to be calculated as [36] ¦ ¥ = ∑ i T ¦ T T (39) and i T and ¦ T are the mole fraction and molecular weight of each constituent gaseous species in the plasma.

Governing Equations -Electric Part
We have already presented in sections 0 and 0 the elementary equations governing the electric-field variables (namely & , 0 , and )) within the MHD channel. Out of these 3 electric fields, & and 0 are mandatory because they form the MHD source terms in the fluid-flow momentum and energy equations. We recall that the magnetic field ! is part of the momentum source term, but its spatial profile should be a given input rather than an unknown to solve for. We propose below 3 different computational approaches to solve for & and 0 ; all of which solve for the scalar potential function ) first.

Solving the Electric Fields: Approach 1
This approach is the most-straightforward, where it solves the tensor-based Poisson equation for the potential scalar ). The equation is repeated below for convenience.
The LHS is discretized implicitly, whereas the source RHS is treated explicitly. When integrating Eq. (17) numerically, we can apply a Dirichlet (or first-type) boundary condition at the anode and cathode, fixing the value of ) at each electrode to the desired value, for example Thus, ) = 0 at the low-voltage anode (as if it is grounded), but ) = §¨x W. at the high-voltage cathode, and §¨x W. is the desired voltage difference across the channel which is also the voltage difference across the external load connected to the electrodes (if voltage drops in the external circuitry is neglected). Even if the anode is not grounded, we can still set its potential as zero because the absolute value of ) is not important, but its spatial variation is what matters.
The boundary condition at the other (electrically nonconducting) channel walls should be a Neumann (or secondtype) boundary condition of the form [37] © U = 1* × ! 2 ªx'U. ⋅ S « where S « is a unit normal pointing perpendicularly away from the boundary. Once the scalar electrical potential function ) is solved for, a discrete gradient operator is applied to obtain the electrostatic field using Eq. (5), and then the modified (with tensorial conductivity) generalized Ohm's law, Eq. (15), is applied to find the spatial distribution of the vector current density 0 at the current time instant.

Solving the Electric Fields: Approach 2
Whereas the 1 st approach can work smoothly, the numerical behavior may deteriorate if the electron mobility 3 4 (thus the Hall parameter ¬) becomes high, where in this case the 4NN becomes less positive-definite. As an alternate (2 nd ) approach, the LHS of Eq. (17) is kept unchanged, but the source RHS will be adapted and expressed in terms of ) and 0 rather than * and ! . To this end, the modified (with tensorial conductivity) generalized Ohm's law is manipulated, and re-written as, As in the 1 st approach, the LHS is treated implicitly, whereas the 2 RHS terms are treated explicitly. It should be mentioned that while the first RHS term is analytically zero by virtue of Eq. (11), it will evaluate to a non-zero value during the calculations, which partly drives the solution of the elliptic equation. Once the scalar electrical potential function ) is solved for, & and 0 are obtained using algebraic equations in the same manner described in the 1 st approach, i.e., using Eq. (5) for & and Eq. (15) for 0 .

Solving the Electric Fields: Approach 3
In the 1 st and 2 nd approach, the solution of ) was based on the tensorial 4NN , which exhibits a poor numerical character for implicit discretization as the electron mobility (thus the Hall parameter ¬ ) increases. In such case, 4NN loses the diagonal dominance and approaches singularity, with its condition number growing to high levels as can be seen from Eq. (C-5) in Appendix C. This may lead to a failure in the computation even with the remedy of the 2 nd approach. A 3 rd approach is presented here, which does not use 4NN altogether. Instead, the scalar is used, forming an alternate scalar-based Poisson equation for ). The derivation of this equation is similar to the one given for Eq. (43), but we start from a simplified version of the original generalized Ohm's law, Eq. (8), where the Hall term is omitted, i.e., Taking the divergence and using the analytical condition that ∇ ⋅ 0 = 0 and the relation & = −∇) , we obtain a simplified Poisson equation for ), having the form Going back to the simplified Ohm's law, Eq. (44), we rewrite it as Then, we take the divergence of both sides, And we use this expression to replace the RHS source in Eq. (45), and again express this source in terms of ) and 0 rather than * and ! , as done in the 2 nd approach. We now obtain the final simplified Poisson equation to be attempted in this 3 rd approach ∇ ⋅ 1σ ∇) -2 = ∇ ⋅ 0 + ∇ ⋅ 1σ ∇) -2 The tilde above ) is used just to emphasize that solving this equation gives a estimate of the electric potential (thus, denoted by ) -), because the Hall effect was not accounted for. The LHS is treated implicitly whereas the RHS is an explicit source. After this step, an estimated electrostatic field is obtained as Then, a series of stationary-iterative correction is performed for the current density 0 using this estimated -& through the scalar-based generalized Ohm's law in Eq. (50). This correction recovers the ignored Hall contribution to 0 . Thus, the following equation is solved iteratively; in every iteration, an updated (more-corrected) value of 0 in the LHS using the old value of 0 in the RHS.
After a pre-specified number of corrections (or a tolerancebased termination criterion, testing that ∇ ⋅ 0 = 0 is satisfied within an acceptable tolerance in the local computational cell), we end up with the final value of 0 . While this approach may be the most tolerant among the presented 3 approaches for the electric field, there is a lagged influence of the Hall effect on the electric field -& . However, as the problem reaches a steady-state condition, the error fades away.

Analysis of Electric Conductivity
The electrical conductivity is a factor affecting directly the amount of electricity that can be extracted from plasma [38]. Careful understanding of its dependence on plasma conditions is thus important. We implemented the electric submodel presented in section 4 using a Python code, and we here examine quantitatively the effect of the plasma-gas composition on the for a typical MHD channel condition having a temperature of 2800 K, a potassium volume fraction of 1%, and a total pressure of 3 atm (both being before the ionization). The estimated electric conductivity is compared for different gaseous species in Table 2 (the plasma-gas is composed only of that species for each case). This exposes the retarding effect of the molecules on the mobility of the free electrons, and thus on the electric conductivity. We augment also the fictitious case of zero collision cross section, which corresponds to the limiting value of . A base case of air (79% N 2 + 21% O 2 , by volume) is also added as a reference. The last two columns represent the estimated if only the electron-neutral collisions are considered (while the scattering by ions and electrons is completely ignored -corresponds to setting m = 0 in Eq. (24)), and the relative erroneous gain in due to this simplification. The species are ordered with an ascending order of . Argon is distinct in its low resistance for electron motion, giving a value very close to the upper bound (the ratio is 0.977). Special attention should be paid to the typical major gases in combustion products, namely CO 2 , H 2 O, and N 2 . CO 2 and N 2 have a similar collision effect, which is far less than the effect of H 2 O. With H 2 O, is only 1/3 of its value in the case of CO 2 . This means that the hydrocarbon fuel type has a notable impact on the MHD channel performance, because the larger carbon content, the larger plasma conductivity and the better power generation, and vice versa. The use of low C/H fuel should be avoided. The base case (air) lies midway in the table. For species with larger , the error in the simplified grows, with a limiting upper bound of 46.4% at this temperature and pressure. However, for the 3 typically major gases in combustion products, the error is limited to 13.1%. This alleviates the effect of the correction in combustion plasma especially that the uncertainties in the model itself are comparable to such error. To further examine the factors affecting , we show in Figure 1 its dependence on the temperature in the range 2200 -3400 K where the gas is air (same base case in Table  2), with either potassium or cesium used as the ionizing alkali metal. This temperature range covers the reasonable operating temperature for an MHD channel. The lower half of the range is relevant to air-fired MHD combustor, where the presence of large amount of diluting air-born N 2 reduces the temperature. On the other hand, the upper half of the temperature range is more relevant to oxygen-fired MHD combustor, where the absence of N 2 raises the combustion temperature. It is evident that the latter combustion technique greatly benefits from the large boost in due to the elevated temperatures. The amplification in (we refer to the corrected value, by default) as a result of using cesium instead of potassium is clear, and is a direct consequence of the lower ionization energy of the former. This amplification decreases from 3.15 at 3400 K to 1.63 at 2200 K. We also compare in Figure 1 the corrected and simplified (i.e., with and without the m correction, respectively). We make the same note made for the data presented for Table 2: the correction effect is pronounced at large . This can be attributed to the larger amount of free electrons, which intensifies the ion-electron and electron-electron scatterings, which are captured by the correction term, reducing the overpredicted simplified .
We also point out that the electric conductivity is very sensitive to the temperature, but the dependence is not accurately approximated as the exponential function as in Eq.
(1), which will under-predict . The presence of the power term in the Saha equation and the nonlinear collision effects complicate the overall dependence, but it is better approximated as a quadratic profile based on our curvefitting analysis.

Conclusion
We presented a set of partial differential equations for describing the temporal and spatial evolution of the threedimensional fluid-flow variables and electric variables of alkali-seeded combustion plasma subject to an applied magnetic field. Derivations and auxiliary expressions are given and a modified form of the generalized Ohm's law that is explicit in the current density vector was derived. The model is supported with a complete submodel for predicting the electric properties of the plasma locally, accounting for scattering due to neutrals, electrons, and ions. This submodel is a prerequisite for applying the presented mathematical model. The fluid-flow equations can be integrated numerically in space and time using one of the several available computational fluid dynamics (CFD) solvers, utilizing the finite volume or the finite difference method. To make the flow solver capable of handling the Faraday induction and Hall effect for a plasma gas, one should augment one source term into each of the momentum and energy equations, and solve for the electrostatic field and the current density. Different computational algorithms were proposed to solve for the electric fields, all of them solve a Poisson equation for the electric scalar potential, but the formulation and solution procedure vary to allow coping with computational instabilities at high electron motilities and Hall parameters; in which case the conductivity tensor of the Poisson equation becomes nearly skew-symmetric. The sensitivity of the electric conductivity, a key parameter in model, to different design choices was discussed in light of the implemented electric-properties submodel. Fuels which are richer in carbon than hydrogen lead to appreciably elevated electric conductivity, thus more electricity extraction. Correction due to electron-and-ion scattering cannot be ignored for conductivity values approximately above 40 S/m.
The mathematical model is a preliminary step toward building a computational electro-fluid dynamics (CEFD) solver to act as a tool for predicting the performance of magnetohydrodynamic channels with the goal of maximizing the electric power output from the plasma through varying one or more of a multitude of geometric and operational parameters.

Plasma
A gaseous conductor (ionized gas) with a sufficient level of charged particles (ions and free electrons) to allow the passage of electric current and the response to magnetic fields upon motion.

Partiallyionized plasma
A class of plasma where only a fraction of the molecules are ionized. Alternative terms include weaklyionized plasma, and low-density plasma. Combustion-gas plasma, which is our interest, is an example. Thus, we focus here on this class of plasma. Fully-ionized plasma A class of plasma where the gases are completely ionized. Examples include solar winds, stellar interiors (e.g., the core of the Sun), and fusion plasma.

Equilibrium plasma
A class of plasma where the temperature of the free electrons is the same as the temperature of the ions. Thus the plasma gas is characterized by a single common temperature value. Due to their smaller mass, electrons are much more mobile than ions (the other current carrier particle in plasma). Thus, the energy of the ohmic heating (Joule heating) of the plasma gas is given to the electrons, which are normally maintained at essentially the same temperature as the ions and neutral particles as a result of the high collision frequencies and energy transfer per collision. Therefore, the free electrons and heavy gas particles (ions and neutrals) are in thermal equilibrium. An alternative term for this plasma class is thermal plasma, and quasi-equilibrium plasma. Examples include combustion-gas plasma [5] solar plasma, and arc discharge (electrical breakdown of a normally-nonconductive gas). We here focus on this class of plasma.

Non-Equilibrium plasma
A class of plasma where the temperature of the free electrons is higher than the temperature of the ions and neutrals. This situation occurs if ionizing current densities flowing through the plasma are sufficiently high, then the electron temperature will be significantly higher than that of the main body of the plasma gas. This is most easily attained in a monatomic gas with a large atom/electron mass ratio (typically argon), because in the absence of molecular vibrational excitation and the large disparity between the electron mass and heavy-particle mass, only a small fraction of the energy difference is exchanged on collision [39]. An alternative term for this plasma class is non-thermal plasma. Examples include aurora borealis, glow discharge (a plasma formed by the passage of electric current through a low-pressure gas), corona discharge (electrical discharge brought on by the ionization of a fluid surrounding a conductor that is electrically charged) [40].

Faraday field
This is an induced electric field (and current density) due to the motion of plasma gas while subject to an applied magnetic field.
Hall field This is another secondary induced electric field (and current density) due to the motion of the electrons (as a consequence of the primary Faraday current) as conductors while subject to the same applied magnetic field that induced the primary Faraday field.

Hall effect
The Hall effect describes the presence of the Hall field, leading to the inclination of the current density vector such that it is no longer collinear with the electrostatic field vector within the plasma.

J \
Ionization energy for an alkaline species (in eV, electron volts) Note: 1 joule is equivalent to 4 -.®n × n¯°± = 6.2415 × 10 s eV [42] 0 Electric-current density (or simply the current density); A/m 2 The current density is calculated here as product of the electron charge, electron number density, and electron drift velocity [28] The minus sign is a consequence of the negative charge of the electron. It is worth mentioning that the general form of Eq. (51) is where S TxU and * -xU are the number density and absolute velocity of the ions, respectively as the other charge-carrying particle in the plasma, with a positive charge. However, since the ions are much heavier than electrons, their velocity is taken to be same as the bulk gas velocity; i.e., * -xU = * . Using a velocity relation (which is explained at a later part in the current section) that * 4 = * + -. , the above equation can be written as [43,44] 0 = (S TxU * -xU − S 4 * − S 4 -. ) Because the plasma is macroscopically neutral [3], the number of electrons is equal to the number of ions, leading to S TxU = S 4 . These two assumptions lead to Eq. (51 Absolute (that is being relative to laboratory/fixed frame) velocity vector of the electrons; m/s -.
Drift velocity of the electrons (relative to the plasma medium); m/s In vacuum and in the absence of an electric field, the electron possesses a random motion with a zero average. In the presence of an electric field affecting an electron in vacuum, a net force acts on the electron accelerating it in the same direction of the electric field (e.g., toward the morepositive electrode in case of an electrostatic field). In the presence of an electric field in a material, the accelerated electron exhibits frequent collisions, leading to a finite average terminal velocity, called the drift velocity, which is used in the definition of the electron mobility. The drift velocity refers to an average velocity because actually the electron does not travel in a straight line, but its path is erratic (with curved segments between two collisions if an applied magnetic field is also present). The drift velocity is normally much smaller than the random velocity of the electrons, whose RMS value is 4,³+& k OE b c + a 6213 √ , where 4,³+& is in m/s and is in kelvins.
From the principles of kinematics [49], we find that * 4 * , -. In ionized gases (plasma) subject to a magnetic field, β is the ratio of the electron gyrofrequency (the angular frequency of the electron circular motion in the plane perpendicular to an applied magnetic field, also called the cyclotron frequency) ¸4 4 + a (in rad/s), to the mean electron collision frequency m 4 (in Hz). It is also the ratio of the electron meanfree-path 4 to the gyroradius (the radius of the circular motion of the electron ¹ 4 , also called the Larmor radius) [50]. This Hall parameter is also the product of the electron mobility 3 4 and the magnitude of the applied magnetic field !.
Moreover, for the simple configuration of axially-moving plasma moving in a horizontal constant-area channel with top anode (low-voltage electrode) and bottom cathode (highvoltage electrode) and applied magnetic field in the binormal direction (perpendicular to both the axial plasma motion and the vertical electrode-normal direction), then the induced electric current density vector 0 within the plasma will be vertically down if the Hall parameter is zero (practically very small) because of having a zero component in the axial direction. However, in general this current density will be inclined with an angle (Hall angle) º p from the vertical due to the presence of an axial component. This angle is the arctan of the Hall parameter ( Figure 3). This feature is mathematically explained in detail in Appendix C. Finally, the Hall parameter is also equal to the product of the gyrofrequency and the mean electron collision time š 4 » a .
With this, the Hall parameter has multiple equivalent mathematical expressions,  3 It is a medium property measuring the capacity of that medium to permit electric field lines. More specifically, it is the ratio of the generated electric displacement field (in C/m 2 ), manifested in slight separation of the positive and negative charges -atomic nuclei and their electrons -thus causing a local electric dipole moment, in response to an electric field (in V/m).
The electric permittivity of vacuum ε 0 (also called the permittivity of free space, and the electric constant) is a universal constant having the approximate value of (8.854 × 10 -12 F/m [51]). It is a material property measuring how quickly an electron moves through that material when pulled by an electric field. More specifically, it is the ratio of the average electron speed, the drift velocity, (in m/s) in response to an applied electric field (in V/m). The electron mobility is related to the magnitude of electron charge , electron mass ˆ4, and mean electron collision frequency m 4 according to the following expression [52]: 3 + :Magnetic permeability (or simply the permeability); henry/m (H/m) = (Wb/A)/m = kg·m/A 2 ·s 2 It is a medium property measuring the capability of that medium to be magnetized. More specifically, it is the ratio of the generated magnetic field (in Wb/m 2 ) in response to a driving magnetizing field (in A/m) produced by electric current flow in a coil of wire. The magnetic permeability of vacuum µ 0 (also called the permeability of free space, and the magnetic constant) is a universal constant having the exact value of (4π × 10 -7 H/m [53]).
The non-zero values of vacuum permittivity and vacuum permeability are consistent with the finiteness of the speed of light in vacuum c 0 , where: n‹ † n 3 n = 1. These non-zero values are attributed to the magnetization and the polarization of continuously appearing and disappearing fermion pairs, and the vacuum is hypothetically considered a medium filled with continuously appearing and disappearing charged fermion pairs [54].
: ; :Net electric-charge density due to ions and electrons (or simply the charge density); C/m 3 . For plasma, : ; = 0 because the plasma is globally neutral.
:Plasma electric conductivity (or simply the conductivity); (A/m 2 )/(V/m) = A/V·m = siemens/m (S/m) = 1/ohm·m (1/Ω·m) = A 2 -s 3 /kg-m 3 It is a material property measuring the capability of that material to conduct electricity. More specifically, it is the ratio of the generated current density in the material (in A/m) in response to a driving electric field (in V/m). The plasma electric conductivity is the product of the magnitude of electron charge , electron number density S 4 , and electron mobility 3 4 , = S 4 3 4 (57) m 4 :Mean frequency of the electron collision with other particles (neutrals, ions, electrons); 1/s = Hz relatively large mass of these heavy particles compared to the electrons. The last two forces combine together as a magnetic-field force due to the absolute electron velocity: -* 4 × ! Applying the dynamic force balance equation with safely neglected inertial term (due to the small mass of the electron [58]), we obtain: which is the generalized Ohm's law given in Eq. (8).

Appendix C: Effective-Conductivity Tensor (and Hall Angle)
We showed that the generalized Ohm's law can be expressed as We here aim at giving the analytical expressions for the elements of 4NN = σ L M , in terms of the electron mobility 3 4 and the Cartesian components of the applied magnetic field vector ! = È! É , ! Ê , ! Ë Ì .
The coefficient matrix L has the following elements: L = A 1 3 4 ! Ë −3 4 ! Ê −3 4 ! Ë 1 3 4 ! É 3 4 ! Ê −3 4 ! É 1 F (C-1) To simplify the expression for L M , let , ¶, · denote 3 4 ! É , 3 4 ! Ê , 3 4 ! Ë , respectively. Then We found that 4NN = σ L M has the following analytical expression: x + 1 x y − · x z + ¶ x y + · y + 1 y z − x z − ¶ y z + z + 1 We point that the elements of both L and L M are dimensionless. In the special, but not trivial, case of having the magnetic field acting purely in the lateral direction, ! = Ò0, 0, ! Ë Ó , Eq. (C-3) reduces to 4NN  This model case demonstrates the unfavorable effects of the Hall parameter, where it (1) causes a parasite Hall current density 0 É , and (2) reduces the magnitude of the useful Faraday current density à0 Ê à. It also explains the presence of the Hall angle º p for the current density vector 0 (measured from the vertical y-axis) in this case, as given in Figure 3, with º p = tan M (¬).

Appendix D: Analytical Fits for Collision Cross-Sections
In the following expressions, a prime in a symbol emphasizes that its unit is not according to the SI system of units. The electron-molecule collision cross-section for each gaseous species ( q &h4 T4& \ in cm 2 ) is expressed as a function of the electron energy * \ (in eV), calculated in turn as a function of the absolute temperature (in kelvins) as follows: is the mean magnitude of electron random velocity, expressed in cm/s. 1. For argon (Ar), we derived an eighth-degree polynomial fit given below q á³ \ = 10 M ® exp(a n â + a â + a â + a OE â OE + a l â l + a X â X + a ® â ® + a • â • + a s â s )