CFD Simulation of Two-phase and Three-phase Flows in Internal-loop Airlift

(1)

LAPPEENRANTA UNIVERSITY OF TECHNOLOGY Faculty of Technology

Department of Mathematics and Physics

CFD Simulation of Two-phase and

Three-phase Flows in Internal-loop Airlift Reactors

The topic of this Master’s thesis was approved by the departmental council of the De- partment of Mathematics and Physics on27^th May, 2010.

Supervisors: Professor Heikki Haario D.Sc. Arto Laari

Examiners: Professor Heikki Haario D.Sc. Arto Laari

In Lappeenranta August 24, 2010.

Giteshkumar N Patel

Teknologiapuistonkatu 2 B 28 53850 Lappeenranta

Phone: +358466477168

Email: Giteshkumar.Patel@lut.fi & gitesh.maths@gmail.com

(2)

Abstract

Lappeenranta University of Technology Department of Mathematics and Physics

Giteshkumar N Patel

CFD Simulation of Two-phase and Three-phase Flows in Internal-loop Airlift Reactors

Master’s thesis 2010

73 pages, 48 figures, 7 tables

Key words: Airlift reactors, Gas holdup, Scale effect, CFD, Two-phase flow, Three-phase flow, Solids distribution

Airlift reactors are pneumatically agitated reactors that have been widely used in chemical, petrochemical, and bioprocess industries, such as fermentation and wastewater treatment. Computational Fluid Dynamics (CFD) has become more popular approach for design, scale-up and performance evaluation of such reactors. In the present work numerical simulations for internal-loop airlift reactors were performed using the transient Eulerian model with CFD package, ANSYS Fluent 12.1. The turbulence in the liquid phase is described using the κ − model. Global hydrodynamic parameters like gas holdup, gas velocity and liquid velocity have been investigated for a range of superficial gas velocities, both with 2D and 3D simulations. Moreover, the study of geometry and scale influence on the reactor have been considered. The results suggest that both, geometry and scale have significant effects on the hydrodynamic parameters, which may have substantial effects on the reactor performance. Grid refinement and time-step size effect have been discussed.

Numerical calculations with gas-liquid-solid three-phase flow system have been carried out to investigate the effect of solid loading, solid particle size and solid density on the hydrodynamic characteristics of internal loop airlift reactor with different superficial gas velocities. It was observed that averaged gas holdup is significantly decreased with increasing slurry concentration. Simulations show that the riser gas holdup decreases with increase in solid particle diameter. In addition, it was found that the averaged solid holdup increases in the riser section with the increase of solid density. These produced results reveal that CFD have excellent potential to simulate two-phase and three-phase flow system.

(3)

Acknowledgements

First of all, I wish to convey my gratitude to the Department of Mathematics for the scholarship during my studies at Lappeenranta University of Technology.

I would like to express my sincerest gratitude to Professor Heikki Haario, for his supervision, advice, and guidance from the very early stage of this thesis as well as giving me all facilities through out the work.

I gratefully acknowledge to D.Sc. Arto laari for his valuable advice and comments, supervision and crucial contribution, which made him a backbone of this thesis. It has been a great pleasure working with and learning from him.

I am extending my heartfelt and very special gratitude to Yogi for her contributions, patience, understanding and love. She provided constant inspiration, constructive ideas and strength to hurdle all the obstacles in the completion this work. Thanks for pushing me to be the best, and giving me a shoulder when times were hard. May GOD bless her in all her endeavors because without her extraordinary support, completion of this study would not have been possible.

I am also grateful to all my friends and colleagues at the Department of Mathematics for my memorial stay in Lappeenranta.

Last but not the least, my sincere acknowledgments are addressed to my family and the one above all of us, the omnipresent God, for answering my prayers and for giving me the strength, thank you so much Dear Lord. I am grateful to my dad, Narayanbhai, for giving me the life I ever dreamed and his constant inspiration and guidance kept me focused and motivated. I can’t express my gratitude for my mom, Sitaben, in words, whose unconditional love has been my greatest strength. The constant love and support of my elder brother, Rakeshbhai, his wife, Jayshree, and their child Krrish are sincerely acknowledged. I convey special acknowledgement to my elder sister, Surekhaben, and her family. I would like to dedicate this thesis to my family and Yogi, the most important person in my life.

Lappeenranta, August 24, 2010

Giteshkumar N Patel

(4)

List of Tables

1 Details of computational grids for the Configurations I, II and III . . . 40

2 Two-phase simulation settings . . . 42

3 Phase properties used in two-phase simulation . . . 43

4 Details of grids used in mesh-independence studies . . . 53

5 Computational time used by each grid for 2D and 3D domains. . . 56

6 Simulation setup for three phase flow calculation. . . 61

7 Phases properties used for three phase flow simulation. . . 61

(7)

List of Figures

1 Multiphase flow regimes . . . 4

2 Flow patterns of multiphase flow . . . 6

3 Flow patterns of gas-liquid flow in vertical tubes . . . 7

4 Flow patterns in horizontal tubes . . . 9

5 Flow regime map for vertical gas-liquid flow . . . 10

6 Flow regime map for the horizontal flow of an air/water mixture in a 5.1cm diameter pipe (Weisman (1983)) . . . 11

7 Overview of CFD . . . 15

8 Discretisation of flow in CFD . . . 16

9 Interface calculation . . . 23

10 Different types of ALRs [10]. . . 37

11 Schematic overview of computational geometry for Configurations I, II and III. . . 39

12 (a) Boundary conditions applied on computational domain, (b) computational grid, and (c) closer view of the bottom region of reactor. . . 41

13 Transient approach to steady-state of the gas and liquid velocities at the centre of the reactor. . . 43

14 Contours of (a) gas velocity, and (b) liquid velocity of Configuration-I whenUg =0.02 m/s. . . 44

15 Gas fraction contours from 2D and 3D cases whenUg= 0.02m/s. . . 45

16 (a) Liquid velocity vector of 2D case and (b) liquid iso-surfaces from 3D simulation. . . 45

17 Gas velocity distribution in the riser tube at selected axial location. . . 46

18 Liquid velocity distribution at selected axial location. . . 46

19 Gas holdup distribution in the riser tube at selected axial location. . . 47

(8)

20 Radial distribution of gas velocity from 2D and 3D simulations at selected axial location in the riser tube. . . 48 21 Radial distribution of liquid velocity from 2D and 3D simulations at se-

lected axial location. . . 48 22 Radial distribution of gas holdup from 2D and 3D simulations at selected

axial location in the riser tube. . . 48 23 Contours of (a) gas velocity, (b) liquid velocity, and (c) gas holdup for

Configuration-II. . . 49 24 Contours of (a) gas velocity, (b) liquid velocity, and (c) gas holdup for

