Contents

Chapter 1
Interface Tracking Methods

1.1 The general approach

Interface Tracking Methods (ITM) are used for free-surface and interfacial flows which involve two or more immiscible fluids (e.g., liquid-liquid or liquid-gas systems) separated by sharp interfaces which evolve in time. Interface tracking/capturing schemes are methods that are able to locate the interface by capturing the interface by keeping track of the evolution of an appropriate field.

The formulation is based on solving a single set of transport equations for the whole computational domain and treating the two phases as a single fluid with variable material properties. Changes in these properties are taken into account by advecting the phase-indicator function χ(x,t):

(1.1)

where the interface velocity u^I reduces in isothermal flow conditions (no phase change) to the local fluid velocity u.

Mass exchange terms and capillary forces are represented as momentum and energy sources, described as delta functions or smoothed gradients at or across the interface. The mass, momentum and energy balance equations for incompressible interfacial flow are written as

1.2 The VOF methodology

The VOF approach relies on the definition of the liquid volume-fraction field occupied by one of the phases within volume V , and is conventionally denoted in a discrete form by F_ijk:

(1.9)

In the VOF context, Eq. (1.1), therefore, represents the evolution of the liquid volume-fraction, identifying flow regions containing pure liquid (where F_ijk = 1) from pure gas flow regions (where F_ijk = 0). Interfacial cells are such that 0 < F_ijk < 1. The VOF method does not amount solely to the solution of:

(1.10)

It requires accurate algorithms for advecting the volume fraction function so as to preserve the thickness of the interface. In order to keep the interface sharp a special advection scheme referred to as the Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) (Ubbink, 1997; Ubbink & Issa, 1999) is used.

1.3 The Level Set method

1.3.1 The scheme

The Level Set approach (Osher & Sethian, 1988) consists in solving a hyperbolic equation:

(1.11)

to track the interface on a fixed Eulerian grid, where ϕ(x,t) is a smooth signed-distance function referring to the shortest distance to the front. Negative values correspond to one of the fluids and positive values to the other. The exact location of the interface corresponds to the zero level of ϕ. To update material properties like density, viscosity and thermal conductivity, a Heaviside function H(ϕ) is introduced and defined by

(1.12)

This implies that ϕ_ijk is linked in discrete form to the liquid volume fraction field F_ijk through

(1.13)

A modified Heaviside function denoted by H(ϕ) is employed to smooth the physical properties across an interface thickness of 2δ.

(1.14)

where H^k(ϕ) is the modified Heaviside function defined by

(1.15)

The same modified Heaviside function is concurrently employed to determine the surface tension.

1.3.2 Re-distancing scheme

As numerical errors cause the contours of the Level Set field to deform as the phase moves, a re-distancing algorithm is required to regularise the function: advecting the initial distance function ϕ(x,0) will not be maintained as such. An extra re-distancing algorithm preserving |∇ϕ| = 1 around the zero level of ϕ is required.

The scheme of Sussman et al. (1994) consists of advecting an intermediate equation for a correction field d(x,t) during a period of intermediate time :

(1.16)

with sgn(ϕ) = 2H(x) - 1 denotes the signum function, d₀(x,t) = ϕ(x,t), and

(1.17)

The last term in the above expression is added to keep the bubble volume constant during the reinitialisation procedure. The equation is iterated until |∇d| = 1. The corrected ϕ–field is then obtained by setting ϕ(x,t) = d(x,ε), where ε is the time elapsed for Eq. (1.16) to reach convergence. In TransAT, Eq. (1.16) is solved after each advection step of the Level Set equation using the non-oscillatory (WENO) 3^rd order scheme.

1.3.3 Fast-marching re-distancing

This is an extreme case of the narrow band method, in that the band in question is taken as a single cell. Of interest in Level Set methods is to preserve the normal distribution of the iso-contours of function ϕ(x) through the Eikonal equation:

(1.18)

with boundary conditions ϕ(x) = g(x),x ∈ Γ ⊂ Rⁿ. In the simplest case where the speed function f(x) is a constant, usually equal to one, then the solution of the Eikonal equation above will represent the shortest distance from a point x to the zero distance curve, usually given by Γ, that is if g(x) = 0. To construct an approximation to the Eikonal equation for a grid point x using the surrounding points x_i, x_n, first define the unit directional vector pointing from point x_i to point x as P_i = (x-x_i)∕||x-x_i||. Then let v_i(x) be the value of the approximate directional derivative in direction P_i and point x.

We then solve a quadratic equation for the distance ϕ(x) in the form

(1.19)

where Q is defined as Q = (PP^T )^-1 and the vectors a and b by their respective constituents a_i and b_i. This equation is then solved to find two roots of which the largest one is taken as a possible new distance for ϕ(x). This new distance can only be accepted if the update comes from within the polygon spanned by x,x_i,c…,x_n and if the calculated distance value is smaller than the old one. Thus we only update if ϕ(x) < ϕ^old(x).

1.3.4 Thin-film boundary conditions

For flows with very low Capillary numbers, there exists a thin water film between the gas phase and the wall. A thin-film boundary condition can then be defined. The viscosity at the film–boundary wall is set so as to satisfy the shear according to Couette flow. For a liquid film condition we get (see Fig. 1.1),

(1.20)

