Single-bubble EHD behavior into water two-phase flow under electric-field stress and gravitational acceleration using PFM

In this study, single-bubble electro-hydrodynamic effects on the two-phase laminar flow of water under electric field stress are investigated using numerical modeling. A 2D axisymmetric model is also developed to study the growth and departure of a single bubble. The phase-field method is applied to track the interphase between liquid and gas. The growth of the attached vapor bubble nucleus to a superheat at 7.0 °C and 8.5 °C are evaluated with 50° and 90° contact angles. The results show that the enhancement of the contact angle changes the velocity and temperature fields around the bubble. It is observed that the growing size and base of the bubble is increased with increasing the wall superheat, but the bubble departure diameter and time are decreased. The electric field results in raising the number of detached bubbles from the superheat at a certain time interval but decreasing the bubbles departure size. Additionally, the formation of stretched bubbles enhances the rate of heat flux and there is a non-linear relationship between the applied voltage and heat flux.


INTRODUCTION
Boiling heat transfer, like nucleate boiling, is one of the efficient types of heat transfer 1 . High heat fluxes can be achieved at low superheats in the nucleate boiling. The thermal function of heat exchangers can be improved by using active and passive heat transfer enhancement techniques 2 . The electric field is one of the active methods to improve heat transfer 3 . In recent years, the use of the electric field in biphasic flows has been increased. The study of the dynamics of electrically charged fluids is known as electrohydrodynamics (EHD). In EHD, the fluid flow equations are evaluated under the influence of an electric field 4 . The electric field can be used to study the hydrodynamics of the bubble and understand the dynamic behavior of bubble deformation and decomposition. It is important for its application in many industries and for the synthesis of materials. Many researchers have studied the effect of the electric field on bubble dynamics [5][6][7][8][9][10] . Tomar et al. simulated the heat and mass transfer rates in the film's boiling region and under the influence of a uniform electric field. A combination of the volume of fluid method (VOF) and the level set method (LSM) was used for the simulation. It was found that mass and heat transfer increases with increasing electric field intensity 11 . Gambhire and Thaokar compared the oscillations of the interphase of the two fluids under the influence of a non-uniform electric field for two semiconductor and conductive models. The results showed that controlling interphase instability was better when the electric field was oscillating. It was also shown that oscillating electric fields could control the instabilities at the interphase by increasing the frequency of the applied voltage, for the steady-state system at zero frequency 12 . Ishimoto et al. studied the effect of the nonuniform magnetic field on the behavior of bubbles in magnetic fluid experimentally. It was concluded that it is important to consider the behavior of the bubbles in such cases in many engineering applications such as modern energy conversion systems, but it is difficult to see the bubbles in such currents 13 . The effect of a uniform magnetic field on the air and vapor bubbles in a magnetic field was experimentally investigated. It was found that applying a magnetic field to the bubble flow of the thermal pipes can increase the heat transfer rate and allows it to be controlled. However, it was reported that the observation of the air and vapor bubbles in the liquid under the influence of the magnetic field is difficult 14 . The boiling of a single bubble in a fluid of water under normal and low gravity was examined and it was concluded that the lifetime of the bubble could be adjusted based on the water temperature 15 . As can be seen, many research works have been done experimentally [16][17][18][19] and it has seen a number of limitations and difficulties. Therefore, it will be highly valuable to investigate these systems using modeling and simulation methods 20,21 . For this reason, this research study was conducted using modeling and simulation.
In present work, comprehensive numerical modeling and simulation using the phase-field method were developed to study single-bubble electro-hydrodynamic behavior in the two-phase flow of water under electric field stress. Also, the bubble hydrodynamic was evaluated without electric field stress or gravitational acceleration to see how the electric field can affect the bubble hydrodynamic. The bubble motion and its deformation were studied at operating conditions of laminar flow, incompressible, and unsteady two-phase fluid flow under gravitational and electric fields. Furthermore, a 2D axisymmetric model was developed to investigate the growth and departure of a single bubble. The developed model was verified with experimental data in the absence of an electric field in terms of bubble departure diameter. The results of the current research study can provide valuable insights into relevant experimental research in the field of aerospace research and terrestrial processes.

