Cite this Article
Zhu Chang-Sheng, Han Dan, Xu Sheng. Phase field simulation of single bubble behavior under an electric field*

Project supported by the National Natural Science Foundation of China (Grant Nos. 51661020, 11504149, and 11364024), the Postdoctoral Science Foundation of China (Grant No. 2014M560371), and the Funds for Distinguished Young Scientists of Lanzhou University of Technology, China (Grant No. J201304).

. Chinese Physics B, 2018, 27(9): 094704  
Permissions
Phase field simulation of single bubble behavior under an electric field*

Project supported by the National Natural Science Foundation of China (Grant Nos. 51661020, 11504149, and 11364024), the Postdoctoral Science Foundation of China (Grant No. 2014M560371), and the Funds for Distinguished Young Scientists of Lanzhou University of Technology, China (Grant No. J201304).

Zhu Chang-Sheng1, 2, †, Han Dan1, Xu Sheng1
School of Computer and Communication, Lanzhou University of Technology, Lanzhou 730050, China
State Key Laboratory of Gansu Advanced Processing and Recycling of Non-Ferrous Metal, Lanzhou University of Technology, Lanzhou 730050, China

 

† Corresponding author. E-mail: zhucs_2008@163.com

Project supported by the National Natural Science Foundation of China (Grant Nos. 51661020, 11504149, and 11364024), the Postdoctoral Science Foundation of China (Grant No. 2014M560371), and the Funds for Distinguished Young Scientists of Lanzhou University of Technology, China (Grant No. J201304).

Abstract

Based on the Cahn–Hilliard phase field model, a three-dimensional multiple-field coupling model for simulating the motion characteristics of a rising bubble in a liquid is established in a gas–liquid two-phase flow. The gas–liquid interface motion is simulated by using a phase-field method, and the effect of the electric field intensity on bubble dynamics is studied without electric field, or with vertical electric field or horizontal electric field. Through the coupling effect of electric field and flow field, the deformation of a single rising bubble and the formation of wake vortices under the action of gravity and electric field force are studied in detail. The correctness of the results is verified by mass conservation, and the influences of different electric field directions and different voltages on the movement of bubbles in liquid are considered. The results show that the ratio of the length to axis is proportional to the strength of the electric field when the air bubble is stretched into an ellipsoid along the electric field line under the action of electrostatic gravity and surface tension. In addition, the bubble rising speed is affected by the electric field, the vertical electric field accelerates the bubble rise, and the horizontal direction slows it down.

PACS: 47.11.-j;47.55.dd
Keyword:electro-hydro dynamics;numerical simulation;phase field method;bubble
1. Introduction

In nature and engineering applications, bubble rising in a viscous fluid subjected to an imposed electric field is a common and abundant gas–liquid two-phase hydrodynamic phenomenon, playing a crucial role in boiling heat transfer and gas–liquid flow.[1–4] So far, the motion characteristics of the rising bubble without electric field have been extensively studied experimentally, theoretically, and numerically. A number of numerical methods have been proposed for complex gas–liquid interface motion deformation and topological changes. The existing methods can be divided into two major categories: interface tracking[5] and interface capture.[6]

The interface tracking method is based on Lagrange’s theory and uses dynamic grid technology to determine the phase interface by tracking the mesh on the interface. This method has more grid reconstruction work and a large amount of computation. Based on Euler’s theory, the interface capture method uses a fixed grid to calculate the twophase liquid, and a scalar indication function is used to identify and reconstruct the phase interface at each step in time. More mature methods of interface capture include the widespread use of volume of fluid (VOF) method, the level set method, and phase field methods. The VOF method[7] uses the volume of fluid to determine the free surface to track the fluid change, but the interface reconstruction is very complicated. The level set method[8] determines the fluid boundary location by defining a distance function, which can easily track the topological structure of the object. The phase field method[9,10] is the most promising method of solving multiphase flow problems. It uses the phase field variable to obtain the information about the interface layer and does not directly track the change of the interface between two fluids. The evolution of phase variable is controlled by the 4order partial differential Cahn–Hilliard equation. The surface tension is equivalent to the product of the gradient of the field variable and the chemical potential, and it is added into the Navier–Stokes equation as a volume force.