Configuration-III. . . 50 25 (a) Gas velocity and (b) liquid velocity distributions of Configuration-I. . 50 26 (a) Gas velocity and (b) liquid velocity distributions of Configuration-II. . 51 27 (a) Gas velocity and (b) liquid velocity distributions of Configuration-III. 51 28 Gas holdup distributions in (a) Configuration-I and (b) Configuration-II. . 52 29 Comparison between average gas holdup in the riser of Configuration-I

and Configuration-II with range ofUg. . . 52 30 Instantaneous snapshots of gas holdup contours with (a) Grid A, (b) Grid

B, (c) Grid C, and (d) Grid D whenU_g =0.02 m/s. . . 54 31 Gas velocity distribution with 2D and 3D simulations for different grids. . 55 32 Gas holdup distribution with 2D and 3D simulations for different grids. . 55 33 Liquid velocity distribution with 2D and 3D simulations for different grids. 55 34 Profiles of gas hold produced from each selected grid and time-step size. . 57 35 Profiles of liquid velocity produced from each selected grid and time-step

size. . . 58 36 Contours of velocity magnitude for (a) liquid, (b) gas, and (c) solid of

three-phase flow system. . . 62 37 Instantaneous snap shot of solid concentration in the reactor when reactor

is initially loaded with5%of solids. . . 63 38 Average gas holdup in the riser with different superficial gas velocity. . . . 64

(9)

39 Average solids holdup in the riser with different superficial gas velocity. . . 64 40 Radial distribution of solids volume fraction at the center of the reactor

with different solids loading conditions. . . 65 41 Averaged gas holdup in the riser with different superficial gas velocity. . . 65 42 Average solids holdup in the riser with different superficial gas velocity. . . 66 43 Radial distribution of solids volume fraction at the center of the reactor

with different solid particle size. . . 67 44 Instantaneous snap shot of solids concentration in the reactor withρs =

2000kg/m³. . . 68 45 Instantaneous snap shot of solids concentration in the reactor withρs =

3000kg/m³. . . 68 46 Average solids holdup in the riser with different superficial gas velocity. . . 69 47 Average gas holdup in the riser section with different superficial gas velocity. 69

(10)

Nomenclature

A cross sectional area [m²]

b_d bubble diameter [m]

CD drag coefficient dimensionless

C_l lift coefficient dimensionless

e coefficient of restitution dimensionless

f force [N]

g gravitational acceleration [m/s²]

h enthalpy J/kg

K interphase exchange momentum [N/m³]

m mass [kg]

N number of nodes

p pressure [n/m²]

Q volumetric flow rate of the phase [m³/s]

R_e Reynolds number dimensionless

t time [s]

u velocity [m/s]

U_g superficial gas velocity [m/s]

V volume [m³]

(11)

Greek Letters

α phase volume fraction dimensionless

dissipation rate of turbulent kinetic energy [m²/s³]

ε holdup dimensionless

κ turbulent kinetic energy [m²/s²]

µ viscosity [kg/ms]

µt turbulent viscosity [kg/ms]

ρ density [kg/m³]

σ surface tension [N/m]

τk viscous stress tensor [kg/ms²]

∇ gradient operator

Abbreviations

ALR Airlift Reactos

CFD Computational Fluid Dynamics FDM Finite Difference Method FEM Finite Element Method FVM Finite Volume Method PDE Partial Differential Equation VOF Volume of Fluid Method

(12)

Subscript or Superscript

c referring continuous phase d referring dispersed phase D drag

G referring to gas m referring mixture M T mass transfer L referring to liquid LF referring to lift force

LU B referring to lubrication force S referring to solid

l, s referring to fluid phase and solid phase respectively p, q referring to phase p and q respectively

T D turbulent dispersion V M virtual mass

(13)

1 Introduction

The topic of the thesis is CFD simulation of two-phase and three-phase flows in internal- loop airlift reactors. The term airlift reactor (ALR) covers a wide range of gas-liquid or gas-liquid-solid pneumatic contacting devices that are characterized by fluid circulation in a defined cyclic pattern through channels built specifically for this purpose [10]. Airlift reactors are characterized by three distinct regions namely riser, gas-liquid-separator and downcomer [11]. There are two types of ALR, ‘internal’ and ‘external’ loop ones. Internal loop reactors consist of concentric tubes or split vessels, in which a part of the gas is entrained into the downcomer, whereas external loop reactors are two conduits connected at the top and the bottom, in which little or no gas recirculates into the downcomer. The part in which the sparger is located is called the riser, and the other one, the downcomer.

Liquid circulation is induced by injecting gas at the bottom of the riser, thus creating a density differences between riser and downcomer [12]. In such reactors, the required interaction is provided by the density differences between the gas and the liquid. This fact enables extremely interaction and complicated gas-liquid or gas-liquid-solid reactions to take place in such reactors. The airlift is popular as a gas-liquid contactor because it can handle large quantities of liquid and gas on a continuous basis. In addition, it has no moving parts, requires limited amount of energy for its operation and exhibits good mass and heat transfer characteristics [13]. Due to this, ALR’s are finding increasing applications in chemical industry, biochemical fermentation and biological wastewater treatment processes [16][17].

Several publications have established the potential of computational fluid dynamics for describing the hydrodynamics of bubble columns [14][15]. The design, scale-up, and performance evaluation of such reactors all require extensive and accurate information about the gas-liquid flow dynamics, particularly as CFD has become more popular approach [29] . An important advantage of the CFD approach is that column geometry and scale effects are automatically accounted for. CFD can be used to simulate and optimize mixing, gas hold-up and mass transfer coefficients and distribution of the phases [18].

In addition, scale and operating conditions of reactors have significant effects on global hydrodynamic parameters such as: gas holdup, liquid circulation velocity, overall mass transfer rate in gas-liquid flow. In gas-liquid-solid flow, solid loading, solid particle size and solid density significantly affect the hydrodynamics characteristics of internal loop air lift reactor (Snape et al 1995, Vial et al 2000).

(14)

1.1 Objectives of the thesis

The objectives of the thesis have been set as follows:

• To simulate the gas-liquid two-phase flow with internal-loop airlift reactor application by using Eulerian model in CFD. To investigate the 2D and 3D simulations with various operating conditions (e.g., superficial gas velocity) and to check the effect on global hydrodynamic parameters like gas holdup, gas velocity and liquid velocity. To study the influence of geometry and scale on the reactor hydrodynamics. Furthermore, grid and time-step sensitivity are performed to investigate the influence of grid refinement and time-step size effect respectively on the results.

