Phase Averaging and Closures

Figure 1. EQUATION_DISPLAY

ϕ = \bar{ϕ} + ϕ' = {〈 ϕ 〉}_{k} + ϕ "

(2523)

where:

$\bar{ϕ}$ is the Reynolds average of the instantaneous variable $ϕ$ .
${〈 ϕ 〉}_{k}$ is the phase average of the instantaneous variable $ϕ$ formed by a Reynolds average that is weighted with the volume fraction for phase $k$ :

Figure 2. EQUATION_DISPLAY

{〈 ϕ 〉}_{k} = \frac{\overline{α_{k} ϕ}}{\overline{α_{k}}}

(2524)

where $α_{k}$ is the instantaneous volume fraction of phase $k$ .

Elsewhere in the User Guide, $α_{k}$ is used for the Reynolds average volume fraction, $\overline{α_{k}}$ .

The following Eulerian multiphase turbulence equations are derived in terms of phase-averaged turbulence stresses, kinetic energy, and dissipation [468]:

Figure 3. EQUATION_DISPLAY

R_{k} = {〈 {v_{k}}^{"} {v_{k}}^{"} 〉}_{k}

(2525)

Figure 4. EQUATION_DISPLAY

k_{k} = \frac{1}{2} {〈 {v_{k}}^{"} \cdot {v_{k}}^{"} 〉}_{k}

(2526)

Figure 5. EQUATION_DISPLAY

ε_{k} = \frac{μ_{k}}{ρ_{k}} {〈 2 {S_{k}}^{"} \cdot {S_{k}}^{"} 〉}_{k}

(2527)

Figure 6. EQUATION_DISPLAY

