Drag

For general applications, the inter-phase drag force is computed as a function of the drag coefficient.

For a continuous-dispersed phase interaction, the force acting on the dispersed phase $i$ due to the drag of phase $j$ is implemented as:

Figure 1. EQUATION_DISPLAY

{F_{i j}}^{D} = A_{D} v_{r}

(1927)

where $A_{D}$ is the linearized drag coefficient and $v_{r}$ is the relative velocity between the phases $i$ and $j$ . $A_{D}$ acts as a linear multiplier of the relative velocity that is defined as:

Figure 2. EQUATION_DISPLAY

v_{r} = v_{j} - v_{i}

(1928)

For a multiple flow regime phase interaction, the total drag force is calculated as:

Figure 3. EQUATION_DISPLAY

F_{ij}^{D} = \sum_{t = f r, i r, sr}^{} W_{t} A_{D, t} v_{r}

(1929)

where $W_{t}$ is the weight function.

The calculation of the linearized drag coefficient is analogous to the continuous-dispersed phase interaction. For the first dispersed regime, the primary phase is considered as the continuous phase and the secondary phase as the dispersed phase. For the second dispersed regime, the secondary phase is considered as the continuous phase and the primary phase as the dispersed phase.

The linearized drag coefficient is related to the standard engineering definition of the drag coefficient $C_{D}$ for particles by:

Figure 4. EQUATION_DISPLAY

A_{D} = C_{D} \frac{1}{2} ρ_{c} | v_{r} | (\frac{a_{cd}}{4})

(1930)

where $a_{c d}$ is the interfacial area density. The factor $\frac{a_{c d}}{4}$ represents the projected area of the equivalent spherical particle. Corrections for the shape of non-spherical particles can be absorbed into the model for the drag coefficient $C_{D}$ .

The standard drag coefficient $C_{D}$ is composed of two factors:

Figure 5. EQUATION_DISPLAY

C_{D} = {f_{D} C}_{D∞}

(1931)

where $C_{D∞}$ is the single-particle drag coefficient, which is determined from the measurement of terminal velocity of a particle in an unbounded stagnant continuous phase. $f_{D}$ is a drag correction factor for the shape of non-spherical particles, which is determined from the measurement of terminal velocity for multiple particles of the same size traveling together at a given concentration. This correction can cover a number of physical effects that can be difficult to separate in practice, such as hindering due to increased turbulence interaction, or co-operative motion due to swarming or coalescence.

Force Balance Under Terminal Velocity Conditions

Under terminal velocity conditions, the pressure buoyancy and drag forces are in equilibrium and the two-phase momentum conservation equations for the continuous and dispersed phases reduce to the result:

Figure 6. EQUATION_DISPLAY

V^{2} = \frac{8 α_{c} α_{d}}{a} \frac{Δ ρ}{ρ_{c}} \frac{g}{C_{D}}

(1932)

where:

$Δ ρ = | ρ_{c} - ρ_{d} |$

$a$ is the interfacial area density for the two-phase interaction.

There are two alternative standard models for interfacial area density in Simcenter STAR-CCM+, which leads to:

Figure 7. EQUATION_DISPLAY