• By performing CFD simulation, to simulate gas-liquid-solid three-phase flow in internal-loop airlift reactor system and to investigate the effects of superficial gas velocity and particle loading on the gas holdup in the riser section. To study the distribution of solid phase in the reactor with different operating conditions. To determine the effect of solid particle diameter and solid density on averaged gas holdup with the range of superficial gas velocity.

1.2 Thesis structure

This chapter consists of the objectives and the methodology of the thesis work. The reader is then introduced to the main content of the following chapters.

Chapter 2 gives the brief details of multiphase flow. It describes the division of multiphase flow regimes. The details information of two-phase flow are provided in it. Some fundamental terms are described there. At the end, CFD is introduced concisely with discretisation in solving fluid flow.

Chapter 3 demonstrates the numerical approaches for multiphase flow in CFD. The different models for solving multiphase flow are explained. The last part contains the details of turbulence models for multiphase flow.

Chapter 4 displays the details about the interphase momentum transform. Some information for interphase exchange coefficients are derived in it.

Chapter 5 presents work done on the two-phase flow with internal-loop airlift reactors by using Eulerian model with different operating conditions. It shows the comparison between 2D vs. 3D results of global hydrodynamics parameters. Scale effect are explained in it. Grid and time-step size effect cases are discussed at the end.

Chapter 6 contains the three-phase flow work. In this chapter, the effects of the superficial

(15)

gas velocity and the particle loading on the gas holdup in the riser are discussed. Solids distribution in the reactor are studied there. Some results about the effect of solid particle diameter on gas holdup with the range of superficial gas velocity are explained.

Chapter 7 describes conclusions of the whole thesis. This focuses on the objectives of the work and how they are achieved throughout the thesis.

(16)

2 Topologies of multiphase flow

A phase is simply one of the states of matter and it can be either a gas, a liquid, or a solid [19]. Multiphase flow is the simultaneous flow of several phases. Two-phase flow is the simplest case of multiphase flow [20]. From a practical engineering point of view one of the major design difficulties in dealing with multiphase flow is that the mass, momentum, and energy transfer rates and processes can be quite sensitive to the geometric distribution or topology of the components within the flow [8]. An appropriate starting point is a phenomenological description of the geometric distributions or flow patterns that are observed in common multiphase flows. First, we will discuss about the flow regimes of multiphase flow.

2.1 Multiphase flow regimes

The definition of the flow regime is a description of the morphological arrangement of the components, or flow pattern [19]. It is important to appreciate that different flow regimes occur at different fluid flow rates and differences also occur for different materials.

Multiphase flow regimes can be grouped into four categories: gas-liquid or liquid-liquid flows; gas-solid flows; liquid-solid flows; and three-phase flows [21]. The following Figure 1 represent a schematic diagram of the multiphase flow regimes.

Figure 1: Multiphase flow regimes

(17)

Gas-liquid or liquid-liquid flows

The combined flow of gas and liquid is called a gas-liquid flow. Bubbly flow, slug slow, droplet flow and stratified/free-surface flow, etc. are the flow regimes types of gas-liquid flow. Later on we will get brief details about it.

Gas-solid flows

The simultaneous flow of gas and solid is a gas-solid flow. The following regimes can be separated as:

• Fluidized bed: Fluidized bed is a vertical cylinder mechanism. It contains the solid particles (homogeneous material or heterogeneous material). In fluidized bed a gas is introduced through a distributor. The gas rising through the bed suspends the particles. Depending on the gas flow rate, bubbles appear and rise through the bed. It can be seen at CFB (Circulating Fluidized Bed) or BFB (Bubbling Fluidized Bed) boiler combustion technology.

• Particle-laden flow: This type of flow regime contains discrete particles in a continuous gas. Particle-laden flow examples include cyclone separators, air classifiers, dust collectors, and dust-laden environmental flows [21].

• Pneumatic transport: This is a flow regime that depends on factors such as solid loading, Reynolds numbers, and particle properties. Typical patterns are dune flow, slug flow, and homogeneous flow. Pneumatic transport examples include transport of cement, grains, and metal powders [21].

Liquid-solid flows

The mixture of solid particles in liquid form a liquid-solid flow. Examples of it are slurry flow or sedimentation.

• Slurry flow: This flow is the transport of particles in liquids. The fundamental behavior of liquid-solid flows varies with the properties of the solid particles relative to those of the liquid. Slurry flow examples include slurry transport and mineral processing.

• Hydrotransport: This describes densely-distributed solid particles in a continuous liquid. Hydrotransport examples include mineral processing and biomedical and physio-chemical fluid systems.

(18)

Figure 2: Flow patterns of multiphase flow

• Sedimentation: This describes a tall column initially containing a uniform dispersed mixture of particles. At the top of the column, a clear interface will be appeared.

At the bottom, the particles will slow down and form a sludge layer and in the middle a constant settling zone will exist. Sedimentation examples include mineral processing [21].

Figure 2 displays examples of multiphase flow with various flow patterns.

Three-phase flows

Three-phase flows are combinations of the other flow regimes. It means a combination of gas-liquid-solid or two solid phases and one gas phase, etc. This types of flow can be seen at petroleum refinery, in chemical separation technology or in combustion. In the present work, we will concentrate on both two-phase flow and three-phase flow. The following section describes the fundamental things about gas-liquid flow.

(19)

2.2 Gas-Liquid two phase flow regimes

The flow of gas-liquid mixtures in pipes and other items of process equipment is common and extremely important. It is important to appreciate that different flow regimes occurs at different gas and liquid flow rates and differences also occur for different materials [22].

In this section we will discuss about flow patterns at different position of tubes.

2.2.1 Flow patterns in vertical tubes

When the pipe is oriented vertically, the regimes of gas/liquid flow are illustrated in Figure 3 (see, for example, Hewitt and Hall Taylor 1970, Butterworth and Hewitt 1977, Hewitt 1982, Whalley 1987). The flow patterns are transformed from bubbly flow to disperse flow gradually where ’G’ refers for gas phase and ’L’ for liquid. The sequence shown is that which would normally be seen as the ratio of gas to liquid flow rates is increased.

In thebubbly regime, there is a distribution of dispersed bubbles throughout the continuous liquid phase. The bubbles may vary widely in size and shape but they are typically nearly spherical and are much smaller than the diameter of the tube itself. As the gas flow rate increases, the average bubble size increases [22].

Figure 3: Flow patterns of gas-liquid flow in vertical tubes

The next regime occurs when the gas flow rate is increased to the point when many bubbles collide and coalesce to produce slugs of gas or to form larger bubbles. Generally it is known asslug flow. The gas slugs have spherical noses and occupy almost the entire

(20)

cross section of the tube, being separated from the wall by a thin liquid film. They are commonly referred to as Taylor bubbles [19]. Between slugs of gas there are slugs of liquid in which there may be small bubbles entrained in the wakes of the gas slugs.