The interface velocity (tangential to the wall) u^I is then given as,

(1.21)

The effective viscosity is then set according to the following equation,

(1.22)

Or simplifying Eq. (1.22) we obtain,

(1.23)

1.3.5 Marangoni effect

Marangoni force can be modelled, which will make the surface tension coefficient depend on temperature

(1.24)

where, σ₀ is the surface tension coefficient at temperature T₀, and T_crit is the temperature of the material at which surface tension is zero.

The thermo-capillary (Marangoni) contribution to the Navier-Stokes equation is the projection tangential to the interface of the shear stress generated by the temperature dependence of the surface tension and writes:

(1.25)

1.3.6 Dynamic contact angle treatment

The treatment of wetting dynamics is based on the physical forces associated with triple lines (Friess & Lakehal, 2004; Narayanan & Lakehal, 2006). An additional triple-line contact force (c_i) is included in the momentum equation; this extended momentum equation provides a physically adequate description of wetting dynamics, eliminating the need for any particular boundary condition specifying the contact angle. The triple-line force (c_i) is based on a consideration of interfacial free energy. Accordingly, it contains only two parameters: the interfacial tension between the fluids γ and the equilibrium contact angle θ_eq.

(1.26)

where θ_dy is the instantaneous dynamic contact angle, δ^t is a smooth Dirac delta function vanishing everywhere except near the triple lines, with the property that for any volume V , the integral ∫ _V δ^tdV is equal to the length of the triple line segment contained in V , and b is the unit vector normal to the triple line and parallel to the wall surface. The triple-line force, is obtained by considerations similar to the derivation of Youngs Law and is referred to as the unbalanced Young force (De Gennes, 1985).

In addition, the shear stress is assumed to have a logarithmic profile as it nears the contact line (Qian et al., 2004). At distances less than a slip length from the triple line, full slip is assumed. In the finite volume containing the triple line, the integrated shear stress is applied. The slip length is taken to be a small value ≈ 10^-9. This formulation resolves the stress singularity that arises due to the no–slip boundary condition.

1.4 Phase change modelling

In the presence of phase change, equation (1.11) is rewritten

(1.27)

where ṁ^′′ stands for the mass transfer per unit area per unit time. A new source term also appear in the velocity divergence equation (1.2), which becomes

(1.28)

Finally, a new source term is added to the energy equation (1.4)

(1.29)

It now remains to model the mass transfer term ṁ^′′.

1.4.1 Direct phase change

When using direct phase change, mass transfer is given by the formula

(1.30)

(1.31)

with the latent heat.

1.4.2 Correction of saturation temperature

The saturation temperature can be locally adjusted to take into account the effect of curvature, using the equation

(1.32)

where, F(ρ) = ρ_G^-1 + ρ_L^-1, κ is the local curvature, σ is the local surface tension coefficient, and L is the latent heat of phase change.

1.4.3 Correlation-based models

Different models are based on global correlation. In this case, local ṁ do not vary in space and time, but the global value will change because of the changes of interfacial area.They are all based on the same expression

(1.33)

where Nu_L is the liquid Nusselt number, is the latent heat and T_sat is the saturation temperature. T_sc and L_L,ṁ are given parameters. Nu_L can be given by different correlations.

• In the Hiemenz correlation (Davis & Yadigaroglu, 2004), it is written

  (1.34)

  where

  (1.35)

  is the liquid Reynolds number, and u_L,ṁ is a correlation velocity.
  

• Other correlations have also been proposed by Lim et al. (Lim et al., 1984):

  (1.36)

  for smooth interface, and

  (1.37)

  for wavy interface. In these expressions,

  (1.38)

  is the liquid Reynolds number.

1.4.4 Interface Mass Transfer (IMT)

This section describes the interfacial mass transfer models that have been implemented in TransAT (Lakehal & Labois, 2011a). For the sake of simplicity, all of them will be written in the form

(1.39)

where ṁ is the mass transfer rate, h is the heat transfer coefficient, u′ is a characteristic velocity of the flow, Pr is the Prandtl number, Sc is the Schmidt number and Re_t a turbulent Reynolds number.

• The surface-divergence model (Banerjee et al., 2004; Lakehal et al., 2008b) can be written:

  (1.40)

  where C = 0.35 if Sc (or Pr) < 4, 0.45 otherwise. This model is adapted for low Reynolds number. It is also implemented with modifications, for the small eddy context

  (1.41)

  This model is adapted for very high turbulent Reynolds numbers. Finally, an adaptive surface divergence version has been implemented

  (1.42)

  with

  (1.43)

  where Re_t,int is the turbulent Reynolds number at the interface. This model is adapted for high Reynolds numbers.

• The small-eddy model (Banerjee et al., 1968; Theophanous et al., 1976) can be written

  (1.44)

• The surface renewal model (Lakehal et al., 2008a) links the mass transfer with the interfacial frictional Reynolds number Re_⋆. It reads

  (1.45)

  where the following expression is used for C

  (1.46)

• The Hughes-Duffey model (Hughes & Duffey, 1991) reads

  (1.47)

• The last model from (Coste, 2004) is inspired by the renewal model and reads

  (1.48)

A turbulent Reynolds number is needed for the interfacial mass transfer model. Several choices are available in TransAT:

