A model for predicting bubble velocity in yield stress fluid at low Reynolds number

Chemical Engineering Science 201 (2019) 325–338

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces

A model for predicting bubble velocity in yield stress fluid at low Reynolds number Zhiyuan Wang ⇑, Wenqiang Lou, Baojiang Sun, Shaowei Pan, Xinxin Zhao, Hui Liu Offshore Petroleum Engineering Research Center, School of Petroleum Engineering, China University of Petroleum (East China), Qingdao 266580, China

h i g h l i g h t s
A new model was proposed to predict the bubble rise velocity in the yield stress fluid.
The energy conservation of rising bubble was considered in the model.
A model for predicting the velocity field around the rising bubble was established.
The model can predict the velocity of bubbles with a Reynolds number less than 0.4.

a r t i c l e

i n f o

Article history: Received 7 October 2018 Received in revised form 16 February 2019 Accepted 25 February 2019 Available online 6 March 2019 Keywords: Yield stress fluid Rising velocity Low Reynolds number Predicting model Velocity field

a b s t r a c t Currently, studies on bubble velocity in yield stress fluids with low Reynolds numbers are lacking. In most engineering applications, such value is often calculated using predictive models developed for viscous fluids. However, these models are only applicable to yield stress fluids in limited conditions because the bubble rise velocity is affected by the yield stress in these fluids significantly. In this study, a velocity field model around a bubble is first developed by performing an infinitesimal analysis of the velocity field in the equatorial plane of the bubble and combining the continuity equation with the constitutive equation of yield stress fluids. Next, a new model is proposed for predicting bubble velocity in yield stress fluids at a low transport rate by integrating the model of velocity field around the bubble with energy conservation equations. Finally, experimental studies are performed to analyze the bubble transport behavior in yield stress fluids at low Reynolds numbers. The experimental results showed that under a low Reynolds number, the bubble rise velocity is dominated by the yield stress of the fluid. Hence, we further validated the proposed bubble velocity prediction model with experimental measurements, which match closely with each other. We conclude that for a Reynolds number less than 0.4, the model developed in this study can be used to predict the bubble velocity in yield stress fluids with high accuracy. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Understanding the rising behavior of bubbles in static, nonNewtonian fluids has been a popular research topic; such phenomenon is typically encountered in many different engineering applications such as food processing, oil drilling, and fabrication of pharmaceuticals (Premlata et al., 2017; Lopez et al., 2018). Bubbles are not desired in many engineering applications. For example, in deep water drilling, gas intruding the wellbore and rising to the wellhead when the well is closed in typhoon weather will increase the risk of well opening (Sun et al., 2015). In cementing operations, gas intrusion and flow destroy the cementing effect ⇑ Corresponding authors. E-mail address: [email protected] (Z. Wang). https://doi.org/10.1016/j.ces.2019.02.035 0009-2509/Ó 2019 Elsevier Ltd. All rights reserved.

(Pinto et al., 2012). Therefore, it is very important to understand the flow dynamics and predict the rising velocity of bubbles, which is of great significance to improve the safety of engineering. The rising behavior of bubbles is affected by the properties of the surrounding fluids significantly. Because the viscosity of nonNewtonian fluids is dependent on the shear rate and yield stress (Beris et al., 1985), the rising velocity of bubbles in a nonNewtonian fluid can be distinctly different from that in a Newtonian fluid. In yield stress fluids, the bubble must overcome the yield stress to move. Otherwise, it will remain suspended in the fluid (Mirzaagha et al., 2017; Terasaka and Shibata, 2003). The fluids surrounding the bubble will be affected by its rising motion (Ohta et al., 2011); however, such effect is only limited to a certain region owing to the yield stress. The boundary of such region is known as the yield surface (Samson et al., 2017; Zenit and Feng,

326

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

Nomenclature Bn C C1–C4 dp/dz f Fr g hr hh i I1 I2 j k m n N p QL Qb r

Bingham number integration constant integration constant pressure gradient in the fluid shooting speed of the camera (s1) Froude number gravitational acceleration (m/s2) calculation step of radius integral calculation step of angle integral grid position index bubble wall yield surface bubble position index consistency factor (Pasn) total number of bubble velocities obtained from the image power law index number of grids pressure (Pa) flow rate of the fluid (m3) gas flow rate (m3) radius of an arbitrary point (m)

2018; Dubash and Frigaard, 2007; Tsamopoulos et al., 2008). The region sandwiched between the yield surface and bubble wall, known as the shear region, will undergo shear deformation owing to the impact from the bubble. The location of the yield surface is a critical factor for analyzing the bubble rising velocity, as it determines the flow field distribution around the bubble and affects the bubble rising velocity substantially. Dubash and Frigaard (2007), in performing the energy conservation analysis of the rising bubble, proposed the width of the shear region to be equal to the equivalent radius of the bubble, and the flow field follows a linear distribution in the shear region. Tsamopoulos et al. (2008) derived a numerical solution of the yield surface and velocity field around the bubble in a yield stress fluid using a mixed finite element method, and the Papanastasiou constitutive model. Mougin et al. (2012) measured the fluid velocity near a gas bubble in the BF Goodrich C940 solution using the PIV(Particle Image Velocimetry) technique and validated part of the numerical results derived by Tsamopoulos et al. (2008). Jianying Zhang (2012) simulated the velocity field around a gas bubble in a yield stress fluid using an augmented Lagrangian method. He later captured the yield surface around the bubble using an enhanced mesh adaptive technique (Zhang, 2018). Dimakopoulos et al. (2013) analyzed the yield surface and flow field around a gas bubble in a Herschel–Bulkley fluid by combining the augmented Lagrangian method with the Papanastasiou constitutive model. He proposed an empirical correlation between the yield surface around the bubble and the Bingham number (Bn = sy/qgRb) based on the calculation results. Duangsuwan et al. (2009, 2010, 2011) studied the shape and dynamics of gas and alcohol gas-liquid composite drops in vegetable oil in detail, analyzed the flow field distribution around the composite drops, and discussed the influence of temperature on the dynamics of composite drops. Yu and Fan (2010) based on the two-phase flow lattice Boltzmann model of multirelaxationtime interaction-potential, the bubble shape and rising velocity under low viscosity and high surface tension were simulated and predicted, and the validity of MRT-LBM in two-phase flow was verified. Hitherto, most studies on the velocity field around a gas bubble in a yield stress fluid were performed based on experiments or numerical simulations where the derivation of the velocity field was typically cumbersome and time consuming.

Rb Re Rsc Rx Ub Ub, 0

Uz Vb Wb Wy Wv z

c lc s sy e1 e2 q qg

j

radius of the bubble (m) Reynolds number radius of the yield surface (m) point with a zero velocity derivative bubble velocity (m/s) bubble velocity obtained between two adjacent frames (m/s) velocity of the fluid particle (m/s) bubble volume (m3) work performed by the bubble buoyancy force (W) work performed by the bubble yield stress (W) work performed by the bubble viscous force (W) position of the centroid of the bubble (m) shear rate (s1) characteristic viscosity of the yield stress fluid (Pa/s) shear stress (Pa) yield stress (Pa) continuity equation error limit energy conservation error limit liquid density (kg/m3) gas density (kg/m3)