Increasing the velocity of the flow, the structure of the flow becomes unstable with fluid traveling up and down in an oscillatory fashion but with a net upward flow. At higher flow rates, a slug flow pattern is destroyed and a chaotic type of flow, generally known aschurn flow. Over most of the cross section there is a churning motion of irregularly shaped portions of gas and liquid.

Further increase in the gas flow rate causes a degree of separation of the phases, the liquid flowing mainly on the wall of the tube and the gas in the core area. This flow regimes is known asannular flow [23]. The increment in the liquid flow rate, the concentration of liquid is increased in the gas core. It is known as wispy-annular flow. The main differences between the wispy-annular and the annular flow regimes are that in the former the entrained liquid is present as relatively large drops and the liquid film contains gas bubbles, while in the annular flow regime the entrained droplets do not coalesce to form larger drops.

2.2.2 Flow patterns in horizontal tubes

Consider a cocurrent gas-liquid flow in horizontal pipes, which displays similar flow patterns to those for vertical flow. However, asymmetry is caused by the effect of gravity, which is most significant at low flow rates. The sequence of flow regimes in horizontal tubes as identified by Alves (1954) is shown in Figure 4. In the figure, the ’G’ is gas phase and ’L’ for liquid phase.

In the bubbly regime the bubbles are confined to a region near the top of the pipe. On increasing the gas flow rate, the bubbles become larger and coalesce to form long bubbles which is known as the plug flow regime. If the gas flow rate as increased then the gas plugs join and form a continuous gas layer in the upper part of the pipe. This type of flow, in which the interface between the gas and the liquid is smooth, is known as the stratified flow regime. Where the gas has lower viscosity and lower density, it will flow faster than the liquid.

As the gas flow rate is increased further, the interfacial shear stress becomes sufficient to generate waves on the surface of the liquid producing the wavy flow regime. If the gas flow rate continues to rise, the waves which travel in the direction of flow, grow until their crests approach the top of the pipe and, as the gas breaks through, liquid is distributed over the wall of the pipe. This is known as theslug flow regime.

At higher gas flow rates anannular flow regime is found as in the vertical flow. At very

(21)

Figure 4: Flow patterns in horizontal tubes

high flow rates the liquid film may be very thin, the majority of the liquid being dispersed as droplets in the gas core. This type of flow may be called as the spray or mist flow regime. It may be noted that similar flow regimes can be seen with immiscible liquid systems. If the densities of the two liquids are close the flow regimes for the horizontal flow will more nearly resemble those of the vertical flow.

2.3 Examples of flow regime maps

The prediction of the flow regime in gas-liquid two-phase flow is rather uncertain partly because the transitions between the flow regimes are gradual and the classification of a particular flow is subjective [22]. In addition, there are many industrial processes in which the mass quality is a key flow parameter and therefore mass flux maps are often preferred.

There are various flow regime maps in the literature which gives the idea of the flow patterns and behavior, two of which are given in the following Figures 5 and 6. For vertical flow of low pressure air-water and high pressure steam-water mixtures, Hewitt and Roberts (1969) have determined a flow regime map shown in Figure 5. Here,jG and jL denote the volumetric fluxes of the gas and liquid.

(22)

Figure 5: Flow regime map for vertical gas-liquid flow For the gas phase the volumetric flux is given by,

jG= QG

A (1)

and for the liquid phase is given by,

j_L= Q_L

A (2)

whereQ_G and Q_L are the volumetric flow rates of the gas and the liquid, and A is the cross-sectional area. Here in Figure 5 the axes represent the superficial momentum fluxes of the gas and liquid. In addition to allowing the flow regime for a specified combination of gas and liquid flow rates to be determined, the diagram shows how changes of operating conditions change the flow regime.

In particular it can be seen that the sequence of flow regimes described above is produced by increasing the gas momentum flux and/or reducing the liquid momentum flux. If the momentum flux of gas and liquid are relatively low then the flow regime can be churn and for higher rate of momentum flux, it will be wispy annular.

Figure 6 displays the flow regime map for horizontal gas-liquid flow patterns. It was given by Weisman (1983). It shows the occurrence of different flow regimes for the flow of an air/water mixture in a horizontal, 5.1cm diameter pipe. HereG_L andG_G, denote

(23)

Figure 6: Flow regime map for the horizontal flow of an air/water mixture in a 5.1cm diameter pipe (Weisman (1983))

the superficial mass fluxes of the liquid and the gas. For the liquid phase,

GL= M_L

A =jLρL (3)

and for gas phase it is given by,

G_G= M_G

A =j_Gρ_G (4)

WherejLandjGdenote the volumetric fluxes of the liquid and gas respectively. Density of liquid phase is ρL and for gas isρG.

2.4 Basic concepts of multiphase flow

This section introduces a few definitions that are fundamental to multiphase flows. For convenience, the term discrete or dispersed phase will be used for the particles, droplets, or bubbles, while carrier or continuous phase will be used for the carrier fluid. The dispersed phase is not materially connected.

(24)

Dispersed Phase and Separated Flows

Dispersed phase flows are flows in which one phase consists of discrete elements, such as droplets in a gas or bubbles in a liquid. The discrete elements are not connected. Where in a separated flow, the two phases are separated by a line of contact. An annular flow is a separated flow in which there is a liquid layer on the pipe wall and a gaseous core.

In other words, in a separated flow one can pass from one point to another in the same phase while remaining in the same medium.

Volume fraction and Densities

The volume fraction of the dispersed phase is defined as

αd= lim

δV→V⁰

δVd

δV (5)

where δV_d is the volume of the dispersed phase in volume δV. The volume δV⁰ is the limiting volume that ensures a stationary average. Unlike a continuum, the volume fraction cannot be defined at a point. Equivalently, the volume fraction of the continuous phase is

αc= lim

δV→V⁰

δVc

δV (6)

where δVc is the volume of the continuous phase in volume. This volume fraction is sometimes referred to as the void fraction and in the chemical engineering literature, the volume fraction of the dispersed phase is often referred to as holdup [24]. By definition, the sum of the volume fractions must be unity, i.e.,

αd+αc= 1 (7)

The bulk density (or apparent density) of the dispersed phase is the mass of the dispersed phase per unit volume of mixture or, in terms of a limit, is defined as

ρ_d= lim

δV→V⁰

δM_d

δV (8)

whereδM_dis the mass of the dispersed phase. The bulk density is related to the material densityρ_d by

(25)

ρ_d=α_dρ_d (9) The sum of the bulk densities of the dispersed and continuous phases is the mixture density

ρ_d+ρ_c=ρm (10)

Superficial and Phase Velocities