• The usual Reynolds number can be written Re_t = . It is taken in the gas core flow. A length scale must be associated to this Reynolds number, it is computed as the weighted-average (in the gas side) of ε, with the characteristic velocity being u_∞ = min(|u|,C_μ^1∕4k^1∕2)
• As an alternative to Re_t, another Reynolds number based on turbulence is available, R_y = . It is also evaluated in the gas core flow.
• The shear Reynolds number reads Re_⋆ = with u_⋆ = the shear velocity at the interface, ϕ the distance of the point to the interface. u_⋆ and ϕ are taken in the gas near the interface, where u_⋆ is at its maximum. u_⋆ is also taken as the associated characteristic velocity.

Because they need the variables k and ε to be evaluated, Re_T and R_y can be used only when the k -ε turbulence model is enabled. Otherwise, Re_⋆ together with the surface renewal model will be used. The following table gives the available interfacial mass transfer models available for the different Reynolds numbers.

Chapter 2
Phase-Average Multi-phase Flow Modelling

2.1 Ensemble-averaged Two-phase Flow Models

We consider a mixture of two incompressible phases. Thus, the following equations are solved (Drew & Passman, 1998):

(2.1)

(2.2)

where the subscript G denotes the gas phase, m denotes a mixture property, and T stands for mixture turbulent property. In the following, the subscript L will also be introduced for the liquid phase. α_G is the gas volume fraction and ρ^m = α_Gρ_G + α_Lρ_L is the density of the mixture. V ^m is the velocity of the mixture, p^m the mixture pressure, and T^m is the mixture temperature. S_αG and S_mom^drift stand for the source terms due to mass transfer and drift velocities in the mass and momentum equations, respectively.

The phasic ensemble averaged velocities split into mixture and drift velocities as follows,

(2.3)

The mixture velocity is defined as follows:

(2.4)

From Eq. (2.4) it follows that,

(2.5)

and the phasic drift velocities are related to each other by,

(2.6)

where Y _k = α_kρ_k∕ρ^m is the mass fraction of the phase k. The slip velocity defined as V _si = V _Gi -V _Li can be further related to the gas drift velocity as,

(2.7)

2.1.1 Homogeneous Model

In the homogeneous model, it is considered that both phases have the same velocity (V_k^d = 0). The equations can then be simplified to,

(2.8)

Only the mixture equations are solved. The divergence of the mixture velocity is zero.

(2.9)

2.1.2 Algebraic Slip Model

In this case, an algebraic model is used to estimate the relative velocity between the two phases or between the mixture and the gas phase. Thus, several drift velocity models have been implemented. With this notation, the algebraic slip model (ASM) can be written as follows,

(2.10)

(2.11)

(2.13)

Models for gas drift velocity V _G^d must now be derived.

2.1.3 Algebraic slip model with drag forces

Ishii & Mishima (1984) have written the drag-induced momentum transfer term as

(2.14)

where C_D is the drag coefficient and R_b the bubble radius.

Then, Manninen & Taivassalo (1996) derived from the conservation equations the following expression, neglecting diffusion terms due to turbulent fluctuations and assuming the local equilibrium.

(2.15)

From these equations, an expression for the gas drift velocity can be found, assuming that all bubbles have the same radius R_b.

(2.16)

where Re_b = is the bubble Reynolds number. Introducing the friction factor , the drag coefficient can be written

(2.17)

and the drift velocity is rewritten

(2.18)

Friction factor remains to be modelled. In the Stokes Regime, it is simply

(2.19)

This model is adapted for very small bubbles Reynolds number, usually less than 2.

Schiller & Naumann (1933) proposed the friction factor coefficient as:

(2.20)

Clift (cited in Crowe et al. (1998)) developed a model that is adapted to bubble Reynolds numbers up to 1000. In this model, friction factor is written

(2.21)

The model derived by Tomiyama (1998a) also takes into account the shape of the bubbles, using the Eötvös number Eo = , where ||g|| is the gravitational acceleration and σ the surface tension. The friction factor is then given (for highly contaminated system here):

(2.22)

2.1.4 Algebraic slip model with drag and lift forces

The expression of the lift force is given by Drew & Passman (1998)

(2.23)

Momentum transfer term 2.14 can then be rewritten as

(2.24)

and solving this equation will give the gas drift velocity. In this expression, lift coefficient is still unknown. It can be taken as a constant, with values of (Drew & Passman, 1998) or (Drew & Lahey, 1987) proposed in the literature for bubbly flows. Tomiyama (1998b) derived a more complex model, in accordance with deformable bubbles correlation. It evaluates the long axis d_H of a deformable bubble, and proposes to compute the lift coefficient based on this axis. The model can then be written

(2.28)

2.1.5 Wall lubrication force

When the lift force is used an unrealistic accumulation of small bubbles at the wall could appear. To counterbalance this effect an additional wall lubrication force is typically included. The expression of the wall lubrication force is given by Rzehak et al. (2012):

(2.29)

where n_W is the unit vector perpendicular to the wall, pointing into the fluid. The second part of the force was added to keep in account the fact that the lubrication force acts only close to the wall in a exponentially way and not at the center of the vessel. Momentum transfer term 2.14 can then be rewritten as:

(2.30)

(2.31)

Based on the observation of single bubble trajectories, Tomiyama et al. (1995) proposed this formulation:

(2.32)

with the coefficient C_W that depends on the Eotvos number. According to the Tomiyama formulation, the coefficient C_W of the Eq. 2.32 is:

(2.33)

Another correlation of the C_W , depending on the Eotvos number, was proposed by Hosokawa et al. (2002) as:

(2.34)

Frank et al. (2004) generalized the Tomiyama model to produce the Frank Wall Lubrication force model, that is given by:

(2.35)

where C_W is the coefficient defined in 2.33, that preserves the same dependence on Eotvos number as the Tomiyama model.

2.1.6 Turbulent dispersion

The drift velocity models of the previous sections do not take into account the effect of the turbulence. In this aim, Ishii (1975) proposed a term to be added to the gas drift velocity:

(2.36)

where ν_T is the turbulent kinematic viscosity and Pr_T is the turbulent Schmidt number.

2.2 Phase change modelling

Phase change phenomenon is divided into phase change happening in the bulk and phase change happening near a wall (heated or cooled). The latter is called heterogeneous phase change.

2.2.1 Phase change in the bulk

Phase change happening in the bulk of the fluid (away from the flow boundaries) can be phenomenologically split into different phase change modes, such as,

• Dispersed phase change
• Interfacial phase change
• Cavitation
• Homogeneous condensation

where the first two are thermally limited phase change phenomena and the last two are inertial phase change phenomena where the spontaneous phase change is related to rapid pressure changes. Currently, homogeneous condensation models are not available in TransAT.

Dispersed-phase phase change