Recently, the movement of bubble under the action of electric field has gradually become a research focus. Marco et al.[11] studied the influences of the electric and gravity field on the separation and movement of bubble experimentally. It was found that the bubble did not fall off at low gas flow rate without electric field; at higher gas flows, even in the absence of buoyancy, the dynamic effect was enough to cause bubble separation. Peng et al.[12] solved the continuity equation of incompressible fluid by adding the momentum equation and VOF equation of electric field force, and obtained the process of detaching single bubble from the wall under the electric field. It was found that the bubble detachment time decreased in the electric field and elongated along the direction of the electric field. Andalib et al.[13] studied the effects of uniform and non-uniform DC electric fields on the motion and deformation of a single bubble in a dielectric liquid medium. It was found that in the case of thick bubble, the electro dynamic force is negligible compared with the high hydrodynamic instability, and the electro hydro dynamic (EHD) flow plays an important role in controlling the motion. Based on the assumption of peristaltic flow, Chen et al.[14] calculated the velocity field distribution inside and outside the bubble under the action of electric field, and compared it with the velocity field without electric field but with incoming flow.

In general, most of the previous researches focused largely on the effect of EHD on solid wall separated bubbles, but less attention was paid to the effect of the electric field intensity and direction on the rising bubble.

In this paper, a three-dimensional (3D) phase field model is established to study the dynamic characteristics of a single bubble in gravitational field. The floating process of bubble under the action of electric field is simulated by the coupling of flow field and electric field to study the shape and dynamic evolution of the rising bubble in the electric field.

2. Mathematical model
2.1. Control equation

The compressibility of the two phases is neglected, so the mass conservation of the whole simulation domain, including the bulk phase and its interface, can be expressed as follows:

The fluid motion is described by Navier–Stokes equations[15–17] for incompressible fluid under the action of gravitational force, electric field force, and surface tension. According to the Newton’s second law of motion, the momentum conservation equation is

where u is the velocity vector of the fluid, ρ (ϕ) is the fluid density, t is the time, p is the pressure, μ (ϕ) is the dynamic viscosity coefficient, ρ (ϕ)g is gravity, Fst is the source of surface tension, and F is the external force, which is referred to here in this paper as the electric field force.

2.2. Phase field equation

In the study of ascending movement of dispersed phase bubble, the interface change between dispersed phase and continuous phase is the key to research. In this paper, the phase field method is used to simulate and track the movement of the gas–liquid interface. The phase field method is used to describe the dynamics of gas–liquid twophase flow by Cahn–Hilliard equation,[18,19] the equation is given as follows:

In the formulas, the φ variable[20] is the auxiliary variable of the phase field, the dimensionless phase field variable ϕ changes from −1 to 1, ∂f/∂φ is the external free energy derivative, which is generally 0, γ is the mobility parameter, whose value determines the time factor of the Cahn–Hilliard equation diffusion with a given value of 1, ξ is the parameter that determines the thickness of the diffusion interface. When ξ → 0, ξ, the mixed energy density λ, and the surface tension coefficient σ satisfy the following relationship:[21]

2.3. Electric field equation

For the two-phase flow interface, the electric force needs to be specified. The bubble under the action of the electrostatic field F can be given by the divergence of the Maxwell stress tensor T as[22]

The Maxwell stress tensor T due to the difference in electrical property at the interface of the two phases can be expressed as[23,24]
where ε0 is the dielectric constant of vacuum and equals 8.8542 × 10−12 F/m, ε(ϕ) is the fluid relative permittivity, E is the electric field strength, and I is the second-order unit tensor. In the solution domain, the Maxwell stress tensor is 0 except for at the other positions in the two-phase interface.

Therefore, the electric field force Fcan be expressed as

On the right-hand side of Eq. (8), the first term is the electric field applied to the surface charge volume force. The second term is the dielectrophoresis force of the electric field. This term exists only in non-uniform medium, which is related to its corresponding spatial position of field strength and the spatial variation of the dielectric constant. The third term is the force of electrostriction, which is caused by the change of the medium density. If the fluid is incompressible, it can be ignored.[1]

The relationship between the electric field strength E and the potential V is

Since the dynamic current generated by the dynamics is small, the effect of the magnetic effect can be neglected. Therefore, the electric field can be regarded as a non-rotating field. Then the potential of any point can be solved according to the Laplace equation[25,26] with a variable permittivity

