Particle Equations of Motion

Drag Force

The drag force in Eqn. (2957) is defined as:

where:

• $C_d$ is the drag coefficient of the particle.
• $\rho$ the density of the continuous phase.
• ${\mathbf{\text{v}}}_s=v-v_p$ is the particle slip velocity with $v$ being the instantaneous velocity of the continuous phase.
• $A_p$ is the projected area of the particle.

Eqn. (2962) can also be written as:

where ${\tau }_v$ is the momentum relaxation time-scale:

Drag Coefficient

The drag coefficient $C_d$ in Eqn. (2962) is a function of the small-scale flow features around the individual particles.

These features are impractical to resolve spatially, and so the usual practice is to obtain the drag coefficient from correlations, typically derived from experiment or theoretical studies. These correlations depend on the nature of the dispersed phase. Droplets, bubbles, and solid particles have different correlations. The shape of the particle, the presence of inter-phase mass and energy transfer, and so on, affect these correlations.

Simcenter STAR-CCM+ provides the following methods for defining the drag coefficient:

• Schiller-Naumann correlation
• Liu dynamic drag coefficient
• Haider and Levenspiel

The following drag force coefficients are available with the DEM drag force model:

• Clift
• Di Felice
• Gidaspow

Schiller-Naumann Correlation

The Schiller-Naumann correlation [696] is suitable for spherical solid particles, liquid droplets, and small-diameter bubbles. It is formulated as:

$$C_d=\{\begin{array}{cc}\frac{24}{{\text{Re}}_p}(1+0.15{\text{Re}}_p^{0.687})\hfill & {\text{Re}}_p\le {10}^3\\ 0.44\hfill & {\text{Re}}_p>{10}^3\end{array}$$

(2965)

where ${\text{Re}}_p$ is the particle Reynolds number that is defined as:

$${\text{Re}}_p\equiv \frac{\rho \left|{\mathbf{\text{v}}}_s\right|D_p}{\mu }$$

(2966)

where $D_p$ is the particle diameter and $\mu$ is the dynamic viscosity. This correlation is available only when the continuous phase is viscous.

Liu Dynamic Drag Coefficient

The Liu dynamic drag coefficient [673] accounts for the dependence of the drag of a liquid droplet on its distortion under the action of aerodynamic forces. As a basis, it uses the following expression for the drag coefficient of an undistorted sphere:

$$C_{d,sphere}=\{\begin{array}{cc}\frac{24}{{\text{Re}}_p}(1+\frac{1}{6}{\text{Re}}_p^{2/3})\hfill & {\text{Re}}_p\le 1000\\ 0.424\hfill & {\text{Re}}_p>1000\end{array}$$

(2967)

As the distortion of the droplet increases, its shape is assumed to become a disk whose axis is aligned with the slip velocity. This shape increases the drag on the droplet. The Liu drag coefficient models this effect by noting that the high Reynolds number limit of the drag coefficient of a disk is 1.54. It then assumes that the disk drag is 1.54/0.424 higher than the sphere drag at all Reynolds numbers, and that the drag of intermediate shapes can be interpolated between those two extremes. So:

$$C_d=C_{d,sphere}(1+2.632y)$$

(2968)

The interpolation factor $y$ is 0 for a sphere and 1 for a disk. It is identified as the TAB distortion, Eqn. (3110), and hence the Liu drag coefficient requires the TAB distortion model to calculate this quantity.

Clift Method

This method is available only with cylindrical particles. The drag coefficient $C_d$ depends on the particle Reynolds number and the aspect ratio of the cylinder. The Clift drag coefficient is given as:

• For long cylindrical particles (aspect ratio > 1):

  $$\begin{array}{c}\begin{array}{c}\begin{array}{c}{\text{Re}}_p\le 5:C_d=9.689{{\text{Re}}_p}^{-0.78}(1+0.147{{\text{Re}}_p}^{0.82})\end{array}\\ 5<{\text{Re}}_p\le 40:C_d=9.689{{\text{Re}}_p}^{-0.78}(1+0.227{{\text{Re}}_p}^{0.55})\end{array}\\ \begin{array}{c}40<{\text{Re}}_p\le 400:C_d=9.689{{\text{Re}}_p}^{-0.78}(1+0.0838{{\text{Re}}_p}^{0.82})\\ 400<{\text{Re}}_p:C_d=0.99{\gamma }^{-0.12}{E}^{-0.08}\end{array}\end{array}$$

  (2969)

  where $\gamma ={\rho }_p/\rho$ is the ratio of particle density to fluid density and $E=h/d$ is the aspect ratio.

• For flattened cylindrical particles (aspect ratio < 1):

  $$\begin{array}{c}\begin{array}{c}\begin{array}{c}\begin{array}{c}{\text{Re}}_p\le 0.01:C_d=\frac{64}{\pi {\text{Re}}_p}(1+\frac{{\text{Re}}_p}{2\pi })\\ \\ 0.01<{\text{Re}}_p\le 1.5:C_d=\frac{64}{\pi {\text{Re}}_p}(1+{10}^x)\end{array}\\ (x=-0.883+0.906\text{log}\left({\text{Re}}_p\right)-0.025{\left(\text{log}{\text{Re}}_p\right)}^2)\end{array}\\ 1.5<{\text{Re}}_p\le 133:C_d=\frac{64}{\pi {\text{Re}}_p}(1+0.138{\text{Re}}_p^{0.792})\\ \end{array}\\ 133<{\text{Re}}_p:C_d=1.17\end{array}$$

  (2970)

Di Felice Drag Coefficient Method

This drag method introduces an extra term in the fluid drag force expression (Eqn. (2962)) to account for the presence of other particles around a particle. The Di Felice drag coefficient is given as:

$$C_d={(0.63+\frac{4.80}{\sqrt{{ϵ}_i{\text{Re}}_p}})}^2{ϵ}^{2-\zeta }$$

(2971)

where ${\mathrm{Re}}_p$ defines the particle Reynolds Number, ${ϵ}_i$ is the void fraction around a particle, and:

$$\zeta =3.7-0.65\text{exp}\left[-0.5{\left(1.5-\text{log}[ϵ{\text{Re}}_p]\right)}^2\right]$$

(2972)

The term ${ϵ}^{2-\zeta }$ accounts for the effect of enhanced drag on a particle, due to the presence of other particles around it.

Gidaspow Drag Coefficient Method

The Gidaspow model is useful in systems with regions of high particle density, for example, in fluidized beds. The model covers a wide range of void fractions by combining the Wen-Yu and the Ergun equations [723].

Simcenter STAR-CCM+ incorporates the equations in drag coefficient form as follows:

$$C_d=\frac{4}{3}(150\frac{1-ϵ}{ϵ{\mathrm{Re}}_p}+1.75)\text{ if }ϵ<{ϵ}_{min}$$

(2973)

Otherwise:

$$C_d=\begin{array}{c}\\ \end{array}\{\begin{array}{cc}\frac{24(1+0.15\cdot {\mathrm{Re}}_p^{0.687})}{ϵ{\text{Re}}_p}{ϵ}^{1-w}& \text{if}ϵ{\text{Re}}_p\le 1000\\ \hfill 0.44\hfill & \text{if}ϵ{\text{Re}}_p>1000\end{array}$$

(2974)

where:

• $ϵ$ is the void fraction.
• ${ϵ}_{min}$ is the cutoff void fraction.
• $w$ is the Wen-Yu exponent.

EMMS Drag Coefficient Method

The EMMS model is useful in systems with regions of high particle density, for example, in fluidized beds. The model is suitable for simulations where the influence of heterogeneity from particle clusters becomes significant. The EMMS drag method models interactions within the fluidized bed by considering the meso-scale effects and an energy minimization approach for particle interchange for various interaction scales.

For a heterogeneous fluidized bed, the EMMS drag coefficient ([568]) is:

$$C_d=\frac{3}{4}\frac{{\alpha }_pϵ{\rho }_g}{d_p}{C}_{D0}\left|{\mathbf{\text{v}}}_p-\mathbf{\text{v}}\right|{ϵ}^{-2.65}{H}_D$$

(2975)

where:

• $d_p$ is the particle diameter.
• ${\mathbf{\text{v}}}_p$ and ${\mathbf{\text{v}}}_f$ are the particle velocity and the fluid velocity, respectively.
• ${\alpha }_p$ and $ϵ$ are the volume fraction of the particle and the void fraction, respectively.
• $C_{D0}$ is the drag coefficient of an individual particle that is calculated from the correlation of Schiller and Naumann in Eqn. (2965).
• $H_D=f\left({\alpha }_g\right)$ is the heterogeneity index which represents the ratio of the mean drag force that is experienced by particles in a heterogeneous system to the drag force in a homogeneous system. In Simcenter STAR-CCM+, the Lu, Wang and Li correlation ([509]) is used to define the heterogeneity index for both Geldart A and B particle groups and bubbly, fast, and pneumatic fluidization regimes. The flow regimes are identified based on the gas velocity and the resulting fluidization characteristics in the system ([469]).The Geldart classification system groups particles based on their size, density, and fluidization characteristics.

Haider and Levenspiel Drag Coefficient Method

The drag coefficient $C_d$ depends on the particle Reynolds number and the particle sphericity:

$$C_d=\frac{24}{{\text{Re}}_p}(1+A{\text{Re}}_p^B)+\frac{C}{(1+\frac{D}{{\text{Re}}_p})}$$

(2976)

where $\varphi$ is the particle sphericity and:

$$\begin{array}{ccc}A& =& 8.1716e^{-4.0665\varphi }\hfill \\ B& =& 0.0964+0.5565\varphi \hfill \\ C& =& 73.690e^{-5.0746\varphi }\hfill \\ D& =& 5.3780e^{6.2122\varphi }\hfill \end{array}$$

(2977)

Pressure Gradient Force

The pressure gradient force is as defined as:

where $V_p$ is the volume of the particle and $\nabla p_{static}$ is the gradient of the static pressure in the continuous phase.

Virtual Mass Force

The virtual mass force [646] is:

where $C_{vm}$ is the virtual mass coefficient. The default value of 0.5 for this coefficient is for a sphere in a uniform, inviscid, incompressible flow [675]. The operator D/Dt is the material (substantive) derivative.

Drag Torque

For DEM particles, the drag torque $M_b$ (see Eqn. (2960)) reduces the difference in the rotational velocities between a particle and the fluid in which it is immersed.

The drag torque $M_b$ applied to a DEM particle is defined as [112]:

where:

• $C_R$ is the rotational drag coefficient.
• $\mathbf{\text{Ω}}$ is the relative angular velocity of the particle to the fluid.

where:

• $\mathbf{\text{v}}$ is the fluid velocity.
• ${\mathbf{\text{ω}}}_p$ is the angular velocity of the particle.

The rotational drag coefficient $C_R$ in Eqn. (2980) is defined as [112]:

${\text{Re}}_R$ is the rotational Reynolds number, defined as:

Although Sommerfeld restricts the rotational Reynolds number to values below 1000, Simcenter STAR-CCM+ places no upper limit on it.

Lift Forces

Lift forces in DEM simulations can arise from particle spin, particle shear, or both. In this context, lift forces refer to mean forces normal to the particle velocity—they are not necessarily forces in the upward direction. Lift forces in statistical Lagrangian simulations can arise from particle shear only.

Particle Spin Lift Force

This force applies to a spinning particle moving relative to a fluid. The force is given by:

$${\mathbf{\text{F}}}_{LR}=\frac{\rho \pi }{8}{D}_p^2{C}_{LR}\left|{\mathbf{\text{v}}}_s\right|\frac{\mathbf{\text{Ω}}\times {\mathbf{\text{v}}}_s}{\left|\mathbf{\text{Ω}}\right|}$$

(2984)

Sommerfeld [112] gives the coefficient of rotational lift $C_{LR}$ as:

$$C_{LR}=0.45+\left(\frac{{\text{Re}}_R}{{\text{Re}}_p}-0.45\right)e^{-0.05684{{\text{Re}}_R}^{0.4}{{\text{Re}}_p}^{0.3}}$$

(2985)

Sommerfeld reports this formula for $C_{LR}$ as valid for particle Reynolds numbers below 140. You can also define $C_{LR}$ using a field function.

Particle Shear Lift Force

This force applies to a particle moving relative to a fluid where there is a velocity gradient in the fluid orthogonal to the relative motion.

Saffman [733] gives the shear lift force as:

$${\mathbf{\text{F}}}_{LS}=C_{LS}\frac{\rho \pi }{8}{D}_p^3({\mathbf{\text{v}}}_S\times \mathbf{\text{ω}})$$

(2986)

where:

• $C_{LS}$ is the shear lift coefficient.
• $\rho$ is the particle density.
• ${\mathbf{v}}_S=\mathbf{v}-{\mathbf{v}}_p$ is the slip velocity, the difference between the fluid velocity and the particle velocity.
• $\mathbf{\text{ω}}=\nabla \times \mathbf{\text{v}}$ , the curl of the fluid velocity.

You can define $C_{LS}$ as a field function or use either of two published definitions:

• Saffman [733] provides a definition that recovers the original Saffman asymptotic solution for low Reynolds numbers:

  $$C_{LS}=\frac{4.1126}{{{\text{Re}}_S}^{0.5}}$$

  (2987)
• Sommerfeld [112] provides a definition for a broader range of Reynolds numbers:

  $$C_{LS}=\frac{4.1126}{{{\text{Re}}_S}^{0.5}}f{\text{Re}}_p,{\text{Re}}_S$$

  (2988)

  where:

  $$f{\text{Re}}_p,{\text{Re}}_S=\{\begin{array}{cc}\left(1-0.3314{\beta }^{0.5}\right)e^{-0.1{\text{Re}}_p}+0.3314{\beta }^{0.5},& ({\text{Re}}_p\le 40)\\ 0.0524{\left(\beta {\text{Re}}_p\right)}^{0.5},& ({\text{Re}}_p>40)\end{array}$$

  (2989)

  and:

  $$\beta =0.5\frac{{\text{Re}}_S}{{\text{Re}}_p}$$

  (2990)

  where:

  ${\text{Re}}_S$ is the Reynolds number for shear flow:

  $${\text{Re}}_S=\frac{\rho D_p^2\left|\mathbf{\text{ω}}\right|}{\mu }$$

Gravity Force

The gravity force is defined as:

where $\mathbf{\text{g}}$ is the gravitational acceleration vector.

Forces in Moving Reference Frames

where $\omega$ is the angular velocity vector of the rotating reference frame and $r$ is the distance vector to the axis of rotation.

Coulomb Force

Particles can enter the domain bearing electrical charges or, under some circumstances, can acquire electric charges after entering, resulting in an additional force acting on the particles. Charged particles appear in many applications, such as electrostatic spray painting and electrostatic precipitation in air pollution control.

A charged particle with a charge $q$ in an electric field $E$ experiences a force. This force is called Coulomb force and is given by:

User-Defined Body Force

The user-defined body force in Eqn. (2958) is defined as:

where ${\mathbf{\text{f}}}_u$ is a user-specified body force per unit volume.

Particle Position and Velocity

The variables describing the state of the particle change as a function of time and space due to the forces enumerated in the previous sections. In the Lagrangian framework, the state is evolved on a particle by particle basis while the Eulerian fields are frozen. The particle position and particle velocity are obtained by numerically integrating over time the particle momentum equation Eqn. (2956) and the position equation:

where $r_p\left(t\right)$ is the instantaneous position vector. A particle trajectory is the locus of points obtained by the integration of the position equation.

In the statistical Lagrangian approach, parcels of individual particles are tracked. The term parcel is used to denote a set of particles which all have the same state and evolve in the same way. A single integration of the state is applied to all the particles in the parcel, and any exchange with the Eulerian phase is multiplied by the number of particles that the parcel represents.

In subsequent sections, a quantity which is identical for a parcel and its particles, such as velocity in the preceding paragraph, is described as a “particle” quantity and given the subscript p. The term “parcel” and the subscript $\pi$ are reserved for quantities which are exclusively defined on parcels.

Continuous phase quantities reported by the parcel are cell-centered values, however the continuous phase velocity is interpolated to the parcel position. Cell velocity is the only continuous phase value that is interpolated.