This mode of phase change is defined as phase change that happens under thermal non-equilibrium under conditions where the phases are well mixed. The phase change rate is limited by the temperature difference (T - T_sat and the interfacial area density. The model available in TransAT is based on the Ranz-Marshall Nusselt number model (Ranz & Marshall, 1952). The mass transfer rate (kg∕m³) is given as,

Interfacial phase change

The main difference between the interfacial and the dispersed phase change models is the fact that for conditions where the phases are separated, one can obtain an estimate for the real interfacial area using the gradient of the volume fraction. Interfacial phase change has been split into two cases in TransAT,

1.

Interfacial thermal diffusion phase change model, and

2.

Interfacial shear phase change model.

The interfacial shear phase change model is based on the proposal by Banerjee (Lakehal & Labois, 2011b) where,

Substitution the expression for Nu into the heat flux equation we obtain,

(2.42)

where C is the only model parameter. Note that this closed form expression is a special case due the fact that the exponent of the Reynolds number in the Nusselt number expression is 1.

In case the velocities are zero, the interfacial shear phase change model will predict zero mass transfer rate. The phase change in this case can occur due to thermal non-equilibrium which propagates by diffusion. Two models are available for interfacial diffusion phase change, viz. the Schrage model and the relaxation model.

For the Schrage model the heat transfer rate is given as,

(2.43)

where, A_c is the accommodation coefficient that is a material constant.

For the relaxation model (which is only valid for very diffusive problems) the heat transfer rate is given as,

(2.44)

When the dispersed and interfacial models are both active, then the two models are blended based on the magnitude of ∇α, where the interfacial model is chosen if the gradient is significant and vice-versa.

Cavitation models

We start with the condition of vapour bubbles in liquid of radius R_B with a nucleate density per unit volume of liquid (n₀). If we consider a volume Ω = Ω_V + Ω_L, then,

The bubble radius can be expressed using Eq. (2.45) as,

(2.47)

The other starting point for cavitation models is the simplified form of the Rayleigh-Plesset equation for bubble growth given as,

(2.48)

The radius grows if p < p_V and shrinks if p > p_V . An expression for dα_V ∕dt using Eq. (2.45) is obtained by assuming that n₀ does not change with time,

(2.49)

The definition of mixture density (ρ = ρ_L + α_V (ρ_V - ρ_L)) is used to simplify the mixture continuity equation to obtain,

(2.50)

under the assumption of incompressibility. Further in the absence of slip between the vapour and liquid phases (homogeneous flow), we can write the above equation as,

(2.51)

Under the assumptions of incompressibility and homogeneous flow the following expression for the divergence of the velocity field is obtained,

(2.52)

where the subscript e denotes evaporation. Using Eqs. (2.51,2.52,2.49) we obtain an expression for the evaporation rate as,

(2.53)

Sauer cavitation model

The model of Yuan et al. (2001) is derived by plugging Eqs. (2.48,2.47) in to Eq. (2.53) to obtain,

Singhal cavitation model

Eq. (9) in Singhal et al. (2002) was found to be not correct. This means that Eq. (10) in Singhal et al. (2002) should actually be,

(2.59)

using the notation of Singhal et al. (2002). In fact, with this correction the basic unadulterated expression for the mass transfer rate (Singhal et al. (2002):Eq. (13)) is the same as Eq. (2.53)

In this model the closure for R_B is provided using critical Weber number considerations based on stable bubble size under aerodynamic drag.

(2.60)

where D_d∞ is the bubble diameter at the end of the fragmentation process (stable diameter). Typically a value between 5 < We_d∞ < 20 is quoted in the literature (Kolev, 2002). By rearranging equation 2.60 for D_d∞, we get

(2.61)

Substituting Eq. (2.61) into Eq. (2.53) we obtain

(2.62)

In order to derive the exact model given in Singhal et al. (2002), we rearrange the equation as follows,

(2.63)

Expressing the term in the square brackets as CV _ch∕σ, where C is a (dimensional with units of velocity) empirical constant, and V _ch is a characteristic velocity, we get,

(2.64)

The final form of the model is obtained by using different coefficients for the evaporation and condensation rates. In addition, the volume fraction of vapour is removed from the evaporation rate equation and the volume fraction of of liquid is removed from the condensation rate equation. This is necessary since for a pure liquid α_V is initially zero thereby making the evaporation rate zero. Introducing these changes, we get

Implementation in TransAT

The evaporative mass transfer rate is written using an indicator function due to the change in coefficients and the appropriate volume fraction as follows.

(2.73)

where the indicator function I_e is 1 if p < p_V and 0 other wise and can be expressed as

(2.74)

2.2.2 Wall Heat Flux Modeling during Subcooled Nucleate Boiling

General Formulation

For the liquid phase, the heat flux from the wall to the liquid are transferred through the following means: heat convection to the liquid phase away from the nucleation sites, microlayer evaporation and quenching. In nucleate subcooled boiling, the wall heat is utilised in forming bubbles and heating the liquid. According to Kurul & Podowski (1991), the total convective heat flux from the wall is the sum of three models:

(2.75)

where q_1ϕ^′′ is the single-phase convective heat flux, q_e^′′ is the heat flux associated with phase change and q_Q^′′ is the quenching flux, during the dwell time.

Single-Phase Heat Flux

In the portion of the wall unaffected by the bubbles, the heat transfer from the wall to the liquid can be determined by:

(2.76)

Now, A_1ϕ^′′ is the fraction of the wall unaffected by the nucleation sites, ρ_l is the liquid density, c_pl is the liquid specific heat, T_w is the wall temperature, T_l is the temperature of the liquid. St_P is the Stanton Number calculated from heat transfer correlation in terms of the local liquid velocity and Prandtl Number. St_P could be written as,

(2.77)

where Nu, the Nusselt Number, is provided by the Ranz-Marshall correlation.

(2.78)

or from the correlation found in Vyskocil & Macek (2008)

(2.79)

Re, Pr and Pe are the Reynolds, Prandtl and Peclet numbers respectively.

Evaporative Heat Flux

Taking into account the effect of area overlapping between the neighboring sites, the the evaporative heat flux , q_e^′′ is given by E. Krepper & Egorov (2007), becomes

(2.80)

where ρ_v is the vapor density, h_fg is the latent heat of vaporisation, a_b² is the bubble influence factor having a value of 2. The area fraction under the influence of bubbles, A_2ϕ^′′ is given by

(2.81)

where d_det is the bubble diameter at departure, which depends on the local liquid subcooling:

(2.82)

Here, ΔT_sub is the difference between the saturation temperature and the liquid temperature.

(2.83)

In Eq. 2.80, f_det is the frequency of nucleation. In general, A simpler estimation of the bubble departure frequency as ratio of the terminal rise velocity over the bubble departure diameter is:

(2.84)

where g is the acceleration due to gravity. N^′′ is the nucleation site density

(2.85)

Quenching Heat Flux

(2.86)

where h_Q is the quenching heat transfer, given by

(2.87)

where, k_l is the liquid thermal conductivity. Here, t_w is the waiting time between the departure and the inception of the next bubble. It is assumed that the waiting time takes 80% of the bubble departure period. Hence,

(2.88)

To close the model, Eq. 2.75 can be solved iteratively for the local wall temperature T_w

2.3 Ensemble-averaged N-phase flow

2.3.1 N-phase model

A mixture of N phases can be simulated with the following set of equations. Definitions are given below

2.3.2 Turbulence equations

For modelling the effect of turbulence for the N-phase flow, only the turbulent kinetic energy k and dissipation rate ε of the mixture are solved for. Also, it is assumed that all the constants of the k - ε model are the standard values for single phase flow. Equations for the turbulent variables can then be written

2.4 Mixture viscosity

2.4.1 Suspension viscosity model

In a suspension, the flow is disturbed by the presence of a dispersed phase. If the concentration of the dispersed phase increases, the particle or droplet will feel an increased resistance of the flow, which results in a higher drag coefficient. This can be taken into account by modifying the mixture viscosity (Ishii & Zuber, 1979; Hallanger et al., 1995; Manninen & Taivassalo, 1996). In this section, subscript c denotes continuous phase, while d denotes dispersed phase.

2.4.2 Suspended drops or bubbles

For a mixture including suspended drops or particles, Ishii & Zuber (1979) developed an apparent viscosity expression, which reads

(2.100)

2.4.3 Solid suspensions

Several correlations are available for suspensions of solid particles. Ishii & Zuber (1979) propose the following expression

(2.101)

with α_CP is the close packing volume fraction, usually taken between 0.62 and 0.65 for solid particles.

Mooney (cited in Manninen & Taivassalo (1996)) proposed the following correlation

(2.102)

Finally, Mills (1985) derived the following equation for a suspension of hard spheres of equal sizes

(2.103)

with α_CP = .

Chapter 3
Particle Tracking Models

3.1 Mathematical formulation

3.1.1 Momentum equations

The equation governing the particle motion is the well known example given by Maxey & Riley (1983). The equation is valid for a rigid sphere in non-uniform flow such that the sphere is isolated and far from any boundary, the Reynolds number for the relative motion between the particle and the fluid is small, and the particle size is much smaller than the smallest flow length scale. The equation is given as follows,

Other forms of the forces and additional forces can be used under different circumstances.

Drag force

In the Stokes regime, drag force applied to the particle is written

(3.2)

For a higher Reynolds number based on particles, the Clift-Gauvin correlation can be used, where

(3.3)

with a friction factor f being equal to

(3.4)

with the particle Reynolds being defined by

(3.5)

Lift force

The standard small particle (under the Stokes flow, Re_p << 1 condition) does not include a lift force under the assumption of uniform flow upstream of the particle. However under non-uniform fluid flow upstream of a particle a lift force exists as proposed by Saffman,

(3.6)

Auton (1987) derived the lift force acting on a sphere in shear flow under inviscid flow conditions as,

(3.7)

For small particles, both of these lift forces are small when compared to the drag and buoyancy forces. The lift force is also small if the dispersed phase is much heavier than the fluid (solid particles or droplets in gas). However, for bubbly flow the lift force could be significant.

Magnus Lift

The lift force due to particle rotation (Magnus force) is given for low Reynolds number by Rubinow & Keller (1961) as,

(3.8)

where ω_p is the particle rotation, and ∇× is the fluid rotation at the particle position. Rewriting it in a general form with a lift coefficient C_L,Mag we get,

(3.9)

where A is the projected area of the particle, and ω_r = ω_p -∇× is the relative rotation between the particle and the fluid phase. The lift coefficient for low Reynolds number is then given as

(3.10)

Lift coefficient of Oesterl & Dinh (1998)

(3.11)

where γ, which is called the spin rate, is given as

(3.12)

where Re_p is the particle Reynolds number.

Two-way Coupling

When two-way coupling is activated the forces acting on the particles also enter the Navier-Stokes equations by means of an additional fluid-particle interaction force F^fp. For example the incompressible Navier-Stokes equations become

(3.15)

where n is the particle index, N_p is the total number of particles, V _p is the particle volume, ρ_f is the fluid density, V _m is the volume of the finite-volume cell, f_i is the component of the fluid-particle interaction force in direction i on a single particle located at xⁿ and W is the projection weight of the force onto the cell node m located at x^m.

3.1.2 Particle Angular momentum equation

(3.16)

where I_p is the moment of inertia, and T_i is the torque acting on the particle. The torque acting on a spherical particle in Stokes flow is given as

(3.17)

For higher Reynolds numbers (up to 1000)

(3.18)

where Re_ω = ρ_fω_pd_p²∕4μ_f.

3.1.3 Temperature equation

Each particle is assumed to be represented by a uniform temperature T_p. This assumption is valid under low particle Biot (Bi_p) number conditions, where Bi_p = Nu_pλ∕λ_p. The rate of change of particle temperature is given as:

(3.19)

where, Nu_p is the particle Nusselt number given as

(3.20)

τ_p^T is the particle thermal response time defined as

(3.21)

T_f(x_p) is the fluid temperature at the particle position, C_p is the particle heat capacity, and λ is the fluid thermal conductivity.

3.2 Particle Dispersion

This section describes the capabilities of TransAT for particle dispersion modelling using Langevin equation. The particle dispersion model is used in combination with the Reynolds-Averaged turbulence modelling approach.

3.2.1 Langevin model

The particle fluctuation variance (σ_p) is defined as a function of the fluid fluctuation variance (σ_f) and the particle response time as:

(3.22)

where τ_f = MIN(τ_e,τ_r) is the effective eddy time scale, τ_p = ρ_pd_p²∕18μ is the particle inertial response time, K = 1 + 0.15Re_p^2∕3 is the drag correction for finite Reynolds number, and σ_f = . The model constants C_τ and C_σ correspond to the input variables tauL_constant and sigmaF_constant, respectively and are assigned default values of 0.11 and 1.69 ≈ 0.5, respectively (Lakehal et al., 1995).

In τ_f, τ_e = C_τk∕ϵ is the local eddy life time in an Eulerian sense, and τ_r = L_e∕|u_p -u_f| is the time taken by the particle to travel through the eddy (residence time), where L_e = C_μ^3∕4k^3∕2∕ϵ is the characteristic length scale of an eddy.

The particle velocity fluctuation equation is written as follows:

(3.23)

where τ′_pi is the particle response time to turbulent fluctuations and ζ(t) is a random function with a Gaussian probability density distribution and zero mean. The turbulent particle response time is given as:

(3.24)

where u_ri is the terminal velocity due to action of body forces such as gravity.

3.3 Dense Particle Modelling

3.3.1 Eulerian-Lagrangian dense particle model

The Eulerian-Lagrangian formulation for dense particle flows featuring non-negligible volume fractions (α_P > 1%) in incompressible flow conditions is implemented in the CMFD code TransAT as follows (mass and momentum equations for the fluid phase and Lagrangian particle equation of motion):

(3.25)

(3.26)

(3.27)

(3.28)

(3.29)

where α is the in-cell volume of fluid (α + α_P = 1), u is the velocity of the carrier phase, u_P is the particle velocity, u[xp] is the fluid velocity interpolated on to the particle position, Re_P is the particle Reynolds number, ΔU_f-P is the relative velocity between the particle and fluid at the particle position, Π is the sum of viscous stress σ and pressure p, τ is the turbulent stress tensor (depending whether RANS, V-LES or LES is employed to deal with turbulence). Sources terms in (3.26) denote body forces, F_b,and the rate of momentum exchange per volume between the fluid and particle phases, F_P . Various drag correlations are available in the literature to account for higher volume fractions of particles; also the fluid viscosity can be modified based on the volume fraction of particles. In (3.27), F_coll denotes the inter-particle stress. The momentum equation (3.26) presented here does not neglect viscous and turbulent diffusion mechanisms in the fluid phase. The inter-phase drag model in (3.27) is set according to Gydaspow (1989)

3.3.2 Four-way coupling modelling

The fluid-independent collision force F_coll is made dependent on the gradient of the inter-particle stress, π, using

(3.30)

Collisions between particles are estimated by the isotropic part of the inter-particle stress (its off-diagonal elements are neglected.). π is modelled as a continuum stress (Harris & Crighton, 1994), viz.

(3.31)

In TransAT, the particle field is predicted in a Lagrangian way first, then high-order accurate interpolations are resorted to map the Eulerian field (to estimate π), then back again to the Lagrangian system to determine F_coll). The constant P_s has units of pressure, and α_c is the particle volume fraction at close packing (typically specified to be 0.6), and constant β is set according to F.M. Auzerais & Russel (1988). The ε is a small number in the order of 10^-7.

