S-Gamma Phase Interactions

The S-Gamma equations (for S 0 and S 2 ) are solved with source terms added to model breakup, coalescence of bubbles and other particle types, and their entrainment in a presence of large scale interface.

Both these processes conserve the total volume of the dispersed phase. No extra term is required in the volume fraction equation which is used together with S 0 and S 2 to give the mean diameter.

Pre-Integrated S-Gamma

This model evaluates the integrals using an analytical method.

Breakup

The source terms for breakup are found by integrating the effect of a bubble or droplet of size d breaking up over all diameters. The general form for the breakup is from Lo and Zhang [507]:

Figure 1. EQUATION_DISPLAY
d d t S γ = d l o d h i n P ( d ) ( N ( d ) 1 - γ / 3 - 1 ) τ ( d ) d γ d ( d )
(2198)

Where:

  • N is the number of fragments, or child bubbles, produced when a bubble breaks down.
  • τ is the timescale for breakup.
  • d l o and d h i are the upper and lower limits of a particular breakup regime.

The implementation in Simcenter STAR-CCM+ considers both viscous and inertial breakup regimes. For laminar flow, the lower limit for viscous breakup is the critical diameter d c r and the upper limit is infinite. For turbulent flow, the upper limit is the Kolmogorov length scale L k ; if d c r is greater than L k then viscous breakup plays no role. For inertial breakup the integral is performed from the greater of L k and d c r to infinity. By writing N ( d ) and τ ( d ) in the form:

Figure 2. EQUATION_DISPLAY
N ( d ) = N' d b
(2199)
Figure 3. EQUATION_DISPLAY
τ ( d ) = τ′ d a
(2200)

The integral can be written in the form:

Figure 4. EQUATION_DISPLAY
d d t S γ = N' 1 - γ / 3 τ′ d l o d h i n P ( d ) d δ 1 d ( d ) - 1 τ′ d l o d h i n P ( d ) d δ 2 d ( d )
(2201)

Where δ 1 = γ + b - b γ / 3 - a and δ 2 = γ - a . This form of the equation is from Hill [474].

Viscous Breakup

For viscous breakup, the critical diameter is given as [507], [474]:

Figure 5. EQUATION_DISPLAY
d c r = 2 σ Ω c r μ c γ ˙
(2202)

Where:

  • μ c is the viscosity of the continuous phase.
  • σ is the surface tension.
  • γ ˙ is the local shear rate.
  • Ω c r is the critical capillary number [474].

The fragmentation number is assumed to be binary [507], that is:

Figure 6. EQUATION_DISPLAY
N ( d ) = 2
(2203)

The breakup timescale is from [507]:

Figure 7. EQUATION_DISPLAY
τ ( d ) = μ c d σ f ( λ )
(2204)

Where:

  • μ c is the continuous phase dynamic viscosity.
  • σ is the surface tension.
  • log f ( λ ) = p 0 + p 1 log ( λ ) + p 2 ( log ( λ ) ) 2 and λ = μ d / μ c and the constants are given in [474].
Inertial Breakup

For inertial breakup, the critical bubble diameter and related values are calculated as follows [507], [474]:

Figure 8. EQUATION_DISPLAY
d c r = ( 1 + C α α d ) ( 2 σ W e c r ρ c ) 3 / 5 ε - 2 / 5
(2205)
Figure 9. EQUATION_DISPLAY
τ ( d ) = 2 π k b r ( 3 ρ d + 2 ρ c 192 σ ) 1 / 2 d 3 / 2
(2206)
Figure 10. EQUATION_DISPLAY
N ( d ) = 2
(2207)

Where:

  • ρ d and ρ c are the dispersed and continuous phase densities.
  • ε is the continuous phase turbulent dissipation rate.
  • W e c r is the critical Weber number which you define with a default value of 0.25.
  • C α and k b r are constants with default values of 4.6 and 0.2, respectively.
