Particle Induced Turbulence Source

$k - ε$ Turbulence

Troshko Hassan

The Troshko-Hassan particle-induced turbulence source specification method begins by assuming that the whole of the work that is done by the drag force is an unconditionally positive production of continuous-phase pseudo-turbulence. Another feature of this model is that it assumes that this strong turbulence source is dissipated locally using a particle relaxation time scale.

The Troshko and Hassan model describes bubble induced turbulence effects [559] by providing source terms to the continuous K-Epsilon model. It uses the Virtual Mass phase interaction model. If the Virtual Mass phase interaction model is not selected, then the virtual mass coefficient defaults to the value of a spherical bubble:

Figure 1. EQUATION_DISPLAY

C_{V M} = \frac{1}{2}

(2471)

The source term is derived by considering that the correlation between velocity fluctuation at the interface and interfacial force density is exactly equal to the work of interfacial force density per unit time. For bubbly flows, the drag force dominates the interfacial force. The original form from Troshko and Hassan [559] for defining the is:

Figure 2. EQUATION_DISPLAY

S_{k_{c}} = v_{r} \cdot F_{i j}^{D} = \frac{3}{4} \frac{C_{D}}{l_{cd}} α_{d} ρ_{c} {| v_{r} |}^{3}

(2472)

where:

$C_{D}$ is the bubble drag coefficient and $v_{r}$ is the mean slip velocity.

In Simcenter STAR-CCM+, Eqn. (2469) is rewritten in terms of the linearized drag coefficient, $A_{D}$ , to account for multi-particle effects:

Figure 3. EQUATION_DISPLAY

S_{k_{c}} = A_{D} {| v_{r} |}^{2}

(2473)

The dissipation term uses the energy source term that decays with a characteristic time, the Bubble Pseudo-Turbulence Dissipation Relaxation (BPTDR) time, $τ_{b}$ , and is scaled with a calibration constant $C_{3}$ :

Figure 4. EQUATION_DISPLAY

S_{ε_{c}} = \frac{C_{3} S_{k_{c}}}{τ_{b}}

(2474)

For the data that was analyzed, Troshko and Hassan found $C_{3} = 0.45$ . The bubble pseudo-turbulence dissipation relaxation time is:

Figure 5. EQUATION_DISPLAY

τ_{b} = \frac{2 C_{V M} d}{3 C_{D} | v_{r} |}

(2475)

Assuming that $ρ_{d} « ρ_{c}$ for bubbly flow and again substituting $C_{D}$ with $A_{D}$ gives:

Figure 6. EQUATION_DISPLAY

τ_{b} = \frac{ρ_{c} C_{V M} α_{d}}{2 A_{D}}

(2476 2491)

where $C_{V M}$ is the virtual mass coefficient.

If this model is used with the Reynolds Stress equations, its contribution is assumed to be isotropic:

Figure 7. EQUATION_DISPLAY

S_{R c} = \frac{2}{3} S_{k c} I

(2477)

Gosman

Reynolds averaging of continuous-phase velocity fluctuation times the drag source term from the fluctuating part of the continuous-phase momentum equation results in the following source term for the continuous-phase turbulent kinetic energy:

Figure 8. EQUATION_DISPLAY

S_{k c} = {A_{c d}}^{D} (q_{c d} - {2 q}_{c}^{2} - v_{T D} \cdot v_{r})

(2478)

where:

$q_{c d}$ is the dispersed-phase average of the product of dispersed and continuous phase velocity fluctuations.
$q_{c}$ is half the dispersed-phase average of the product of continuous and continuous phase velocity fluctuations.

Also keep in mind that:

Only the contribution of the drag force has been considered in this derivation.
The closure that is used for the drag force is good under Stokesian conditions but only approximate if ${A_{c d}}^{D}$ represents a non-Stokesian drag law.

Gosman et al. [468] describe one model for this drag-related source term, which is essentially:

Figure 9. EQUATION_DISPLAY

S_{k c} = A_{c d} (C_{t} - 1) {2 K}_{c} - {F_{c d}}^{T D} \cdot v_{r}

(2479)

and can be arrived at by approximating:

Figure 10. EQUATION_DISPLAY

