Electromagnetohydrodynamic Effects on Steam Bubble Formation in Vertical Heated Upward Flow

In this paper, the modeling of a steady state two phase flow heated through a vertical upward flow under electro-magneto-hydro-dynamic forces is presented. The thermal non-equilibrium, non-homogeneous, two-phase flow model consisting of mass, momentum and energy conservation in each phase has been adjusted for subcooled inlet conditions close to saturation. The P-1 approximation, viscous dissipation and Joule heating are included in the energy equations. It was seen that the Lorentz force can decrease and postpone the bubble generation, as well as affect the slip velocity, flow forces, viscous dissipation and Joule heating. Furthermore, two correlations for the slip velocity under magnetohydrodynamic (MHD) forces are presented. As shown, skin friction and Joule heating increase with the magnetic field strength.


Introduction
Two-phase flow in complex systems, such as blood flow in the body [1][2][3][4], bubbly flows [5] and nano-bubble mixtures through a microchannel [6], are in the interest of scientist and engineers.Specially, the two-phase flow of liquid water and vapor liquid occurring in a pipe is important in many applications, such as fluid engineering mechanics, multiphase heat exchanger design [5], solar heating systems [7], special heat exchangers [8], flow separation devices in central heating systems [9], the MHD generator [10], the MHD propulsor [11], etc.Based on the electrical conductivity of the water and vapor, they can be controlled and measured by electric and magnetic fields.An electromagnetic flow meter can measure air-water accurately at low void fractions (<30%).
The problem of two-phase magnetohydrodynamic flow and heat transfer in a vertical channel is of interest for many researchers experimentally [12] and numerically [13].The fluids in the two phases are usually assumed to be immiscible, incompressible and electrically conducting.Further, the two fluids are assumed to have different viscosities, thermal and electrical conductivities.The transport properties of the two fluids were taken to be a function of pressure and temperature.A review by Branover of two-phase MHD energy conversion systems used in the Ericson cycle mentioned that the kinetic energy of the liquid flow is converted to electric power via a near-isothermal gas expansion process [14].A study by Inoue et al. [15] illustrated that the ratio of the two-phase heat transfer coefficient with a magnetic field with that without a magnetic field increased almost proportionally with the magnetic field strength.There has been some research on the sub-cooled boiling water at the low-pressure condition in the horizontal [16] and vertical boiling channel [7][8][9][10][11][12][13][14][15][16][17].Furthermore, there is research on the computational fluid dynamics(CFD) simulation of the upward sub-cooled boiling flow of other fluids, such as refrigerant-113, in the literature [18,19].The problem of two-phase magnetohydrodynamic flow and heat transfer in a horizontal channel has been investigated analytically [13].As shown, the effect of increasing the Hartmann number is to accelerate the velocity and to increase the temperature field.There have also been numerical [20][21][22] and experimental [23] studies on the Hartmann flow of a conducting fluid in the channel between two horizontal insulating plates showing that significant increases can be obtained in the flow-rate of the conducting fluid for suitable ratios of the depths and viscosities of the two fluids.The pressure head of the micropump in a 2-mm diameter tube was measured under the applied electric and magnetic field (up to 60 Volts and 0.44 Tesla) until the bubbles generated by the electrolysis of the conducting liquid were measured [23].Furthermore, the effect of slag layers on the heat transfer characteristics of a coal-fired MHD generator at fixed fluid/gas ratios has been studied analytically [10,24].Shoukri et al. [25][26][27][28] by applying 10 4 V to R134a through the 1-cm diameter tube with the mass fluxes of the order 10 2 kg/m 2 •s, it was found that when the Masuda number (or dielectric Rayleigh number) is of the same order of magnitude as the square of the liquid Reynolds number, hence electrohydrodynamics (EHD) interfacial forces are important in the pressure drop.The use of electromagnetic force on the two-phase flows has been used to affect the onset of bubble generation [29], the momentum transfer between the phases [30], the control of the other properties of flows and the creation of a different flow structure [31].
Motivated by the above-mentioned literature, the purpose of the present study was to develop a numerical tool with a non-equilibrium model to analyze and predict the steady-state two-phase flow boiling, which will remain accurate (in comparison with the experimental results) for a large range of sub-cooling inlet conditions.To do this, first the proper prediction of the void fraction profile with experimental results was performed.Then, we investigated the effect of electromagnetic forces on the thermo-physical properties of flows.