As a bubble rises in the fluid, the work performed by the buoyance force is dissipated by the work performed by the inertial force, yield stress, viscous force, and surface tension of the bubble (Chhabra, 2006). The work performed by the yield stress and viscous force are dissipative works induced by the shear deformation of the fluids surrounding the bubble (Papanastasiou et al., 1999; White and Corfield, 2006). These works are closely related to the location of the yield surface and flow field around the bubbles. Dubash and Frigaard (2007) made two assumptions when performing the energy conservation calculation of the bubble: (1) the width of the shear region is equal to the equivalent radius of the bubble; (2) the velocity of the fluid particles in the shear region contains only the circumferential component and follows a linear distribution. The ratio between the work performed by friction force and buoyancy force, which should ideally be equal to 1, was found to be 0.3228 from experimental measurements. Such deviation is caused by two factors: (1) an error exists that is associated with the range of the shear region in the assumption; (2) a linear velocity distribution fails to consider the shear thinning behavior of yield stress fluids. Sikorski et al. (2009) validated Dubash’s model (Dubash and Frigaard, 2007) experimentally and concluded that the work performed by inertial force is negligible at a low Reynolds number. Currently, studies exploring the rising behavior of bubble in a yield stress fluid are scarce. Combining experimental measurements with an empirical correlation (Astarita and Apuzzo, 1965) and drag coefficient (Dijkhuizen et al., 2010; Lopez et al., 2018) is the primary method for predicting the bubble velocity in a yield stress fluid. A model for predicting bubble velocity in a yield stress fluid at a low Reynolds number has not been reported yet, as such value is often calculated based on predictive models developed for viscous fluids (Yin et al., 2016). Studies on bubble velocity prediction model in viscous fluids are relatively mature. Stokes (Hutter and Wang, 2016) calculated the velocity of a viscous fluid pass around a spherical object under laminar conditions and obtained a velocity prediction model for spherical bubbles in laminar flow (Re < 1). Peebles (1953) performed the force analysis of a gas bubble in a vertical tube of small diameter. By combining such analysis with bubble transport experimental data, they developed a bubble final velocity prediction model based on the fluid surface tension,

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

327

fluid density, and bubble diameter. Their model agrees well with Stokes’ model in laminar flow. By analyzing the rising velocity of a bubble in a power-law fluid, Rodrigue et al. (1999) proposed a drag coefficient at a low Reynolds number. His model is of a relatively simple structure, and is often the primary choice for analysis at low Reynolds numbers. However, the bubble experiences a strong resistance from the yield stress when rising in a low Reynolds number flow. Therefore, the bubble rising velocity prediction model developed for Newtonian and viscous fluids can only be used in limited conditions for yield stress fluids. In summary, the flow field distribution around a bubble is determined by the yield stress at a low Reynolds number. The yield dissipation in this case contributes to a majority fraction of total energy dissipation, and therefore cannot be neglected when analyzing the rising velocity of a bubble. The current empirical models for a low Reynolds number flow fail to consider the impact of yield stress on the rheology of the fluid, and therefore do not fit well with yield stress fluids. In the present work, we perform the following mechanistic studies to explore the bubble rising velocity. In Section 2, a model for describing the velocity field around a bubble is proposed and combined with the energy conservation equation to yield the bubble velocity prediction model. In Section 3, we discuss the experimental setup and properties of the experimental fluids, and further analyze the experimental results. In Section 4, we verify the model with the experimental results and determine the application range of the model. Fig. 1. Schematic of the geometric model for bubble transporting in a yield stress fluid. h is the rotation angle of the plumb to an arbitrary position in clockwise direction.

2. Bubble velocity prediction model The yield surface affects the bubble dynamics. The fluid beneath the yield surface is under the effect of the bubble and undergoes shear deformation. The fluid beyond the yield surface, however, is not affected by the bubble and remained still. In this study, the bubble and the yield surface are assumed to be of spherical shape. The equivalent radius of the bubble is expressed by

Rb ¼ ð3V b =4pÞ1=3

ð1Þ

where Vb is the volume of the bubble. The geometric configuration of the bubble and yield surface is shown in Fig. 1 (Tsamopoulos et al., 2008; Dimakopoulos et al., 2013; Lopez et al., 2018). Within the shear region, the mechanical energy of the flowing fluid is converted into internal energy irreversibly by shear deformation. As the bubble rises with constant velocity, the work performed by the bubble buoyancy force is dissipated by the shear deformation of the fluid, and an energy conservation equation can be constructed based on this feature. The flow field distribution in the shear region can be further derived by performing an infinitesimal analysis of the fluid in the equatorial plane of the bubble. Finally, the bubble velocity prediction model in a yield stress fluid can be obtained by substituting the velocity field of the fluid in the shear region into the energy conservation equation.