2.4. Fluid properties parameters

In this numerical simulation, there are corresponding flow parameters that can characterize the two-phase mixture. For the average flow density ρ (ϕ), viscosity coefficient μ(ϕ), and relative permittivity ε(ϕ) in the preceding equation, the simulation formulas are as follows:

where ρi, μi, and εi (i = 1,2) represent the density, viscosity, and relative permittivity of the two phases (air and fluid), respectively.

2.5. Dimensionless parameters

The following dimensionless parameters are used to describe the motion characteristics of the bubble: the Eotvos number E0, the Morton number M0, the Reynolds number Re, and the Weber number We,

where g is the gravitational acceleration, d is the bubble diameter, M0 represents the characteristics of the external fluid, and E0 characterizes the ratio of the float lift to the surface tension. These two numbers are controllable. And Re and We describe the bubble rise rate and bubble characteristics, respectively. Re represents the relative size of the viscous force and the inertial force in the fluid motion, We denotes the relative size of the inertial action and the surface tension, which are available only under the premise that the velocity is known.

3. Initial conditions
3.1. Geometric motion model

The rising bubble in the liquid is simulated and the whole system is considered to be affected by external electric fields (horizontal electric field EX and vertical electric field EZ). The motion area of the bubble is shown in Fig. 1.

Fig. 1. (color online) Schematic diagram of the movement area.
3.2. Initial condition setting

The main geometric parameters of the simulation are shown in Fig. 2. A 120 mm × 120 mm × 280 mm 3D gas–liquid simulation model is established, in which the bubble radius is 30 mm, the water depth is 250 mm, and the initial velocity of the bubble is 0 m/s.

Fig. 2. (color online) Schematic diagram of two-dimensional geometrical model in numerical simulation.

The flow field of the bubble is studied as a low Reynolds number flow field with a density ratio between two phases of 1000 and a viscosity ratio of 100. The main physical parameters used in the simulation are listed in Table 1.

Table 1.

Material properties in the model.

.
4. Results and analysis
4.1. Rising bubble without electric field

Figure 3 shows a dynamic change diagram of the bubble rising process as predicted by the numerical phase field model. The bubble is ideally spherical at t = 0 s. Then the bubble starts to rise gradually, the bottom of the bubble becomes flat and concave, and the indentation of the bottom of the bubble turns more and more serious over time. The shape of the bubble continues to change into a ball-cap shape. The simulation results in Fig. 3 are the same as the observations in Ref. [17].

Fig. 3. (color online) Bubble rising without electric field: (a) t = 0 s, (b) t = 0.08 s, (c) t = 0.12 s, (d) t = 0.18 s.

Figure 4 shows the flow field and pressure distribution in the computing domain at t = 0.12 s. It can be seen from the figure that due to the effect of gravity, the hydrostatic pressure distribution gradually increases from top to bottom. The pressure distribution in the bubble is greater than that of the external liquid because of the effect of the surface tension.

Fig. 4. (color online) Distributions of stream of flow field and pressure at t = 0.12 s.

In order to better observe the direction and change of velocity caused by bubble motion, the velocity vector around the magnified bubble is shown in Fig. 5. The arrow shows the direction of the velocity vector and the length of the arrow shows the magnitude of the velocity vector. The flow line of the bubble and its surrounding are obtained by the average rising velocity of the bubble, and the whole flow field is symmetrical. As the bubble density is less than that of the liquid, the bubble rises gradually from the rest and drives the liquid movement under the action of buoyancy, gravity, and other forces. The flow around the bubble is relatively dense, and a backflow zone is produced at the end of the bubble, with a vertically upward component of velocity.

Fig. 5. (color online) The velocity vector at t = 0.12 s.

The conservation of mass in the process of bubble movement is an important basis for studying the reliability and accuracy of numerical results. Since there is no reaction nor flow passing through the boundary, the total mass of each fluid should be constant. Figure 6 shows the total mass of the bubble as a function of time during the bubble movement. It can be seen from the figure that the bubble mass loss is less than 1% during the movement of the bubble, which indicates that the conservation of the bubble mass holds during simulation and also proves that the model is reliable and accurate.

Fig. 6. (color online) Total mass of the bubble as a function of time.
4.2. Effect of electric field intensity on bubble dynamics
4.2.1. Vertical electric field