Problem Definition
In general, the one-dimensional codes consume less computational time to perform a simulation as compared to 3D approaches.Consequently, they are preferred to investigate two-phase system behavior when there are no swirling flows, separated flows and flows due to complex geometries where averaging is inappropriate.Figure 1 shows the sectional view of the vertical tube and its boundary condition.This problem is the same as the Case 3 experiment from [29], where the vertical concentric annular is a part of the low pressure water boiling loop.The fluid entered the channel with a mass flux of 283.1 kg/m 2 •s and 19.7 • C sub-cooling.The outer tube had a 25.4-mm inner diameter, while the outside diameter of the inner tube was 12.7 mm with a length of 114.6 cm.The inner tube of the middle section of the pipe, from 34 cm from the beginning to 64.6 cm, was heated uniformly.The other surfaces of the channel were adiabatic.Therefore, the vapor formed in the boiling section with a 478.5 kW/m 2 wall heat flux was condensed in the adiabatic downstream section.The surfaces were impermeable, and there was unidirectional flow in the y-direction between the impermeable boundaries.The pressure at the outlet was held constant at a value of 1.23 bar.At y = 0, the inlet temperature (T in ) and the velocity have been assumed constant (U in ), and the velocity and temperature are assumed as the outflow boundary condition at the outlet.Furthermore, in order to align the Lorentz force J × B through the flow line, the electrical field E should be perpendicular to the magnetic field B. The amount of the current density corresponding to the motion of the charge continuum is calculated from the ohm law ( J = σ E ( E + V × B)).The various electromagnetic conditions are presented in Table 1.

Mathematical Model
When the length of the pipe is much larger than the diameter, the one-dimensional method to analyze the fluid flow leads to acceptable results [29].The general case of the mathematical formulation of the proposed model involves 14 unknowns: seven variables of the flow (pressure, velocity, temperature of each phase and the void fraction), three interfacial transfer terms (mass, momentum and energy) and four transfer terms at the wall assumed impermeable (transfer of the momentum and the energy of each phase at the wall).By taking into account the assumption of 1D two-phase flow, the local balance equations averaged over the cross-section of the pipe lead to a system of six partial differential equations.The complete formulation of the model has up-to-now been quasi and not feasible due to the limited and practical knowledge of the transfer processes in multiphase systems.Some approximations must be made in order to reduce the size of the problem.

Conservation of Mass: Continuity Equation
While the flow does not include external mass sources or sinks, there is the mass transfer due to evaporation or condensation between phases throughout the process.Hence, there is liquid generation, which is equal to the negative of the vapor generation.For a two-phase fluid, the conservation of mass for each phase is satisfied if [32]: where the α i for the gas is the common α and for the fluid is 1 − α.This equation is also called the continuity equation.The left term represents the rate of mass of each phase, and the right term represents the integral of interfacial mass transfer from the inlet to the current location.Here, the steady state solution is required.However, for the numerical method, the transient formulation is applied until the variation in the variables approaches zero.The mass generation at the vapor phase can occur in the bulk of the fluid or near the walls.For the simplicity of the problem, just the vapor generation at the bulk of the fluid is assumed, so [33,34]: The wall vapor condensation rate can be modeled by considering the bubbles' coalescence and breakup and the turbulence effects [35,36], but because of the smallness of the hydraulic diameter and for the sake of simplicity, here it is assumed equal to:

Conservation of Momentum
The conservation of momentum equations for each phase are written under these assumptions: • The Reynolds stresses, viscous stresses and co-variance terms are neglected.
• The surface tension forces are negligible.Therefore, the same pressure for both of the phases and interfacial are assumed.• The interfacial momentum storage is neglected.
• The interface force terms consist of both pressure and viscous stresses.
• The constant area is assumed.
• The virtual mass formulation is used (for the more stability in the numerical solution, the continuity equation at each phase is multiplied by a percentage of the corresponding phase velocity, and the consequential equation is subtracted from the related momentum equation).
The momentum equation for the each phase is [37]: In some references, the total interfacial force considered in the present study includes the effects of the drag force, lift force, wall lubrication force and turbulent dispersion force [38].However, in the present study, a more simplified force derived from [32][33][34][35][36][37][38][39] was used for the wall friction coefficient.The interfacial drag force (between phases) is [40]: The inertial force term in the interfacial momentum exchange model has been combined into a single expressions, which is continuous for all flow regimes.Thus, no flow regime maps are required [41].
In wall forces terms, the boundary layer constitutes the major resistance to heat transfer, and a variety of external forces, centrifugal, vibrational and electrostatic, have been used to reduce its thickness and to promote drainage.

Conservation of Energy
The Reynolds heat flux, the covariance terms, the interfacial energy storage and the internal heat transfer are neglected in the conservation of energy equations.The thermal energy equation for each phase is [32]: The interphase energy transfer, the Joule heating, the required pumping power to conquer the gravitational force and the heat loss of viscous dissipation are the components of the volumetric heating of the energy equation whose calculation is straight forward.However, the calculation of the radiative term at the simplest assumptions (surface-to-surface radiation model) needs a two-dimensional domain.Thermal radiative heat transfer is important in a wide range of applications from molecular dynamics [42] to macro-scale problems, such as boundary layer flow [43] and channel flows [44].The problem of radiative transfer in a gray medium (or on a spectral basis) is determining the radiative intensity with three space coordinate variables and two direction coordinate variables.The method of spherical harmonics is an approximate solution of a truncating of a series of simultaneous partial differential equations.The highest value, in P-1 approximation, is one (the approximations of odd order are never used).The water layer with an order of centimeter thickness is not the optically-thick medium (high amount of photon absorbing), for which the Rosseland model could be applied.Therefore, in this research, the P-1 approximation for the radiation model is assumed for the transport of thermal radiation heat flux through the material [45].
The heat sources or sinks due to radiation are calculated using the equation: where: When the wall emissivity is equal to one, the flux of the radiation, q r,w , at the walls, caused by the incident radiation, H w , is given as: For the wall temperature condition of the energy equation, it should be noted that the heat flux at the wall has different components (evaporation, quenching and convection).The wall temperature can be found from the heat balance at the wall [5].The interfacial energy exchange rate represents the rate of energy transferred from one phase to the other.This transfer can be due to either conduction, which is a function of the temperature distribution, or mass transfer.If one considers the interface to be infinitesimally thin and at saturated conditions, then the energy transfer can be modeled.Defining the energy transfer as positive when the vapor receives the energy, the energy transfer rate may be written as in Tu et al. [33,38] q" = 0.52nd 3 In sub-cooled and saturated boiling conditions, the interface-to-vapor energy transfer can be appropriately modeled by considering the vapor to be at saturated conditions.In order to maintain the vapor at saturated conditions when the bulk liquid is sub-cooled, a relatively high rate of vapor-to-interface heat transfer is required.The force fields, such as magnetic forces, electrostatic fields and ultrasonic vibrations, could be used to improve the critical heat flux and increase boiling heat transfer rates.