RESULTS AND DISCUSSION
To validate the model, the process was simulated without an electric-field and obtained results were compared with the experimental data from the published literature 22 . Then, the 1 Isfahan University of Technology, Department of Chemical Engineering, 84156-83111 Isfahan, Iran. 2 Institute of Research and Development, Duy Tan University, Da Nang 550000, Viet Nam. 3 The Faculty of Environment and Chemical Engineering, Duy Tan University, Da Nang 550000, Viet Nam. ✉ email: mahdighadiri@duytan.edu.vn effect of the electric-field on the desired parameters was investigated. Modeling results were consistent with the results of experimental research on the evaporation of water in the presence of an electric-field. Figure 1A and B compare the bubble growth time for the experimental data and modeling values when the fluid comes into contact with the superheat with a temperature of 7°C and 8.5°C 22 . As can be observed, there is a great agreement between the experimental data and modeling values. Also, the modeling results are more accurate compared to the previous numerical work 22 . From Fig. 1A, the highest error for the bubble growth time was obtained at 11%. In terms of the bubble departure diameter, the maximum error was found to be 2.6%. However, in the superheat with the temperature of 8.5°C, the error values were slightly higher and those were calculated 16 and 8% for the growth time and the departure diameter, respectively. It can be seen in Fig. 1A and B, the obtained diameter in the simulation is less than the experimental and previous numerical values.
Furthermore, in Fig. 1A and B, mesh size sensitivity analysis was investigated. As can be seen, the change in the meshing from 9632 to 11086 does not affect the results and the error is less than 2%. More information about meshing can be found in Supplementary Section 2 (Supplementary Table 2 and Fig. 2).
The change in the contact angle as a function of the interphase velocity for a single bubble in the current work and the experimental data 23 from the literature was shown in Fig. 2 (A). From Fig. 2 (A), there is a good agreement between modeling values and the experimental data. Furthermore, when the shape of the bubbles is compared with the empirical results, it can be seen that the shapes of the two bubbles match well for a superheat with temperature and contact angle of 8.5°C and 50°( Fig. 2 (B) 22 . Also, velocity and temperature fields were given in Supplementary Section 4 ( Supplementary Fig. 4).
Electrohydrodynamic force can change the shape of the bubble. This force is the divergence of the Maxwell stress tensor. From Eqs. (1) to (4) 24 , the Maxwell stress tensor consists of two positive and negative terms.  In these equations, E is the electric field, B is the magnetic field, S is a closed surface, ε is the absolute permittivity of space between charges, ε 0 is vacuum permittivity, and ε r is the relative permittivity or dielectric constant of the medium. Figure 3 shows the first, the second, and complete terms of the radial component of the Maxwell stress tensor in the sections of (A), (B), and (C), respectively. Negative values indicate the force to the left and positive values indicate the force to the right. As shown in section (A) of Fig. 3A, the force from the first term is applied more to the middle part of the bubble and it causes a transition of the middle part of the bubble from the liquid bulk, to the central axis of the bubble. The second part of the Maxwell stress tensor is zero, except in the small regions of the top and bottom corners of the bubble (shown in Fig. 3B). This force direction is to the right, but it is smaller than the first term. Based on Fig. 3C, there is a force in the middle of the vapor bubble in the left direction but the force , the x-axis is the bubble diameter and y-axis is the bubble height and the unit of the x-axis and y-axis is meter. The unit of rainbow colorbars is N (newton). Arrows in G and H are radial and vertical EHD forces respectively. direction on the top and bottom of the bubble is upward and downward respectively. Furthermore, the amount of force in the left direction is dominant in the system, therefore, the size of the bubble was decreased with increasing the amount of voltage. From Fig. 3D and F, it was observed that the vertical forces are mostly upward and it can lead to the acceleration of bubble departure from the surface. On the other hand, these forces in the inner corners of the bubble in the upper and lower parts also cause stretching of the bubble shape. The direction of radial and forces were provided in Fig. 3G and H for a clear understanding of effects of different terms in the Maxwell stress tensor.
The effect of electric field on the growing bubble base and height was shown in Fig. 4A and B respectively. In Fig. 4A, the empty symbols on the x-axis are the bubble departure time at different applied voltage. From Fig. 4 (A), the bubble base increased and reached a peak as a function of time, then, it was seen decreasing in the bubble base until its departure from the surface. In the absence of electric field, the maximum bubble base of 1 mm was obtained at time of 0.019 s and the bubble departure from the surface happens at 0.04 s. Applying 4000 V electric field deceased the bubble departure time to 0.022 s. There was an increase in the maximum base from 1 mm to 1.18 mm and 1.20 mm when 1000 V and 2000 V voltages were applied respectively but the amount of time needed for reaching the maximum base was not changed. Moreover, the maximum base was found to be 1 mm when 3000 V voltage was applied but time was decreased from 0.019 s to 0.016 s. Applying 4000 V voltage decreased the maximum value of base and time to 0.64 mm and 0.016 s, respectively. The empty symbols on the horizontal X-axis are the bubble departure time. In addition, the bubble diameter was decreased by applying 3000 V and 4000 V voltage. As it was explained in Fig. 3H, the direction of the vertical forces is upward and it resulted in the decrease in the bubble departure time from the surface (Fig. 4A). Furthermore, the upward and downward vertical forces in the bubble internal corners at the top and bottom of the bubble led to the increase in the growing bubble height when the bubble diameter was the same. There was a bit increase in the bubble height at voltage of 1000 V. The same trend was reported by Zu and Zhang works 25,26 . But, with decreasing the bubble size at a higher voltage value (Fig. 4B), it can be seen a bit decrease in the bubble height at the voltage of 3000 V, and the bubble height was decreased 0.5 mm by applying 4000 V voltage in comparison with the system without any electric field.
The bubble departure diameter can be affected by different forces including buoyancy force, surface tension, and applied electric force. Therefore, electric force can change the bubble size. Electrohydrodynamic force effect on the bubble diameter at the time of departure for the superheat with temperatures of 7°C and 8.5°C was presented in Fig. 4C. The bubble departure diameter was decreased from 2.9 mm to 1.0 mm and 3.1 mm to 1.7 mm when the voltage was increased from 0 to 4000 V for the superheat with temperatures of 7°C and 8.5°C. As can be seen, the higher electric force was required to significantly decrease the departure diameter of the bubble. This means that a larger inward force is needed to overcome the expansion pressure and the volumetric growth of the bubble which enters from the vapor phase to the liquid phase. The decrease in the bubble diameter could be attributed to the existence of radial electrohydrodynamic force in the middle of the bubble (Fig. 3G).  Figure 4D shows the change in the bubble velocity as a function of the applied voltage. In the present study, the bubble velocity was measured at 1 cm from the bottom of the container. It was observed that the bubble velocity was increased from 250 mm/s to about 360 mm/s when the applied voltage was increased from 0 to 4000 V. This result is consistent with the experimental work of Kweon and Kim 27 .
After the bubble growth and its departure from the hot surface, the cooler saturated liquid fills the position of the separated bubble from the hot surface, then, a certain time is required for heating of liquid and nucleation of the bubble. The growth time (t g ) and waiting time (t w ) can be used for the calculation of the average bubble frequency (f) for a bubble using the following equation: 27 The change in growth time, waiting time, and the average bubble frequency as a function of the applied voltage is shown in Fig. 5A. As can be seen, the growth time is always higher than the waiting time for all voltages and the case without electrical force. Kweon et al. 27 reported the same behavior for water but based on the fluid type it is possible to have a waiting time higher than the growth time. There was a decrease in the growth and waiting times with increasing applied voltage and subsequently the amount of the average frequency was increased with the enhancement of applied voltage. It is clear based on Eq. (5) and the same trend was reported by Kweon et al. 27 Moreover, it was observed that the waiting time was not changed when the applied voltage was increased from 2000 V to 4000 V. It is because there is a direct relationship between waiting time and thermal boundary layer. The voltage higher than 2000 V was not able to change the thermal boundary layer. Therefore, it was remained constant between 2000 V and 4000 V.
The fluid velocity is zero on the surface of the hot wall, therefore, there are two mechanisms for the heat flux through the wall include conductive heat flux and the heat flux due to phase change from liquid to vapor and it is maximum at the triple point on the wall. The heat flux can change with the change in the vapor nucleation, growth, and its departure from the surface. The change in heat flux as a function of time (900 ms) for the superheat with the temperature of 7°C (9 th bubbles) and 8.5°C (20 th bubbles) was shown in Fig. 5B. The vapor phase thermal conductivity coefficient is too lower than the liquid phase. Therefore, thermal conductivity was too much lower where the vapor was formed. The heat flux was decreased during bubble growing and increasing of bubble base, but it was then increased with decreasing bubble base and beginning of the bubble departure step. It can be seen a slight decrease in heat flux when the process goes from one bubble to the next one. In addition, it was observed that there are small peaks between large peaks when the superheat temperature is 7.5°C. It is because the waiting time is longer for this case and vapor nucleation occurs after a certain time interval when a bubble departure from the surface.
Using dimensionless number is more appropriate to compare different samples or parameters. Nusselt dimensionless number is  Nu ¼ l 0 q kΔT sup (6) where q and k are heat flux and thermal conductive coefficient. The difference between Nusselt numbers at different applying voltage is lower than the difference between heat flux. It is because ΔT sup is the denominator in the fraction (Fig. 5C). As it can be seen, the distance between the Nusselt diagrams is less than the distance between the graphs of the two superheated fluxes (Fig. 5B), because although the heat flux is higher at 8.5°C the lower ΔT sup at 7°C reduces the difference between the two Nusselt number values. Figure 6A and B show the effects of electric field and time on the heat flux passing through hot surfaces with temperatures of 7°C and 8.5°C during a certain time, respectively. The contact angle was 50 degrees. It was observed that increasing applying voltage increased the heat transfer rate and the number of bubbles which departure the surface. For example, the number of bubbles was increased from 3 to 9 when the applied voltage was enhanced from 0 to 4000 V for the case with a superheat temperature of 7°C. As it was explained, applying voltage decreases growth and waiting times. Therefore, the heat transfer rate through superheats increases with increasing of the average bubble frequency and bubble departure. Also, Electrohydrodynamic force can change streamlines close to the superheat and bubble. All these changes can affect the temperature field in the system and subsequently improve the heat transfer rate. In addition, it was found that increasing applying voltage was decreased the thermal boundary layer and this is also increased the heat flux rate.
From the comparison of Fig. 6A (7.0°C) and 6B (8.5°C), it was found that increase in the heat flux is higher for the system without any electric field as the heat flux increased from 26000 W/ m 2 to 38000 W/m 2 while when 4000 V voltage was applied there was an only 7000 W/m 2 increase in the heat flux (from 48000 W/ m 2 to 55000 W/m 2 ). Therefore, the impact of superheat temperature was higher in the absence of an electric field.
The Nusselt number as a function of time and voltage is shown in Fig. 6C. The maximum Nusselt number was increased from 16 to about 29 with increasing applying voltage from 0 to 4000 V.
The distribution of electrical potential depends on the shape and the position of the vapor bubbles in the computational domain. Figure 7 shows the distribution of dimensionless electrical potential in the computational domain. Electrical permeability is different for the liquid and gas phases which is clearly shown in Fig. 7. The difference in the electrical permeability of liquid and vapor phases could be attributed to the different electrical volume force gradian in each phase. In fact, the bubble acts as a barrier for the passing of the magnetic field due to its lower magnetic permeability. On the other hand, the fluid flow tends to remove this barrier for the passing of the electrical field. The applied force is only effective on the interphase of two phases. Therefore, it can be called magnetic surface tension. The electrical field is a straight line for one phase. But, it can be seen a change in the electrical field, when it penetrates from the bubble to the liquid phase.
Further results and discussion are provided in Supplementary  Section 4 (Supplementary Figs. 4 to 8)

Governing equations
The system geometry and fluid properties were provided in Supplementary Section 1 (Supplementary Fig. 1 and Table 1). The main equations for the model building are conservation equations for the mass and momentum in the incompressible state. The velocity and pressure fields for the liquid phase were modeled by the Navier-Stokes equations, according to Eqs. (7) and (8) 28 .
∇:u l ¼ 0 where ρ l , u l , η l are the fluid density, fluid velocity, and fluid dynamic viscosity, respectively, and l subtitle represents the liquid phase. For the vapor phase, a weak form of Navier Stokes equations was used, according to Eqs. (9) and (10) 28 : where u and η are the fluid bulk velocity and the dynamic viscosity of the fluid bulk, respectively, u v is the fluid velocity of the vapor phase, and ρ v is the fluid density of the vapor phase. The energy equation was solved only for the vapor phase based on Eq. (11) 29 : where C p is the heat capacity of the vapor phase and k v is the heat transfer coefficient of the vapor phase. The interphase temperature of the vaporliquid was considered to be equal to the saturation temperature, therefore, the heat transfer equation was solved only for the vapor phase. The Poisson equation (Eq. (12)) was also solved as the electric field equation.
This function is Poisson's equation for the dielectric materials 30 .

Domain equations
In the phase-field method, instead of direct tracking of the interphase between the two phases, the intermediate layer is obtained by the phasefield variable. It is assumed that the state of the system is described at any time by the phase-field variable ϕ, which is a function of the position vector. The Cahn-Hilliard equation is used for describing the dynamics of the interphase in two-phase flow as follows (shown in Eq. (13)) 31 : where G (Pa) is the chemical potential and γ (m 3 .s/kg) is the mobility parameter. The free energy is a function of the dimensionless parameter of the phase-field is presented based on Eq. (14): where ε is the value of the interfacial thickness and f tot (J/m 3 ) is the total free energy density of the system. The free energy density relating to a mixture of isotherm consisting of two insoluble fluids is a summation of mixing energy and elastic energy. The mixing energy was estimated by the Ginzburg-Landau Eq. (15): 31 The dimensionless parameter of the phase-field ; is determined in such a way that the volume fraction of the fluid components is ð1 þ ;Þ=2 and ð1 À ;Þ=2. The symbol λ N ð Þ is the mixing energy density. This parameter is related to Eq. (13 in SIF file), the surface tension coefficient, σ N=m ð Þ, and the thickness of the interphase. The degree of mobility determines the time scale of the diffusion of the Cahn-Hilliard, and it should be large enough to have a constant interfacial thickness, and also it should be small enough to avoid the over-damping of the convective terms. The degree of mobility is related to the thickness of the interphase using the mobility tuning parameter χ m:s=kg ð Þwhich is determined according to γ ¼ χε 2 . The value of the chemical potential is obtained by the following Eq. (16) 32 : Hence, the right side of Eq. (12) is for the minimization of the total free energy using the residence time, which this time was controlled by motion parameter γ m 3 s=kg ð Þ . The Cahn-Hilliard function forces the parameter ;, except in very narrow regions of the fluid interphase, to accept two values Fig. 7 Distribution of dimensionless electric potential for the superheat with temperature of 8.5°C and contact angle of 50°. The unit of x-axis and y-axis is meter. of 1 or −1 and by breaking the fourth-order equation into the following two second-order equations Eqs. (17) and (18): The term of free energy can sometimes include other sources. These sources can be considered by modifying the Eq. (18) to the following equation: where f ext is the determined free energy by the user (J/m 3 ). The external free energy is mostly zero 31 .

Phase-field method in solution of two-phase problem
To solve the two-phase problem of the boiling process with the phase-field method (PFM), firstly, the boundary and initial conditions and fluid properties in both phases and the motion parameter of the phase-field model must be determined. Then, the variables of the laminar flow model of this method include the dimensionless parameter of phase (;), and the auxiliary variable ψ in the computational domain must be initialized. In the next step, the transient analysis type is selected so that the main solution of the model is performed with all variables of flow equation including (p) pressure, the dimensionless variable of phase (;), auxiliary variable (ψ), (u) and (ν) velocities, energy equation including temperature variable (T) and electrostatic equation including variable of potential difference (V). In this step, a number of terms is added to the equations as a source to consider the effect of phase change. The following terms are added to the equations of momentum, energy, and phase, Eqs. (20), (21) and (22) respectively: ∂; ∂t þ u: where δ is the length of the interphase between two phases based on the Eq. (23): In this way, the variables of the model are obtained over time, and then the physical and thermal properties of the fluid are calculated. The volume fraction amount of vapor and liquid is calculated using the phase-field variable given in Eq. (24): Min and max operators were used to measure the volume fraction between the lower and upper limits of 0 and 1. Thus, mixed properties such as specific heat density, thermal conductivity coefficient, and electrical permeability are obtained using the volumetric component and the following Eq. (25): where A can be any of the properties mentioned, and the subtitles l and g are related to the liquid and gas phases, respectively. Surface tension force was applied as a volume force, shown in Eq. (26): where G is obtained by Eq. (16). The value of mean curvature is also obtained by Eq. (27)

Boundary conditions
The boundary conditions for the boiling model are complex. The interphase velocity of the two phases is not necessarily equal to the velocity of the liquid phase or the velocity of the vapor phase, shown in Eq. (28): In Eq. (28), n is the unit normal vector and its direction is from the liquid phase to the vapor phase. Also, _ m is the rate of vapor production. The boundary conditions of the vapor phase in the interphase were considered as follows (based on Eq. (29)): n:ρ v :u v ¼ _ m 1 À ρ v ρ l þ n:ρ v :u l ð Þ (29) If no phase change occurs, u v ¼ u l ¼ u int , and the velocity of the solution when passing through an interphase of the two phases is continuous. In addition, if the density of the vapor and liquid phases is equal, the first term on the right side of the Eq. (29) will be zero and the velocity of the two phases is equal. When the phase change occurs and the density of the two phases is not equal, the first term on the right side of the equation has a value and a stream flows into the interphase of the two phases. The second term on the right gives the velocity in the outward direction of the interphase, therefore discontinuity is created in the field of velocity passing through the interphase. The velocity of the liquid in the outer direction of the vapor phase is greater than the vapor phase. The mass flux that leaves the interphase can be written as follows: 34 where M w is the molecular mass of the vapor and ΔH l:v is the vapor enthalpy. This equation is approximated by ignoring the kinetic energy and work done by viscous forces. In the interphase of the two phases, three forces affect the liquid phase, so the boundary conditions resulting from the balance of these forces for the liquid phase can be expressed as follows: 35 n: Àp l I þ η l ∇u l þ ∇u l ð Þ T h i ¼ _ m u l À u v ð Þþn: The first term on the right, the equation is related to the force that separates the vapor from the liquid. This force is always negligible. The second term is the sum of compressive and viscosity forces which is applied to the liquid phase from the vapor phase. The mass flux that leaves