Collision/Contact Pressure Models

Snider Model

The collisional pressure model of Snider (2001) is defined as,

(3.32)

where P_p, β are model constants and ε = 10^-7 is a small number that prevents the sign of the collision pressure from changing for α_p > α_p,cp. Typical values for the model constants are given in Table 3.1.

Jackson Model

The frictional collisional pressure of Johnson & Jackson (1987) is defined as,

(3.33)

where typical values for the constants are given in Table 3.2. The Jackson model is envisaged as an add-on to kinetic models for granular flows and therefore has a cut-off particle volume fraction below which it is set to zero. The Jackson model becomes qualitatively similar to the Snider model (Snider, 2001), if α_p,min = 0.

In order to apply the Jackson model in Lagrangian particle tracking, the change in sign of the denominator for α_p > α_p,cp must be disallowed. Therefore, the Jackson model is implemented in TransAT in a way similar to the Snider model.

(3.34)

3.4 Particle-boundary collision

The modelling of particle-boundary collision in TransAT is based on the hard sphere model described in Crowe et al. (1998). This model is applied to interactions of particles with solid wall boundaries, embedded objects and symmetry planes.

3.4.1 Hard sphere model

The hard sphere model for particle-wall collision is based on the integrated form of the momentum equations of classical dynamics, where the deformation of the particle does not appear explicitly in the formulation.
The translational velocity _p and angular velocity _p after the collision is calculated from the corresponding quantities before collision, the coefficient of restitution e and kinetic (sliding) friction f.
The deformation of the particle upon collision is assumed to be negligible. The particles are allowed to slide along the wall, where Coulomb’s friction law is applied during this process. It is assumed that once a particle has stopped sliding there is no further sliding.