The effect of the vertical electric field (the electric field applied in the Z direction) on the rising bubble is studied. The bubble of radius R = 30 mm under each of three different voltages (250000 V, 350000 V, and 500000 V) is simulated by the combination of buoyancy, electrodynamics, and gravity in fluid, and the bubble deformation and velocity field are numerical simulated.

In order to quantitatively study the effect of electric field on the rising bubble, figure 7 shows the effects of the vertical electric field intensity on the top position of the bubble under three different voltages. It can be seen that in the vertical electric field, the position of the bubble rising at the same time is higher than that in the case of non-electric field. As the intensity of the vertical electric field increases, the position of the bubble rises: the bigger the intensity of the vertical electric field is, the higher position the bubble reaches. This shows that the vertical electric field helps to accelerate bubble rising in hydrostatic water.

Fig. 7. (color online) Time-dependent bubble top positions under different vertical electric fields.

We define the aspect ratio of the bubble as α = W/H as shown in Fig. 8, where W is the long axis (width) of the bubble, and H is the short axis (height) of the bubble.

Fig. 8. (color online) Illustration of bubble measurement dimensions.

Figure 9 shows the bubble deformation under vertical electric field voltage of 0 V, 250000 V, 350000 V, and 500000 V at t = 0.2 s, and the bubble shape parameters at the same time are shown in Table 2. As shown in the table, the length of the horizontal axis of the bubble gradually decreases, the length of the vertical axis increases, and the horizontal-to-vertical deformation ratio decreases. The table shows that the electric field is directly proportional to the tensile strength of the bubble and the applied electric field strength.

Fig. 9. (color online) Bubble deformation diagram at t = 0.2 s under different vertical electric field voltages of (a) 0 V, (b) 250000 V, (c) 350000 V, and (d) 500000 V.
Fig. 10. (color online) Time-dependent bubble top positions under horizontal electric fields.
Table 2.

Bubble shape parameters at t = 0.2 s under different vertical electric fields.

.
4.2.2. Horizontal electric field

The dynamic process of bubble rising under the horizontal electric field (the electric field applied in the X direction) is studied. The voltage parameters are set to be 250000 V, 350000 V, and 400000 V, respectively.

In order to quantitatively study the effect of electric field on rising bubble, as previously described, figure 1 shows the effect of the horizontal electric field intensity change on the top position of the bubble. Obviously, the effect of the horizontal electric field is opposite to that of the vertical electric field, which slows the rise of the bubble. At the same time, the rising position of the bubble decreases with the increase of the voltage.

In order to better compare the influence of electric field direction on rising bubble, the time taken for the bubble to rise to the top of the region is shown in Fig. 11. From the figure, we can clearly see that the bubble rises faster under the vertical electric field, taking a shorter time to reach the top of the bubble; in the case of horizontal electric field, the bubble rises more slowly and takes longer time to reach the top.

Fig. 11. (color online) Time taken for the bubble to rise to the top.

Figure 12 illustrates the deformation of bubbles for each of the horizontal electric field voltages of 0 V, 250000 V, 350000 V, and 400000 V, the bubble shape parameters at the same time are given in Table 3. Compared with the case of a vertical electric field, the length of the horizontal axis of a bubble increases with the increase of the electric field intensity in the horizontal electric field, and the length of the vertical axis decreases, and the ratio of the horizontal to vertical deformation increases. A large value of W/H corresponds to the high resistance, which hinders the rising of the bubbles. The bubble directly deforms into an oblate ellipsoid. The reason is that the horizontal electric field tends to stretch the bubble into an oblate shape, which is beneficial to the action of the flow jet.

Fig. 12. (color online) Bubble deformation diagram at t = 0.2 s under horizontal electric field.
Table 3.

Bubble shape parameters at t = 0.2 s under horizontal electric field.

.

We know that during the rising of the bubble, the excessive movement of the lower gas will squeeze the upper gas and cause it to move to the left and right sides to form two arc-shaped bulges. The bubble forms a crescent arc as a whole. This study finds that when the voltage intensity is particularly high, the tail of the bubble is dragged to form a wake. As time goes on, the left and right two gas snare bags will eventually split out of the mainbody air bubble, becoming two separate small bubbles, while the remaining mainbody part of the bubble will continue to move upwards at a certain speed as shown in Fig. 13. Figure14 shows the shape of the bubble at U = 500000 V and t = 0.12 s: when the bubble buoyancy force is almost negligible compared with the electric field force, the bubble does not rise and is deformed in place due to the magnitude of the electric field force, and it is elongated in the direction of the field.