For multiphase flow, the superficial velocity of each phase is the mass flow rateM˙of that phase divided by the cross sectional area Aand material density or in other words, the superficial velocity is nothing but the velocity of a fluid in a pipe, conduit, column etc in the absence of packing or obstruction. Like in packed columns the actual velocity of the fluid through it is actually the volumetric flow rate divided by the cross sectional area.

so the velocity achieved by the same fluid in the same column in absence of the packing is called superficial velocity. The superficial velocity for the dispersed phase is

U_d= M˙_d

ρ_dA (11)

The phase velocityuis the actual velocity of the phase. The superficial velocity and the phase velocity are related by the volume fraction

U_d=α_du_d (12)

In addition, the same relations hold for the carrier phase.

Quality, Concentration, and Loading

Another parameter important to the definition of dispersed-phase flows is the dispersed- phase mass concentration, which is the ratio of the mass of the dispersed phase to that of the continuous phase in a mixture as,

C = ρ_d

ρ_c (13)

(26)

This parameter will sometimes be referred to as the particle or droplet mass ratio. Some- times dispersed-phase volume fraction is designated as concentration. The quality of a liquid–vapor mixture where the liquid in the dispersed phase is

x= ρ_d

ρ_m (14)

Another term in common use in multiphase flows is loading, which is the ratio of mass flux of the dispersed phase to that of the continuous phase:

z= m˙d

˙

m_c (15)

However, loading has also been used to denote concentration.

2.5 Computational Fluid Dynamics

Fluid dynamics is the science of fluid motion. Fluid flow is commonly studied in one of three ways:

• Experimental fluid dynamics

• Theoretical fluid dynamics

• Numerically: computational fluid dynamics (CFD)

CFD is one of the branches of fluid mechanics [25]. It is the science of predicting fluid flow, heat transfer, mass transfer, chemical reactions, and related phenomena by solving the mathematical equations which govern these processes using numerical methods and algorithms. In order to provide easy access to their solving power all commercial CFD packages include sophisticated user interfaces to input problem parameters and to exam- ine the results. Hence all CFD codes contains three main elements: (i) a pre-processor, (ii) a solver and (iii) a post-processor [26].

• Pre-processor: pre-processing consists of the input of a flow problem to a CFD program by means of an operator-friendly interface and the subsequent transfor- mation of this input into a form suitable for use by the solver [26]. The region of fluid to be analysed is called the computational domain and it is made up of a number of discrete elements called the mesh (or grid). The users need to define the properties of fluid acting on the domain before the analysis is begun; these include external constraints or boundary conditions, like pressure and velocity to implement realistic situations.

(27)

• Solver: a program that calculates the solution of the CFD problem. Here the governing equations are solved. This is usually done iteratively to compute the flow parameters of the fluid as time elapses. Convergence is important to produce an accurate solution of the partial differential equations.

• Post-processor: used to visualize and quantitatively process the results from the solver. In a contemporary CFD package, the analysed flow phenomena can be presented in vector plots or contour plots to display the trends of velocity, pressure, kinetic energy and other properties of the flow.

Figure 7: Overview of CFD

When solving fluid flow problems numerically, the surfaces, boundaries and spaces around and between the boundaries of the computational domain have to be represented in a form usable by computer. This can be achieved by some arrangement of regularly and irregularly spaced nodes around the computational domain known as the mesh. Basically, the mesh breaks up the computational domain spatially; so that calculations can be carried out at regular intervals to simulate the passage of time, as numerical solutions can give answers only at discrete points in the domain at a specified time. The process

(28)

of transforming the continuous fluid flow problem into discrete numerical data which are then solved by the computer is known as discretisation. Generally, there are three major parts of discretisation in solving fluid flow

• Equation discretisation

• Spatial discretisation

• Temporal discretisation

The following figure shows the procedures of performing discretisation on a typical fluid flow problem.

Figure 8: Discretisation of flow in CFD

Equation discretisation

As mentioned above, the governing equations consist of partial differential equations.

Equation discretisation is the translation of the governing equations into a numerical analogue that can be solved by computer. In CFD, equation discretisation is usually performed by using the finite difference method (FDM), the finite element method (FEM) or the finite volume method (FVM) [26].

The FDM employs the concept of Taylor expansion to solve the second order partial differential equations (PDE) in the governing equations of fluid flow. This method is straightforward, in which the derivatives of the PDE are written in discrete quantities of variables resulting in simultaneous algebraic equations with unknowns defined at the nodes of the mesh. FDM is famous for its simplicity and ease in obtaining higher order accuracy discretisation. However, FDM only applies to simple geometries because it

(29)

employs a structured mesh.

Unlike FDM, unstructured mesh is usually used in FEM. The computational domain is subdivided into a finite number of elements. Within each element, a certain number of nodes are defined where numerical values of the unknowns are determined. In FEM, the discretisation is based on an integral formulation obtained using the method of weighted residuals, which approximates the solutions to a set of partial differential equations using interpolation functions. FEM is famous for its application around complex geometries because of the application of unstructured grid. But numerically, it requires higher computer power compared to FDM. So the finite volume method (FVM), which is mathematically similar to FEM in certain applications, but requires less computer power, is the next consideration in CFD applications.

In FVM, the computational domain is separated into a finite number of elements known as control volumes. The governing equations of fluid flow are integrated and solved iteratively based on the conservation laws on each control volume. The discretisation process results in a set of algebraic equations that resolve the variables at a specified finite number of points within the control volumes using an integration method. Through the integration on the control volumes, the flow around the domain can be fully modelled.

FVM can be used both for the structured and unstructured meshes. Since this method involves direct integration, it is more efficient and easier to program in terms of CFD code development. Hence, FVM is more common in recent CFD applications compared to FEM and FDM.

Spatial discretisation

Spatial discretisation divides the computational domain into small sub-domains making up the mesh. The fluid flow is described mathematically by specifying its velocity at all points in space and time. All meshes in CFD comprise nodes at which flow parameters are resolved. Three main types of meshes commonly used in computational modelling are structured, unstructured and multi-block structured mesh.

A structured mesh is built on a coordinate system, which is common in bodies with a simple geometry such as square or rectangular sections. However, a structured mesh performs badly when the geometry is complex, which is quite common in industrial applications. In the view of this, unstructured meshes were introduced.

In an unstructured mesh, the nodes can be placed accordingly within the computational domain depending on the shape of the body, such that different kinds of complex compu-

(30)

tational boundaries and geometries can be simulated. An unstructured mesh works well around complex geometries but this requires more elements for refinement compared to a structured mesh on the same geometry, leading to higher computing cost. To com- pensate between computing cost and flexibility, we turn our attention to the multi-block structured mesh.