The coefficient of restitution e is regarded as a property of the material and is the ratio of the normal component of impulse in the compression period to that in the recovery period:

(3.35)

Coulomb’s law states that the friction force is the product of the normal force and the coefficient of kinetic friction f. The friction force results in a change of the tangential component of impulse during the sliding period according to

(3.36)

with corresponding to the unit vector in the direction of tangential velocity.

Note that the particle collision model is formulated for stationary boundaries. If the boundaries are moving, the particle tangential velocities are transformed into a new, moving frame of reference, which has the same velocity as the moving boundary, before the model equations are applied. Afterwards they are transformed back to the original frame of reference.

3.4.2 Collision with wall boundary

Depending on whether the particle stops sliding during the collision process or continues throughout, the translational and angular velocity of the particle after collision are calculated as follows (assuming a plane wall with a normal vector in the positive y-direction):
If v_y⁰ |v_p| < - 2_____ 7f(e+1), then

If - 2_____ 7f(e+1) < v_y⁰ |v_p| < 0, then

The superscript 0 denotes the initial state before collision and ε_x,z corresponds to the components of the tangential velocity unit vector in the respective direction, i.e. they satisfy the following equation.

(3.49)

3.4.3 Collision with symmetry boundary

When particles collide with a symmetry boundary the collision model coefficients of restitution and friction take the following values

In this case the denominator of the condition for distinction becomes zero and the following relation holds for all non-zero particle velocities.

(3.51)

The corresponding equations for the translational and rotational velocities after collision simplify to

i.e. only the sign of the wall normal velocity component is changed, the remaining translational velocity components and the rotational velocities do not change, thus the total translational and rotational momentum is conserved.

Chapter 4
Hydrates Modelling

Hydrates modelling is available in TransAT in conjunction with the N-phase mixture model (See Chapter Phase-Average Mixture Models).

4.1 Hydrate Formation Kinetics

Two models have been implemented for hydrate formation: finite–rate kinetics model (Vysniauskas & Bishnoi, 1983) and instantaneous formation models as formation kinetics can be considered to be negligible compared to the other time scales of the flow.