Numerical Procedure
In this study, the system of Equations ( 1)-( 8) is solved by the power law-implicit scheme in finite-volume formulation with a staggered grid for the mass, energy and momentum.The above set of equations is solved with a semi-implicit numerical method as discussed in [46].The conservation equations are first approximated by a linearized set of finite difference equations.The density and internal energy in these equations are then eliminated in favor of the pressure and temperatures by using the equations of state.The momentum equations are then used to derive a relationship between velocity and pressure, which can be used to eliminate the velocity in the mass and energy equations.As a result of the difficulty of the development of the gas volume fraction and pressure of the system and velocity and internal energy distributions of the each phase, the resolution of the local and instantaneous balance equations remains very difficult.These operators introduce some additional transfer terms related to the constitutive laws and to the boundary conditions.The iteration process is terminated when the following condition is satisfied: where φ stands for either P, α, U g , U g , v g and v f and m denotes the iteration step.This method converges on the pressure distribution with the user controlling the convergence criteria and the number of iterations (both inner and Newton).If convergence is not attained in the specified number of iterations, the code has the capability of reducing the mesh size.The advantage of this feature is that if the mesh size is small enough, the method will always converge.Hence, this solution method is extremely reliable for any type of flow conditions.At each time step, the Q r is evaluated by a two-dimensional model with the known temperatures of the previous iteration.The index of the refraction of water and the characteristic penetration depth of the collimated radiation in the water, which is the inverse of the spectral absorption coefficient of water (l = 1/α) in the wavelength range from 0.2 µm to 7 µm, are shown in Figure 2 [47][48][49][50].The absorption and scattering of the radiation by spherical particles like bubbles that are greater than the radiation wavelength and MHD electromagnetic wavelength can be calculated by the Mie theory.Thus, in the current case, the method published in [51,52] was used to calculate the effective absorption and scattering coefficient of radiation for P-1 approximation in the system of semi-transparent water containing vapor spheres.

