Effect of the wall presence on the bubble interfacial forces in a shear flow field

Accepted Manuscript

Effect of the Wall Presence on the Bubble Interfacial Forces in a Shear Flow Field Jinyong Feng , Igor A. Bolotnov PII: DOI: Reference:

S0301-9322(17)30494-9 10.1016/j.ijmultiphaseflow.2017.10.004 IJMF 2658

To appear in:

International Journal of Multiphase Flow

Received date: Revised date: Accepted date:

13 July 2017 18 September 2017 2 October 2017

Please cite this article as: Jinyong Feng , Igor A. Bolotnov , Effect of the Wall Presence on the Bubble Interfacial Forces in a Shear Flow Field, International Journal of Multiphase Flow (2017), doi: 10.1016/j.ijmultiphaseflow.2017.10.004

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights 

Effect of wall distance, bubble deformation and relative velocity on the interfacial forces are evaluated with level-set interface tracking method Full flow field is resolved without the need for interfacial closure laws: the

CR IP T



interfacial force coefficients are evaluated based on the momentum balance

The sign change of net transverse lift force is observed at certain conditions



A bubble migration map based on the simulation results is created

AC

CE

PT

ED

M

AN US



ACCEPTED MANUSCRIPT

Effect of the Wall Presence on the Bubble Interfacial Forces in a Shear Flow Field

CR IP T

Jinyong Feng, Igor A. Bolotnov*

Department of Nuclear Engineering, North Carolina State University, Raleigh, NC 27695, USA

AN US

[email protected], [email protected]

Abstract

Understanding the interactions between the bubbles and solid walls is important for

M

predicting how the bubbly flow characteristics, such as void fraction and statistical distribution of bubble sizes, are affected by the presence of solid boundaries. In the present

ED

research, the numerical simulations with interface tracking method are performed to evaluate the interfacial forces‟ closures in shear flow field under well-controlled conditions.

PT

The proportional-integral-derivative bubble controller is adopted to keep a single bubble at

CE

fixed position near the wall by balancing the forces in all three directions. This allows the extraction of the interfacial forces and the calculation of the corresponding coefficients.

AC

Validation tests have previously been conducted by the authors to compare the interfacial closures‟ results against the experimental data. In this work, the influence of wall distance, relative velocity and bubble deformation on the interfacial forces are investigated. The

*

Corresponding author. Tel: +1 518 542 8939 E-mail address: [email protected]

ACCEPTED MANUSCRIPT

results show that as the bubble approaches the wall, it experiences a transition from attractive force to repulsive force normal to the wall. The bubble Reynolds number also has significant impact on the interfacial forces: the net lateral force changes direction at

CR IP T

𝑅𝑒𝑏 = 35. The coupled phenomenon of bubble deformation and wall effect is investigated as well. The net lift sign change is observed at 𝐸𝑜 = 2.3 when the dimensionless wall distance (the ratio between the distance from the bubble center to the wall and the bubble diameter) is equal to 1. The presented studies complement the state-of-the-art knowledge

AN US

and those are aimed to help improve the closure laws for interfacial forces and contribute to the multiphase computational fluid dynamics (M-CFD) closure model development.

Keywords

AC

CE

PT

ED

M

Interfacial forces; wall distance; deformable bubble; shear flow

ACCEPTED MANUSCRIPT

1. Introduction The two-phase bubbly flows have wide applications in the engineering fields due to the outstanding characteristics of the momentum and heat transfer efficiency. The bubble rise

CR IP T

velocity and bubble volume fraction (or void fraction) distribution is influenced by the presence of solid boundaries. Therefore, it is important to understand the bubble migration behavior near the wall and the interfacial forces acting on the bubble. However, the physics behind the bubbly flows behavior near the wall can be very complex due to the various

AN US

parameters and phenomena involved, such as surface contamination (Takemura and Magnaudet, 2003), bubble deformation, relative velocity and wall distance. Numerous experimental and analytical work (Saffman G., 1965; McLaughlin, 1991; Vasseur and Cox, 1977; Cox and Hsu, 1977; Clift et al., 2005; Ambari et al., 1983; Drew and Lahey, 1987a)

M

are reported on the interfacial forces acting on the solid sphere, while only a few studies

ED

analyzed the bubbles‟ motion near the wall. Antal et al. (Antal et al., 1991) proposed the first wall force model based on the theoretical derivation for the flow past two cylinders.

PT

The void fraction distribution and liquid velocity profile agreed quite well with the trends found experimentally (Nakoryakov et al., 1986; Gorelik et al., 1987), while the prediction

CE

of zero void fraction positions are different. Takemura et al. (Takemura et al., 2002; Takemura and Magnaudet, 2003) experimentally investigated the spherical bubble rising

AC

near a plane vertical wall in a quiescent liquid with different bubble Reynolds numbers, 𝑅𝑒𝑏 . They found that at 𝑅𝑒𝑏

35, the net transverse lift force acting on the bubble was

directed away from the wall and it could be solved analytically based on the Oseen approximation (Vasseur and Cox, 1977). However, for higher 𝑅𝑒𝑏 , the bubble tended to migrate toward the wall. They explained that the distribution of the vorticity around the

ACCEPTED MANUSCRIPT

bubble was responsible for the weak interactions observed for Reynolds numbers of (10). For the high Reynolds numbers where the inertia effects are dominant, the interaction with the wall results in an attractive lift force. At high 𝑅𝑒𝑏 , the bouncing process of the bubbles

CR IP T

moving sufficiently close to the wall may be due to the competition between these two mechanisms.

Thanks to the increasing computer power and advanced algorithms, numerical simulation of multiphase flows, especially computational fluid dynamics (CFD) approach,

AN US

emerges as an effective tool for exploring the bubbly flow. Due to the wide range of spatial and temporal scales in the industrial size systems, it is virtually impossible to capture all the details of the flow field with the current available computational resources. Depending on the scales, three approaches are mostly used to simulate bubbly flows: the Eulerian-

M

Eulerian (E-E) approach, the Eulerian-Lagrangian (E-L) approach and the interface

ED

tracking methods (ITM) which is often coupled with direct numerical simulation (DNS) or large eddy simulation (LES) for turbulence modeling. DNS/ITM typically adopts the single

PT

fluid formulation where the two-phase flow is simulated by solving a single set of NavierStokes equations with spatially variable fluid properties. The three approaches have their

CE

own advantages and disadvantages and their specific range of applicability. In the E-E

AC

approach, which is also referred as the two-fluid model (Ishii and Mishima, 1984; Podowski, 2009; Lahey and Drew, 2001), both phases are treated as continuum fluids. The ensemble-averaged mass and momentum conservation equations are used to describe the time-dependent motion of both phases. In case of adiabatic two-phase flow without phase change, the closure terms are the interfacial forces terms which are generally considered to include the drag force (Tomiyama et al., 1998; Ishii and Chawla, 1979), lift force

ACCEPTED MANUSCRIPT

(Tomiyama et al., 2002; Drew and Lahey, 1987b), wall force (Antal et al., 1991; Tomiyama et al., 1995), and turbulent dispersion force (Podowski, 2008; Lopez de Bertodano, Martin A, 1998; Lahey et al., 1993). At this level of E-E continuum models, bubbles lose their

CR IP T

discrete identity, which enables the simulation of relatively large systems and the study of various flow regimes, like bubbly flow, slug flow, churn-turbulent and annular flow, in a single channel. In the E-L approach, each bubble is separately tracked, while the liquid phase is treated as a continuum. The interaction between the bubbles and the liquid is

AN US

accounted through a source term (Tomiyama et al., 1998; Hosokawa et al., 2002) in the momentum equations. Since the E-E approach treats the gas bubble or particles as a continuum phase and the E-L approach treats the gas bubble as a non-deformable small spherical particle, both cannot be used for describing a deformable bubble behavior.

M

DNS is a useful tool to improve the current understanding of the local and instantaneous

ED

properties of the flow fields. The flow field is obtained by solving the governing equations numerically on sufficiently fine grids so that all the flow details in the turbulent flow or the

PT

bubble-induced turbulence in the laminar flow are fully resolved. DNS coupled with interface tracking method (ITM) has emerged as a reliable tool for engineers and scientists,

CE