4.1.1 Finite Rate Kinetics

A finite-rate kinetics model for the formation of hydrates (based on interfacial super-saturation) was implemented in TransAT. It is called the Bishnoi model. This model (Vysniauskas & Bishnoi, 1983) is as follows:

(4.1)

with the default values used for the constants listed in Table 4.1.

---

Constant	Value
A	4.55 × 10-26
ΔEa	106204 J∕gmol
γ	2.986
a	0.078Kb
b	2.411

Table 4.1: Hydrates kinetics model (Vysniauskas & Bishnoi, 1983)

---

4.1.2 Instantaneous Kinetics

The instantaneous formation model has been implemented in TransAT as follows.

• Similar to the finite-rate model, the temperature in each cell is compared to the T_eq for formation of hydrates. Several curves are available in TransAT for the following materials:
  • Methane
  • Ethane
  • Propane
  • Isobutane

  If the temperature in the computational cell is higher than T_eq then no hydrates are formed.

• Note that if the hydrate formation rate is high, then locally more heat will be released resulting in a higher temperature. If the temperature goes higher than T_eq then hydrate formation is inhibited. The iterations per time step in the implicit method allow one to obtain the converged solution of the correct temperature and reaction rate.
• The number of moles of hydrate forming phase (hf) in a computational cell is given as:

  (4.2)

  where, Mw_hf is the molecular weight of hydrate forming phase, M_hf is the number of moles of the hydrate forming phase in the cell, and Ω_i is the volume of the computational cell.

• The number of moles of water in the computational cell is given as:

  (4.3)

  where Mw_w is the molecular weight of water, and M_w is the number of moles of water available for hydrate formation.

• The rate of formation of hydrates ṙ is then given as:

  (4.4)

  where n_hy is the hydration number of the hydrate forming phase. The reaction rate is calculated such that all the hydrate forming phase present in a cell is consumed within one time step (provided that the temperature is less than T_eq).

• The hydrate formation rate ṙ is updated with a relaxation factor β_hy in every iteration for a given time step as follows.

  (4.5)

  where (k + 1) denotes the current iteration, and (k) denotes the previous iteration. Relaxation is a measure to prevent oscillations in the reaction rate which in turn prevents the convergence of the volume fraction equations. A relaxation factor of 0.4 is used in TransAT.

4.2 Hydrate Dissociation and Melting

4.2.1 Instantaneous Hydrates Dissociation

An instantaneous hydrates dissociation model has been implemented in TransAT.

• The temperature in each cell is compared to the T_eq for formation of hydrates. The same curve as for hydrate formation is used. If the temperature in the computational cell is lower than T_eq by ΔT_sh then hydrates can form, where ΔT_sh is the equilibrium range allowed.
• If the temperature in the computational cell is higher than T_eq then hydrates dissociate.
• Note that if the hydrate dissociation rate is high, then locally more heat will be absorbed resulting in a lower temperature. If the temperature goes lower than T_eq then hydrate dissociation is inhibited. The iterations per time step in the implicit method, allow one to obtain the converged solution of the correct temperature and reaction rate.
• The number of moles of hydrates (M_hy) in a computational cell is given as:

  (4.6)

  where, Mw_hy is the molecular weight of the hydrate phase, M_hy is the number of moles of the hydrate forming phase in the cell, and Ω_i is the volume of the computational cell. The rate of dissociation of hydrates ṙ is then given as:

  (4.7)

  The reaction rate is calculated such that all the hydrate phase present in a cell is consumed within one time step, provided that the temperature is higher than T_eq + ΔT_sh.

• The hydrate dissociation rate ṙ is updated with a relaxation factor β_hy in every iteration for a given time step as follows.

  (4.8)

  where (k + 1) denotes the current iteration, and (k) denotes the previous iteration. Relaxation is a measure to prevent oscillations in the reaction rate A relaxation factor of 0.4 is used in TransAT.

4.2.2 Heat Absorption

The formation of hydrates is exothermic. The dissociation of hydrates into water and gas phases is endothermic and results in absorption of heat locally. The enthalpy of dissociation is taken to be the same as the enthalpy of formation.

The heat absorbed is included in the mixture temperature equation as a sink term given by,

(4.9)

where H_hy,diss is the enthalpy of hydrate dissociation, and ṙ_diss is the rate of dissociation of hydrates. With the dissociation and the related heat absorption, the hydrate volume fraction at any location is a non-linear dynamic balance between possible formation and dissociation (with a gap of T_sh - T_sc between the formation and dissociation limits for numerical stability).

4.3 Hydrate Rheology

A special hydrate rheology model has been implemented in TransAT.

First, the continuous phase is defined depending on the Water Cut. If it is more than 50%, water is the continuous phase. Otherwise, oil is the continuous phase.

The volume fraction of the dispersed phase is then calculated as follows

(4.10)

where α_{water in hydrates} = Hα_hydrates, with H the hydration number. The viscosity of the slurry is then calculated using Mills law (Mills, 1985):

(4.11)

with

(4.12)

where τ denotes the instantaneous shear stress. For numerical reasons, in the implementation τ₀ cannot be smaller than 10^-6 Pa and α_eff must not be larger than 0.997α_max.

Bibliography