Results and Discussion
As stated in the problem definition, steady state one-dimensional case experiment results performed on the geometry of the vertical concentric annular tube (Figure 1) by Shoukri [29] where sub-cooled water is entering upward were compared to the results of the numerical calculations of the model presented here.This experiment simulated a low pressurized tube with water that was heated at the middle section of the tube.The void fraction predictions of the code are illustrated in Figure 3a.As can be seen, the experiments show that the boiling occurred earlier, and the code accurately predicted sub-cooled boiling.Additionally, the void fraction distributions in Figure 3b are significantly different, which is also due to the differences in the applied MHD conditions in Table 1.
Figure 3a shows the comparison of the void fraction distribution results obtained by numerical modeling with the experimental results of the Shoukri experiment.Figure 3a shows the predicted gas void fraction profile in the y direction (continuous line) and the experimental data (squares).As seen, the values of void fraction obtained by the numerical simulation are in good agreement with the experimental results data, and the maximum difference between the present study and experimental study of Shoukri et al. is about 5%. Figure 3b examines the effect of various parameters and models on the numerical prediction.Figure 3 shows the variation of the void fraction along the y direction.The void fraction prediction is the most important characteristic of a low pressure boiling problem.For example, all of the numerical predictions of the [29] are about the void fractions.(b) effects of the electric field and the magnetic field (see Table 1).
After an initial adiabatic axial distance where distilled water passed through the vertical concentric annular test section, bubbles are generated in the middle section because of the exposure to the heat from the electrically-heated inner tube.This heated section was followed by another adiabatic section.Accordingly, uniformly, heat was generated in the middle section where the sub-cooled boiling took place.As shown, the density of bubbles was maximum near the end of the middle section walls due to the convection phenomena, and the bubbles continued to exist to the end of channel.The peak local void fraction was observed in the vicinity of the heated surface because of the large number of bubbles generated from the active nucleation sites on the heated surface.The vapor formed in the boiling section was condensed in the adiabatic downstream section.With increasing the length of the channel void fraction, the the void fraction is decreased, due to the diffusion of the bubbles.Because the volume flow was less than 0.3 (α ≤ 0.3), the bubbly-slug transition did not happened.As shown in Figure 3b, increasing the positive Lorentz force caused the void fraction to decrease.Although approximately the same heat flux conditions were maintained along the middle axial (along the length) section of the pipe, the void profile, as shown in Figure 3b, indicates a significant drop.This void drop is due to the pressurization effect of the Lorentz force.Furthermore, the bubble generation drop along the length of the pipe was small, but it had a significant effect on the mixture properties.One effect of an electromagnetohydrodynamics field applied across a vapor-liquid interface is to produce a force on the liquid phase tending to accelerate the fluid towards the solid surface.This results in the vapor film in the film boiling region being destabilized and, hence, increased wetting of the heater surface.Heat transfer coefficients in the nucleate boiling region are also increased by the electrical field.The profiles of the friction coefficient and friction as shown in Figures 4 and 5 clearly indicate a significant variation along the length of the pipe.It can be noticed that the friction of water has less sensitivity to the Lorentz force than the gas component.By increasing the Lorentz force, the gas friction coefficient dropped from 1.6 to about a 0.1 value before exiting the pipe (see Figure 4b).Generally, the friction coefficient profile of both phases increases along the pipe length after the heating region, as illustrated in Figure 4.It is clear that when the sub-cooled water at the pipe entrance partially changes to vapor, after undergoing the wall heat flux, the generated bubbles stay in the saturation state until exiting the pipe, and both velocities increase.Furthermore, by increasing the positive Lorentz force, the wall friction increases at both phases.The pressure profile along the pipe in each case is not shown in the figures, because those profiles are approximately linear (except the steep slope at the end of pipe).The net forces to fluid flow through the pipe for Cases A-G are 0.01, 1.82, 3.64, 18.2, 54.6, 127.4 and 362.544N, respectively.Therefore, the net force on the fluid dramatically increased.The drift flux model can be used with or without reference to any particular flow regime.The drift flux approach, therefore, satisfactorily accounts for the influence of mass velocity on the void fraction, as seen in the separated flow model, to provide the required relationship between the void fraction and the independent flow parameters.As stated by Wallis [53] and used in some numerical justification of the experimental results [29], the following equation is used for slip velocity in drift flux models: where the value 1.18 was reported as 1.41 [54] and 1.53 [29] in other references.However, using the Equation (11) in the current study produced great errors.As the buoyancy force plays an important role in the void fraction development, the MHD forces can affect the body force on both phases.Therefore, the following correlation: is suggested to give the proper estimation for calculating slip velocities with MHD flows at low Lorentz forces (see Cases A-D in Figure 6).For Cases E-G, which have the high Lorentz forces, another correlation is suggested: Slip ratio correlations, which are applicable to the current regime, which is beyond the bubbly flow, are simple and useful to describe the balance between two phase momentum.If the engineer uses the slip ratio on the simple annular flow theory or homogeneous theory, the calculations leads to approximately equal friction pressure gradients.As shown in Figure 6, while at low Lorentz forces (Cases A-D), increasing the positive Lorentz forces increases the slip velocity, but the opposite trend is observed at high Lorentz forces (Cases E-G).In the current study, instead of using the mixture momentum equation, the distinct equations were used for each phase because of the important effect of the Lorentz force on the fluid flow.Therefore, the calculated slip velocity is precise and can be used to evaluate Equation (11).In Figure 6, the gas velocity is calculated by the use of Equations ( 12)- (13) and is presented by a prime sign on each case letter.As revealed, the suggested Equations ( 12)- (13) give an acceptable value for the gas velocities in MHD regimes.This is because both liquid-extraction and electro-convection phenomena generate a secondary fluid motion in the liquid that increases as the axial momentum decreases.The prediction of the slip ratio by using Smith's method as a function of local vapor quality leads to great values of errors, as the local vapor quality is small.The most important non-dimensional parameter of the flow, i.e., the Reynolds number, is plotted in Figure 7.The Reynolds number of the gas phase as a function of axial location and Lorentz forces is illustrated in Figure 7.As demonstrated, an increase in the Lorentz forces increases the gas phase Reynolds number.Because of the effects of MHD and drag forces, the vapor velocity increases rapidly along the channel.Furthermore, the maximum Reynolds numbers occurred at the end of the heated region.By increasing the positive Lorentz forces, the vapor phase Reynolds number at low Lorentz forces (Cases A-D) increases and decreases at high Lorentz forces (Cases E-G).Furthermore, all cases take place in the laminar regime.After the heated section, the vapor phase Reynolds number increases, while the liquid phase Reynolds number remains constant.The liquid phase Reynolds number shows the same trend with lesser variation between various cases and fluctuates about the 10 4 .Therefore, ignoring the turbulence contribution [38] is justified, and it is not necessary to equip the formulation with a turbulence model.Furthermore, the maximum vapor phase Reynolds number occurred between Cases C and D. It is recommended that the pressure drop in bubbly flow be calculated using the separated flow model with the values of the void fraction calculated from the slip relationships.The wall shear force in bubbly flow is usually small, and the frictional pressure gradient can be evaluated using a drift flux model with little error.The liquid MHD forces are approximately constant throughout the pipe, but the volumetric vapor MHD forces differ, as revealed in Figure 8. Generally, increasing the Lorentz forces liquid MHD forces to increase (A-F), except in Case G, where they decrease because of the lower void fraction value.The liquid and vapor viscous energy dissipations are plotted in Figure 9.For the liquid viscous energy dissipation, there is less variation.Increasing the positive Lorentz forces increases the liquid viscous energy dissipation for the low Lorentz forces (Cases A-D) and decreases them for the high Lorentz forces (Cases E-G).On the other hand, the vapor viscous energy dissipation at low Lorentz forces is one order of magnitude lesser than the high Lorentz forces and in Case E has the maximum heat loss.The reason is that by decreasing the ratio of gas to liquid, the electromagnetohydrodynamic force affects more parts of the liquid, and the resultant perturbation has more effects in rising heat transfer.The incident radiation for Case G and the volumetric contribution of thermal radiation (W/m 3 ) along the pipeline for various MHD forces are plotted in the Figure 10.As shown, the incident radiation cannot completely be considered one-dimensional, and its maximum occurs around the end of the heated region.In the two adiabatic sections, the temperatures of the fluid and gas are approximately equal, and so, the volumetric contribution of thermal radiation is not significant in those regions.