Kocamustafaogullari
The Kocamustafaogullari model, is based on the stochastic secondary droplet (SSD) model. Kocamustafaogullari [488] correlated the experimental data for droplet size distribution in annular flows by using Levich [501] formulations of the equation of motion for a droplet in a turbulent flow. Eqn. (2255) is averaged over all possible realizations of the slip velocity Eqn. (2261), giving:
Figure 11. EQUATION_DISPLAY
B r e a k u p R a t e = B 1 { ( 8 3 π ρ g ρ l w t 2 D p 2 ) 1 / 2 ( 1 + β W e c r ) exp ( β W e c r ) , U = 0 ρ g ρ l w D p , U W e c r
(2208)

where:

  • W e c r is the critical Weber number. The default value is 12.
  • B 1 is a parameter regulating the speed of breakup. Large values signify that breakup occurs later. The default value is 2 3 .
  • ρ l is the density of the liquid.
  • ρ g is the density of the gas.
  • D p is the droplet diameter
  • w t is the turbulent slip velocity of the droplet through the continuous phase.
This method is the same as the corresponding method in the AMUSIG model. For more details on the formulation of the method, see AMUSIG Breakup and Coalescence.
Effect of the Liquid Viscosity

The effect of the liquid viscosity on the breakup rate is characterized by the Ohnesorge number [470]:

Figure 12. EQUATION_DISPLAY
Oh = μ l ρ l D p σ
(2209)

where μ l is the liquid viscosity and σ is the surface tension of the liquid.

To account for the enhanced stability of more viscous droplets, a correlation for the critical Weber and Ohnesorge numbers is defined as:

Figure 13. EQUATION_DISPLAY
W e c r = W e c r | O h = 0 ( 1 + a O h b )
(2210)

where a and b are constants, derived by Gelfand in [464].

Coalescence

Performing the double integral required for an exact evaluation of the coalescence term is not practicable. A simplified source term is developed based on an equivalent diameter d e q for each moment, which here is taken to be equal to d 3 γ .

The source term that is given in [507] is:

Figure 14. EQUATION_DISPLAY
d d t S γ = F c l ( 2 γ / 3 - 2 ) ( 6 α π ) 2 k c o l l u r e l P c l ( d e q ) d e q γ - 4
(2211)

Where:

  • F c l is the calibration coefficient.
  • k c o l l is a collision rate.
  • u r e l is the typical phase velocity difference.
  • P c l ( d e q ) is the probability of collision leading to coalescence.

To split the source term into a collision rate (in collisions per unit volume per unit time) and a coalescence probability, either of which can be replaced with functions you define. This formulation has been rearranged slightly to give:

Figure 15. EQUATION_DISPLAY
t S γ = F c l ( 2 γ / 3 - 2 ) K c o l l P c l ( d e q ) d e q γ
(2212)

Where K c o l l is the collision rate and is:

Figure 16. EQUATION_DISPLAY
K c o l l = k c o l l u r e l d e q 2 ( 6 α π d e q 3 ) 2
(2213)

For a turbulent flow, viscous and inertial effects lead to coalescence and a source term of the form that is shown above is evaluated for each. By successively setting γ to 0 and 2, the source terms for the zeroth and second moment S-Gamma transport equations are obtained.

Viscous Coalescence

For viscous coalescence, the collision rate is [507]:

Figure 17. EQUATION_DISPLAY
K c o l l , v = ( 8 π 3 ) 1 / 2 ( γ ˙ d e q ) d e q 2 ( 6 α π d e q 3 ) 2
(2214)

With γ ˙ being the shear rate. For a laminar flow, the shear rate is calculated using the velocity field. For a turbulent field, the shear rate is γ ˙ = ε ρ c / μ c [505] which uses the turbulent dissipation rate, density, and viscosity from the continuous phase. The probability is from [507]:

Figure 18. EQUATION_DISPLAY
P v ( d e q ) = exp ( - t d γ ˙ )
(2215)

Physically, this is comparing the time that the bubbles interact with the time required for the film initially to separate and then drain away. Four models for drainage modes are from [505]. Starting from the model that gives the slowest drainage time (and hence the lowest coalescence rate) toward the fastest, these models are called:

  • Fully Immobile:
Figure 19. EQUATION_DISPLAY
t d = 3 μ c d e q 2 F i 64 π σ 2 h c r 2
(2216)
  • Partially Mobile Interface (Short Collision time):
Figure 20. EQUATION_DISPLAY
t d = 3 2 ( F i d e q 2 μ d ρ d 32 π σ 2 h c r )
(2217)
  • Partially Mobile Interface (Quasi-Steady state):
Figure 21. EQUATION_DISPLAY
t d = π μ d F i 2 h c r ( d e q 4 π σ ) 3 / 2
(2218)
  • Fully Mobile:
Figure 22. EQUATION_DISPLAY
t d = 3 μ c d e q 4 σ ln ( d e q 8 h c r )
(2219)

Where F i = ( 3 / 2 ) π μ c γ ˙ d e q 2 .

The critical film thickness is:

Figure 23. EQUATION_DISPLAY
h c r = ( A H d e q 24 π σ ) 1 / 3
(2220)

Where:

  • A H is the Hamaker constant and its default value is 5 × 10 - 21 .
  • σ is the surface tension.

Reference [507] provides all equations except the fully immobile expression.

Inertial Coalescence

For inertial coalescence, the collision rate and probability are from [507]:

Figure 24. EQUATION_DISPLAY
K c o l l , i = ( 2 π 15 ) 1 / 2 ( ε d e q ) 1 / 3 d e q 2 ( 6 α π d e q 3 ) 2
(2221)
Figure 25. EQUATION_DISPLAY
P ( d e q ) = Φ m a x π ( 1 - k c l , 2 2 ( W e - W e 0 ) 2 Φ m a x 2 ) 1 / 2
(2222)

Φ m a x , the maximum phase difference is:

Figure 26. EQUATION_DISPLAY
Φ m a x = ( 2 h 0 2 ρ c σ ) / ( W e 0 μ d 2 d e q )
(2223)

Also:

Figure 27. EQUATION_DISPLAY
W e 0 = 0.8 W e c r
(2224)
Figure 28. EQUATION_DISPLAY
h 0 = 8.3 h c r
(2225)

You can set W e c r and k c l , 2 , the latter having a default value of 12.7.

Discrete Quadrature S-Gamma

This model evaluates the integrals associated with the breakup and coalescence using an adaptive discrete quadrature method. It also models the creation of new dispersed phase due to entrainment of one phase into another phase in a multiple flow regime scenario.

Coulaloglou and Tsouris Coalescence Efficiency

The Coulaloglou and Tsouris coalescence efficiency is given by (Tsouris and Tavlarides [560]):

Figure 29. EQUATION_DISPLAY
λ = exp ( C ρ c μ c ϵ c σ 2 ( 1 + α d ) 3 d e f f 4 )
(2226)

where:

  • C [ m 2 ] is a dimensional calibration constant

  • ρ c is the density of the continuous phase

  • μ c is the dynamic viscosity of the continuous phase

  • ϵ c is the turbulence dissipation rate of the continuous phase

  • σ is the surface tension

  • α d is the volume fraction of the dispersed phase.

O'Rourke Coalescence Efficiency

The O'Rourke coalescence efficiency model ([524]) detects collisions of the liquid droplets in a continuous gas phase. The colliding liquid droplets undergo a chaotic motion due to the turbulent fluctuations of the continuous phase. This method is the same as the corresponding method in the AMUSIG model. For more details on the formulation of the method, see O'Rourke Coalescence Efficiency.

Entrainment

Entrainment of a phase into another phase is observed physically in both natural processes and engineering applications. The S-Gamma entrainment implementation in Simcenter STAR-CCM+ allows specification of particles (bubbles/droplets) creation rate and diameter. Two methods for predicting gas entrainment rate are available. With these methods, the entrainment process is accounted for by identifying the high turbulent dissipation rate cells at the large scale interface and when the turbulent dissipation reaches a specified critical value, bubble creation is permitted.

The large scale interface cells are identified using the Large Interface Marker Band, and cells with entrainment probability are those where D k d t > 0 . For these cells, the turbulent dissipation rate ϵ is compared against the critical turbulent dissipation rate number, and entrainment is permitted when ϵ ϵ c r .

The most widely used critical dissipation rate method is the Yu et al. ([576]) methods, which gives:

Figure 30. EQUATION_DISPLAY
ϵ c r = C ϵ g 5 4 ( σ ρ d ) 1 4
(2227)
and:
Figure 31. EQUATION_DISPLAY
C ϵ = 2 3 2 W e l 1 4 ( F r u 2 ) 5 4
(2228)

where:

  • g is gravity.
  • W e l is the lower bound of the Weber number.
  • F r u is the upper bound of the Froude number.
Yu (Scale Separation) Entrainment Rate

The Yu (Scale Separation) entrainment rate model ([576]) detects creation and entrainment of bubbles in a dispersed phase. This method assumes that for the critical dissipation rate, upper r u and lower r l bounds of bubble radius are equal r u ( ϵ c r ) = r l ( ϵ c r ) , where:

Figure 32. EQUATION_DISPLAY
r u = 4 ( F r u 2 g ) 3 ϵ 2
(2229)
Figure 33. EQUATION_DISPLAY
r l = 2 8 5 ( σ W e l ρ ) 3 5 ϵ 2 5
(2230)

Then the critical bubble radius r c is assumed as the radius across which the gravity effects dominate the entrainment and vice-versa. When gravity and surface tension force are in balance, the corresponding bubble radius is given as:

Figure 34. EQUATION_DISPLAY
r o = r c = 0.5 ( σ ρ g ) 1 2
(2231)

While r u and r l will vary from cell to cell , r o does not, and the relation r l r o r u will hold true in all candidate cells.

The size of the bubble implemented in Simcenter STAR-CCM+ uses a limited adaptation of the Hinze scale ([576]), which assumes that the slope changed (gravity effects) dominate the entrainment. The Hinze estimated bubble diameter is given as:

Figure 35. EQUATION_DISPLAY
d e = 2 r
(2232)
for:
Figure 36. EQUATION_DISPLAY
r = r l r H r u
(2233)

Yu et al. suggest that to entrain a bubble, the potential as well as surface energy of the interface need to be overcome, where bubble creation rate per unit volume (/m3/s) is defined as:

Figure 37. EQUATION_DISPLAY
B ( d e ) = c b ρ α ϵ 2 3 d e 2 3 t t u r b ( π 3 ρ l d e 4 + π d e 2 σ )
(2234)
where:
  • t t u r b is the turbulent time scale.
  • c b is the entrainment calibration factor.

This equation ensures that number density varies at a different exponent when entrainment is dominated by gravity or surface tension based on bubble size.

Ma (Bubble Surface Energy) Entrainment Rate

The Ma (Bubble Surface Energy) entrainment rate model ([576]) detects creation and entrainment of bubbles in a dispersed phase. This method assumes that to entrain a bubble surface tension needs to be overcome and the bubble creation at the large interface is linearly proportional to ϵ . Ma et al. defines the bubble creation rate per unit volume (/m3/s) as:

Figure 38. EQUATION_DISPLAY
B ( d e ) = c b π ( ρ σ ) α l d e 2 ϵ
(2235)
For the S-Gamma entrainment model, S 0 and S 2 are solved with an added source term added for the rate of bubble creation.
Figure 39. EQUATION_DISPLAY
B S 0 e = B ( d e )
(2236)
Figure 40. EQUATION_DISPLAY
B S 2 e = ρ 2 3 g B ( d e ) d e 2
(2237)