{S_{k}}^{"} = \frac{1}{2} {\nabla {v_{k}}^{"} + {(\nabla {v_{k}}^{"})}^{T}}

(2528)

where $v_{k}$ is the instantaneous velocity of phase $k$ .

Elsewhere in the User Guide, $v_{k}$ is used for the phase average velocity, ${〈 v_{k} 〉}_{k}$ .

Boussinesq Closure

Using the eddy-diffusivity or Boussinesq hypothesis, the Reynolds flux of any scalar $φ$ due to turbulent velocity fluctuation ${v'}_{c}$ in the continuous phase can be modeled as:

Figure 7. EQUATION_DISPLAY

\overline{ϕ' {v'}_{c}} = \overline{ϕ' v_{c}} = \overline{ϕ {v'}_{c}} \approx - \frac{{ν_{c}}^{t}}{σ_{ϕ}} \nabla \bar{ϕ}

(2529)

where $σ_{ϕ}$ is the ratio of turbulent diffusivities of momentum and $ϕ$ .

Tchen Closures

Thai-Van et al. [555] summarises closures for in-phase and cross-phase turbulent diffusivities that are based on the work of by Simonin and others (for example see Deutch and Simonin [449]) in extending the results of Tchen theory for small particles in steady homogeneous turbulence (for example see Hinze [475]).

These closures are used in Simcenter STAR-CCM+ for the Tchen options for turbulent dispersion force, for dispersed-phase turbulence, and for particle induced turbulence. They provide models for the following turbulent diffusivities:

Continuous-phase turbulent diffusion:

Figure 8. EQUATION_DISPLAY

{ν_{c}}^{t} = \frac{2}{3} {q_{c}}^{2} τ_{T}

(2530)

Interphase turbulent dispersion:

Figure 9. EQUATION_DISPLAY

D_{c d} = τ_{I} \frac{1}{3} q_{c d} I

(2531)

Dispersed phase turbulent diffusion:

Figure 10. EQUATION_DISPLAY

{ν_{d}}^{t} = τ_{I} \frac{1}{3} q_{c d} + \frac{1}{2} τ_{R} \frac{2}{3} {q_{d}}^{2}

(2532)

constructed using these correlations:

${q_{c}}^{2} = \frac{1}{2} {〈 {v_{c}}^{"} \cdot {v_{c}}^{"} 〉}_{c}$  
$q_{c d} = {〈 {v_{c}}^{"} \cdot {v_{d}}^{"} 〉}_{d}$  
${q_{d}}^{2} = \frac{1}{2} {〈 {v_{d}}^{"} \cdot {v_{d}}^{"} 〉}_{d}$

and these time scales:

$τ_{T}$ is the eddy turnover time.
$τ_{I}$ is the particle-eddy interaction time.
$τ_{R}$ is the particle relaxation time.

According to Thai-Van et al., Eqn. (2532) was originally derived for turbulent gas flows laden with heavy particles. They commented that an eddy-viscosity closure for a bubbly phase could also be derived but had a negligible effect on predictions.

If a turbulence model is used for the continuous phase, then ${q_{c}}^{2}$ is readily available as the continuous-phase turbulent kinetic energy $k_{c}$ , and ${ν_{c}}^{t}$ corresponds to the turbulent diffusivity provided by a Standard K-Epsilon turbulence model.

From Thai-Van et al., the other correlations are estimated using:

Figure 11. EQUATION_DISPLAY

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

(2533)

Figure 12. EQUATION_DISPLAY

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

(2534)

Figure 13. EQUATION_DISPLAY

τ_{I} = \frac{τ_{T}}{σ_{0}} \frac{1}{\sqrt{1 + C_{β} ξ^{2}}}

(2535)

Figure 14. EQUATION_DISPLAY

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

(2536)

Figure 15. EQUATION_DISPLAY

τ_{R} = τ_{D} (1 + \frac{ρ_{c}}{ρ_{d}} C_{V M})

(2537)

where:

$τ_{D}$ is a basic particle timescale. It is used here more generally to represent the inverse of the linearized drag coefficient. For more information on this coefficient, see Contribution of Drag to Turbulent Dispersion.
$b$ is the ratio of the coefficients of the continuous/dispersed phase acceleration terms in the particle equation of motion:

Figure 16. EQUATION_DISPLAY

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

(2538)

$C_{V M}$ is a virtual mass coefficient for particle response, with a default value of 0.5 if a virtual mass force model has not been specified for the mean flow.
$η$ is the particle-eddy interaction time, which is scaled by particle relaxation time:

Figure 17. EQUATION_DISPLAY

η = \frac{τ_{I}}{τ_{R}}

(2539)

$C_{β}$ is a calibration coefficient for the cross-trajectories effect. Following Deutsch and Simonin [449], choose it as 0.45 for the direction parallel to the mean relative velocity, and as 1.8 for orthogonal directions. The default value is 1.8.
$ξ$ is the particle slip velocity, scaled by turbulent fluctuation velocity:

Figure 18. EQUATION_DISPLAY

ξ = \frac{| v_{r} |}{\sqrt{\frac{2}{3} k_{c}}}

(2540)

Although Thai-Van et al. include a correction for fluctuations in this slip velocity estimate, this correction is neglected in Simcenter STAR-CCM+.

Contribution of Drag to Turbulent Dispersion

Under Stokesian (viscous drag) conditions, the instantaneous drag force on a solid spherical particle of fixed diameter $d$ and density $ρ_{d}$ , due to velocity relative to a fluid of viscosity $μ_{c}$ , is:

Figure 19. EQUATION_DISPLAY

\frac{ρ_{d} α_{d}}{τ_{D}} (v_{d} - v_{c})

(2541)

Figure 20. EQUATION_DISPLAY

τ_{D} = \frac{ρ_{d} d^{2}}{18 μ_{c}} = \frac{d^{2}}{18 ν_{c}} (\frac{ρ_{d}}{ρ_{c}})

(2542)

The Reynolds average of the instantaneous drag force can be written as the sum of a mean drag force term and a turbulent dispersion term:

Figure 21. EQUATION_DISPLAY

\bar{\frac{ρ_{d} α_{d}}{τ_{D}} (v_{d} - v_{c})} = \frac{ρ_{d}}{τ_{D}} \overline{α_{d}} (v_{r} + v_{T D})

(2543)

Figure 22. EQUATION_DISPLAY

v_{r} = {〈 v_{d} 〉}_{d} - {〈 v_{c} 〉}_{c}

(2544)

Figure 23. EQUATION_DISPLAY

v_{T D} = {〈 v_{c} 〉}_{c} - {〈 v_{c} 〉}_{d}

(2545)

Using the Boussinesq closure to model the last term:

Figure 24. EQUATION_DISPLAY

v_{T D} = \frac{\bar{α_{c} v_{c}}}{\overline{α_{c}}} - \frac{\bar{α_{d} v_{c}}}{\overline{α_{d}}} = \frac{\bar{α_{c} v_{c}'}}{\overline{α_{c}}} - \frac{\bar{α_{d} v_{c}'}}{\overline{α_{d}}} \approx - \frac{{ν_{c}}^{t}}{σ_{α}} (\frac{\nabla \overline{α_{c}}}{\overline{α_{c}}} - \frac{\nabla \overline{α_{d}}}{\overline{α_{d}}})

(2546)

For non-Stokesian drag laws, and for other interfacial area models, the result of Reynolds averaging no longer has a straightforward closure, but the above result is still used as an approximation.

For small particles, it is possible to compensate partly for the non-Stokesian drag laws. This compensation is done by including a fluctuating contribution to the mean particle relative velocity used in the drag coefficient calculation (for example see Thai-Van et al.). However, this practice is not followed here, since more generally it is difficult to say whether particle-eddy interaction reduces or enhances mean slip [431].

For the general case, replace the Stokesian coefficient

Figure 25. EQUATION_DISPLAY

\frac{ρ_{d}}{τ_{D}} \overline{α_{d}}

(2547)

with the mean linearized drag coefficient ${A_{c d}}^{D}$ evaluated using a mean slip velocity $v_{r}$ and mean interfacial area $a_{c d}^{'''}$ :

Figure 26. EQUATION_DISPLAY

{A_{c d}}^{D} = \frac{a_{c d}^{'''}}{8} ρ_{c} C_{D} | v_{r} |

(2548)

to give these mean drag and turbulent dispersion forces:

Figure 27. EQUATION_DISPLAY

{F_{c d}}^{D} = {A_{c d}}^{D} v_{r}

(2549)

Figure 28. EQUATION_DISPLAY

{F_{c d}}^{T D} = {A_{c d}}^{D} v_{T D}

(2550)

This form of turbulent dispersion force and its closure has been known for some time. For example, see Deutsch and Simonin [449] or Gosman et al. [468]. It has also recently been reconciled with alternative approaches in Drew [452] and Burns et al. [437].

Contribution of Virtual Mass to Turbulent Dispersion

The instantaneous virtual mass force reads:

Figure 29. EQUATION_DISPLAY

ρ_{c} C_{V M} α_{d} (a_{d} - a_{c})

(2551)

where $a$ is the acceleration of the corresponding phase. Averaging the above equation yields:

Figure 30. EQUATION_DISPLAY

\begin{matrix} ρ_{c} C_{V M} \overline{α_{d}} ({〈 a_{d} 〉}_{d} - {〈 a_{c} 〉}_{d}) = \\ ρ_{c} C_{V M} \overline{α_{d}} ({〈 a_{d} 〉}_{d} - {〈 a_{c} 〉}_{c}) + ρ_{c} C_{V M} \overline{α_{d}} (\frac{〈 α_{c}' a_{c} 〉}{\overline{α_{c}}} - \frac{〈 α_{d}' a_{c} 〉}{\overline{α_{d}}}) \end{matrix}

(2552)

The first term is for the inertia of the mean flow. The second one contributes to the turbulent dispersion force. It can be modeled using the Boussinesq closure, assuming particle relaxation time $τ_{R}$ as an estimate of the acceleration time scale:

Figure 31. EQUATION_DISPLAY

\frac{〈 α_{c}' a_{c} 〉}{\overline{α_{c}}} - \frac{〈 α_{d}' a_{c} 〉}{\overline{α_{d}}} \approx - \frac{{ν_{c}}^{t}}{τ_{R} σ_{α}} (\frac{\nabla \overline{α_{c}}}{\overline{α_{c}}} - \frac{\nabla \overline{α_{d}}}{\overline{α_{d}}})

(2553)

Figure 32. EQUATION_DISPLAY

τ_{R} = \frac{ρ_{d} \overline{α_{d}}}{{A_{c d}}^{D}} (1 + \frac{ρ_{c}}{ρ_{d}} C_{V M})

(2554)

Combining the contributions from drag and virtual mass terms:

Figure 33. EQUATION_DISPLAY

\begin{matrix} {A_{c d}}^{D} v_{T D} - ρ_{c} C_{V M} \overline{α_{d}} \frac{{ν_{c}}^{t}}{τ_{R} σ_{α}} (\frac{\nabla \overline{α_{c}}}{\overline{α_{c}}} - \frac{\nabla \overline{α_{d}}}{\overline{α_{d}}})) = \\ (1 + \frac{ρ_{c} C_{V M}}{ρ_{d} + ρ_{c} C_{V M}}) {A_{c d}}^{D} v_{T D} \end{matrix}