(

c¼0 ; s 6 sy i s ¼ kjcj þ sy jcj1  c ; s > sy h

n1

ð2Þ

where c is the shear rate, s is the shear stress, k is the consistency factor, and n is the power law index. When n = 1, Eq. (2) becomes the constitutive equation for the Bingham plastic. When the bubble moves at a low velocity in a viscous fluid, the inertial effect is negligible compared to the viscous effect during the bubble rise (White and Corfield, 2006). When the bubble moves at a low velocity in a yield stress fluid, the work performed by the inertial force is negligible as it is far less than the work performed by the viscous force and yield stress (Sikorski et al., 2009). Unlike the zigzag rising path in a viscous fluid (Shu and Yang, 2013; Tripathi et al., 2015), the rise of a bubble in a yield stress fluid follows almost a straight line once the bubble velocity becomes stable. Furthermore, the shape of the bubble remains stable during the rise. For a stable spherical bubble, the work performed by the surface tension during the rise is zero (Dubash and Frigaard, 2004, 2007). Therefore, as the bubble moves at a low velocity, the energy dissipation by surface tension and inertial force can be neglected, and the energy conservation equation can be simplified into Eq. (3) as

2.1. Energy conservation equation

Wb ¼ Wv þ Wy

During the rise of the bubble, the work performed by different forces satisfies the energy conservation principle. In other words, the work performed by the buoyance force is dissipated by the work performed by the inertial force, yield stress, viscous force, and surface tension of the bubble (Dubash and Frigaard, 2007; Sikorski et al., 2009). The viscous force and yield stress of fluid around bubbles are determined by the constitutive equation of fluid. The constitutive equation describing the nature of yield stress fluids was first proposed by Bingham (1922) and later expanded by Herschel and Bulkley (1926). The constitutive relation of yield stress fluids is given by Eq. (2).

where Wb, Wy, and Wv are the work performed by the bubble buoyancy force, yield stress, and viscous force, respectively. The buoyancy force is the only driving force for the bubble rising in the fluid. The relationship between the work performed by the bubble buoyancy force, bubble volume Vb, and bubble velocity Ub is given by Eq. (4) (Dubash and Frigaard, 2007) as

  W b ¼ q  qg gV b U b  qgV b U b

ð3Þ

ð4Þ

where q is the liquid density, qg is the gas density, g is the gravitational acceleration, Vb is the bubble volume, and Ub is the bubble

328

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

velocity. Under normal temperature and pressure conditions, we have qg/q  1. Therefore, the density difference between the liquid and gas is replaced by the liquid density q for calculating the work performed by the buoyancy force. Work dissipation is caused by the shear deformation of the yield stress fluid, which can be expressed by sc (White and Corfield, 2006; Papanastasiou et al., 1999). Based on the constitutive Eq. (2) of a yield stress fluid, the work dissipation can be divided into the work dissipated by yield stress and that dissipated by the viscous force. Within the shear region that is located between the bubble wall Rb and yield surface Rsc, the viscous dissipation is proportional to kcn+1 while the yield stress dissipation is proportional to syc. Based on the bubble geometric model shown in Fig. 1 and the coordinate system, the viscous dissipation and yield stress dissipation equations during the bubble rise can be constructed as shown in Eqs. (5) and (6):

Wv ¼

Wy ¼

Z pZ 0

Rb

Z pZ 0

Rsc

Rsc

Rb

2pr2 sinhk

 nþ1 dU drdh dr

2pr 2 sinhsy

  dU drdh dr

ð5Þ

ð6Þ

where Rb is the radius of the bubble, Rsc is the radius of the yield surface, r is the radius of an arbitrary point in the shear region, and dU/dr is the velocity gradient of the flow field in the shear region (c = dU/dr in the shear region). Based on Eqs. (5) and (6), we found that the velocity field of the fluid in the shear region is the key for calculating the yield stress and viscous dissipation. 2.2. Velocity field around the bubble The velocity distribution of fluid particles in the shear region is the basis for calculating the energy dissipation, and an important factor for analyzing the rising behavior of the bubble. The velocity field in the shear region around the bubble in a yield stress fluid is affected by the bubble size, yield stress of the fluid, and shear thinning behavior. Such characteristic is distinctive from that of a Newtonian fluid. A stable rising bubble exhibits axial symmetry (Premlata et al., 2017; Zhang, 2018; Tripathi et al., 2015) A similar axial symmetric feature is also found in the yield surface (Samson et al., 2017; Zhang, 2018). Therefore, concentric and cylindrical differential elements with a height of dh can be obtained in the bubble and its

surrounding shear region along the equatorial plane, as shown in Fig. 2b. The velocity vector of fluid particles around the gas bubble can be decomposed into the radial and circumference components (Papanastasiou et al., 1999; White and Corfield, 2006). The velocity vector of fluid particles in the shear region on the equatorial plane contains only the circumferential component while the radial velocity component is zero (Tsamopoulos et al., 2008). The velocity vector of fluids in the cylindrical differential element contains only the vertical component. Therefore, the fluid flow in the shear region on the equatorial plane can be simplified as the flow of yield stress fluid in a concentric ring. A stable rising bubble will create a displacing effect to the fluid beneath the yield surface. Such displacing effect can be considered as an extra force acting on the fluid. The stable flow of the fluid in the differential element on the equatorial plane in the shear region creates a stable pressure gradient. As the bubble rises stably, the flow pattern along the equatorial plane is as shown in Fig. 2c. The boundary of the differential element I1 (bubble wall) will continue to rise at a velocity of Ub while the boundary I2 (yield surface) remains still. A stable positive pressure gradient (dp/dz > 0) exists in the flow field between I1 and I2. The governing equations of the fluid flow on the equatorial plane in the shear region consist of the continuity and momentum equations, which are given by Eqs. (7) and (8) (Liu and Zhu, 2010) as

rU ¼0

q

@U þ ðU  rÞU ¼ rp þ r  s @t

ð7Þ ð8Þ

where q, U, p, and s are the density, velocity, pressure, and deviatoric stress tensor of the fluid, respectively. In a symmetric cylindrical coordinate system (r, h, z), the flow field on the equatorial plane contains the velocity component Uz along the z-axis direction (Filip and David, 2003). Therefore, the continuity and momentum equations are automatically satisfied. The axial momentum equation can be simplified into Eq. (9) as

dðsrz rÞ dp ¼ rdr dz

ð9Þ

where srz is the only nonzero component of the deviatoric stress tensor s and dp/dz is the pressure gradient in the fluid. Integrating over the momentum, Eq. (9) yields the general solution to the stress equation of the fluid on the equatorial plane:

Fig. 2. Schematic of flow field in a differential element on the equatorial plane of the bubble (h = 90°). I1 is the bubble wall and I2 is the yield surface.

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

srz ¼

r dp C þ 2 dz r

ð10Þ

where C is the integration constant. As shown in the flow pattern in Fig. 2b, it is extremely difficult to obtain the analytical solution to the velocity distribution of a Herschel–Bulkley fluid. In this study, we derived the solution using a power law index of n = 1. The constitutive equation of the yield stress fluid with n = 1 is given by Eq. (11) as



srz ¼ k þ jdUszy=drj dU z dr



dU z dr

for jsrz j > sy

¼ 0 for jsrz j 6 sy

ð11Þ

where dUz/dr is the velocity gradient of the flow field on the equatorial plane. A general solution of the velocity distribution Uz on the equatorial plane can be obtained by combining Eqs. (10) and (11) as

Uz ¼

  r 2 dp C 1 sy dU z r þ C2 þ lnr  sgn 4k dz k k dr

ð12Þ

where C1 and C2 are the integration constants, and Uz is the velocity of fluid particles on the equatorial plane. The fluid is not completely slippery on the bubble wall (Tsamopoulos et al., 2008; Dimakopoulos et al., 2013). Therefore, a point with zero velocity gradient Rx must exist in the region between I1 and I2, as shown in Fig. 2c (Liu and Zhu, 2010). The entire velocity field is divided into two segments by Rx. The analytical expression of the velocity distribution in the shear region is given by Eq. (13)

( Uz ¼

s dp þ Ck1 lnr þ ky r þ C 2 4k dz sy C3 r 2 dp þ k lnr  k r þ C 4 4k dz r2

Rb 6 r 6 Rx Rx 6 r 6 Rsc

ð13Þ

where C3 and C4 are integration constants and Rx is the point with a zero velocity derivative. 2.3. Bubble velocity prediction model A set of equations describing the bubble motion in the fluid is constructed by combining the velocity field around the gas bubble and the energy conservation equation. The bubble velocity under different conditions can be obtained by solving these equations. The introduction of proper boundary conditions is key for calculating the velocity field around the bubble. When a bubble is rising in the fluid, the boundary condition at the bubble wall I1 is vastly different from that imposed on a solid sphere. Traditionally, a no-slip boundary condition is used in studies on the transport of solid sphere and internal flows (Beris et al., 1985; Liu and Zhu, 2010; Ahonguio et al., 2014). However, a zero shear boundary condition is imposed on the bubble wall for analyzing bubble motions (Van Hout et al., 2002; Tsamopoulos et al., 2008; Dimakopoulos et al., 2013). Tsamopoulos et al. (2008) proposed the zero shear boundary condition when analyzing the yield surface of a bubble. He applied this boundary condition in numerical simulations and obtained the flow field around the bubble at different Bingham numbers. Using the velocities of fluid particles at the bubble wall from his simulation results, we obtained the relationship between the velocity of the fluid particle and the Bingham number Bn via linear regression analysis. Eq. (14) is used to describe the zero shear boundary condition at the interface I1.

U 0z ¼ U b  ð0:2242  4:1138Bn Þ 0

ð14Þ

where Uz is the velocity of the fluid particle at interface I1. The relative position of the yield surface I2 and bubble wall I1 determines the fluid velocity distribution in the shear region, the viscous force experienced by the bubble during rising, and the magnitude of the work performed by the yield stress.

329

Dimakopoulos et al. (2013) obtained the radius of yield surface Rsc with different Bingham numbers by numerical simulations. Using the literature data from his study, we obtained the relationship between the radius of the yield surface Rsc and the Bingham number Bn by the least-squares fitting method, as shown in Eq. (15):

Rsc ¼ 0:6465Rb  B0:591 n

ð15Þ

In the solution with the same rheological parameters, the bubble radius increases, the Bingham number Bn decreases, and the effect of inertia term on the rising velocity of bubbles increases. In this paper, the influence of inertia term is neglected, and the model is only applicable to low-speed migration. Therefore, the applicability of boundary condition Eqs. (14) and (15) needs to be limited. Eqs. (14) and (15) According to the results of Tsamopoulos et al. (2008), the boundary conditions have higher accuracy under the condition of Bn > 0.03. The higher the number of Bn, the lower the rising velocity of bubbles, the stronger the applicability of the model. Other boundary conditions used for solving Eq. (13) include the following: the velocity is zero at the yield surface Rsc (in order to simplify the calculation, the boundary of neglecting dUz/dr = 0 is neglected), the velocity is continuous at Rx, the derivative of the velocity is zero at Rx (i.e., the shear stress equals to the yield stress sy). These three boundary conditions are described by Eqs. (16), (17), and (18) as

U z jr¼Rsc ¼ 0

ð16Þ

U z jr¼Rx ¼ U z jr¼Rþx

ð17Þ

        sjr¼Rx  ¼ sjr¼Rþx  ¼ sy

ð18Þ

When the bubble enters the yield stress fluid domain, the flow rate Q at any cross section in the fluid domain should be equivalent if no mass transfer occurs at the surface of the gas bubble. When a single bubble reaches a steady rising state in the yield stress fluid, continuous gas phases do not occur, and the local flow rate should be zero in the fluid domain. With any cross section of a stable rising bubble, the gas flow rate Qb at the cross section induced by the rise of the bubble should be equal to the flow rate QL generated by the flow of the yield stress fluid in the shear region (Van Hout et al., 2002; Sousa et al., 2005; Nogueira et al., 2006). As shown in Fig. 2b, a continuity equation is constructed for the differential element on the equatorial plane of the bubble. The integral form of the continuity equation is given by Eq. (19):

pR2b  U b ¼

Z

Rsc Rb

2pr  U z dr

ð19Þ

We substituted the aforementioned boundary conditions (14)– (18) into the velocity field Eq. (13) and combined it with the continuity Eq. (19) to obtain a velocity field model consisting of six equations and seven unknowns: C1, C2, C3, C4, Rx, Ub, and dp/dz. Apparently, the velocity field model cannot be solved by itself because the number of unknowns exceeds the number of equations. These equations must be combined with the energy conservation equation to yield a closed-form equation system. The velocity of the fluid particles in the shear region is mandatory for calculating the energy dissipation during the rise of the bubble. In this study, the same assumptions proposed by Dubash and Frigaard (2007) were adopted in our analysis, where the radial velocity components are neglected owing to their slight impact on the energy dissipation. Therefore, the energy conservation calculation is performed based on the circumferential velocity component U. The circumferential velocity of any fluid particle near the bubble

330

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

is related to the velocity on the equatorial plane. Such relation is given by Eq. (20) (Tsamopoulos et al., 2008):

U ¼ U z  sinh

ð20Þ

where U is the circumferential velocity component of an arbitrary fluid particle around the bubble. The energy dissipation of any fluid particle in the shear region can now be solved based on the velocity distribution on the equatorial plane. The relationship between the energy conservation equation and the velocity field on the equatorial plane of the bubble is given by Eq. (21):

qgV b U b ¼

 nþ1 dU z drdh dr 0 Rb   Z p Z Rsc dU z 2 drdh þ 2pr 2 sin hsy dr 0 Rb Z pZ

Rsc

2pr2 sin

nþ2

hk

ð21Þ

The bubble velocity model is composed of Eqs. (13), (19), and (21). Substituting the boundary condition Eqs. (14)–(18) into the bubble velocity prediction model yields a set of equations with seven unknowns and seven equations. The detailed equations for solving the model are shown in Eq. (22). Eqs. (22a)–(22e) is a velocity field equation after introducing boundary conditions. Eq. (22f) is a continuity equation and Eq. (22g) is an energy conservation equation. With known bubble dimensions and fluid properties, the bubble rise velocity can be obtained by solving the equation sets.

Rx dp C 1 þ ¼ sy 2 dz Rx

ð22aÞ

Rx dp C 3 þ ¼ sy 2 dz Rx

ð22bÞ

2sy Rx þ ðC 1  C 3 ÞlnRx þ C 2  C 4 ¼ 0 k

ð22cÞ

R2b dp C 1 sy þ lnRb þ Rb þ C 2  U b  ð0:2242  4:1138Bn Þ ¼ 0 4k dz k k ð22dÞ 

0:6465Rb  B0:591 n

2

 dp C 1  þ ln 0:6465Rb  B0:591 n dz 4k k  sy  0:591  0:6465Rb  Bn þ C2 ¼ 0 k

p  Ub ¼ R2b

Z

Rsc

Rb

2pr  U z dr

 nþ1 dU z nþ2 qgV b U b ¼ 2pr2 sin hk drdh dr 0 Rb   Z p Z Rsc dU z 2 drdh 2pr 2 sin hsy þ dr 0 Rb Z pZ

ð22eÞ

3.1. Experimental apparatus Gas transport experiments were performed in a transparent plexiglass sink. A square sink with a side length of 40 cm and height of 50 cm was used. Fig. 3 shows the experimental apparatus used in this study. The ratio between the bubble diameter and sink height was less than 1/20 in our experiments. The effect of sink wall on the bubble rising behavior was neglected. Compared to a cylindrical tube, a square sink resolves the dimension distortion issue when recording the experiments with a camera. In this case, the volume of the droplet recorded during the experiment was closer to its actual value. Air was pressurized using a hand pump and injected into the square sink through a 3-mm pipeline. Bubbles formed near the nozzle located at the center bottom of the square sink. To ensure a stable gas injection process and reduce the pressure instability induced by the friction of the pipeline, a gas reservoir was implemented beneath the nozzle for stabilizing the operating pressure (Terasaka and Tsuge, 2001). The volume of the bubble was controlled by the injection speed from the hand pump and the opening degree of the valve during the experiment. The rising phenomena of the bubble were recorded by a highspeed Olympus i-SPEED TR camera during the experiment. The backlight was placed on the other side of the camera. The highspeed camera was affixed at 15 cm above the top edge of the sink for obtaining the size and final velocity of the gas bubble through images. We confirmed in our preliminary experiments that the bubble shape and velocity were stabilized at such height. The entire experiment was performed at the room temperature of 23 ± 2 °C. In the process of shooting, 200 frames per second were selected, and the corresponding resolution was 1280  1024 pixels. The experimental image was processed by Image-Pro Plus software, and the geometry and dynamics of the bubble were obtained. This information includes the position of the centroid, the cross-section area, the height and width of the bubble, the relative error of data is less than 0.5%. Bubble rising velocity is obtained by data processing obtained by software.

U b;j ¼

ð23Þ

where f is the shooting speed of the camera, z is the position of the centroid of the bubble, and Ub, j is the bubble velocity obtained between two adjacent frames. The average bubble velocity obtained by image processing is taken as the final bubble velocity.

Ub ¼ ð22fÞ

zjþ1  zj 1=f

m 1 X U b;j m j¼1

ð24Þ

where, m is the total number of bubble velocities obtained from the image.

Rsc

3.2. Fluid preparation and rheology

ð22gÞ

3. Experiment and result analysis Compared to viscous fluids, the presence of yield stress induces changes in the friction force caused by the bubble when rising, and affects the bubble rise velocity substantially (Santos and Azar, 1997). When the bubble is moving at a low velocity, the energy dissipation is dominated by the yield stress (Sikorski et al., 2009) that determines the rise velocity of the bubble. Therefore, it is necessary to perform bubble transport experiments in a yield stress fluid to analyze the impact of yield stress on bubble velocity.

Carbopol is a highly molecular polymer that forms a stable, transparent, and thixotropic colloid when dispersed in water and neutralized (Mitsoulis and Tsamopoulos, 2017; Curran et al., 2002; Piau, 2007). The yield stress of Carbopol gel can be tuned by controlling the polymer concentration and pH value. When the bubble is moving at a low velocity, the surrounding fluids are at low shear rate conditions. In this case, the elastic behavior of Carbopol gel can be neglected (Sikorski et al., 2009), and such material can be approximated as an ideal yield stress fluid. Carbopol 940 was used in this study. The gel of Carbopol 940 is a highly transparent fluid ideal for capturing videos. First, 80 L deionized water was added into the container. Subsequently, a fixed amount of powder was added slowly to water while the stirrer

331

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

3.3. Experimental results Fig. 5 shows the relationship between the bubble rise velocity and bubble volume in different yield stress fluids. As shown in the figure, increasing the yield stress in the solution results in a gradual decrease in the bubble rise velocity with the same volume.

100

100

10

10

1

1

viscosity (Pa·s)

operated at 200–400 rpm continuously. Subsequently, Carbopol powders were added to the container, and the mixture was further stirred for 6–8 h for dispersing the powder uniformly in water for hydrolysis. After a complete hydrolysis, the stirrer was shut down and the solution was placed in static for 24 h to ensure the collapse of foam and release of gas. At this stage, the pH of the solution was approximately 3. The solution although viscous lacked yield stress. Neutralizing the solution will enable it to become more viscous and exhibit yield stress. When performing the neutralization, the stirrer continued to operate at a low speed of 50–100 rpm (to prevent the formation of gas bubble), and saturated NaOH solution at room temperature was added slowly into the solution along the stirring impeller. After the solution attained the yield stress, the rotation speed of the stirrer was increased slightly with the continuous addition of NaOH solution until the pH of the solution reached a stable value of 6.5–7. The shear deformation of the Carbopol solution exhibits a historical effect. After one bubble rises in the solution, a path will remain that affects the rising behavior of subsequent gas bubbles (Mougin et al., 2012). After transferring the Carbopol solution into the experiment sink slowly, the solution was stirred gently to remove any historical effect. After completing the transport experiment with a single bubble, the solution was stirred again to remove the path generated by the previous bubble during rising. Carbopol solutions with five different mass concentrations (0.10%, 0.105%, 0.11%, 0.115%, and 0.12%) were used in the experiments. The rheological curves of the solutions were determined using the Physica MCR 301 rheometer at 25 °C. The stress and viscosity of each test sample under different shear rates are summarized in Fig. 4. After fitting the rheological parameters of the solutions with the constitutive equations of the yield stress fluid, we found that Carbopol solutions satisfied the constitutive equation of the Herschel–Bulkley fluid excellently. Table 1 shows the detailed parameters of the five solutions obtained by curve fitting.

Shear stress (Pa)

Fig. 3. Schematic diagram of the experimental setup.

Shear stress Viscosity 0.1 0.1

1

10 Shear tate(s-1)

100

0.1 1000

Fig. 4. Flow curves for the 0.105% Carbopol solutions. The Herschel-Bulkley model fits to the stress data are also shown (dotted line).

When the bubble is moving at a low velocity, the bubble transport is affected by the yield stress of the fluid substantially, where increasing the yield stress creates a stronger friction force to the bubble. Therefore, the bubble rise velocity reduces with increasing volume. The existence of yield stress also induces a significant distinction in bubble dynamics in a yield stress fluid compared to that in other non-Newtonian fluids. When the buoyancy of the bubble is not sufficiently strong to overcome the yield stress, the bubble will remain floating in the solution (Dubash and Frigaard, 2004, 2007; Sikorski et al., 2009; Dimakopoulos et al., 2013). As shown in Fig. 5, the projection of the bubble transport velocity on the horizontal axis is the critical volume of a suspended bubble, which increases gradually with increasing yield stress. Based on the experimental data, the minimum and maximum critical volumes of a suspended bubble are approximately 2.38  107 m3(0.10%) and 5.79  107 m3(0.12%), respectively. Within a yield stress fluid, the bubble dynamics is affected by the viscous force, yield stress, and inertial force simultaneously (Dhiman et al., 2007; Islam et al., 2015). The impact from each factor on the bubble transport can be analyzed using dimensionless

332

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

Table 1 Rheology of the Carbopol solutions. Solution

Mass percentage (%)

pH

r(N/m)

sy(Pa)

k(Pasn)

n

S1 S2 S3 S4 S5

0.10 0.105 0.11 0.115 0.12

6.65 7.15 6.55 6.56 6.55

0.0662 0.0658 0.0660 0.0657 0.0663

5.69 5.96 6.38 7.20 7.73

5.00 6.92 9.21 9.72 10.49

0.3623 0.3331 0.3127 0.3108 0.3184

0.09

0.5

S1 S2 S3 S4 S5

0.07

Bn S1 S2 S3 S4 S5

0.4

0.3

0.16

0.12

Re

U/(m/s)

0.06

Re

Bn

0.08

0.2

0.05 0.04

0.2

0.08

0.1

0.04

0.03 0.02 0.01

0 0.003

0 0

5

10 15 V/( 10-7m3)

20

25

Fig. 5. Relationship between terminal velocity and bubble volume.

parameters that primarily include the Reynolds number Re, Bingham number Bn, and Froude number Fr. These parameters are defined by the following equations:

Re ¼

2qU b Rb

lc

ð25Þ

sy qgRb

ð26Þ

Ub Fr ¼ pffiffiffiffiffiffiffiffi gRb

ð27Þ

Bn ¼

In the equations above, lc is the characteristic viscosity of the yield stress fluid (Lopez et al., 2018), which is defined by Eq. (28) as

lc ¼

sy þ kcn Ub ; c¼ c Rb

ð28Þ

During the experiment, it is not possible to control the value of a certain dimensionless parameter directly because their values vary with different bubble volume. The effects of dimensionless parameters on the bubble dynamics are shown in Figs. 6–8. Fig. 6 shows how the Reynolds number and Bingham number change with different bubble radii. For all bubble transport experiments performed in this study, the Reynolds number was always less than 0.5. In addition, increasing the equivalent radius of the bubble results in a higher Reynolds number and a smaller Bingham number. This suggests that the yield stress and viscosity dominate the bubble dynamics when the bubble volume is small. By increasing the bubble volume, the impact from yield stress decreases slowly while the effect of inertia becomes increasingly important. In addition, from the relationship between Bingham number and bubble radius in Fig. 6, it can be seen that all Bingham numbers satisfy

0.004

0.005

0.006 R/(m)

0.007

0.008

0 0.009

Fig. 6. Reynolds and Bingham numbers vs equivalent radius.

Bn > 0.08 in the experimental data in this paper. At this time, the Bingham numbers of all bubbles can satisfy the Bn > 0.03 applicable condition of Eqs. (14) and (15). The rising speed of bubbles can be predicted by the model in this paper. The relationship between the Reynolds number and Bingham number is shown in Fig. 7(a). For the same solution, a higher Bingham number results in a smaller Reynolds number and a slower bubble velocity. When the Reynolds number approaches zero, the bubble velocity becomes almost zero. In this case, the Bingham number is maximized, which is known as the critical Bingham number of bubble transport. According to the data shown in Fig. 7, the critical Bingham number of a suspended bubble is 0.153. If Bn < 0.153, the bubble can move in the fluid. In reality, the critical Bingham number will be even higher than such value. This is because apart from the yield stress, the bubble is affected by the surface tension in our experiments, thus preventing any bubble smaller than the critical size to escape from the nozzle (Lopez et al., 2018). The effects of the inertial term and yield stress term on the bubble dynamics are illustrated in Fig. 7(b). The results show that the Froude number Fr decreases with increasing Bingham number Bn. For the same solution, a smaller bubble volume results in a larger Bingham number and a smaller Froude number. This indicates that the effect of yield stress on the bubble velocity is more significant than the effect of the inertial term. By increasing the bubble volume, the effect of yield stress on bubble transport becomes weaker while the impact from the inertial term becomes stronger. 4. Model verification and application 4.1. Comparisons with experimental result Simpson formula is used to solve the integral equation to solve the bubble velocity prediction model with seven variables, we

333

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

0.08

0.5

0.07

S1 S2 S3 S4 S5

0.05

Re

0.3

0.06

U/(m/s)

0.4

S1-Experiment S3-Experiment S5-Experiment S1-Model S3-Model S5-Model

Re=0.41 Re=0.40 Re=0.4

0.04 0.03

0.2

0.02

0.1

0 0.08

0.01 0 0.003

0.1

0.12

0.14

0.004

0.005

0.16

Bn

0.006 R/(m)

0.007

0.008

(a)

(a) 0.08

0.3 0.07

S1 S2 S3 S4 S5

Re=0.43

Re=0.43

U/(m/s)

0.05

Fr

0.2

0.06

S2-Experiment S4-Experiment S2-Model S4-Model

0.04 0.03

0.1

0.02 0.01

0 0.08

0 0.003

0.1

0.12 Bn

0.14

0.16

0.004

0.005

0.006 R/(m)

0.007

0.008

(b)

(b) Fig. 7. (a) Relationship between Reynolds numbers and Bingham numbers. (b) Relationship between Froude numbers and Bingham numbers.

obtained the relationship between bubble rise velocity and equivalent radius in different yield stress fluids. Using the rheological parameters of the testing fluid shown in Table 1 and the velocity prediction model, we predicted the bubble velocity in different yield stress fluids. The relationship between bubble velocity and equivalent radius in solutions S1–S5 is shown in Fig. 8. The dashed lines represent the bubble velocity curves obtained from simulation calculations. Comparing the values in Fig. 8(a) and (b), we found that the error between the simulation and experimental results is less than 20% when the Reynolds number is less than 0.4. Therefore, our proposed model can provide a decent prediction of the bubble velocity. Lopez et al.(2018) performed systematic experimental studies of bubble rising in yield stress fluids. Their experimental conditions are similar to ours described herein. In their study, the Reynolds number was less than 0.5, and the ratio between the bubble diameter and length of experimental sink was less than 1/10. Therefore,

Fig. 8. Comparison between the proposed model predictions and experimental results. (a) Data sets S1, S3 and S5; (b) Data sets S2 and S4.

the effect of sink wall on the bubble rising behavior can be neglected. We validated two test groups in their study that were associated with a concentration of 0.12% and 0.15%, respectively. Because our model is highly sensitive to rheological parameters, it is necessary to first determine these parameters of the target solution. As the rheological curves of the testing fluids described in the literature study contains certain discrete points, we reperformed a fitting of the rheological parameters after removing two abnormal points measured under a low rotation speed. The rheological parameters are summarized in Table 2. Based on the rheological parameters of the testing fluids listed in Table 2, we predicted the bubble rise velocities. The comparisons between the prediction results from our model and the experimental measurements are shown in Fig. 9. The data shown in the figure indicate that the model proposed herein can provide a decent prediction of the bubble rise velocity for Re < 0.4. The shear region is the key factor affecting the distribution of velocity field around bubbles. Yield stress sy is the key factor

334

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

tion on the application range of the proposed model. For Re < 0.4, the model developed herein can provide an accurate prediction of the bubble rise velocity in yield stress fluids. For Re > 0.4, however, the error generated from the simulation exceeds 20%. Therefore, the effect of the inertial term must be considered. In addition, increasing the yield stress can result in a higher aspect ratio of the bubble (Samson et al., 2017; Lopez et al., 2018), as shown in Fig. 10, thereby changing the flow field around the bubble. When the yield stress sy exceeds 15.8 Pa, the aspect ratio of the bubble can exceed 2.5 (Lopez et al., 2018). In this case, the assumption of a spherical bubble shape in our model is no longer valid, thereby degrading the accuracy of our predictive model.

0.09 0.08 0.07

S6-Lopez’s experiment S7-Lopez’s experiment S6-Model S7-Model Tsamopoulos' result Present model

Re=0.41 Re=0.46

U/(m/s)

0.06 0.05 0.04 0.03 0.02

4.2. Model application

0.01 0 0.0015

0.0035

0.0055

0.0075 R/(m)

0.0095

0.0115

Fig. 9. Comparison between the proposed model predictions and Lopez’s experimental results (Lopez et al., 2018) and Tsamopoulos’ result (2008). In Tsamopoulos’ results, the yield stress sy = 3 Pa, k = 0.0601 Pas, n = 1.

Table 2 Rheological parameters of the testing fluids used in Lopez’s experiments (Lopez et al., 2018). Solution

Mass percentage (%)

pH

sy

k

n

S6 S7

0.12 0.15

– –

11.6 15.8

3.505 6.137

0.4158 0.3690

affecting the position of yield interface. Power exponent n does not affect the size of shear region, so power exponent n has relatively small influence on the velocity field around bubbles. In Fig. 9, the model in this paper is in good agreement with Tsamopoulos’ results (2008). It can be seen that the calculation results of this model are less sensitive to the power exponent and can be applied to the prediction of bubble velocity in n  1 yield stress fluid. The energy conservation equation used in the bubble velocity prediction model is developed assuming the energy dissipation is composed primarily of viscous dissipation and yield dissipation. The effect of inertia on the bubble rise is neglected in this equation. Because yield stress fluids exhibit a relatively large characteristic viscosity, the bubble rise velocity in yield stress fluids is generally higher under the same Reynolds number. With increasing bubble rise velocity, the increasing amount of work performed by bubble buoyancy force is now dissipated by inertial force (Sikorski et al., 2009). This induces greater error in the energy conservation equation. The energy dissipation by the inertial term imposes a limita-

The low velocity bubble transport in yield stress fluids is a typical phenomenon in many different engineering scenarios. However, the severity of such issue is typically neglected. For example, when predicting the bubble rise velocity in a drilling fluid during oil drilling, a bubble velocity prediction model developed for viscous fluids is often used (Yin et al., 2016). However, owing to the difference in the rheological properties of the operating fluid, models developed for viscous fluids have a limited application range and cannot be used to provide the accurate prediction of bubble rise velocity in the yield stress drilling fluid. The Peebles model (Peebles, 1953), Stokes model (Hutter and Wang, 2016), and Rodrigue small Reynolds number method (Rodrigue et al., 1999) are the most typically used velocity prediction models for viscous fluids in laminar flow. When applying these models to predict the bubble velocity in yield stress fluids, the viscosity of the model must be corrected prior to the calculation. The characteristic viscosity of yield stress fluid is expressed by Eq. (28). Therefore, bubble velocity calculation models compatible with yield stress fluids can be obtained by substituting the characteristic viscosity into the velocity prediction models developed for viscous fluids. We calculated the bubble velocity by substituting the rheological properties of yield stress fluids in the corrected Stokes model (Hutter and Wang, 2016) and the Rodrigue small Reynolds number model (Rodrigue et al., 1999). Fig. 11 summarizes the calculation results obtained using the model developed herein, the Stokes model (Hutter and Wang, 2016), and the Rodrigue small Reynolds number model (Rodrigue et al., 1999). The comparison of the results shows that the bubble velocity prediction model developed for viscous fluids can only be used in limited scenarios and cannot provide the accurate prediction of the bubble rise velocity in a yield stress fluid. Yield stress in the fluid is the key factor restricting the applicability of classical models. Because the fluid exhibits a higher characteristic viscosity under the effect of yield stress, the bubbles are affected by a higher equivalent friction

Fig. 10. Schematic diagram of bubble aspect ratio in different yield stress fluids. a: the width of the bubble a is w = 5 mm, the length h = 5.7 mm, and the aspect ratio h/ w = 1.14; b: the width of the bubble a is w = 8.9 mm, the length h = 11.9 mm, and the aspect ratio h/w = 1.34; c: the width of the bubble a is w = 12 mm, the length h = 18.8 mm, and the aspect ratio h/w = 1.57; d: the width of the bubble a is w = 15.3 mm, the length h = 30.7 mm, and the aspect ratio h/w = 2.01.

335

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

0.09

1.6 Present experiment

0.08 0.07

Stokes's model Rodrigue's model

Bn=0.14

Present model

1.4

0.06

Bn=0.15

1.3

0.05

U/Ub

U/(m/s)

Bn=0.13

1.5

Re=0.40

0.04

1.2

0.03

1.1

0.02

1

0.01 0 0.004

0.005

0.006 R/(m)

0.007

0.008

0.9 1

1.5

(a)

2

2.5 r/Rb

3

3.5

4

Fig. 12. Velocity field distribution around the bubble.

0.09 0.08 0.07

Lopez ’s experiment Stokes's model Rodrigue's model Present model

shear region increases monotonically. The rate at which the maximum velocity increases is also larger at higher Bingham numbers. This is because the distance between the yield surface and the bubble wall becomes smaller with increasing Bingham number. Therefore, the flow region of the fluid on the equatorial plane becomes smaller. Once the yield surface coincides with the bubble wall completely, the flow interface of the fluid becomes zero and the bubble will be suspended in the yield stress fluid (Dimakopoulos et al., 2013; Zhang, 2018).

Re=0.41

U/(m/s)

0.06 0.05 0.04 0.03 0.02

5. Conclusions

0.01

The primary difference between yield stress fluids and other non-Newtonian fluids is the yield stress in the fluid. The relative magnitude of the yield stress and the bubble buoyancy force determines whether the bubble could move in the fluid. When the bubble moved in the fluid, the yield stress affected the position of the yield surface around the bubble, thus further affecting the flow field surrounding the bubble. Therefore, the yield stress of the fluid affected the bubble rise velocity significantly, and is a distinct feature from other non-Newtonian fluids. Based on the energy conservation equation and considering the rheological properties of yield stress fluids, a bubble velocity prediction model at a low Reynolds number was developed in this study. Compared to the bubble velocity model developed for viscous fluids, the model proposed in this study considered the energy dissipation by the yield stress in the fluid. Therefore, the impact of yield stress on the bubble rise velocity was described more accurately. Compared with the traditional energy conservation analysis of bubbles, our model considered both the impact of bubble yield surface and the rheological parameters of the fluid on the flow field distribution in the shear region. Thus, the energy dissipation calculated from our model was more accurate compared to that calculated based on a linear velocity field. Bubble transport experiments were performed in this study. The results indicated that the bubble rise velocity decreased with increasing Bingham number. For Bn > 0.153, the bubble remained suspended in the testing fluid. The dimensionless analysis revealed that when the droplet was moving at a low velocity, both the Reynolds number and Froude number were relatively small. In such condition, the inertia of the fluid exhibited a slight impact on the

0 0.008

0.009

0.01 R/(m)

0.011

0.012

(b) Fig. 11. Comparison of bubble velocity between the proposed model, Stokes’s model (Hutter and Wang, 2016) and Rodrigue’s model (Rodrigue et al., 1999) predictions (a) Present experimental data(S3); (b) Lopez’s experimental data(S7).

force during the rise. Therefore, the velocity prediction model developed for viscous fluids generally yields a smaller result. A bubble immersed in a yield stress fluid is surrounded by a velocity field restricted by the yield stress. This distinctive velocity field affects the bubble transport, suspension, and deformation behaviors significantly (Dubash and Frigaard, 2007; Dimakopoulos et al., 2013). Therefore, it is a critical factor for analyzing bubble dynamics. When constructing the bubble velocity prediction model, we propose a new method to solve the velocity field around the bubble. When deriving the bubble velocity by coupling multiple variables, a velocity field satisfying the continuity equation can be obtained. Fig. 12 shows the distribution curve of the relative velocity of each fluid particle in the velocity field. The velocity distribution curve shows that the fluid particles beyond the yield surface are relatively still, while the fluid particles near the bubble wall are affected by the bubble transport. With increasing Bingham number, Bn, the maximum velocity of the fluid particles in the

336

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

bubble rise velocity, while the yield stress was crucial in the bubble rise velocity. Comparing the results obtained from the model and experiments, we found that the model developed in this study could provide an accurate prediction of bubble rise velocity in yield stress fluids for Re < 0.4. The model proposed in this study addressed the issue of poor compatibility at low Reynolds numbers for bubble velocity models developed for viscous fluids. Specifically, our model could provide a more accurate prediction of bubble velocity. In addition, when constructing the bubble velocity model, we proposed a submodel for solving the velocity field around the bubble. Such submodel could be used to analyze bubble transport and the suspension behavior in yield stress fluids. The yield stress caused the fluid to exhibit a higher characteristic viscosity compared to other non-Newtonian fluids. This feature resulted in a relatively lower Reynolds number for bubbles moving in the yield stress fluid. Considering that the bubble transport at low Reynolds numbers in yield stress fluids is a typical issue in many engineering applications, the model developed herein could be used widely in a variety of engineering practices. Declaration of interest The authors declared that there is no conflict of interest. Acknowledgments The work was supported by the National Natural Science Foundation-Outstanding Youth Foundation (51622405), the Shandong Natural Science funds for Distinguished Young Scholar (JQ201716), the Changjiang Scholars Program (Q2016135), the Construction Project of Taishan Scholars, and the Program for Changjiang Scholars and Innovative Research Team in University (IRT_14R58). Appendix A. Equation discretization Simpson’s formula is used to solve the integral equation. The computational intervals [Rb, Rsc] are equally divided. The meshing is shown in Fig. A.1. The approximate interpolation calculation is carried out for each component in each grid.

Fig. A1. Schematic diagram of mesh generation in shear region.

Continuity equation numerical integral equation:

Ql ¼

N  hr X pðri þ riþ1 Þ U z;i þ 4U z;iþ1=2 þ U z;iþ1 6 i¼0

ð29Þ

Viscous force energy dissipation numerical integral equation:

Wv ¼

nþ1 N r þ r 2 dU dU z;iþ1=2 dU z;iþ1 hr hh X z;i pk i iþ1 þ4 þ 18 i¼0 2 dri dr iþ1=2 dr iþ1 

N   X nþ2 nþ2 nþ2 sin hi þ 4sin hiþ1=2 þ sin hiþ1

ð30Þ

i¼0

The numerical integral equation for energy dissipation of yield stress:

Wy ¼

 N r þ r 2 dU dU z;iþ1=2 dU z;iþ1 hr hh X z;i psy i iþ1 þ4 þ 18 i¼0 2 dr i dr iþ1=2 dr iþ1 

N   X 2 2 2 sin hi þ 4sin hiþ1=2 þ sin hiþ1

ð31Þ

i¼0

where hr is the calculation step of radius integral, hh is the calculation step of angle integral, computation by the following formula:

8 r dp C 1 1 sy > > i þ þ Rb 6 r 6 Rx dU z;i < 2 dz k ri k ¼ > r i dp C 3 1 sy dr i > :  þ Rx < r 6 Rsc 2 dz k ri k hr ¼

Rsc  Rb p ; hh ¼ N N

ð32Þ

ð33Þ

Appendix B. Solution procedure The bubble velocity model consists of the energy conservation equation, continuity equation, and velocity field equation around the bubble. All seven equations in the model are interrelated with each other. Therefore, it is challenging to derive an analytical solution to the equation set. Simpson quadrature formula is used to solve the equation by numerical integration. The process of solution is as follows (Fig. B.1). (1) The rheological parameters of bubble equivalent radius Rb and yield stress fluid are determined: k, sy, n. Set the initial trial value of bubble rising velocity Ub and the initial trial value of velocity extreme point position Rx in shear zone. The shear area is meshed and the number of meshes N is determined. (2) According to the basic parameters in step (1), the boundary 0 condition Uz and yield interface position Rsc are calculated by using Eqs. (14) and (15). (3) The equations of velocity field in shear region are estab0 lished by Uz and Rsc. The data of velocity field in shear region corresponding to current Rx value and all unknowns in the equation are solved. (4) Combining with the velocity field equation in shear region, the flow QL of fluid at equatorial plane is calculated by Simpson formula. (5) If the continuity equation satisfies the error of e1, then Rx is the exact value. Otherwise, Rx is iterated and steps (3) and (4) are repeated until the velocity field satisfying the continuity equation is solved. (6) The velocity field around bubbles is brought into the energy Eqs. (4)–(6), and the work values of buoyancy, yield stress and viscous force during bubbles rising are solved by Simpson formula.

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

337

Fig. B1. Flow chart for bubble velocity prediction. In this paper, N = 105, radius iteration step 4R = hr, velocity iteration step 4U = 0.1 mm/s, continuity equation error limit is e1 = 1%, energy conservation error limit is e2 = 5%.

(7) The energy conservation equation is judged. If the continuity equation satisfies the error of e2, Ub is the exact value. Otherwise, the bubble rising velocity Ub is iterated and the steps (2)–(6) are repeated until the energy conservation equation is satisfied and the bubble rising velocity is output.

References Ahonguio, F., Jossic, L., Magnin, A., 2014. Influence of surface properties on the flow of a yield stress fluid around spheres. J. Nonnewton. Fluid Mech. 206, 57–70. Astarita, G., Apuzzo, G., 1965. Motion of gas bubbles in non-Newtonian liquids. AIChE J. 11 (5), 815–820. Beris, A.N., Tsamopoulos, J.A., Armstrong, R.C., Brown, R.A., 1985. Creeping motion of a sphere through a Bingham plastic. J. Fluid Mech. 158, 219–244. Bingham, E.C., 1922. Fluidity and Plasticity, vol. 2. McGraw-Hill. Chhabra, R.P., 2006. Bubbles, Drops, and Particles in Non-Newtonian Fluids. CRC Press.

Curran, S.J., Hayes, R.E., Afacan, A., Williams, M.C., Tanguy, P.A., 2002. Properties of carbopol solutions as models for yield-stress fluids. J. Food Sci. 67 (1), 176–180. Dhiman, A.K., Chhabra, R.P., Eswaran, V., 2007. Heat transfer to power-law fluids from a heated square cylinder. Numer. Heat Transf., Part A: Appl. 52 (2), 185– 201. Dijkhuizen, W., Roghair, I., Annaland, M.V.S., Kuipers, J.A.M., 2010. DNS of gas bubbles behaviour using an improved 3D front tracking model—drag force on isolated bubbles and comparison with experiments. Chem. Eng. Sci. 65 (4), 1415–1426. Dimakopoulos, Y., Pavlidis, M., Tsamopoulos, J., 2013. Steady bubble rise in Herschel-Bulkley fluids and comparison of predictions via the augmented Lagrangian method with those via the Papanastasiou model. J. Nonnewton. Fluid Mech. 200, 34–51. Duangsuwan, W., Tüzün, U., Sermon, P.A., 2009. Configurations and dynamics of single air/alcohol gas–liquid compound drops in vegetable oil. Chem. Eng. Sci. 64 (13), 3147–3158. Duangsuwan, W., Tüzün, U., Sermon, P.A., 2010. Feasibility of N2/sunflower oil compound drop formation in methanol induced by bubble train. AIChE J. 56 (12), 3274–3278. Duangsuwan, W., Tüzün, U., Sermon, P.A., 2011. The dynamics of single air bubbles and alcohol drops in sunflower oil at various temperatures. AIChE J. 57 (4), 897– 910.

338

Z. Wang et al. / Chemical Engineering Science 201 (2019) 325–338

Dubash, N., Frigaard, I., 2004. Conditions for static bubbles in viscoplastic fluids. Phys. Fluids 16 (12), 4319–4330. Dubash, N., Frigaard, I.A., 2007. Propagation and stopping of air bubbles in Carbopol solutions. J. Nonnewton. Fluid Mech. 142 (1–3), 123–134. Filip, P., David, J., 2003. Axial Couette-Poiseuille flow of power-law viscoplastic fluids in concentric annuli. J. Petrol. Sci. Eng. 40 (3–4), 111–119. Herschel, W.H., Bulkley, R., 1926. Konsistenzmessungen von gummibenzollösungen. Colloid Polym. Sci. 39 (4), 291–300. Hutter, K., Wang, Y., 2016. Creeping motion around spheres at rest in a newtonian fluid. In: Fluid and Thermodynamics. Springer, Cham, pp. 1–45. Islam, M.T., Ganesan, P., Cheng, J., 2015. A pair of bubbles’ rising dynamics in a xanthan gum solution: a CFD study. RSC Adv. 5 (11), 7819–7831. Liu, Y.Q., Zhu, K.Q., 2010. Axial Couette-Poiseuille flow of Bingham fluids through concentric annuli. J. Nonnewton. Fluid Mech. 165 (21–22), 1494–1504. Lopez, W.F., Naccache, M.F., de Souza Mendes, P.R., 2018. Rising bubbles in yield stress materials. J. Rheol. 62 (1), 209–219. Mirzaagha, S., Pasquino, R., Iuliano, E., D’Avino, G., Zonfrilli, F., Guida, V., Grizzuti, N., 2017. The rising motion of spheres in structured fluids with yield stress. Phys. Fluids 29 (9), 093101. Mitsoulis, E., Tsamopoulos, J., 2017. Numerical simulations of complex yield-stress fluid flows. Rheol. Acta 56 (3), 231–258. Mougin, N., Magnin, A., Piau, J.M., 2012. The significant influence of internal stresses on the dynamics of bubbles in a yield stress fluid. J. Nonnewton. Fluid Mech. 171, 42–55. Nogueira, S., Riethmuler, M.L., Campos, J.B.L.M., Pinto, A.M.F.R., 2006. Flow in the nose region and annular film around a Taylor bubble rising through vertical columns of stagnant and flowing Newtonian liquids. Chem. Eng. Sci. 61 (2), 845–857. Ohta, M., Kikuchi, D., Yoshida, Y., Sussman, M., 2011. Robust numerical analysis of the dynamic bubble formation process in a viscous liquid. Int. J. Multiph. Flow 37 (9), 1059–1071. Papanastasiou, T., Georgiou, G., Alexandrou, A.N., 1999. Viscous Fluid Flow. CRC Press. Peebles, F.N., 1953. Studies on the motion of gas bubbles in liquid. Chem. Eng. Prog. 49 (2), 88–97. Piau, J.M., 2007. Carbopol gels: elastoviscoplastic and slippery glasses made of individual swollen sponges: meso-and macroscopic properties, constitutive equations and scaling laws. J. Nonnewton. Fluid Mech. 144 (1), 1–29. Pinto, G.H.V.P., Martins, A.L., Rocha, J.M.S., Martinelli, A.E., 2012. New methodology for gas migration prediction before oilwell cementing. Braz. J. Petrol. Gas 6 (2). Premlata, A.R., Tripathi, M.K., Karri, B., Sahu, K.C., 2017. Dynamics of an air bubble rising in a non-Newtonian liquid in the axisymmetric regime. J. Nonnewton. Fluid Mech. 239, 53–61.

Rodrigue, D., De Kee, D., Man Fong, C.F.C., 1999. A note on the drag coefficient of a single gas bubble in a power-law fluid. Can. J. Chem. Eng. 77 (4), 766–768. Samson, G., Phelipot-Mardelé, A., Lanos, C., Pierre, A., 2017. Quasi-static bubble in a yield stress fluid: elasto-plastic model. Rheol. Acta 56 (5), 431–443. Santos, O.L., Azar, J.J., 1997. A study on gas migration in stagnant non-Newtonian fluids. Latin American and Caribbean Petroleum Engineering Conference. Society of Petroleum Engineers. Sikorski, D., Tabuteau, H., de Bruyn, J.R., 2009. Motion and shape of bubbles rising through a yield-stress fluid. J. Nonnewton. Fluid Mech. 159 (1–3), 10–16. Shu, S., Yang, N., 2013. Direct numerical simulation of bubble dynamics using phase-field model and lattice Boltzmann method. Ind. Eng. Chem. Res. 52 (33), 11391–11403. Sousa, R.G., Riethmuller, M.L., Pinto, A.M.F.R., Campos, J.B.L.M., 2005. Flow around individual Taylor bubbles rising in stagnant CMC solutions: PIV measurements. Chem. Eng. Sci. 60 (7), 1859–1873. Sun, B., Guo, Y., Wang, Z., Yang, X., Gong, P., Wang, J., Wang, N., 2015. Experimental study on the drag coefficient of single bubbles rising in static non-Newtonian fluids in wellbore. J. Nat. Gas Sci. Eng. 26, 867–872. Terasaka, K., Tsuge, H., 2001. Bubble formation at a nozzle submerged in viscous liquids having yield stress. Chem. Eng. Sci. 56 (10), 3237–3245. Terasaka, K., Shibata, H., 2003. Oxygen transfer in viscous non-Newtonian liquids having yield stress in bubble columns. Chem. Eng. Sci. 58 (23–24), 5331–5337. Tripathi, M.K., Sahu, K.C., Karapetsas, G., Matar, O.K., 2015. Bubble rise dynamics in a viscoplastic material. J. Nonnewton. Fluid Mech. 222, 217–226. Tsamopoulos, J., Dimakopoulos, Y., Chatzidai, N., Karapetsas, G., Pavlidis, M., 2008. Steady bubble rise and deformation in Newtonian and viscoplastic fluids and conditions for bubble entrapment. J. Fluid Mech. 601, 123–164. Van Hout, R., Gulitski, A., Barnea, D., Shemer, L., 2002. Experimental investigation of the velocity field induced by a Taylor bubble rising in stagnant water. Int. J. Multiph. Flow 28 (4), 579–596. White, F.M., Corfield, I., 2006. Viscous Fluid Flow. McGraw-Hill, New York. Yin, J.Y., Zhou, Y.C., Zhang, H., Jang, H.W., Zhu, L., 2016. Bubble Rise Velocity in Herchel-Bulkley Fluid o f Drilling Anm ilus. J. Southwest Petrol. Univ. (Sci. Technol. Ed.) 38 (3), 135–143. Yu, Z., Fan, L.S., 2010. Multirelaxation-time interaction-potential-based lattice Boltzmann model for two-phase flow. Phys. Rev. E 82 (4), 046708. Zenit, R., Feng, J.J., 2018. Hydrodynamic interactions among bubbles, drops, and particles in non-newtonian liquids. Annu. Rev. Fluid Mech., 50 Zhang, J., 2012. An augmented Lagrangian approach to simulating yield stress fluid flows around a spherical gas bubble. Int. J. Numer. Meth. Fluids 69 (3), 731–746. Zhang, J., 2018. An enhanced mesh refinement strategy to capture yield surfaces around a spherical gas bubble. GSTF J. Math., Stat. Oper. Res. (JMSOR) 2 (2).