providing valuable information on the temporal and spatial distribution of key flow variables in the two-phase flow fields (Feng and Bolotnov, 2017a; Fang et al., 2016; Lu and

AC

Tryggvason, 2008; Legendre et al., 2006). Although existing CFD codes can be applied at DNS scale with confidence to solve a variety of single-phase flow problems, considerable research efforts are still necessary to investigate the gas-liquid two-phase flows with the same level of confidence. Detailed knowledge on the behavior of single bubble or droplet in complex flow fields is limited. For example, even the behavior of a single air bubble

ACCEPTED MANUSCRIPT

rising in quiescent water is not yet completely understood: not only the physical properties like the density, viscosity, and surface tension (Tomiyama et al., 2002; Dijkhuizen et al., 2010) affect the behavior of the bubbles, but also small amounts of surface impurities

CR IP T

(Fukuta et al., 2008). The information gained by DNS can be employed to improve our understanding on the single bubble behavior and developing closure relations in the framework of E-E computations which allow the simulations of much larger engineering problems.

AN US

Zeng et al. (Zeng et al., 2005; Zeng et al., 2009) carried out DNS for a rigid sphere moving parallel to a plane wall in a quiescent fluid for 𝑅𝑒𝑏

3

using a spectral element

method. They observed that the lift coefficient steadily decreases on increasing Reynolds number up to about 100. With increasing Reynolds number, from 𝑅𝑒𝑏

, the lift coefficient decays with the wall distance while the decay rate steadily

M

𝑅𝑒

to about

𝑅𝑒𝑏

ED

increases. However, the lift coefficient increases dramatically for increasing 𝑅𝑒𝑏 when . Zaruba et al. (Zaruba et al., 2007) numerically investigated the motion of

PT

single bubble rising in an upward shear liquid flow near a vertical wall. Bubbles were found to slide along the wall when their diameter was small. Bubbles could also experience

CE

multiple collisions with the wall at certain experimental parameters. Vhora et al. (Vhora et

AC

al., 2013) investigated the wall effect on a single bubble using DNS coupled with interface tracking method. A spherical bubble was initially released with different inclination angles and migrated towards the wall. Due to the force balance between the wall repulsive force and the buoyancy force, the bubble reached a steady state and stayed at a fixed position away from the wall. The bubble‟ trajectories agreed with the experimental results. Lucas et

ACCEPTED MANUSCRIPT

al. (Lucas et al., 2007) performed the full 3D simulations using the commercial code CFX for simplified monodisperse cases. In the validation of different models for the interfacial forces, a set of Tomiyama lift (Tomiyama et al., 2002) and wall force (Tomiyama, 1998)

the best agreement with the experimental data.

CR IP T

and Favre averaged turbulent dispersion force (Moraga et al., 2003) was found to provide

In the present work, the open-source flow solver, called PHASTA (Jansen, 1999; Jansen, 1993), is used to evaluate the interfacial force closures in laminar flow field under well-

AN US

controlled conditions. The proportional-integral-derivative (PID) bubble controller (Thomas et al., 2015) is adopted to keep a single bubble at fixed position by balancing the forces in all three directions under various laminar flow scenarios. In this setup, the interfacial forces can be extracted based on the force-balance equation and the

M

corresponding coefficients in the closure law expressions can be obtained. Validation tests

ED

have been conducted by the authors before (Feng and Bolotnov, 2017b; Thomas et al., 2015) to compare the drag and lift coefficients‟ results against the experimental data under

PT

laminar flow conditions. In this paper, the research scope is significantly expanded to

CE

include new scenarios and study the influence of wall presence on the interfacial forces.

2. Interfacial forces models

AC

Note, that this discussion is introduced to show how the presented results can inform the

model coefficients for the lower fidelity approach. The presented research does not utilize highly empirical lift and drag coefficients to obtain the results. For the E-E approach, the equations of the two-fluid model are derived by ensembleaveraging the local instantaneous equations (Drew A. and Passman L., 1998). Two sets of

ACCEPTED MANUSCRIPT

balance equations for mass and momentum are obtained. Ignoring the interfacial mass transfer term (and thus assuming adiabatic flow), the mass and momentum conservation equations have the following form:

(

𝑘 𝑘 ⃗⃗⃗⃗𝑘 )

where the indices

(

𝑘 𝑘)

𝑘 𝑘 ⃗⃗⃗⃗𝑘 ⃗⃗⃗⃗𝑘 )

(

=

𝑘

𝑘 ⃗⃗⃗⃗𝑘 )

(

𝑘

refers to the phase ( for liquid,

and stress tensor of phase , respectively;

𝑘 𝑘)

𝑘 𝑘

for gas); ⃗⃗⃗⃗𝑘 and

⃗⃗⃗⃗⃗𝑘

(2)

are the velocity

stands for the volume fraction of phase .

AN US

𝑘

(1)

=

CR IP T