(2555)

For bubbly flows ( $ρ_{c} ≫ ρ_{d}$ ), including a contribution from the virtual mass can double the overall turbulent dispersion coefficient.

While virtual mass is always important in modeling particle response in a turbulent flow, the modeling of the mean virtual mass force is optional, as not all applications include an accelerating mean flow.

Particle Induced Turbulence

This section outlines the derivation of the source term for the continuous-phase fluctuation kinetic energy equation. Out of all possible interphase forces, the drag force is assumed to make the dominant contribution to this term.

Phase averaging of the kinetic energy of the instantaneous flow gives the following composition for the kinetic energy of the velocity fluctuations in the continuous-phase:

Figure 34. EQUATION_DISPLAY

k_{c} = \frac{1}{2} {〈 {v_{c}}^{"} \cdot {v_{c}}^{"} 〉}_{c} = \frac{1}{2} {{〈 v_{c} \cdot v_{c} 〉}_{c} - {〈 v_{c} 〉}_{c} \cdot {〈 v_{c} 〉}_{c}}

(2556)

where:

$\frac{1}{2} {〈 {v_{c}}^{"} \cdot {v_{c}}^{"} 〉}_{c}$ is the mean kinetic energy held in the instantaneous flow. The balance equation can be derived by Reynolds averaging the scalar product of the instantaneous phase velocity $v_{c}$ with the instantaneous phase momentum equation.
$\frac{1}{2} {〈 v_{c} 〉}_{c} \cdot {〈 v_{c} 〉}_{c}$ is the kinetic energy that is held in the mean flow. The equation can be derived from the scalar product of the phase mean velocity ${〈 v_{c} 〉}_{c}$ with the Reynolds averaged phase momentum equation.