In a multi-block structured mesh, the computational domain is subdivided into different blocks, which consists of a structured mesh. A multi-block structured mesh is more complicated to generate compared to a structured and an unstructured mesh but it combines the advantages of both. It is more computer efficient than an unstructured mesh and yet provides ease of control in specifying refinement needed along certain surfaces or walls, especially for meshing around complex geometries.

Temporal discretisation

The third category of discretisation is the temporal or time discretisation. Generally, temporal discretisation splits the time in the continuous flow into discrete time steps.

In transient or time-dependent formulations, we have an additional time variable t in the governing equations compared to the steady state analysis. This leads to a system of partial differential equations in time, which comprise unknowns at a given time as a function of the variables of the previous time step. Thus, unsteady simulation normally requires longer computational time compared to a steady case due to the additional step between the equation and spatial discretisation.

Either explicit or implicit method can be used for unsteady time-dependent calculation.

In an explicit calculation, a forward difference in time is taken when calculating the time tⁿ⁺¹ by using the previous time step value (n denotes state at time t and n+ 1 at time t+ ∆t ). An explicit method is straight forward, but each time step has to be kept to a minimum to maintain computation stability and convergence. On the other hand, implicit method computes values of time step tⁿ⁺¹ at the same time level in a simulation at different nodes based on a backward difference method. This results in a larger system of linear equations where unknown values at time steptⁿ⁺¹ have to be solved simultaneously. The principal advantage of implicit schemes compared to explicit ones is that significantly larger time steps can be used, whilst maintaining the stability of the time integration process . A smaller time step ∆t in an explicit method implies longer computational running time but it is relatively more accurate [9].

(31)

3 Approaches for numerical calculations of multiphase flow

A multiphase flow system consists of a number of single phase regions which are bounded by moving interfaces. In principle, a multiphase flow model could be formulated in terms of the local instant variables relating to each phase and matching boundary conditions at all phase interfaces. It is very complicate to obtain solution of multiphase system directly or in other words it is almost impossible to solve directly. As a starting point for derivation of macroscopic equations which replace the local instant description of each phase by a collective description of the phases.

For the formulation of the multiphase flow, averaging procedures can be classified into three main groups (Ishii 1975), the Boltzmann averaging, the Lagrangian averaging and the Eulerian averaging. These groups can be further divided into sub-groups based on the variable with which a mathematical operator or averaging is defined. Here we will discuss about two numerical approaches for solving multiphase flows in CFD.

3.1 Eulerian - Lagragian approach

This approach is applicable to continuous-dispersed systems and is often referred to as a discrete particle model or particle transport model. The primary phase is continuous and is composed of a gas or a liquid. The secondary phase is discrete and can be composed of particles, drops or bubbles.

In the Eulerian–Lagrangian (E–L) approach, the continuous phase is treated in an Eu- lerian framework (using averaged equations) [27]. Its continuous-phase flow field is computed by solving the Navier-Stokes equations. The dispersed phase is represented by tracking a small number of representative particle streams. For each particle stream, ordinary differential equations representing mass, momentum and energy transfer are solved to compute its state and location. The two phase are coupled by inclusion of appropriate interaction terms in the continuous-phase equations.

In this approach the volume displaced by the dispersed phase is not taken into account.

So, this approach is applicable for low-volume fractions of the dispersed phase. This approach is applicable for situations in which the discrete phase is injected as a continuous stream into the continuous phase. A force balance equation based on Newton’s second law of motion is solved to compute the trajectory of the discrete phase.

The Eulerian-Lagrangian approach is suitable to unit operations in which the volume fraction of the dispersed phase is small, such in spray dryers, coal and liquid fuel combustion, and some particle-laden flow [21]. This approach provides complete information on the behavior and residence time of individual particles. Interaction of individual particle

(32)

streams with turbulent eddied and solid surfaces such as walls can be modeled.

3.2 Eulerian - Eulerian approach

Eulerian-Eulerian approach is the most general approach for solving multiphase flows. It is based on the principle of interpenetrating continua , where each phase is governed by the Navier-Stokes equations [21]. The phases share the same volume and penetrate each other in space and exchange mass, momentum and energy. Each phase is described by its physical properties and its own velocity, pressure, concentration and temperature field.

The interphase transfer between phases is computed using empirical closure relations.

The Eulerian-Eulerian approach is applicable for continuous-dispersed and continuous- continuous systems.

For continuous-dispersed systems, the velocity of each phase is computed using the Navier-Stokes equations. The dispersed phase can be in the form of particles, drops or bubbles. The forces acting on the dispersed phase are modeled using empirical corre- lations and are included as part of the interphase transfer terms. In addition, drag, lift, gravity, buoyancy and virtual-mass effects are some of the forces that might be acting on the dispersed phase. These forces are computed for an individual particle and then scaled by the local volume fraction to account for multiple particles.

There are three different Euler-Euler multiphase models available:

• Volume of fluid method (VOF)

• Mixture model

• Eulerian model

3.2.1 Volume of fluid method

In computational fluid dynamics, the Volume of fluid method is one of the most well known methods for volume tracking and locating the free surface. The motion of all phases is modelled by solving a single set of transport equations with appropriate jump boundary conditions at the interface [1]. It can model two or more immiscible fluids by solving a single set of momentum equations and tracking the volume fraction of each of the fluids throughout the domain. Typical applications include the motion of large bubbles in a liquid, the motion of liquid after a dam break, the prediction of jet breakup, and the steady or transient tracking of any liquid-gas interface.

In general, the steady or transient VOF formulation relies on the fact that two or more

(33)

fluids (or phases) are not interpenetrating. During the numerical calculation in each control volume, the sum of the volume fractions of all phases remains to unity. In addition, the fields for all properties and variables are shared by the phases and represent volume-averaged values, as long as the volume fraction of each of the phases is known at each location. Thus in any given cell, the properties and variables are either purely representative of one of the phases, or representative of a mixture of the phases, depending upon the volume fraction values. Also in other words, if p^th fluid’s volume fraction in the cell is denoted as α_p , then the following three different conditions are possible:

• α_p = 0: The computational cell is empty of the p^th fluid.

• αp = 1: The computational cell is full of thep^th fluid.

• 0< α <1: The computational cell contains interface betweenp^th fluid and one or more other fluids available.

Volume fraction equation

In VOF method, the interface(s) tracking between the phases is established by getting the solution of a continuity equation for the volume fraction of one (or more) of the phases [1] . For thep^th phase, this equation has the following form:

1 ρ_p





∂

∂t(αpρp) +∇ ·αpρp−→υp =S+

n

X

q=1

( ˙mqp−m˙pq)



 (16)

Whereρp is the density of the p^th fluid. Alsom˙pq is the mass transfer from phase p to phase q and m˙qp is the mass transfer from phase q to phase p. This volume fraction equation will be solved for the secondary phase. It will not be solved for the primary phase. The primary-phase volume fraction will be calculated based on the following constraint:

n

X

p=1

α_p = 1 (17)

The volume fraction equation may be solved either through implicit or explicit time discretization scheme.

(34)

Momentum equation

In the VOF method, a single set of momentum equation is solved throughout the whole computational domain. The resulting velocity field is shared among the phases. Where the momentum equation is dependent on the volume fractions of all the phases through the densityρ and viscosityµ. The momentum equation is as follow:

∂

∂t(ρ−→υ) +∇ ·(ρ−→υ−→υ) =−∇p+∇ ·

µ(∇−→υ +∇−→υ^T)

+ρ−→g +−→

F (18)

Surface tension

The VOF method can also include the effects of surface tension along the interface between each pair of phases. The model can be augmented by the additional specification of the contact angles between the phases and the walls. The continuum surface tension force (CSF) of Brackbill et al. [2] have been widely used to model surface tension in multiphase flow in volume of fluid (VOF), level-set (LS) and front tracking (FT) methods [3]. The solver will include the additional tangential stress terms that aries due to the avariation in surface tension coefficient. The effect of variable surface tension are usually important only in zero/near-zero gravity conditions.

The surface tension is a force, acting only at the surface, that is required to maintain equilibrium in such instances. Surface tension aries as a result of attractive forces between molecules in a fluid. For example, consider an air bubble in water. Within the bubble, the net force on a molecule due to its neighbors is zero. At the surface, the net force is radially inward, and the combined effect of the radial components of force across the entire spherical surface is to make the surface contract, thereby increasing the pressure on the concave side of the surface. It acts to balance the radially inward intermolecular attractive force with the radially outward pressure gradient force across the surface. In regions where two fluids are separated, but one of them is not in the form of spherical bubbles, the surface tension acts to minimize free energy by decreasing the area of the interface [4].

Reconstruction based schemes

For the special interpolation treatment to the computational cells that lie near the interface between two phases, there are two reconstruction based schemes as Geo-Reconstruct and Donor-Acceptor.

Figure 9 shows an actual interface shape along with the interfaces assumed during computation by these two methods.

(35)

Figure 9: Interface calculation

• The Geometric Reconstruction Scheme

Whenever a cell is completely filled with one phase or another, the geometric reconstruction approach, the standard interpolation schemes that are used to obtain the face fluxes. When the cell is near the interface between two phases, the geometric reconstruction scheme is used. The geometric reconstruction scheme represents the interface between fluids using a piecewise-linear approach. In ANSYS FLUENT this scheme is the most accurate and is applicable for general unstructured meshes. It assumes that the interface between two fluids has a linear slope within each cell, and uses this linear shape for calculation of the advection of fluid through the cell faces. We can see in the above figure. The procedure of the geometric reconstruction approach is as follow,

• To calculate the position of the linear interface relative to the center of each partially-filled cell, based on information about the volume fraction and its derivatives in the cell

• To calculate the advecting amount of fluid through each face using the computed linear interface representation and information about the normal and tangential velocity distribution on the face

• To calculate the volume fraction in each cell using the balance of fluxes calculated during the previous step

• The Donor-Acceptor Scheme

When the cell is near the interface between two phases, a “the Donor-Acceptor” scheme is used to determine the amount of fluid advected through the face [4]. This scheme

(36)

identifies one cell as a donor of an amount of fluid from one phase and another (neighbor) cell as the acceptor of that same amount of fluid, and is used to prevent numerical diffusion at the interface. The amount of fluid from one phase that can be convected across a cell boundary is limited by the minimum of two values: the filled volume in the donor cell or the free volume in the acceptor cell. The orientation of the interface is also used in determining the face fluxes. It is either horizontal or vertical, depending on the direction of the volume fraction gradient of the p^th phase within the cell, and that of the neighbor cell that shares the face in question. The flux values are obtained by pure upwinding, pure downwinding, or some has a combination of the both. In addition it depends on the interface’s orientation as well as its motion.

3.2.2 Mixture model

The mixture model is a simplified multiphase model that can be used in different ways.

The mixture model can apply to model multiphase flows where the different phases move at different velocities and also it is applicable to model homogeneous multiphase flow and to calculate non-Newtonian viscosity.

The mixture model can model n phases (fluid or particulate) by solving both the continuity equation and the momentum equation for the mixture, where mixture can be a combination of continuous phase and the dispersed phase. In addition, the mixture model solves the energy equation for the mixture, and the volume fraction equation for the secondary phases, as well as algebraic expressions for the relative velocities (if the phases are moving at different velocities). Also it allows us to select the granular phases and we can calculate the different properties for granular phases. It is applicable in the particle-laden flows with low loading, and bubbly flows where the gas volume fraction remains low, cyclone separators, sedimentation and in liquid-solid flows [21]

Continuity equation for the mixture

First of all, we can write the continuity equation as follow from Ishii (1975),

∂

∂t(α_kρ_k) +∇ ·(α_kρ_k−→υ_k) = Γ_k (19) where αk is the volume fraction of the k^th phase and the term Γk represents the rate of mass generation of phasek at interface. From the above equation(19), we obtain by summing over all phases,

(37)

∂

∂t

n

X

k=1

(α_kρ_k) +∇ ·

n

X

k=1

(α_kρ_k−→υ_k) =

n

X

k=1

Γ_k (20)

Because of the mass conservation the right hand side of the equation (20) must be vanish,

so n

X

k=1

Γk= 0 (21)

From it, we obtain a continuity equation of the mixture as

∂ρm

∂t +∇ ·(ρm−→υm) = 0 (22) Here the mixture density and the mixture velocity are defined as,

ρ_m =

n

X

k=1

α_kρ_k (23)

−

→υm = 1 ρm

n

X

k=1

(αkρk−→υk) =

n

X

k=1

ck−→υk (24)

The mixture velocity um represents the velocity of the mass centre. Also ρm varies although the component densities are constants. The mass fraction of phase k is defines as

ck= αkρk

ρm

(25) equation(22) has the same form as the continuity equation for the single phase flow.

Momentum equation for the mixture

The momentum equation fork^th is written as,

∂

∂t(α_kρ_k−→υ_k) +∇ ·(α_kρ_k−→υ_k−→υ_k) =−α_k∇p_k+∇ ·[α_k(τ_k+τ_{T k})] +α_kρ_kg+M_k (26) whereM_k is the average interfacial momentum source for the phasek. τ_k is the average viscous stress tensor and τ_{T k} is the turbulent stress tensor.

(38)

Now the momentum equation for the mixture follows (26) by summing over the phases

∂

∂t

n

X

k=1

α_kρ_k−→υ_k+∇ ·

n

X

k=1

α_kρ_k−→υ_k−→υ_k