\frac{q_{c d}}{{2 q}_{c}^{2}} \approx \frac{| v_{d}' |}{| v_{c}' |} = C_{t}

(2480)

so the $(C_{t} - 1)$ contribution represents an exchange between cross-phase covariance and continuous-phase turbulent kinetic energy.

Although the final term has the appearance of a work term, in fact it has the opposite of the expected sign. If the relative motion of the phases is in the same direction as the turbulent dispersion force between the phases, the final term can become a sink of continuous phase turbulent kinetic energy.

The corresponding contribution to the Reynolds Stress equations is assumed to be isotropic:

Figure 11. EQUATION_DISPLAY

S_{R c} = \frac{2}{3} S_{k c} I

(2481)

The dominant dissipation time scale is assumed to be $k_{c} / ε_{c}$ as for shear-induced turbulence, but only the $(C_{t} - 1)$ exchange term is allowed to contribute.

Figure 12. EQUATION_DISPLAY

S_{ε c} = {A_{c d}}^{D} (C_{t} - 1) {2 ε}_{c}

(2482)

Tchen

Reynolds averaging of continuous-phase velocity fluctuation times the drag source term from the fluctuating part of the continuous-phase momentum equation results in the following source term for the continuous-phase turbulent kinetic energy:

Figure 13. EQUATION_DISPLAY

S_{k c} = C_{0} A_{D} (q_{c d} - {2 q}_{c}^{2} - v_{T D} \cdot v_{r})

(2483)

where:

$q_{c d}$ is the dispersed-phase average of the product of dispersed and continuous phase velocity fluctuations.
$q_{c}$ is half the dispersed-phase average of the product of continuous and continuous phase velocity fluctuations.
Constant $C_{0}$ is not part of the model derivation, but defaults to 1. You can adjust it to test the sensitivity of the solution to this term.

See Particle Induced Turbulence Models for the origin of this term, but also keep in mind that:

Only the contribution of the drag force has been considered in this derivation.
The closure that is used for the drag force is good under Stokesian conditions but only approximate if ${A_{c d}}^{D}$ represents a non-Stokesian drag law.

Taken together, $q_{c d} - {2 q}_{c}^{2}$ represents an exchange between cross-phase covariance and continuous-phase turbulent kinetic energy.

${A_{c d}}^{D} (v_{T D} \cdot v_{r})$ represents the work that is done by the turbulent dispersion force against mean slip velocity. See Particle Induced Turbulence Models.

The corresponding contribution to the Reynolds Stress equations is assumed to be isotropic:

Figure 14. EQUATION_DISPLAY

S_{R c} = \frac{2}{3} S_{k c} I

(2484)

Following Elghobashi and Abou-Arab [455], the dominant dissipation time scale is assumed to be the same as for shear-induced turbulence:

Figure 15. EQUATION_DISPLAY

S_{ε c} = C_{ε 3} \frac{ε_{c}}{k_{c}} max (S_{k c}, 0)

(2485)

A negative turbulent energy source term is not allowed to contribute to the turbulent dissipation equation. Parameter $C_{ε 3}$ has a default value of 1.44, and you can adjust it independently of constants for other terms.

The correlations $q_{c d}$ and ${q_{c}}^{2}$ are closed using Tchen theory as described in Tchen Closures:

Figure 16. EQUATION_DISPLAY

\begin{matrix} q_{c d} = 2 [\frac{b + η}{1 + η}] k_{c} \\ {q_{c}}^{2} = k_{c} \end{matrix}

(2486)

where:

$η$ is the particle-eddy interaction time, which is scaled by particle relaxation time.
$b$ is the ratio of the coefficients of the continuous/disperse phase acceleration terms in the particle equation of motion:

Figure 17. EQUATION_DISPLAY

b = \frac{1 + C_{V M}}{\frac{ρ_{d}}{ρ_{c}} + C_{V M}}

(2487)

After closure, the transfer term $q_{c d} - {2 q}_{c}^{2}$ becomes proportional to $(b - 1)$ .

So for gas-particle flow, $ρ_{d} » ρ_{c}$ , where coefficient $b$ is small, the transfer term becomes a sink of continuous phase turbulence kinetic energy.