The equation for $k_{c}$ can be derived by subtracting these two equations. Then an instantaneous force term, assumed proportional to slip velocity $v_{d} - v_{c}$ , and to disperse phase volume fraction $α_{d}$ ,

Figure 35. EQUATION_DISPLAY

f_{c} = \frac{ρ_{d} α_{d}}{τ_{D}} (v_{d} - v_{c})

(2557)

makes the following contribution to the equation for $k_{c}$ :

Figure 36. EQUATION_DISPLAY

S_{k c} = \bar{v_{c} \cdot f_{c}} - 〈 v_{c} 〉 \cdot \overline{f_{c}} = \bar{{v_{c}}^{"} \cdot f_{c}}

(2558)

Figure 37. EQUATION_DISPLAY

= \frac{ρ_{d}}{τ_{D}} \overline{α_{d}} {{〈 {v_{c}}^{"} \cdot {v_{d}}^{"} 〉}_{d} - {〈 {v_{c}}^{"} \cdot {v_{c}}^{"} 〉}_{d} + {〈 {v_{c}}^{"} 〉}_{d} \cdot [{〈 v_{d} 〉}_{d} - {〈 v_{c} 〉}_{c}]}

(2559)

Figure 38. EQUATION_DISPLAY

= {A_{c d}}^{D} {{〈 {v_{c}}^{"} \cdot {v_{d}}^{"} 〉}_{d} - {〈 {v_{c}}^{"} \cdot {v_{c}}^{"} 〉}_{d} - v_{T D} \cdot v_{r}}

(2560)

The final results are obtained by using the drag linearization coefficient ${A_{c d}}^{D}$ and dispersion velocity $v_{T D}$ from the earlier section Contribution of Drag to Turbulent Dispersion, together with this identity:

Figure 39. EQUATION_DISPLAY

v_{T D} = {〈 v_{c} 〉}_{c} - {〈 v_{c} 〉}_{d} = \frac{\bar{[{〈 v_{c} 〉}_{c} - v_{c}] α_{d}}}{\overline{α_{d}}} = - {〈 {v_{c}}^{"} 〉}_{d}

(2561)

$S_{k c}$ can now be modeled either with the help of a turbulent response model or with the help of the Tchen closures, when taken together with this further approximation:

Figure 40. EQUATION_DISPLAY

{〈 v_{c} \cdot v_{c} 〉}_{d} \approx {〈 v_{c} \cdot v_{c} 〉}_{c} = {2 k}_{c}

(2562)

The first two terms in Eqn. (2560) represent an equilibrium or transfer between the covariance of dispersed-continuous velocity fluctuation and the turbulent kinetic energy of the continuous phase.

The last term can be viewed as the dot product of turbulent dispersion force with the mean slip velocity (or alternatively, of the mean drag force with the turbulent dispersion velocity). Its role depends on the alignment of these two vectors.

The dot product can often be small or zero in a fully developed and converged pipe flow. In a horizontal pipe flow case, this small value is because the magnitude of the slip velocity is likely to be small under fully developed conditions. In a vertical pipe flow case, it is because although the slip is not small, it is orthogonal to any gradients of volume fraction under fully developed conditions.

However, the term can be positive or negative, and large, under developing conditions or during convergence. It makes a positive contribution to continuous phase turbulent kinetic energy when mean drag force and turbulent dispersion force are acting in opposite directions. For example, it contributes when particles are falling, or when bubbles are rising, against their own concentration gradients.