V^{2} = {\begin{matrix} \begin{matrix} \begin{matrix} \frac{4}{3} \frac{Δ ρ}{ρ_{c}} \frac{g d}{C_{D}} α_{c} for Spherical Particle Interfacial Area \end{matrix} \\ \begin{matrix} \frac{4}{3} \frac{Δ ρ}{ρ_{c}} \frac{g d}{C_{D}} for Symmetric Interfacial Area \end{matrix} \end{matrix} \end{matrix}

(1933)

Eqn. (1933) is a key equation for the calibration of single-particle drag models and multi-particle correction factors from experimental measurements of terminal velocity $V$ as a function of particle size $d$ , and continuous phase volume fraction $α_{c}$ .

The results of such a two-phase calibration are applied in the multiphase context, and in flows where particles are not traveling at terminal velocity. The terminal velocity is replaced with the computed relative velocity $| v_{r} |$ when determining the drag coefficient $C_{D}$ and drag correction factor ${f_{i j}}^{D}$ for each phase interaction $i j$ .

Role of Interfacial Area Density Method in Force Balance

The Spherical Particle Interfacial Area method is geometrically exact for spherical particles of a given diameter. Any corrections for non-spherical shape can be absorbed into the drag coefficient. Regardless of particle shape, the term $α_{c} Δ ρ g d$ in the force balance correctly reflects the change in effective buoyancy force with concentration due to the reducing difference in density between the equivalent-volume spherical particle and the two-phase mixture as $α_{c}$ falls.

The Symmetric Interfacial Area method is used where robustness is required across a wide range of volume fraction, so that the maximum packing of spherical particles can be exceeded. If an additional factor of $1 / α_{c}$ is absorbed into the drag correction factor that is used with the Symmetric Interfacial Area, both models give the same result, up to the maximum packing point. This adjustment is done automatically for the standard drag correction factor models that are implemented in Simcenter STAR-CCM+. However, include the adjustment explicitly when using the Symmetric Interfacial Area Density for drag correction models from field functions or exponents.

Special Considerations for Bubble Drag Coefficients

The drag laws for bubbles which allow for reduced drag due to internal circulation includes Tomiyama/Pure, Tomiyama/Moderately Contaminated, Bozzano-Dente, and the Wang drag law. All of these are implemented in Simcenter STAR-CCM+ as a function of ( $R e$ , $E o$ , $M o$ ).

These laws have been formulated in the literature to reproduce the correct drag coefficient at terminal velocity. However, the laws, can result in unphysically low estimates of drag coefficient for small bubbles traveling at high slip velocity.

In Simcenter STAR-CCM+, this problem is avoided by imposing a lower limit on the bubble drag coefficient. This limit is based on the lowest value of drag coefficient observed for a bubble of any size when it is traveling at its terminal velocity. This point corresponds to the slip velocity just before the bubble starts to distort from its spherical shape. It corresponds to a Weber number of about 2, so this is a physically meaningful constraint to apply at all bubble speeds.

Drag Correction Factors

The objective of drag correction methods is to reproduce experimental terminal velocities correctly as a function of particle concentration.

The traditional approach of Richardson and Zaki [534] is widely applicable for spherical particles of a given size. Drag correction for bubble columns is more complex because the bubble size can change as a function of concentration. In traditional approaches (for example, a model for the churn turbulent regime [481]), bubble drag correction has implicitly included the effects of turbulence interactions and of coalescence and breakup. A more recent trend is to isolate the effect of concentration alone on the terminal velocity of bubbles of a specified size (for example Simmonet et al. [545]).

To consider all these possibilities, define correction factors for particle size, terminal velocity, and drag coefficient. Let:

Figure 8. EQUATION_DISPLAY

f_{d} = \frac{d}{d_{\infty}}

(1934)

Figure 9. EQUATION_DISPLAY

f_{V} = \frac{V}{V_{\infty}}

(1935)

Figure 10. EQUATION_DISPLAY

f_{D} = \frac{C_{D}}{C_{D \infty}}

(1936)

where the subscript $\infty$ indicates conditions for a single particle in an unbounded stagnant fluid.

Setting $α_{c} = 1$ in the force balance Eqn. (1933), gives the force balance for a single particle traveling at terminal velocity in an unbounded medium:

Figure 11. EQUATION_DISPLAY

{V_{\infty}}^{2} = \frac{4}{3} \frac{Δ ρ}{ρ_{c}} \frac{g d_{\infty}}{C_{D \infty}}

(1937)

This result applies for either of two interaction area density models. Then taking the ratio of Eqn. (1934) over Eqn. (1936) and combining with Eqn. (1937) leads to the following relation between the three correction factors:

Figure 12. EQUATION_DISPLAY

{f_{V}}^{2} = {\begin{matrix} \begin{matrix} \frac{α_{c} f_{d}}{f_{D}} for Spherical Particle Interfacial Area \end{matrix} \\ \begin{matrix} \frac{f_{d}}{f_{D}} for Symmetric Interfacial Area \end{matrix} \end{matrix}

(1938)

The Simcenter STAR-CCM+ model does not include the size correction factor $f_{d}$ when computing $f_{V}$ for two reasons:

Modern computational methods focus on drag and drag correction for particles of a specified size, and particle size distribution is considered as a separate model.
Size change is only relevant for drag correction of large bubbles whose drag coefficient is independent of Re.

Working with drag correction for particles of a specified size, the above relationship reduces to:

Figure 13. EQUATION_DISPLAY

{f_{V}}^{2} = {\begin{matrix} \begin{matrix} \frac{α_{c}}{f_{D}} for Spherical Particle Interfacial Area \end{matrix} \\ \begin{matrix} \frac{1}{f_{D}} for Symmetric Interfacial Area \end{matrix} \end{matrix}

(1939)

In Simcenter STAR-CCM+ the relationship is used to:

implement a drag correction factor $f_{D}$ for models that are specified in literature as using a velocity correction factor $f_{V}$ .
obtain a velocity correction factor $f_{V}$ for models that are specified in the literature as using drag correction factor $f_{D}$ . The velocity correction factor is required to evaluate the drag coefficient for a single particle.

This last point is important in reproducing terminal velocity where the drag coefficient $C_{D \infty}$ depends on the particle Reynolds number. Evaluate the single-particle drag coefficient at the single particle terminal velocity $V_{\infty}$ not at the terminal velocity $V$ for which you are trying to estimate the drag coefficient in a multi-particle system. Apply the velocity correction factor $f_{V}$ to the slip Reynolds number:

Figure 14. EQUATION_DISPLAY

C_{D \infty} \equiv C_{D} (R e_{\infty}) = C_{D} (\frac{R e}{f_{V}})

(1940)

Output Variables

To help you follow these drag computations, the following fields are available for inspection for each phase pair in a Simcenter STAR-CCM+ simulation:

The particle Eotvos number, $E o$ . This field only applies where the drag model depends on surface tension.
The particle slip velocity Reynolds number, $R e$
The single particle Reynolds number, $R e_{\infty} = \sqrt{\frac{f_{D}}{α_{c}}} R e$ ,

where $f_{D}$ is the drag correction factor in the multi-particle system, $α_{c}$ is the continuous phase volume-fraction, and $R e$ is the phase-pair Reynolds Number.
The single particle drag coefficient, $C_{D \infty}$
The drag correction factor, $f_{D}$

Fixed Exponent Regimes

Where correction factors can be modeled as the continuous phase volume fraction raised to some constant power, that is:

Figure 15. EQUATION_DISPLAY

f_{V} = \frac{V}{V_{\infty}} = α_{c}^{n}

(1941)

Figure 16. EQUATION_DISPLAY

f_{D} = \frac{C_{D} (R e, α_{c})}{C_{D} (R e_{\infty}, 1)} = α_{c}^{n_{D}}

(1942)

Then Eqn. (1942) leads to a simple relation between these exponents:

Figure 17. EQUATION_DISPLAY

n_{D} = {\begin{matrix} \begin{matrix} 1 - 2 n for Spherical Particle Interfacial Area \end{matrix} \\ \begin{matrix} - 2 n for Symmetric Interfacial Area \end{matrix} \end{matrix}

(1943)

For example, under Stokesian conditions ( ${R e}_{\infty} < 0.1$ ), the constant Richardson Zaki [534] exponent of $n = 4.65$ for velocity correction leads to a drag correction exponent of:

Figure 18. EQUATION_DISPLAY

n_{D} = - 8.3

(1944)

when using the Spherical Particle Interfacial Area, or:

Figure 19. EQUATION_DISPLAY

n_{D} = - 9.3

(1945)

when using the Symmetric Interfacial Area.