On the other hand, for bubbly flow, $ρ_{d} « ρ_{c}$ , where coefficient $b$ tends to 3/2 if $C_{V M}$ is about 1/2, the transfer term becomes a source of continuous phase turbulence kinetic energy.

Generic

This approach combines the other particle-induced turbulence source methods into a generic one and provides the flexibility to modify key source term parameters that cannot be adjusted in the other methods. This adaptability makes the generic approach more appropriate for considering the multiscale effects associated with multiphase turbulence.

The drag source term from the fluctuating part of the continuous-phase momentum equation results in the following source term for the continuous-phase turbulent kinetic energy:

Figure 18. EQUATION_DISPLAY

S_{k c} = C_{g} \frac{ρ_{d}}{τ_{D}} \overline{α_{d}} {| v_{r} |}^{2}

(2488)

In Simcenter STAR-CCM+, Eqn. (2488) is rewritten in terms of the linearized drag coefficient, $A_{D}$ given as:

Figure 19. EQUATION_DISPLAY

S_{k c} = A_{D} {(C_{g} | v_{r} |}^{2} - C_{g}^{'} \min [0, 〈 v_{r} 〉 \cdot v_{T D}])

(2489)

where:

$C_{g}$ is drag force modulation parameter which regulates the interfacial turbulence transfer assuming only a fraction of drag contributes to large-scale turbulence, with the rest dissipating into small-scales.
$C_{g}^{'}$ is a modulation factor for scaling the contribution of the turbulence dispersed force to particle-induced turbulence, by default it is 1.

For the turbulence dissipation term, the dominant dissipation time scale is assumed to be

τ_{ε}

and the dissipation term is modelled as:

Figure 20. EQUATION_DISPLAY

S_{ε_{c}} = \frac{1}{τ_{ε}} S_{k c}

(2490)

The formulation of the Generic method is similar to Tchen, Gosman, and Troshko & Hasan PIT methods, with the flexibility on adjusting the methods and parameters for solving the dissipation time scale $τ_{ε}$ .

For the particle relaxation method, the dissipation term uses the energy source term that decays with a characteristic time, the bubble pseudo-turbulence dissipation relaxation time, similar to the Troshko and Hassan method, given as:

Figure 21. EQUATION_DISPLAY

τ_{c} = \frac{ρ_{c} C_{V M} α_{d}}{2 A_{D}}

(2476 2491)

The vortex turnover method for solving the dissipation time scale is given as:

Figure 22. EQUATION_DISPLAY

τ_{c} = \frac{k_{c}}{{C_{3} ε}_{c}}

(2492)

where

C_{3}

is a calibration constant, with a default value of 1.44.

For the realizable vortex turnover method the dissipation time scale is computed as:

Figure 23. EQUATION_DISPLAY

τ_{c} = \frac{k_{c}}{{C_{3} ε}_{c}} (1 + \sqrt{\frac{{ν ε}_{c}}{k_{c}^{2}}})

(2493)

where $ν$ is the kinematic molecular viscosity.

If the Generic method is used with the Reynolds Stress equations, its contribution is assumed to be anisotropic, given by:

Figure 24. EQUATION_DISPLAY

S_{R c} = A_{D} {C_{g} [A_{t} (〈 v_{r} 〉 〈 v_{r} 〉 + {(〈 v_{r} 〉 〈 v_{r} 〉)}^{T}) + (1 - A_{t}) \frac{2}{3} 〈 v_{r} 〉 \cdot 〈 v_{r} 〉] - C_{g}^{'} [A_{t} (〈 v_{r} 〉 v_{T D} + {(〈 v_{r} 〉 v_{T D})}^{T}) - (1 - A_{t}) \frac{2}{3} 〈 v_{r} 〉 \cdot v_{T D} 〉]}

(2494)

where $A_{t}$ is the isotropy coeﬃcient. The default is $A_{t} = 0.5$ . The turbulent dissipation rate for the Reynolds Stress equations is the same as defined in Eqn. (2489).

Sato

The Sato model is the simplest and earliest model for particle-induced mixing, and is a robust alternative to later source-based models for particle induced turbulence. This model is applicable when modeling a dilute bubbly flow.

In a dilute bubbly flow, the shear-induced turbulence and bubble-induced fluctuations operate at different scales. Figure 11 in Lance and Bataille (1991) [497] gives an example of grid turbulence where bubble-induced fluctuations are dominant but the decay rate of the shear-induced part remains unaffected. This situation is not well-described by the more advanced models that solve transport equations for the combined turbulence energy with one length scale equation. However, the Sato model assumes the shear-induced part is totally independent while the bubble-induced part is always in local equilibrium and does not need a transport equation.

It acts through enhanced effective viscosity of the continuous phase:

Figure 25. EQUATION_DISPLAY

ν_{e f f} = ν + ν_{t u r b} + ν_{S a t o}

(2495)

where:

$ν$ is the continuous phase kinematic viscosity
$ν_{t u r b}$ is the turbulent diffusivity
$ν_{S a t o}$ is the Sato bubble-induced viscosity.

The Sato bubble-induced viscosity is given by:

Figure 26. EQUATION_DISPLAY

ν_{S a t o} = k f_{d} f_{B} α_{g a s} \frac{d_{b}}{2} U_{s l i p}

(2496)

where:

$k$ is a model calibration constant
$f_{d}$ is the van Driest damping factor
$f_{B}$ is a bubble diameter shape correction factor
$α_{g a s}$ is the dispersed phase volume fraction
$d_{b}$ is the bubble diameter
$U_{s l i p}$ is the slip velocity between phases.

The van Driest [562] damping factor is:

Figure 27. EQUATION_DISPLAY

f_{d} = {(1 - \exp (- \frac{y^{+}}{A^{+}}))}^{2}

(2497)

where $A^{+}$ is a dimensionless constant.

The bubble diameter shape correction factor is:

Figure 28. EQUATION_DISPLAY

f_{B} = \begin{matrix} 4 y / d_{b} (1 - y / d_{b}) & y < d_{b} / 2 \\ 1 & y > d_{b} / 2 \end{matrix}

(2498)

where $y$ is the distance from the wall.

This superposition of turbulent and bubble-induced viscosities is partly justified by the measurements of turbulent intensity at different volume fractions (seen in Figure 11 of Lance and Bataille [497]). The effect of volume fraction on the decay of grid-generated turbulence suggests that the effect of bubble-induced fluctuations is purely additive. However, this is not strictly true, particularly in a flow with significant shear production and wall effects. This model should be used with caution, particularly at low Reynolds numbers where the bubble-induced part may be much larger than the shear-induced part.

$k - ω$ Turbulence

The interaction of dispersed phase particles with turbulence in the continuous phase has traditionally been modeled in terms of the turbulence dissipation rate, $ϵ$ . This is because within the inertial subrange between the largest energy-containing eddies and the smallest viscous eddies, $ϵ$ describes the relationship between length scale and velocity scale, which allows the magnitude of velocity fluctuations around a particle of given size to be estimated.

For the $k - ω$ model, these $ϵ$ -based models are reused by computing an $ϵ$ field consistent with the $k - ω$ model. However, different turbulence models can have different interpretations of k and $ϵ$ required to successfully reproduce measured velocity profiles and shear stress, and so may compute different values of k and $ϵ$ , particularly near walls. Therefore it is to be expected that multiphase models that are known to physically depend on $ϵ$ , such as the rate of particle break-up, may need some recalibration when the continuous phase turbulence model is switched between $k - ϵ$ and $k - ω$ .

Alternative models for Particle Induced Turbulence differ in two aspects, which must be chosen and calibrated consistently. The first is the strength and derivation of the source term for turbulent kinetic energy, $S_{k}$ . The second is the choice of the dissipation time scale for the corresponding contribution to the turbulent dissipation equation, $S_{ϵ}$ . The manipulation of these two terms to derive the source $S_{ω}$ for the $ω$ equation is outlined below.

Using rate of change notation, $\dot{k}$ , as a shorthand for the continuous-phase balance equation for k, and so on, the contributions of particle-induced turbulence to the three equations are:

Figure 29. EQUATION_DISPLAY

\dot{k} = S_{k}

(2499)

Figure 30. EQUATION_DISPLAY

\dot{ϵ} = S_{ϵ}

(2500)

Figure 31. EQUATION_DISPLAY

\dot{ω} = S_{ω}

(2501)

Assuming $ϵ = β^{*} ω k$ gives:

Figure 32. EQUATION_DISPLAY

S_{ω} = ω {\frac{\dot{ϵ}}{ϵ} - \frac{\dot{k}}{k}} = \frac{1}{k} {\frac{S_{ϵ}}{β^{*}} - ω S_{k}}

(2502)

This result also leads to a useful insight into the behavior of the turbulent timescale, $\frac{k}{ϵ} = \frac{1}{β^{*} ω}$ , in alternative models whenever the particle-induced turbulence term becomes significant.

Since the Gosman model uses $\frac{k}{ϵ}$ as the dissipation time scale for particle induced turbulence (neglecting the turbulent dispersion dot product term which is small in fully developed flows), it gives a precisely zero contribution to the $ω$ equation:

Figure 33. EQUATION_DISPLAY

S_{ϵ} = \frac{ϵ}{k} S_{k}

(2503)

Figure 34. EQUATION_DISPLAY

S_{ω} = 0

(2504)

The Tchen model starts with the same dissipation timescale, $\frac{k}{ϵ}$ , as Gosman, but with a calibration constant $C_{ϵ 3}$ , so that $S_{ϵ}$ is not immediately in local equilibrium with $\frac{ϵ}{k} S_{k}$ .

Then there is an opportunity for particle-induced turbulence to shorten $\frac{k}{ϵ}$ when using the default calibration $C_{ϵ 3} = 1.44$ and the net source $S_{k}$ happens to be positive, as is likely in fully developed flow, and to lengthen it in other circumstances.

Figure 35. EQUATION_DISPLAY

S_{ϵ} = C_{ϵ 3} \frac{ϵ}{k} \max (S_{k}, 0)

(2505)

Figure 36. EQUATION_DISPLAY

S_{ω} = \frac{ω}{k} {C_{ϵ 3} \max (S_{k}, 0) - S_{k}}

(2506)

The Troshko-Hassan model uses a dissipation time scale $τ_{b}$ proportional to bubble diameter over slip velocity:

Figure 37. EQUATION_DISPLAY

S_{ϵ} = C_{ϵ 3} \frac{S_{k}}{τ_{b}}

(2507)

Figure 38. EQUATION_DISPLAY

S_{ω} = \frac{S_{k}}{k} {\frac{C_{ϵ 3}}{β^{*} τ_{b}} - ω}

(2508)

This model gives a clear bubble-related asymptote for the turbulent time scale wherever the bubble-induced term becomes dominant, for example, far from walls:

Figure 39. EQUATION_DISPLAY

\frac{k}{ϵ} \to \frac{τ_{b}}{C_{ϵ 3}}

(2509)

The production rate for the specific dissipation rate for the Generic method is:

Figure 40. EQUATION_DISPLAY

S_{ω} = C_{ϵ 3} \frac{13}{36} A_{D} (C_{g} {| 〈 v_{r} 〉 |}^{2} - C_{g}^{'} \min [0, 〈 v_{r} 〉 \cdot v_{T D}])

(2510)

Large Eddy Simulation

The Particle Induced Turbulence model contributes to the turbulent viscosity according to:

Figure 41. EQUATION_DISPLAY

μ_{t} = μ_{S G S} + μ_{P I T}

(2511)

Here $μ_{S G S} = ρ_{c} Δ^{2} S$ is the SGS contribution to the turbulent viscosity with $ρ_{c}$ being the density of the continuous phase, $Δ$ the length scale, and $S$ a function of the strain rate tensor. The latter two are defined individually for each of the SGS models as in Eqn. (1388), Eqn. (1392), and Eqn. (1398) for $Δ$ , and Eqn. (1387), Eqn. (1391), and Eqn. (1397) for $S$ .

On the other hand:

Figure 42. EQUATION_DISPLAY

μ_{P I T} = Δ ρ_{c} α_{d} | v_{r} |

(2512)

is defined as in Lakehal [495], where $α_{d}$ is the dispersed phase volume fraction and $v_{r}$ is the phase relative velocity.

Particle Induced Turbulence Source

k−ε Turbulence

k−ω Turbulence

Large Eddy Simulation

$k - ε$ Turbulence

$k - ω$ Turbulence