Conclusions
A numerical model of a steady state two-phase flow heated through a vertical upward flow under electromagnetohydrodynamic forces was performed to predict the effect of MHD on the onset of nucleate boiling; in addition, bubble generation and heat and the momentum characteristics of the sub-cooled boiling were developed.The model is based on a nonequilibrium and nonhomogenous two-phase model.The simulated results are compared to the experiments, and good agreement was observed.Then, the effect of the electrohydrodynamic force on fluid flow was simulated.The following conclusions were drawn from the application of an electromagnetic field: • To get a correct understanding of the mechanism of heat transfer enhancement in the present investigation, the system was studied in a steady state condition.• Positive Lorentz force (in the flow direction) decreases and postpones the bubble generation.
• Negative Lorentz force (opposite of the flow direction) increases the bubble generation.
• Increasing of the positive Lorentz force increases the wall friction coefficient and wall friction force in both phases.
• By the increase of the electromagnetodynamic field, The net force on the fluid dramatically increased.• Slip velocity changes such that it cannot be predicted by the previous relationships, and the new correlations for the slip velocity under MHD forces are presented (Equations ( 12)-( 13)).The drift flux model is valuable only when the drift velocity is significant compared to the total volumetric flux.This limits its usefulness to the bubbly, slug and churn flow patterns.• The vapor phase Reynolds number decreases.This leads to diminishing of the velocity in the churn-turbulent region of bubbly flow.• By the increase of the electric and magnetic field, the vapor phase MHD force and, hence, the viscous dissipation in the fluid and liquid phases increases.• By the increase of the liquid phase and postponement of the bubble generation, the density of the high absorbing medium increases, and the incident radiation and volumetric contribution of thermal radiation increases.radiation heat flu (J/m 2 •s) q" density of heat transfer at the wall (J/m 2 •s)

Figure 1 .
Figure 1.Sectional view of the vertical tube and boundary condition.

Figure 2 .
Figure 2. Index of refraction and penetration depth of water.

Figure 3 .
Figure 3. Void fraction (a) comparison of the current numerical study with the experiment [29];(b) effects of the electric field and the magnetic field (see Table1).

Table 1 .
Various electric field and magnetic field values.
Author Contributions: Mohammad Yaghoub Abdollahzadeh Jamalabadi designed the research with theoretical formulations.Mojtaba Mirzaee, Payam Hooshmand, Hamed Ahmadi, Hassan Kavoosi Balotaki, and HamidReza KhakRah together have critical contribution in revising the manuscript including the preparation of responses and English correction.All authors have read and approved the final manuscript.The authors declare no conflict of interest.