=−

n

X

k=1

α_k∇p_k+∇ ·

n

X

k=1

α_k(τ_k+τ_{T k}) +

n

X

k=1

α_kρ_kg+

n

X

k=1

M_k (27)

By using the definitions (23) and (24) of the mixture densityρ_mand the mixture velocity

−

→υ_m, the second term of equation(27) can be written as

∇ ·

n

X

k=1

α_kρ_k−→υ_k−→υ_k=∇ ·(ρm−→υm−→υm) +∇ ·

n

X

k=1

α_kρ_k−→υ_{M k}−→υ_{M k} (28)

where −→υ_{M k} is the diffusion velocity, i.e.,the velocity of the k^th phase relative to the centre of the mixture mass

−

→υ_{M k}=−→υ_k− −→υm (29) The momentum equation takes the form in terms of the mixture variables as,

∂

∂tρm−→υm+∇ ·(ρm−→υm−→υm) =−∇p_m+∇ ·(τm+τ_{T m}) +∇ ·τ_Dm+ρmg+M_m (30) Whereτ_m,τ_{T m}, and τ_Dm are the stress tensors.

Volume fraction equation for the secondary phases

We can obtain the volume fraction equation for secondary phase p from the continuity equation of the secondary phase pas follow:

∂

∂t(αpρp) +∇ ·αpρp−→υm =−∇ ·αpρp−→υ_dr,p+

n

X

q=1

( ˙mqp−m˙pq) (31)

In the mixture model some terms as relative (slip) velocity and drift velocity, interfacial area concentration, solid pressure and granular properties like collisional viscosity, kinetic viscosity, granular temperature are also important.

(39)

3.2.3 Eulerian model

Eulerian model is the most general model for solving multiphase flows. In the present work, we are using Eulerian model to simulate two-phase and three-phase flow. The Eulerian model is the most complex of the multiphase models. It solves a set of n momentum and continuity equations for each phase. In the next section, there are more details about Eulerian model.

3.3 Eulerian model theory

As we discussed above, in the Eulerian approach all the phases are treated as continuum.

Eulerian model solves continuity, momentum and energy equations for each phase.

3.3.1 Volume fraction equation

The description of multiphase flow as interpenetrating continua incorporates the concept of the phasic volume fractions which is denoted as α_q. The volume fractions represent the space occupied by each phase and the laws of conservation of mass and the conservation of momentum are satisfied by each phase individually. The derivation of the conservation equations can be done by averaging the local instantaneous balance for each of the phases or by using the mixture theory approach.

The volume of the phasep is denoted byV_p and is defined by,

Vp= Z

V

αqdV (32)

Where the total volume fraction in the computational cell remains one as,

n

X

p=1

α_p = 1 (33)

The effective density of the phasep is defined by,

ρbp =αpρp (34)

whereρp is the physical density of the phase p.

The volume fraction equation may be solved by either implicit time discretization or

(40)

explicit time discretization.

3.3.2 Conservation equations

The volume-averaged mass and momentum conservation equations in the Eulerian framework are given by following way.

1. Conservation of mass

The continuity equation for the phase q is written as follow,

∂

∂t(α_qρ_q) +∇ ·(α_qρ_q−→υ_q) =

n

X

p=1

( ˙m_pq−m˙_qp) +S_q (35) In the above equation, the velocity of phase q is given by −→υq and m˙pq is characterized by the mass transfer from thep^th toq^thphase andm˙qp is the mass transfer from phaseqto phasep. The last term of the right hand side in Equation (35),Sq, is the source term. In addition, it can be specified as a constant or by user-define function for each phase.

2. Conservation of momentum

The momentum balance for the phaseq is given by the following,

∂

∂t(αqρq−→υq) +∇ ·(αqρq−→υq−→υq) =−α_q∇p+∇ ·τq+αqρq−→g+

n

X

p=1

(−→

Rpq+ ˙mpq−→υpq−m˙qp−→υqp) + (−→ Fq+−→

F_{lif t,q}+−→

Fvm,q) (36)

Where−→

F_q is an external body force, −→

F_{lif t,q} is a lift force and −→

F_vm,q is a virtual mass force. −→

R_pq is an interaction force between phases and the pressure, which is shared by all phases, is given byp.

−

→υ_pq is the interphase velocity which can be defined as follows.

• Ifm˙pq >0(i.e., phasep mass is being transferred to phaseq),−→υpq =−→υp

• Ifm˙_pq <0(i.e., phaseq mass is being transferred to phasep),−→υ_qp =−→υ_q

• Ifm˙_qp >0,−→υ_qp =−→υ_q

• Ifm˙_qp <0,−→υ_qp =−→υ_p

(41)

In addition, the interaction force between phases −→

Rpq depends on the friction, pressure, cohesion, and other effects. In literature, the simple form of interaction term is given by following way,

n

X

p=1

−

→R_qp =

n

X

p=1

K_pq(−→υ_p− −→υ_q) (37) WhereKpq(=Kqp)is the interphase momentum exchange coefficient.

In Equation (36),τqis the stress-strain tensor of theq^thphase and it can be written as follow,

τ_q=α_qµ_q(∇−→υ_q+∇−→υ^T_q) +α_q(λ_q−2

3µ_q)∇ · −→υ_qI (38) In the above equation,µ_q andλ_q express the shear and bulk viscosity of the phase q.

3. Conservation of energy

The conservation of energy in Eulerian multiphase modeling approach is described as a separated enthalpy equation for each phase as follow,

∂

∂t(αqρqhq)+∇·(α_qρq−→υqhq) =αq

∂pq

∂t +τq :∇−→υq−∇·−→q_q+Sq+

n

X

p=1

(Qpq+ ˙mpqhpq−m˙qphqp) (39) Wherehqis the specific enthalpy of theq^thphase,−→qqis the heat flux,Sqis a source term, which is a source of enthalpy due to chemical reaction or radiation. The term Qpq displays the intensity of heat exchange between the p^th and q^th phases, and hpq is the interphase enthalpy.

3.4 Turbulence models

In the study of the multiphase flow, it is important to include turbulence in it. To solve and to describe the effects of turbulent fluctuations of velocities and scalar quantities of flow, there are various closure models of turbulence are available. In comparison to single-phase flows, the number of terms to be modeled in the momentum equations in multiphase flows is large, and this makes the modeling of turbulence in multiphase simulations extremely complex([4]) .

In the present work, we have used ANSYS FLUENT 12.1. It provides three methods for modeling turbulence in multiphase flows within the context of the κ− models. In addition, there are two turbulence options within the context of the Reynolds stress models (RSM). Theκ−turbulence model has the following options,

CFD Simulation of Two-phase and Three-phase Flows in Internal-loop Airlift