Fig. 13. (color online) Bubble shape at U = 250000 V and t = 0.44 s.
Fig. 14. (color online) Bubbles that do not rise but deform in place at U = 500000 V and t = 0.12 s.
5. Conclusions

It is difficult to solve the problem by using the traditional method, and the computational efficiency is low when the bubble rises in the liquid under the action of electric field. In this paper, the phase-field method is used to simulate the gas–liquid two-phase flow. The coupling effect of electric field and phase field on the bubble deformation and dynamic behavior in static liquid is studied. From the present study, some conclusions can be drawn as follows.

(i) A series of dynamic processes of bubble rising and deforming in liquid is simulated and the external flow field induced by the bubble rising process is obtained.

(ii) The feasibility of the model is verified by the conservation of the bubble mass in the rising process, and the effect is better.

(iii) The presence of a vertical electric field increases the rate at which bubbles rise in the liquid, and causes the bubbles to deform into a flat shape with a smaller width/height ratio, whereas a horizontal electric field is reversed.

(iv) The study of bubble rising under the action of horizontal electric field shows that the bubble does not rise when the bubble rising buoyancy force is almost negligible compared with the electric field force, and it is deformed due to the electric field force at the original place.

Reference
[1] Zhao T Liu Y P Lü F C Liu S 2016 J. Syst. Sl. 12 3081
[2] Zhang J S Lv Q Sun C D Lu D Chen L Y 2000 Acta Phot. Sin. 29 952 in Chinese
[3] Zu Y Q Yan Y Y 2009 Int. J. Heat. Fluid Fl. 30 761
[4] Amaya-Bower L Lee T 2010 Comput. Fluids. 39 1191
[5] Hua J Stene J F Lin P 2008 J. Comput. Phys. 227 3358
[6] Liu Q Shui H S Zhang X Y 2002 Adv. Mech. 32 259
[7] Zhang S J Wu C J Wang H M 2005 Hehai. Univ: Nat. Sci. Ed. 33 534
[8] Safi M A Turek S 2015 Comput. Math. Appl. 70 1290
[9] Wang L L Li G J Tian H 2011 Xi'an. Jiao. Tong. Univ. 45 65
[10] Shen J Yang X 2010 Chin. Ann. Math. 31 743
[11] Marco P D Grassi W Memoli G Takamasa T Tomiyama A Hosokawa S 2003 Int. J. Multiphas. Flow. 29 559
[12] Peng Y Feng F Song Y Z 2008 Chin. J. Chem. Eng. 16 178
[13] Andalib S Hokmabad B V Esmaeilzadeh E 2013 Colloids. Surfaces. Ace. 436 604
[14] Chen F Song Y Z Chen M 2006 J. Therm. Sci. Techno. 5 139
[15] Gurtin M E Polignone D Vinals J 1996 Math. Models. Methods. Appl. Sci. 6 815
[16] Lowengrub J Truskinovsky L 1998 Proc. R. Soc. Lond. 454 2617
[17] Sun T Sun J G Ang X Y Li S S Su X 2016 Int. J. Heat. Fluid Fl. 62 362
[18] Zhou C Yue P Feng J J Olliviergooch C F Hu H H 2006 J. Comput. Phys. 219 47
[19] Ding H Spelt P D M Shu C 2007 J. Comput. Phys. 226 2078
[20] Yue P Feng J J Liu C Shen J 2004 J. Fluid Mech. 515 293
[21] Jacqmin D 1999 J. Comput. Phys. 155 96
[22] Huang X Q 1995 J. Nanjing. Normal: Nat. Sci. Ed. 14 41
[23] Allen P H G Karayiannis T G 1995 Heat Recov. Syst. CHP 15 389
[24] Zaghdoudi M C Lallem M 2000 Int. J. Therm. Sci. 39 39
[25] Wu J Xu S Zhao N 2013 Chin. J. Comput. Phys. 30 1
[26] Teigen K E Munkejord S T 2009 IEEE. T. Dlelect. El. 16 475