(

The momentum exchange term, ⃗⃗⃗⃗⃗𝑘 in Eq. ( 2 ), describing the interfacial forces, needs to close and is given as a superposition of the following terms: ⃗⃗⃗⃗⃗ = ⃗⃗⃗⃗⃗𝑙

⃗⃗⃗⃗⃗⃗⃗⃗ 𝑙

⃗⃗⃗⃗ 𝑙

⃗⃗⃗⃗⃗⃗⃗ 𝑙

⃗⃗⃗⃗⃗⃗ 𝑙

(3)

M

⃗⃗⃗⃗𝑙 =

The interfacial force terms on the right hand side of Eq. ( 3 ) are: the drag force, ⃗⃗⃗⃗⃗𝑙 ; the

ED

⃗⃗⃗⃗ ⃗⃗⃗⃗⃗⃗⃗ virtual mass force, ⃗⃗⃗⃗⃗⃗⃗⃗ 𝑙 ; the lift force, 𝑙 ; the turbulent dispersion force, 𝑙 ; and the wall

PT

force, ⃗⃗⃗⃗⃗⃗ 𝑙 . Currently, it is generally agreed that the drag force is largely predominant over other distributions in the streamwise direction, and the lateral distribution of the bubble is

CE

associated with the lift force. The accuracy of the closure relations for the interfacial forces determines the prediction capabilities of the multiphase computational fluid dynamics

AC

approach for the dispersed two-phase flow. Despite considerable research efforts on the closure modeling (Clift et al., 2005; Ervin

and Tryggvason, 1997; Tomiyama et al., 1998; Hosokawa et al., 2002), accurate prediction of the interfacial forces remain an open question in numerical simulations of bubbly flows. For single bubble, its shape varies with the bubble size, liquid phase flow field, and the

ACCEPTED MANUSCRIPT

physical properties of the system. This turns the interfacial force closures into a complex function of the bubble/liquid physical properties, inflow conditions, wall effect, etc. In this paper, the focus is on the influence of wall presence on the interfacial forces, especially

CR IP T

drag and lift forces. This work also aims to expand the current database for interfacial closure development. The brief literature review of the drag, lift and wall forces are given in the following sections which are helpful to gain a complete understanding of the research

AN US

presented in this paper.

2.1 Drag force

In bubbly flow, the drag force plays a significant role on the flow pattern, especially in the streamwise direction. It determines the gas phase residence time and the terminal

M

velocity of the bubbles. For a single bubble, the drag force acting on the bubble is typically

ED

modelled as

=

is the drag coefficient,

PT

where

𝑙

is the bubble cross-sectional area, 𝑙

(4) is the relative

is the liquid density.

CE

velocity between gas and liquid, and

2

The origin of the drag force is the resistance experienced by a bubble moving in the

AC

liquid. Various drag closure laws have been proposed in the literatures in the last few decades (Hadamard, 1911; Moore, 1963; Tomiyama et al., 1998). One of the most widely accepted correlation is proposed by Tomiyama et al. (Tomiyama et al., 1998). This correlation is based on spherical bubbles in pure systems:

ACCEPTED MANUSCRIPT

=

[

6 ( 𝑅𝑒𝑏

. 5𝑅𝑒𝑏 .

)

48 ] 𝑅𝑒𝑏

(5)

where the pure systems are the two-phase flow systems without surfactant contamination. The existence of surfactant on the bubble surface will change the tangential shear stress on

CR IP T

the bubble surface and cause a variation of surface tension. A fully contaminated bubble behaves like a rigid particle in water. In our studies, we focus on the evaluation of interfacial forces acting on the bubble in pure system.

We developed the proportional-integral-derivative (PID) based controller (Thomas et al.,

AN US

2015) to numerically evaluate lift and drag forces acting on a single bubble in wellcontrolled conditions. The theory and implementation of PID bubble controller is described in section 3.3. Typical experimental studies estimate the interfacial forces by measuring the

M

velocity field and bubble trajectories which may introduce additional measurement uncertainties (Tomiyama et al., 1998; Tomiyama et al., 2002). Based on the extracted

ED

interfacial forces from the numerical simulation, we can precisely calculate the corresponding interfacial closures based on Eqs. ( 4 ) and ( 7 ). Based on the laminar

PT

(Thomas et al., 2015) and turbulent flow (Feng and Bolotnov, 2017a) data from DNS

CE

simulation, the authors proposed a modification of the Tomiyama‟s closure which is given in Eq. ( 6 ) and it has good predictions of drag coefficient for 𝑅𝑒𝑏 up to 900 shown in Fig.

AC

1.

=

[

6 ( 𝑅𝑒𝑏

. 5𝑅𝑒𝑏 .

)

48 ( 𝑅𝑒𝑏

3

𝑅𝑒𝑏 .

)]

(6)

AN US

CR IP T

ACCEPTED MANUSCRIPT

Fig. 1. Drag coefficient comparison between Tomiyama‟s correlation (Tomiyama et al., 1998) and DNS-based correlation (Feng and Bolotnov, 2017b).

M

2.2 Lift force

The lift force on the bubble originates from the velocity gradient of the liquid phase and

ED

relative velocity between the bubble and the liquid phase. It has a significant effect on the radial void fraction distribution in bubbly flows. Bubbles traveling through a shear flow

PT

(such as wall boundary layer) will experience lift force transverse to the direction of

CE

motion. The experimental investigations of Tomiyama et al. (Tomiyama et al., 1995) as well as numerical predictions of Ervin and Tryggvason (Ervin and Tryggvason, 1997) have

AC

shown that the direction of the lift force can change its sign if a substantial deformation of the bubble occurs. The lift force exerted on a single bubble is typically modeled as

where

⃗⃗⃗ =

is the lift coefficient,

vapor and liquid, i.e., ⃗⃗⃗⃗ = ⃗⃗⃗⃗

𝑏

𝑏 ⃗⃗⃗⃗

(

⃗⃗⃗⃗ )

is the bubble volume,

𝑙

(7) is the relative velocity between

⃗⃗⃗𝑙 . In the plane shear flow, the lift force is modeled as

ACCEPTED MANUSCRIPT

⃗⃗⃗ =

𝑏

(8)

𝑙

The lift force governs the transverse migration of bubbles in the liquid. The magnitude and direction of the lift force is closely related to the bubble properties, such as size and

CR IP T

deformability as well as liquid flow characteristics. Earlier experimental research (Zun, 1980; Liu J., 1993; Hibiki and Ishii, 1999) in two-phase upflow made observations that small bubbles tend to move towards the pipe wall which causes a wall-peak bubble distribution, whereas large bubbles are preferably found in the center which results in a

AN US

core-peak bubble distribution. For a spherical bubble in upflow, the lift coefficient

is

positive so that the lift force acts (Eq. ( 7 )) in the direction of decreasing liquid velocity, i.e., in the direction towards the channel wall.

M

One of the important experimental research efforts on the lift force evaluation has been carried out by Tomiyama et al. (Tomiyama et al., 2002). They used a linear shear velocity

ED

field in viscous liquids. The following dimensionless numbers were used to describe the bubble characteristics in the experiment: Morton number ( 𝑜) characterizing the shape of

PT

the bubbles in a surrounding fluid or continuous phase; Eotvos number (𝐸𝑜) defining the

CE

ratio of buoyancy to the surface tension forces, and the velocity gradient

of a simple

shear flow. These are defined as

AC

𝑜=

where

(

)

𝑙

𝑙

𝐸𝑜 =

(

𝑙

)

=|

|

(9)

𝑙

is the bubble diameter,

is the surface tension.

Tomiyama‟s experiments (Tomiyama et al., 2002) varied these parameters within the ranges of

5.5

𝑜

2.8, .39

𝐸𝑜

5.74, and

8.3

. The closure

ACCEPTED MANUSCRIPT

law for the lift coefficient proposed by Tomiyama is based on the modified Eotvos number as follows: ( .288 (𝐸𝑜 )

( . 2 𝑅𝑒𝑏 )

𝐸𝑜 = where (𝐸𝑜 ) is expressed as .

(𝐸𝑜 )) 4

(

5𝐸𝑜

)

.

59𝐸𝑜

𝐸𝑜 𝐸𝑜

4 .7

( 10 )

CR IP T

={

( 11 )

. 2 4𝐸𝑜

.474 .

is the

extended bubble diameter and it is defined as the maximum horizontal dimension of a is the spherical bubble diameter. In Eq. ( 10 ),

AN US

bubble shown in Fig. 2. In Eq. ( 5 ),

Tomiyama et. al. (Tomiyama et al., 2002) used modified Eotvos number, 𝐸𝑜 , for deformable bubble and they replaced spherical bubble diameter with extended bubble diameter,

, in their study. Fig. 2 shows a schematic illustrating the difference between

M

spherical bubble diameter and extended bubble diameter and the spherical and deformable

AC

CE

PT

ED

bubble migration direction in laminar shear upflow.

CR IP T

ACCEPTED MANUSCRIPT

AN US

Fig. 2. Illustration of the migration direction for spherical and deformable bubble. represents the lateral migration velocity of the bubble, and represent the streamwise velocity of the bubble and liquid, respectively. In our previous study (Feng and Bolotnov, 2017b), we conducted the numerical simulation with the same parameters as used in Tomiyama‟s experiment (Tomiyama et al., 2002) and observed the lift force sign change due to the bubble deformation. The transition

M

happened at 𝐸𝑜 = 5 which is close to the Tomiyama‟s prediction of 6. The difference may come from the uncertainties of experimental measurement and impurities of fluids.

ED

The general trend of the lift coefficient obtained from the simulation follows the

PT

Tomiyama‟s experiment data while the estimated lift coefficient from the simulation is

CE

slightly smaller (Feng and Bolotnov, 2017b).

2.3 Wall force

AC

When the bubbles migrate toward the wall, the bubble rise velocity and void fraction

distribution are affected by the solid boundary presence. The wall-induced force, often simply referred to as wall force, makes the interfacial force closures model more complex. The near wall behavior and bubble distribution play an important role in the multiphase

ACCEPTED MANUSCRIPT

flow. Therefore, it is important to understand the motion of bubbles near the wall, and moreover, to understand the hydrodynamic forces acting on the bubbles. Antal et al. (Antal et al., 1991) proposed the wall force model based on a potential flow

⃗⃗⃗⃗⃗⃗ = 2 [

wall,

)]

2

⃗

( 12 )

is the distance between the bubble center and the wall, ⃗ is the unit normal to the and

are the wall force coefficients. =

. 6

.

and 4

are defined as:

( 13 )

AN US

where

(

CR IP T

past two cylinders

= . 47

( 14 )

Tomiyama et al. (Tomiyama et al., 1995) pointed out that the coefficient

in Antal‟s

M

model vanishes for bubbles in two parallel walls and modified the Antal‟s model by 2 . By measuring the trajectories of single

bubble rising in glycerol-water solution ( 𝑜

𝑜=

ED

introducing the second term with respect to

2.8), they confirmed that the value of

PT

was zero and proposed the following wall force model ⃗⃗⃗⃗⃗ =

2

(

2

)

( 15 )

⃗

CE

Under a constant Morton number, the coefficient

should be a function of the Eotvos

AC

number, 𝐸𝑜, and deducted the following empirical correlation of ={

.

( .933𝐸𝑜 7𝐸𝑜 . 4

.79)

𝑜 𝑜 5

𝐸𝑜 𝐸𝑜

. 5 33

( 16 )

Hosokawa et al. (Hosokawa et al., 2002) extended the study of the Tomiyama et al.

(Tomiyama et al., 1995) and claimed that

is not only a function of bubble deformation

ACCEPTED MANUSCRIPT

but also the bubble Reynolds number 𝑅𝑒𝑏 . They proposed the following correlation of for 𝑅𝑒𝑏

6

𝑜

7 𝑅𝑒𝑏 .

. 2 7𝐸𝑜}

under the condition of =

{

2.5 and 2.

𝐸𝑜

. . ( 17 )

CR IP T

It is important to investigate the combined effects of lift and wall forces. Sugioka and Tsukada (Sugioka and Tsukada, 2015) observed that the direction of net transverse lift force near the wall changes at bubble Reynolds number of 30. Takemura and Magnaudet

AN US

(Takemura and Magnaudet, 2003) had similar observation that the lift force acting on clean bubble is directed away from the wall force with 𝑅𝑒𝑏

35 and toward it for higher 𝑅𝑒𝑏

while the contaminated bubbles are always repelled by the wall. Additionally, the current study on the effect of the wall presence on the interfacial forces mainly focuses on the

M

influence of the transverse migration. Limited work (Sugioka and Tsukada, 2015) has been reported about the wall effect on the drag force which will determine the bubble residence

ED

time and terminal velocity. The PID bubble controller enables the extraction of the

PT

transverse force experienced by the bubble which consists of both lift force and wall force. The summed transverse force,

𝑏 ⃗⃗⃗⃗

(

⃗⃗⃗⃗ )

𝑙

= ⃗⃗⃗

⃗⃗⃗⃗⃗

( 18 )

CE

⃗⃗⃗⃗ =

, is expressed as follows:

where

is the net transverse lift coefficient where both wall force and lift force are

AC

considered. Since the value of the lift force is well-known, the analysis of summed transverse force is equivalent to the approach of isolating the wall force. Utilizing the PID bubble controller, we can apply the modifications to the interfacial

force closures. The presented study on the wall effect will help understand the mechanism of the wall force as a function of the wall distance and bubble Reynolds number.

ACCEPTED MANUSCRIPT

3. Numerical approach The simulations are performed with a research CFD-solver, PHASTA (Jansen, 1999). PHASTA is a Parallel, Hierarchic, higher-order (from the 2nd to the 5th order accuracy,

CR IP T

depending on the choice of function), Adaptive, Stabilized (finite element method-based (FEM)) Transient Analysis flow solver for both incompressible and compressible flows. The flow solver has been proven to be an effective tool for a multitude of diverse types of simulations including, Reynolds-Averaged Navier-Stokes (RANS), Large-Eddy Simulation

AN US

(LES), Detached-Eddy Simulation (DES), and DNS. Various benchmarks (Behafarid et al., 2013; Bolotnov et al., 2011; Jansen, 1993; Feng and Bolotnov, 2017a; Fang et al., 2017; Bolotnov et al., 2011; Jansen, 1993; Feng and Bolotnov, 2017a; Li et al., 2017; Bolotnov et

M

al., 2011; Jansen, 1993; Feng and Bolotnov, 2017a) have verified and validated the PHASTA code for the simulation of turbulent flows. PHASTA and its advanced features,

ED

such as anisotropic adaptive algorithms (Sahni et al., 2006) and the LES/DES models (Tejada-Martínez and Jansen, 2006) have been extended to two-phase flow simulation

PT

where the level-set interface tracking method is adopted to track the boundary between two

CE

immiscible fluids (Nagrath et al., 2006; Sethian, 1999; Sussman et al., 1999).

AC

3.1 Governing equations The spatial and temporal discretization for the Incompressible Navier-Stokes (INS)

equations used in PHASTA have been described in (Whiting and Jansen, 2001). The strong form of the INS equations is given by: Continuity

𝑖𝑖

=

( 19

ACCEPTED MANUSCRIPT

) ( 20

where

𝑖𝑡

is the density,

𝑖

viscous stress tensor, and

is the 𝑖

𝑡

𝑗 𝑖𝑗

=

𝑖

𝑖𝑗 𝑗

component of velocity,

is the body force along the

𝑡

𝑖

is the static pressure,

) 𝑖𝑗

is the

CR IP T

Momentum

coordinate, including the surface

tension force and gravitational force. The viscous stress tensor for an incompressible flow of a Newtonian fluid is related to the fluid‟s dynamic viscosity, , and the strain rate tensor, as: 𝑖𝑗

=2

AN US

𝑖𝑗 ,

𝑖𝑗

= (

𝑖𝑗

𝑗 𝑖)

( 21 )

The level-set interface tracking method has the superior capability in capturing interface quantities. Based on the level-set distance field ( ), the curvature of gas-liquid interface

M

can be computed as

( )=

ED

(

|

|

)

( 22 )

𝑡

=

( |

|

)

CE

PT

The volumetric force due to the surface tension is then represented by

where

( 23 )

is the surface tension coefficient. We utilize the Continuum Surface Force (CSF)

AC

model proposed by Brackbill et al. (Brackbill et al., 1992) which is formulated to numerically model the surface tension effects and it is included in 𝑖 .

ACCEPTED MANUSCRIPT

3.2. Level-Set interface tracking method To distinguish between the liquid phase and gas phase, the level set method (Sussman et al., 1999; Sussman and Fatemi, 1999; Sussman et al., 1998; Sethian, 1999) is utilized in this

value of a smooth function, , where

CR IP T

research for interface tracking. This involves modeling the interface as the zero level set is often called the first scalar and it represents the

signed distance from the interface. Thus, the interface is defined by

= . In case of no

phase change in the flow the scalar, , is advected with a moving fluid according to ⃗

=

AN US

=

( 24 )

where ⃗ is the flow velocity vector. The liquid phase is indicated by a positive level set. Fig. 3 shows the well-preserved level-set contours for a near wall bubble at steady state

CE

PT

ED

M

where the thicker contour represents the liquid/vapor interface.

AC

Fig. 3. Level-set contours of a near wall bubble. The thicker white line represents the liquid/gas interface. The thinner white lines represent the level-set contours where the spacing corresponds to the interface half-thickness.

Since evaluating the jump in physical properties using a step change across the interface

is numerically challenging, the properties near the interface are defined using a smoothed Heaviside kernel function, 𝐻𝜀 .

ACCEPTED MANUSCRIPT

𝐻𝜀 ( ) = {

2

[

(

)]

| |

( 25 )

where is the interface half-thickness. The fluid properties are then defined as 𝐻𝜀 ( )

(

𝐻𝜀 ( ))

( 26 )

( )=

𝐻𝜀 ( )

(

𝐻𝜀 ( ))

( 27 )

CR IP T

( )=

Although the presented research focus on the interfacial forces under laminar flow scenarios, the bubble-induced turbulence and local shear in the liquid may distort the level

AN US

set contours. Thus, the level set must be corrected with a re-distancing operation (Thomas et al., 2015; Sussman and Fatemi, 1999; Sussman et al., 1999). Additionally, the level-set interface tracking method is known for mass conservation issues. Previously (Fang, 2016),

M

we developed a methodology to ensure statistically constant void fraction for adiabatic twophase flow simulations. When the expected void fraction is known, a simplified PID

ED

controller is used to compare the current void fraction with the prescribed value, and then manipulate the advection velocity ⃗ in Eq. ( 24 ) accordingly to preserve the total volume

PT

fraction. In all the simulations, the void fraction is preserved and the averaged void fraction

CE

variation for each case,

∆𝛼 𝛼

, is within 0.3% from the target value.

AC

3.3 Proportional-Integral-Derivative bubble controller The proportional-integral-derivative (PID) bubble controller is implemented to the

PHASTA code to balance the interfacial forces, like drag and lift force (Thomas et al., 2015). Drag and lift forces originate from the liquid flow field around the bubble which are not explicitly modeled using theoretical correlations in our simulations. Their effect is

ACCEPTED MANUSCRIPT

directly resolved through high resolution interface tracking simulations and then evaluated using PID approach. The existence of the interfacial forces causes the bubble‟s migration in both streamwise and lateral directions. The PID bubble controller will keep the bubble at

CR IP T

fixed position by introducing artificial control forces and balancing the interfacial forces shown in Fig. 4.

The control forces can be extracted from the simulation and used to compute lift and drag coefficients as has been previously demonstrated (Thomas et al., 2015). In our studies,

AN US

the streamwise control force is applied to balance the drag force which is identical to the buoyant force. In this way, the gravity is a part of our solution. By selecting the appropriate relative velocity, we are able to simulate the realistic experiment under earth gravity conditions and compare the interfacial closures given in Eqs. ( 6 ) and ( 10 ).

M

The solution of the Navier-Stokes equations will provide the bubble‟s instantaneous

ED

location. The PID bubble controller will calculate the difference, i.e., the error, between the proposed bubble position and the instantaneous position. The PID controller consists of

PT

three error feedback mechanisms: the proportional part which represents the current error; the integral part which represents the accumulation of past errors; the derivative part which

CE

represents the prediction of future error. The error output for the bubble control algorithm is treated as the distance between the bubble‟s current position and its initial position. For

AC

each simulated time step, the bubble controller assesses the bubble location and velocity and uses those variables as inputs for control loop feedback. The computed control forces will be explicitly added into the right hand side body force term of Eq. ( 20 ). The mathematical form of the PID-based bubble controller is given below (Thomas et al., 2014)

ACCEPTED MANUSCRIPT

(𝑛) 𝑖

where of the

𝑡

is the

𝑡

=𝑐

̅̅̅̅̅̅̅ (𝑛) 𝑖

𝑐 [

(𝑛) 𝑖

𝑐

(𝑛) 𝑖

(𝑛)

𝑐

𝑖

component of the control force at time ,

component of control force at time n,

position difference between time n and time zero,

(𝑛) 𝑖 (𝑛) 𝑖

is the

average velocity at time n.

̅̅̅̅̅̅̅ (𝑛) 𝑖

𝑡

is the

( 28 (𝑛) 𝑐5 𝑖 ]

)

is historical average

component of the bubble

CR IP T

(𝑛+ ) 𝑖

𝑡

component of the bubble‟s

The application and implementation of PID bubble controller is an extensive learning

AN US

process. Through iterative test-trials, we are able to find the appropriate control coefficients to keep the bubble fixed. Note that as long as the bubble well-controlled, specific coefficients used are not expected to influence the behavior of the overall simulation. Table

M

1 provides the values of the control coefficients. This single set of control coefficients are valid for all the cases we discussed in this paper.

ED

Table 1. Summary of control coefficients used throughout this paper. 𝑐 0.2

𝑐 0.8

𝑐 5.

𝑐 5.

𝑐5 5.

PT

Control coefficients Value

CE

4. Results and discussion

The numerical simulations are performed with the PHASTA solver. In our previous

AC

work (Feng and Bolotnov, 2017b), the interfacial force under the two-phase shear flow scenarios was evaluated with the exactly same setup of Tomiyama‟s experiment (Tomiyama et al., 2002) with .67

𝐸𝑜

= 3.8

and .54

𝐸𝑜

5.96 which corresponds to

.34. Good agreement of the bubble shape between the experiment and the

PHASTA simulation was obtained (Feng and Bolotnov, 2017b). The sign change of lift

ACCEPTED MANUSCRIPT

force happened at 𝐸𝑜 = 5 was observed. For comparison, Tomiyama et al. (Tomiyama et al., 2002) predicted the transition happens at 𝐸𝑜 = 6. The transverse migration direction change is related to the presence of a slanted wake behind the deformable bubble which is

CR IP T

caused by the interaction between the wake and shear field. Those studies (Feng and Bolotnov, 2017b) validated the feasibility of CFD solver coupled with interface tracking method on the prediction of interfacial forces acting on single bubble in laminar flow condition without wall presence. In this paper, we will

AN US

significantly expand the scope of the previous study by analyzing the effect of wall presence, bubble deformation and the bubble Reynolds number over a wide range. Part of the results have been reported in the conference (Feng and Bolotnov, 2017c). Those quantities are varied within the range of .625

2.5 and

.7

𝑅𝑒𝑏

43.3 .

M

represents the dimensionless wall distance between the bubble center and the solid wall

( )=

ED

boundary. The constant shear flow shown in Fig. 4, with velocity profile described by , is applied as the inflow conditions. Constant velocity on the top and

PT

bottom surface in wall-normal direction and symmetric boundary conditions in spanwise

CE

direction are also applied here. The zero pressure condition is applied on the outflow surface. The summary of the case setup about the fluid properties, the geometrical and

AC

mesh configuration is given in Table 2.

ACCEPTED MANUSCRIPT

Table 2. Summary of the case setup

PT

ED

M

AN US

CR IP T

Geometrical configuration Dimensional quantities Non-dimensional quantities Bubble diameter (mm) 3.52 1 Bubble diameter ( ) 3 Simulation domain (mm ) 8 8 8 Simulation domain ( ) 5. 5. 5. Mesh configuration Number of elements across bubble diameter 45 Total mesh size (million) 1.87 Fluid properties Liquid Gas 1154 1.161 Density ( ) 0.019 1.845 10-5 Viscosity ( ) 0.04~0.12 Surface tension ( )

CE

Fig. 4. Initial shear velocity profile and the mesh configuration of the simulation domain. (a) shows the force balance using PID bubble controller where , , , and represent the drag force, lift force, streamwise control force and wall-normal control force, respectively. (b) shows the mesh configuration with refinement region near the bubble.

AC

Fig. 4 provides the mesh configuration of the entire simulation domain with the initial bubble shape and location shown as a white surface. To capture the curvature behavior of a very deformable bubble, the unstructured meshing instead of structured Cartesian mesh approach is used. The advantage of using the unstructured meshing capability (Jansen, 1993) is that a specific refinement region could be defined which enable the increase of the

ACCEPTED MANUSCRIPT

mesh resolution in the region containing the bubble while maintaining the total mesh size within the computational budget. For the wall distance study, the location of the refinement

CR IP T

region is adjusted depending on the bubble position.

4.1. Effect of wall distance

The wall distance plays a key role on the magnitude of wall effect on the interfacial forces. Thus, we separately place a single bubble at various distances from the top wall

setup can be found in Table 3.

AN US

while keeping the same relative velocity between liquid and gas. The details of the case

Table 3. Summary of the wall distance cases setup.

(s-1)

Dimensionless parameters 2.58/1.50 𝐸𝑜 𝑒 42.3 𝑅𝑒𝑏 0.625, 0.75, 0.875, 1, 1.25, 1.5, 2, 2.55

ED

3.8

M

Dimensional parameters 3.52 (mm) 0.198 (m/s)

PT

Fig. 5 shows the evolution of control forces versus time for different wall distances. The magnitude of the control force oscillation is small and periodic which is due to the

CE

implementation of the PID bubble controller. Since the oscillation of the control forces is inherent for the interface tracking method and the PID bubble controller, it is reproducible

AC

with different mesh configuration, timesteps and the domain length. However, the magnitude of the oscillation may be slightly different due to the numerical setup of the PID bubble controller. The existence of top wall increases the drag force and consequently increases the drag coefficient shown in Fig. 6. Thus, the wall presence is always found to increase the rise

ACCEPTED MANUSCRIPT

velocity relative to the liquid streamwise velocity. The drag coefficient for the case where the bubble is away from the wall is compared with the experiment-based correlation (Eq. ( 5 )). The 10% difference is observed because the Tomiyama‟s correlation is initially

CR IP T

designed for the spherical bubble. The high liquid viscosity system we used for our cases makes the bubble more deformable than it is supposed to be at the ordinary fluid properties at 25 oC.

As the bubble approaches the wall, lift force sign change is observed at

about 1.75

AN US

which indicates that the wall force influences the bubble migration direction. The sum of lift force and wall force pushes the bubble away from the wall. Fig. 6 shows that when the wall distance (

) is over 2, the existence of wall will not influence the solution. The

domain length in the transverse direction is over 5D, so the interfacial forces acting on the

M

bubble are not influenced by other surfaces. The choice of the simulation domain will not

ED

affect the results. Fig. 6(b) shows that simulation results generally follow Tomiyama‟s correlation, Eq. ( 15 ), for the small wall distance, ⁄

, which for the large wall

AC

CE

PT

distance the results agree with the Hosokawa‟s correlation, Eq. ( 16 ).

Fig. 5. Evolution of control forces versus different wall distances,

.

CR IP T

ACCEPTED MANUSCRIPT

AN US

Fig. 6. Drag and net lift coefficients versus different wall distances. (a) compares the drag coefficient based on the PHASTA data and Tomiyama‟s correlation (Eq. ( 5 )) and (b) compares the lift coefficient based on the PHASTA data, Hosokawa‟s correlation (Eq. ( 17 )), and Tomiyama‟s correlation (Eq. ( 17 )). In Fig. 7, the pressure field for selected wall distances is compared. Note that it is customary to use zero pressure boundary conditions in simulations since only local pressure gradient is important when solving incompressible Navier-Stokes equations. The pressure

M

field as shown in Fig. 7 is colored by the relative pressure and modified locally around the = .25,

ED

bubble due to its presence. For the bubble located away from the wall, about

the pressure field on the near side and far side of the bubble is quite symmetric. As the

PT

bubble approaches the wall, the interaction between the bubble and the wall causes the

CE

asymmetric pressure distribution which is responsible to the transverse migration sign change. The bubble is not axisymmetric when it is near the wall. It is also observed that the

AC

area of the negative pressure field on the bubble bottom increases as the wall distance decreases. Fig. 7 and Fig. 8 show that the flow field around the bubble, especially the wake region behind the bubble, is fully captured with our modeled domain size. The streamline pattern shown in Fig. 8 demonstrates that a strong vortex is formed below the bubble when the bubble is close to wall,

= .625. Considering the Bernoulli‟s equation (Bernoulli,

ACCEPTED MANUSCRIPT

= 𝑐𝑜

1738) in incompressible flow,

, stronger vortex formed below the

bubble indicates that higher pressure difference exists across the bubble in the wall-normal direction and the bubble experiences a repulsive force from the wall as it approaches the

CR IP T

wall. The “wrinkles” on the bubble surface is caused by the stretch of the level-set field which only occurs for highly deformable bubble. While the influence of the „wrinkles” on the interfacial forces may require further verification, it is expected to be minimal

M

AN US

considering the ratio of the “wrinkles” region to the entire bubble surface.

AC

CE

PT

ED

Fig. 7. Comparison of relative pressure field for different wall distance ( ⁄ = .625 . .25) where the white line represents the pressure equal to -5 Pa.

Fig. 8. Comparison of streamline pattern for different wall distances ( ⁄ = .625 . .25).

In Fig. 7, pressure discontinuities are observed near the wall which is due to the

transition from fine mesh to coarse mesh. An additional test is performed by extending the refinement region to the wall as shown in Fig. 9. Fig. 10 compares the pressure field for the baseline mesh (1.79 million elements) and the extended refinement mesh (2.23 million

ACCEPTED MANUSCRIPT

elements). Finer mesh configuration did generate smooth transition of the pressure field. However, a close comparison of the interfacial forces given in Fig. 11 demonstrates that the influence of the additional refinement near the wall will not influence the solution. The

CR IP T

interfacial forces for these two mesh setups are almost identical however additional 25%

ED

M

AN US

computational cost is needed for the extended refinement case.

AC

CE

PT

Fig. 9. Comparison of the baseline mesh (a) and extended refinement mesh (b).

Fig. 10. Comparison of the pressure field for baseline mesh (a) and extended refinement mesh (b).

CR IP T

ACCEPTED MANUSCRIPT

4.1 Effect of bubble Reynolds number

AN US

Fig. 11. Comparison of the control forces for the baseline mesh (red line) and extended refinement mesh cases (black line).

Experimental results (Takemura and Magnaudet, 2003; Tomiyama et al., 1995) indicated

M

that the wall force coefficient is a function of the bubble Reynolds number 𝑅𝑒𝑏 and the Eotvos number 𝐸𝑜. The study on the bubble Reynolds number effect on the interfacial

ED

force is presented. Since the Takemura and Magnaudet‟s experiment studied the spherical bubble, the spherical bubble shape is chosen in our study to provide consistent comparison. = , is used, but a series of different relative velocities are performed,

PT

For all the case,

CE

corresponding to the set of 𝑅𝑒𝑏 from 10.7 to 42.3. The surface tension value is adjusted for each case to have constant Weber number of 1.50. Weber number,

𝑒=

, is

AC

directly proportional to the relative velocity square and serves as a good measurement for the bubble shape when only the relative velocities are varied (Feng and Bolotnov, 2017a). The details about the case setup can be found in Table 4.

ACCEPTED MANUSCRIPT

Table 4. Summary of the bubble Reynolds number cases setup. Dimensionless parameters 2.58/1.50 𝐸𝑜 𝑒 10.7, 21.4, 32.1, 𝑅𝑒𝑏 37.4, 42.3 1

CR IP T

Dimensional parameters 3.52 (mm) 0.05, 0.1, 0.15, (m/s) 0.175, 0.198 -1 3.8 (s )

Fig. 12 shows the evolution of control forces for different bubble Reynolds numbers. The drag and net lift coefficients versus bubble Reynolds number are given in Fig. 13. In Fig. 13 (a), we compare the drag coefficient obtained from the PHASTA simulation and the

AN US

DNS-based drag coefficient closure given in Eq. ( 6 ) are compared. The discrepancy in high bubble Reynolds number region is attributed to the wall presence. From Fig. 13(b), it is observed that the sign change of lift coefficient occurs around 𝑅𝑒𝑏 of 35. A change of the lift force direction due to the bubble Reynolds number was also experimentally observed by

M

Takemura and Magdaudet (Takemura and Magnaudet, 2003) and the transition happens at

ED

𝑅𝑒𝑏 of 35. Fig. 14 compares the velocity field around the bubble with different relative velocities. By keeping the Weber number constant, the bubble shape is nearly identical for

PT

those cases. At low Reynolds number, the influence of wall presence dominates the bubble migration direction. However, when the inertia effects are dominant (the increasing high

CE

velocity field on the bubble bottom region shown in Fig. 14), the velocity field is

AC

essentially confined in the boundary layer and the wake of the bubble. This work shows that the influence of both wall distance and bubble Reynolds number leads to the bubble bouncing behavior near the wall. However, the bubble bouncing behavior near the wall is a complex phenomenon where other parameters, like bubble deformation, will also play an important role.

CR IP T

ACCEPTED MANUSCRIPT

PT

ED

M

AN US

Fig. 12. Evolution of control forces versus bubble Reynolds number where the case setup follows Table 4.

AC

CE

Fig. 13. Drag and net lift coefficients versus 𝑅𝑒𝑏 . The drag coefficient from PHASTA data (solid line) is compared with the DNS informed correlation (dash line), Eq. ( 6 ).

CR IP T

ACCEPTED MANUSCRIPT

Fig. 14. Comparison of velocity field for different bubble Reynolds numbers. The range of legend for each figure is based on the velocity profile, = 3.8 .

AN US

4.3 Effect of bubble deformation

The lift sign change due to the bubble deformation has been reported at 𝐸𝑜 = 3.4 in the previous study (Feng and Bolotnov, 2017b) without considering wall effect. The interaction between the bubble and the wake region causes the different migration behaviors of the

M

spherical and deformable bubble. In two-phase upflow, spherical bubble tends to migrate

ED

toward the pipe wall whereas the deformable bubble tends to migrate toward the pipe center. As the spherical bubbles accumulate and coalescence near the wall, the large

PT

deformable bubbles are formed and then they migrate toward the channel center. In this section, the bubble is placed at fixed position,

= , with different bubble deformation

CE

levels to analyze the integral effect of bubble deformation and wall distance on the

AC

interfacial forces. The summary of the case setup is given in Table 5. Table 5. Summary of the bubble deformation cases setup.

Dimensional parameters

(mm)

3.52

(m/s) (s-1)

0.198 3.8

Dimensionless parameters 0.79, 1.76, 2.79, 𝐸𝑜 3.46, 5.75 42.3 𝑅𝑒𝑏 1

ACCEPTED MANUSCRIPT

Fig. 15 compares the control forces for different bubble deformation levels. The control forces reach a steady state and the observed oscillation is attributed to the periodic variation of the void fraction. As shown in Fig. 16(a), the drag coefficient increases with the bubble

CR IP T

deformation because the bubble shape is extended in the wall-normal direction which blocks the flow field and introduces higher drag force. Fig. 16(b) compares the net lift coefficient with two wall distances,

=1 and 2.55, respectively. The transverse migration

direction is flipped at 𝐸𝑜 = 2.3 for the near wall case. Without solid wall boundary, the 3.4, will migrate toward the channel wall which hydrodynamically

AN US

spherical bubble, 𝐸𝑜

corresponds the top wall in our uniform shear flow simulations. However, the existence of the wall repulsive force will push the bubble away from the wall due to the unbalanced pressure field shown in Fig. 17. Thus, a spherical bubble, 𝐸𝑜 = 2.3, will experience a

M

change in the lateral force direction near the wall, while this will not happen in the case

ED

without the wall presence. Fig. 17 shows that the pressure field in the front of the deformable bubble is tilted and stronger compared to the spherical bubble‟s. Fig. 18

PT

compares the streamline patterns for different bubble deformation levels. Denser streamlines are formed on the top of spherical bubble, while strong vortex is formed behind

AC

CE

the deformable bubble.

CR IP T

ACCEPTED MANUSCRIPT

M

AN US

Fig. 15. Evolution of control forces versus bubble deformation levels.

AC

CE

PT

ED

Fig. 16. Interfacial closures versus 𝐸𝑜. (a) shows the drag coefficient versus 𝐸𝑜. (b) compares the lift coefficients with the case where the bubble is far away from the wall, = 2.55.

CR IP T

ACCEPTED MANUSCRIPT

AN US

Fig. 17. Comparison of pressure field for different bubble deformation levels. Three pressure contours of 5, 0, and -5 Pa are shown.

Fig. 18. Comparison of streamline pattern for different bubble deformation levels.

M

The bubble transverse migration behavior is a coupled phenomenon of bubble deformation, wall distance and relative velocity. Fig. 19 shows a comprehensive coefficient 𝑒 for wall distance

=

ED

map as a function of 𝑅𝑒𝑏 and

. Additional dataset

(surrounded by dashed circle) for high bubble Reynolds number (𝑅𝑒𝑏

PT

bubble deformation level ( 𝑒

5 ) and high

2.5) have been generated and are presented here to gain a

CE

complete picture. The analytical interpolation generates the following expressions for the

AC

prediction of life coefficient sign

= {

𝑅𝑒𝑏

(3

.6𝑒

)

𝑅𝑒𝑏

(3

.6𝑒

)=

𝑅𝑒𝑏

(3

.6𝑒

)

( 29 )

Sugioka and Tsukada‟s DNS data (Sugioka and Tsukada, 2015) and Takemura et al.‟ experimental data (Takemura and Magnaudet, 2003) for spherical bubble study are also

ACCEPTED MANUSCRIPT

included. Positive and negative lift coefficients are marked with blue and red colors, separately. This map can help us identify the region where the bubble migration direction

M

AN US

CR IP T

switches, thus predicting the bubble behavior.

PT

ED

Fig. 19. Lift coefficient map (𝑅𝑒𝑏 𝑒) at ⁄ = . Red and blue symbols represent negative and positive lift coefficient, respectively. The data sources are from Sugioka and Sugioka and Tsukada‟s DNS data (sphere) (Sugioka and Tsukada, 2015), Takemura et al.‟s experimental data (square) (Takemura and Magnaudet, 2003), and PHASTA‟s DNS data (triangle).

CE

5. Conclusions

In the present research, the interfacial forces are evaluated in laminar shear flow using

AC

interface tracking method. Note that highly empirical lift and drag coefficients are not used to obtain the high-resolution simulation results in the presented research. A set of parametric studies, including the wall proximity effect and local bubble Reynolds number, and bubble deformation are performed under well-controlled conditions. The results show that the wall presence starting to impact the bubble motion when the wall distance (

) is

ACCEPTED MANUSCRIPT

less than 2. As the bubble approaches the wall, the drag force increases which implies the increasing bubble rise velocity compared to liquid. The bubble transverse migration direction changes at

around 1.75, while the transition point is determined by both

CR IP T

bubble Reynolds number, 𝑅𝑒𝑏 = 42.3 , and bubble deformation level, 𝐸𝑜 = 2.58 . The effect of bubble Reynolds number is investigated as well and the lift sign change is observed at 𝑅𝑒𝑏 = 35 which is close to previous experimental observations. As the bubble becomes deformable, the net transverse lift force experiences a sign change at 𝐸𝑜 = 2.3,

migration map (

(𝑅𝑒𝑏

AN US

which is smaller than our previous observation without the wall presence. A novel bubble 𝑒)) is produced based on the results obtained. The results and

methods presented herein help us understand the effect of the wall presence on the bubble migration behavior in two-phase flow and provide the insights on the multiphase flow

M

modelling. Future work will include the coupling of wall presence and turbulent flow on

PT

ED

the interfacial forces closure studies.

6. Acknowledgements

CE

This work was supported by the National Science Foundation (award 133399 under CBET-Fluid Dynamics Program). Computational resources were supported by U.S.

AC

Nuclear Regulatory Commission‟s Faculty Development Program. An award of computer time was provided by the ASCR Leadership Computing Challenge (ALCC) program. This research used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. The

ACCEPTED MANUSCRIPT

solution presented herein made use of the Acusim linear algebra solution library provided by Altair Engineering Inc. and meshing and geometric modeling libraries by Simmetrix Inc.

CR IP T

References Ambari, A., Manuel, B.G., Guyon, E., 1983. Effect of a plane wall on a sphere moving parallel to it. Journal de Physique Lettres 44, 143-146.

Antal, S., Lahey, R., Flaherty, J., 1991. Analysis of phase distribution in fully developed

AN US

laminar bubbly two-phase flow. Int. J. Multiphase Flow 17, 635-652.

Behafarid, F., Shaver, D., Bolotnov, I., Antal, S., Jansen, K., Podowski, M., 2013. Coupled DNS/RANS Simulation of Fission Gas Discharge During Loss-of-Flow Accident in Generation IV Sodium Fast Reactor. Nucl Technol 181, 44-55.

M

Bernoulli, D., 1738. Hydrodynamica.

ED

Bolotnov, I.A., Jansen, K.E., Drew, D.A., Oberai, A.A., Lahey, J.T.,R., 2011. Detached Direct Numerical Simulation of Turbulent Two-phase Bubbly Channel Flow. Int. J.

PT

Multiphase Flow 37, 647-659.

Brackbill, J.U., Kothe, D.B., Zemach, C., 1992. A continuum method for modeling surface

CE

tension. Journal of computational physics 100, 335-354. Clift, R., Grace, J.R., Weber, M.E., 2005. Bubbles, Drops, and Particles. Courier

AC

Corporation.

Cox, R., Hsu, S., 1977. The lateral migration of solid particles in a laminar flow near a plane. Int. J. Multiphase Flow 3, 201-222.

ACCEPTED MANUSCRIPT

Dijkhuizen, W., van Sint Annaland, M., Kuipers, J., 2010. Numerical and experimental investigation of the lift force on single bubbles. Chemical engineering science 65, 12741287.

CR IP T

Drew A., D., Passman L., S., 1998. Theory of Multicomponent Fluids. Springer, New York. Drew, D., Lahey, R., 1987a. The virtual mass and lift force on a sphere in rotating and straining inviscid flow. Int. J. Multiphase Flow 13, 113-121.

Drew, D., Lahey, R., 1987b. The virtual mass and lift force on a sphere in rotating and

AN US

straining inviscid flow. Int. J. Multiphase Flow 13, 113-121.

Ervin, E., Tryggvason, G., 1997. The rise of bubbles in a vertical shear flow. Journal of Fluids Engineering 119, 443-449.

Fang, J., Feng, J., Bolotnov, I.A., 2016. Integral and separate effect simulations of bubbly

M

flows using interface tracking approach. 2016 CFD4NRS-6 Workshop.

Simulations.

ED

Fang, J., 2016. Development of Advanced Analysis Toolkit for Turbulent Bubbly Flow

PT

Fang, J., Rasquin, M., Bolotnov, I.A., 2017. Interface tracking simulations of bubbly flows in PWR relevant geometries. Nucl. Eng. Des. 312, 205-213.

CE

Feng, J., Bolotnov, I.A., 2017a. Evaluation of Bubble-Induced Turbulence Using Direct Numerical Simulation. Int. J. Multiphase Flow 93, 92-107.

AC

Feng, J., Bolotnov, I.A., 2017b. Interfacial force study on a single bubble in laminar and turbulent flows. Nucl. Eng. Des. 313, 345-360.

Feng, J., Bolotnov, I.A., 2017c. Development of interfacial forces closures based on DNS data. 17th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH17) .

ACCEPTED MANUSCRIPT

Fukuta, M., Takagi, S., Matsumoto, Y., 2008. Numerical study on the shear-induced lift force acting on a spherical bubble in aqueous surfactant solutions. Physics of Fluids (1994-present) 20, 040704.

CR IP T

Gorelik, R., Kashinskii, O., Nakoryakov, V., 1987. Study of a downward bubbly flow in a vertical pipe. Journal of Applied Mechanics and Technical Physics 28, 64-67.

Hadamard, J.S., 1911. Mouvement permanent lent d'une sphere liquide et visqueuse dans un liquide visqueux. C. R. Acad. Sci. 152, 1735.

AN US

Hibiki, T., Ishii, M., 1999. Experimental study on interfacial area transport in bubbly twophase flows. Int. J. Heat Mass Transfer 42, 3019-3035.

Hosokawa, S., Tomiyama, A., Misaki, S., Hamada, T., 2002. Lateral migration of single bubbles due to the presence of wall. , 855-860ASME 2002 Joint US-European Fluids

M

Engineering Division Conference.

ED

Ishii, M., Mishima, K., 1984. Two-fluid model and hydrodynamic constitutive relations. Nucl. Eng. Des. 82, 107-126.

PT

Ishii, M., Chawla, T., 1979. Local drag laws in dispersed two-phase flow. Jansen, K.E., 1993. Unstructured grid large eddy simulations of wall bounded flows.

CE

Annual Research Briefs, Center for Turbulence Research, NASA Ames/Stanford University , 151.

AC

Jansen, K.E., 1999. A stabilized finite element method for computing turbulence. Comput. Methods Appl. Mech. Eng. 174, 299-317.

Lahey, R.T., Drew, D.A., 2001. The analysis of two-phase flow and heat transfer using a multidimensional, four field, two-fluid model. Nucl. Eng. Des. 204, 29-44.

ACCEPTED MANUSCRIPT

Lahey, R., de Bertodano, M.L., Jones, O., 1993. Phase distribution in complex geometry conduits. Nucl. Eng. Des. 141, 177-201. Legendre, D., Merle, A., Magnaudet, J., 2006. Wake of a spherical bubble or a solid sphere

CR IP T

set fixed in a turbulent environment. Phys. Fluids 18, 048102. Li, M., Feng, J., Bolotnov A, I., 2017. Heat Transfer Efficiency in Various Prandtl Number Turbulent Flows using DNS approach. 2017 ANS Annual Meeting.

Liu J., T., 1993. Bubble size and entrance length effects on void development in a vertical

AN US

channel. Int.Journal of Multiphase Flow 19, 99-113.

Lopez de Bertodano, Martin A, 1998. Two fluid model for two-phase turbulent jets. Nucl. Eng. Des. 179, 65-74.

Lu, J., Tryggvason, G., 2008. Effect of bubble deformability in turbulent bubbly upflow in

M

a vertical channel. Physics of Fluids (1994-present) 20, 040701.

ED

Lucas, D., Krepper, E., Prasser, H., 2007. Use of models for lift, wall and turbulent dispersion forces acting on bubbles for poly-disperse flows. Chemical Engineering

PT

Science 62, 4146-4157.

McLaughlin, J.B., 1991. Inertial migration of a small sphere in linear shear flows. J. Fluid

CE

Mech. 224, 261-274.

Moore, D., 1963. The boundary layer on a spherical gas bubble. J. Fluid Mech. 16, 161-176.

AC

Moraga, F., Larreteguy, A., Drew, D., Lahey, R., 2003. Assessment of turbulent dispersion models for bubbly flows in the low Stokes number limit. Int. J. Multiphase Flow 29, 655-673.

ACCEPTED MANUSCRIPT

Nagrath, S., Jansen, K., Lahey Jr, R.T., Akhatov, I., 2006. Hydrodynamic simulation of air bubble implosion using a level set approach. Journal of Computational Physics 215, 98132.

CR IP T

Nakoryakov, V., Kashinsky, O., Kozmenko, B., Gorelik, R., 1986. Study of upward bubbly flow at low liquid velocities. Izv.Sib.Otdel.Akad.Nauk SSSR 16, 15-20.

Podowski, M., 2009. On the consistency of mechanistic multidimensional modeling of gas/liquid two-phase flows. Nucl. Eng. Des. 239, 933-940.

AN US

Podowski, M.Z., 2008. Multidimensional modeling of two-phase flow and heat transfer. Int. J. Numer. Methods Heat Fluid Flow 18, 491-513.

Saffman G., P., 1965. The lift on a small sphere in a slow shear flow. J. Fluid Mech. 22, 385-400.

M

Sahni, O., Müller, J., Jansen, K.E., Shephard, M.S., Taylor, C.A., 2006. Efficient

ED

anisotropic adaptive discretization of the cardiovascular system. Comput. Methods Appl. Mech. Eng. 195, 5634-5655.

Press.

PT

Sethian, A.J., 1999. Level Set Methods and Fast Marching Methods. Cambridge University

CE

Sugioka, K., Tsukada, T., 2015. Direct numerical simulations of drag and lift forces acting on a spherical bubble near a plane wall. Int. J. Multiphase Flow 71, 32-37.

AC

Sussman, M., Almgren, A.S., Bell, J.B., Colella, P., Howell, L.H., Welcome, M.L., 1999. An adaptive level set approach for incompressible two-phase flows. Journal of Computational Physics 148, 81-124. Sussman, M., Fatemi, E., Smereka, P., Osher, S., 1998. An improved level set method for incompressible two-phase flows. Comput. Fluids 27, 663-680.

ACCEPTED MANUSCRIPT

Sussman, M., Fatemi, E., 1999. An efficient, interface-preserving level set re-distancing algorithm and its application to interfacial incompressible fluid flow. Siam J.on Scientific Computing 20, 1165-1191.

CR IP T

Takemura, F., Magnaudet, J., 2003. The transverse force on clean and contaminated bubbles rising near a vertical wall at moderate Reynolds number. J. Fluid Mech. 495, 235-253.

Takemura, F., Takagi, S., Magnaudet, J., Matsumoto, Y., 2002. Drag and lift forces on a

AN US

bubble rising near a vertical wall in a viscous liquid. J. Fluid Mech. 461, 277-300.

Tejada-Martínez, A.E., Jansen, K.E., 2006. A parameter-free dynamic subgrid-scale model for large-eddy simulation. Comput. Methods Appl. Mech. Eng. 195, 2919-2938. Thomas, A.M., Fang, J., Bolotnov, I.A., 2014. Estimation of Shear-induced Lift Force in

ED

Hydraulics - 2014 (ATH '14).

M

Laminar and Turbulent Flows. International Topical Meeting on Advances in Thermal

Thomas, A.M., Fang, J., Feng, J., Bolotnov, I.A., 2015. Estimation of shear-induced lift

PT

force in laminar and turbulent flows. Nucl Technol 190, 274-291. Tomiyama, A., Sou, A., Zun, I., Kanami, N., Sakaguchi, T., 1995. Effects of Eötvös

CE

number and dimensionless liquid volumetric flux on lateral motion of a bubble in a laminar duct flow. Advances in multiphase flow 1995, 3-15.

AC

Tomiyama, A., 1998. Struggle with computational bubble dynamics. Multiphase Science and Technology 10, 369-405.

Tomiyama, A., Kataoka, I., Zun, I., Sakaguchi, T., 1998. Drag coefficients of single bubbles under normal and micro gravity conditions. JSME international journal.Series B, fluids and thermal engineering 41, 472-479.

ACCEPTED MANUSCRIPT

Tomiyama, A., Tamai, H., Zun, I., Hosokawa, S., 2002. Transverse migration of single bubbles in simple shear flows. Chemical Engineering Science 57, 1849-1858. Vasseur, P., Cox, R., 1977. The lateral migration of spherical particles sedimenting in a

CR IP T

stagnant bounded fluid. J. Fluid Mech. 80, 561-591. Vhora, S., Cancelos, S., Bolotnov A, I., 2013. Interface tracking study of bubble/wall interaction. 109, 1618-1620ANS Winter Meeting.

Whiting, C.H., Jansen, K.E., 2001. A stabilized finite element method for the

AN US

incompressible Navier-Stokes equations using a hierarchical basis. Int. J. Numer. Methods Fluids 35, 93-116.

Zaruba, A., Lucas, D., Prasser, H., Höhne, T., 2007. Bubble-wall interactions in a vertical

science 62, 1591-1605.

M

gas–liquid flow: Bouncing, sliding and bubble deformations. Chemical engineering

ED

Zeng, L., Balachandar, S., Fischer, P., 2005. Wall-induced forces on a rigid sphere at finite Reynolds number. J. Fluid Mech. 536, 1-25.

PT

Zeng, L., Najjar, F., Balachandar, S., Fischer, P., 2009. Forces on a finite-sized particle located close to a wall in a linear shear flow. Phys. Fluids 21, 033302.

CE

Zun, I., 1980. Transverse migration of bubbles influenced by walls in vertical bubbly flow.

AC

Int. J. Multiphase Flow 6, 583-588.