Altair AcuSolve EDEM Coupling

Altair AcuSolve-EDEM coupling sequence for bi-directional coupled simulations is shown below.

Note: The detailed DEM simulation sequence is not shown here. For detailed information about the DEM simulation sequence, refer to the EDEM help manual.

Note:

The temperature and species equations are only solved when the heat transfer and/or mass transfer physics models are active.
Since the DEM time step is usually multiple orders lower than the CFD time step, the DEM solver loop is repeated multiple times per single CFD time step to ensure that the physical time is synchronized in both the solvers.
For coupled simulations, the data is always exchanged in SI units.
For unidirectional coupling, once the coupling forces are calculated and shared with EDEM, the fluid momentum equation is not updated with the coupling force because the effect of particles on fluid is ignored.

AcuSolve-EDEM coupling uses the Eulerian-Lagrangian approach for modeling fluid-particle flows where the fluid transport equations are solved in an Eulerian framework and the dispersed phase is represented as Lagrangian particles. The fluid phase is solved based on the volume-averaged Navier-Stokes equations and the Discrete Element Method (DEM) is used for computing the motion of the solid phase. This coupling strategy allows you to study the momentum and heat transfer at the individual particle scale.

Governing Equations

The volume-averaged Navier-Stokes equations for CFD-DEM momentum coupling are given by:(1)

\frac{\partial (ε_{f} ρ_{f})}{\partial t} + \nabla \cdot (ρ_{f} ε_{f} u_{f}) = 0

(2)

\frac{\partial}{\partial t} (ε_{f} ρ_{f} v_{f}) + \nabla \cdot (ρ_{f} ε_{f} v_{f} v_{f}) = - ε_{f} \nabla p - ε_{f} \nabla \cdot (τ_{f}) + ρ_{f} ε_{f} g - F_{p}

where,

$ε_{f}$ is the fluid volume fraction or porosity;
$v_{f}$ is the fluid velocity;
$ρ_{f}$ is the fluid density;
$p$ is the fluid pressure;
$τ_{f}$ is the viscous stress tensor;
$g$ is the gravity vector;
$F_{p} = \frac{1}{V_{c e l l}} \sum_{i} {(f^{d} + f^{l})}_{i}$ is the particle-fluid force exchange term;
$V_{c e l l}$ is the volume of fluid cell;
$f^{d}$ is the fluid drag force;
$f^{l}$ is the fluid lift force.

Note: Wherever applicable the following notations are used throughout the document. The subscript

f

denotes a fluid property,

p

denotes a particle property and

i

denotes that the calculation is done for an individual particle.

Drag Models

The drag models available in AcuSolve-EDEM coupling are listed below:

Ergun-Wen Yu:
The Ergun-Wen Yu model, also known as the Gidaspow model, reads as:(3)
$f_{i}^{d} = \{\begin{matrix} \frac{A (1 - ε_{f})}{18 ε_{f}^{2}} + \frac{B}{18 ε_{f}^{2}} {Re}_{i} ε_{f} \leq 0.8 \\ \frac{C_{d}}{24} {Re}_{i} ε_{f}^{- 3.65} ε_{f} > 0.8 \end{matrix}$

where ${Re}_{i}$ is the particle Reynolds number ${Re}_{i} = \frac{ρ_{f} ε_{f} |v_{f} - v_{p}| d_{p}}{μ_{f}}$ .

$v_{p}$ is the particle velocity and $d_{p}$ is the particle’s volume equivalent sphere diameter, A = 150 and B = 1.75. The values of these coefficients can be modified by you while specifying the model inputs.(4)
$C_{d} = \{\begin{matrix} \frac{24}{{Re}_{i}} (1 + 0.15 {Re}_{i}^{0.687}) {Re}_{i} \leq 1000 \\ 0.44 {Re}_{i} > 1000 \end{matrix}$

here $C_{d}$ is the drag coefficient.

The Ergun-Wen-Yu model is one of the most widely used drag models and is recommended for most of the fluid-particle flows since it works well for both dense phase and dilute phase regimes. In this model, the Ergun equation is used for fluid volume fractions less than 0.8 and the Wen-Yu equation for fluid volume fractions greater than 0.8.
DiFelice(5)
$f_{i}^{d} = \frac{C_{d}}{24} {Re}_{i} ε_{f}^{- χ}$

where,(6)
$χ = 3.7 - 0.65 e^{(- 0.5 {(1.5 - \log_{10} {Re}_{i})}^{2})}$
(7)
$C_{d} = {(0.63 + 4.8 {Re}_{i}^{- 0.5})}^{2}$

Unlike the Ergun-Wen-Yu correlation, the Di Felice correlation is a monotonic function of Reynolds number and porosity and does not have the step change in drag force evaluation.
Beetstra(8)
$f_{i}^{d} = (\frac{180 (1 - ε_{f})}{18 ε_{f}^{2}}) + ε_{f}^{2} (1 + 1.5 \sqrt{(1 - ε_{f})}) + (\frac{0.413}{24 ε_{f}^{2}}) (\frac{ε_{f}^{- 1} + 3 (1 - ε_{f}) ε_{f} + 8.4 {Re}_{i}^{- 0.343}}{1 + 10^{3 (1 - ε_{f})} {Re}_{i}^{- (1 + 4 (1 - ε_{f})) / 2}}) {Re}_{i}$
Rong(9)
$f_{i}^{d} = 0.5 C_{d} \frac{π d_{e}^{2}}{4} ρ_{f} |v_{f} - v_{p}| (v_{f} - v_{p}) ε_{f}^{2 - β - λ}$

where,(10)
$β = 2.65 (ε_{f} + 1) - (5.3 - 3.5 ε_{f}) ε_{f}^{2} e^{[- \frac{{(1.5 - \log_{10} {Re}_{i})}^{2}}{2}]}$
(11)
$λ = (1 - ε_{f}) (C - D e^{- 0.5 {(3.5 - \log_{10} {Re}_{i})}^{2}})$
(12)
$C = 39 φ - 20.6$
(13)
$D = 101.8 {(φ - 0.81)}^{2} + 2.4$

$φ$ is the sphericity of the particle, $d_{e}$ is the diameter of volume equivalent sphere.

Since the sphericity of the particle is considered while calculating the drag force, this model is strongly recommended for non-spherical particles compared to the other models available in AcuSolve.
Syamlal-O’Brien
The correlation reads as:(14)
$f_{i}^{d} = \frac{C_{d} {Re}_{i} ε_{f}}{24 v_{r, p}^{2}}$

where,(15)
$C_{d} = {(0.63 + \frac{4.8}{\sqrt{\frac{{Re}_{i}}{v_{r, s}}}})}^{2}$
(16)
$v_{r, s} = (A - 0.06 {Re}_{i} + \sqrt{{(0.06 {Re}_{i})}^{2} + 0.12 (2 B - A) + A^{2}})$
(17)
$A = ε_{f}^{4.14}$
(18)
$B = \{\begin{matrix} 0.8 ε_{f}^{1.28} ε_{f} \leq 0.85 \\ ε_{f}^{2.65} ε_{f} \geq 0.85 \end{matrix}$
Wen-Yu
The Wen-Yu correlation reads as:(19)
$f_{i}^{d} = \frac{C_{d}}{24} {Re}_{i} ε_{f}^{- 3.65}$
(20)
$C_{d} = \{\begin{matrix} \frac{24}{{Re}_{i}} (1 + 0.15 {Re}_{i}^{0.687}) {Re}_{i} \leq 1000 \\ 0.44 {Re}_{i} > 1000 \end{matrix}$
Schiller Nauman
The correlation reads as:(21)
$f_{i}^{d} = \frac{C_{d}}{24} {Re}_{i}$
(22)
$C_{d} = \{\begin{matrix} 24 (1 + 0.15 {Re}_{i}^{0.687}) {Re}_{i} \leq 1000 \\ 0.44 {Re}_{i} > 1000 \end{matrix}$

The drag force calculated does not consider the effect of surrounding particles, that is, volume fraction is not accounted for, and hence this model is strictly valid only for dilute phase flows.

Non-Spherical Drag Coefficient Models

The effect of a particle’s shape can be taken into account by using non-spherical drag coefficient models. There are two types of models available in AcuSolve which are listed below. If the non-spherical drag coefficient model is set to none the particles are assumed to be of spherical shape.

Isometric (Haider Levenspiel)
In this model the drag coefficient is considered to be a function of particle Reynolds number and sphericity. The instantaneous orientation of the particle is not taken into account. This type of model is applicable for particles with shapes closer to a sphere such as rocks, some grains (beans), and when the orientation of the particles is not critical. The user inputs required for this model are particle’s volume and sphericity. The Haider-Levenspiel correlation is given by:(23)
$C_{d}^{n s} = \frac{24}{{Re}_{i}^{}} [1 + A_{1} {({Re}_{i}^{})}^{A_{2}}] + \frac{A_{3}}{1 + \frac{A_{4}}{{Re}_{i}^{}}}$

where,(24)
$A_{1} = e^{(2.3288 - 6.4581 φ_{i} + 2.4486 φ_{i}^{2})}$
(25)
$A_{2} = 0.0964 + 0.5565 φ_{i}$
(26)
$A_{3} = 73.69 e^{(- 5.0748 φ_{i})}$
(27)
$A_{4} = 5.378 e^{(6.2122 φ_{i})}$
Non-spherical (Ganser and Holzer-Sommerfeld)
The Ganser and Holzer-Sommerfeld models consider both the shape and orientation of the particle. Since the orientation of the particles is also considered, this model is applicable to particle shapes such as disk, ellipsoid and elongated cylinder. The user inputs for these models are volume and aspect ratio of the particles.

The Ganser correlation is given by:(28)
$\frac{C_{d}^{n s}}{k_{2}} = \frac{24}{k_{1} k_{2} {Re}_{i}^{}} [1 + 0.1118 {(k_{1} k_{2} {Re}_{i}^{})}^{0.65657}] + \frac{0.4305}{1 + \frac{3305}{k_{1} k_{2} {Re}_{i}^{}}}$

Here $k_{1}$ and $k_{2}$ are Stokes and Newton shape factors respectively.

The Holzer-Sommerfeld correlation is given by:(29)
$C_{d}^{n s} = \frac{8}{{Re}_{i}^{}} \frac{1}{\sqrt{φ_{i}^{⊥}}} + \frac{16}{{Re}_{i}^{}} \frac{1}{\sqrt{φ_{i}}} + \frac{3}{\sqrt{{Re}_{i}^{}}} \frac{1}{φ_{i}^{\frac{3}{4}}} + 0.42 \times 10^{0.4 {(- \log_{10} φ_{i})}^{2}} \frac{1}{φ_{i}^{⊥}}$

where, $φ^{⊥}$ is the crosswise sphericity which is defined as the ratio of projected area of the volume equivalent sphere to the projected area of the particle perpendicular to the flow.

Lift Models

Generally the lift force acts in a direction normal to the relative motion of the fluid and particle. The two components of the lift force considered are Saffman force and Magnus force. The Saffman lift force is due to the pressure gradient on a non-rotating particle in the presence of a non-uniform shear velocity field while the Magnus lift force is due to the particle rotation in a uniform flow. Unlike spherical particles, the behavior of non-spherical particles in turbulent flows is much more complicated and the lift force acting on them can no longer be neglected. As the particle’s principal axis becomes inclined with the flow direction, the effect of lift force on the particle motion becomes significant.

There are three lift models available in AcuSolve:

Saffman-Magnus
This model is for spherical particles and hence the orientation is neglected whereas the last two models take the particle orientation into account while calculating the lift forces.

The correlation for the Saffman force is given by:(30)
${\vec{f}}_{i}^{S a f f m a n} = SLC \times C_{l s} w_{i} d_{i}^{2} \sqrt{μ_{f} ρ_{f}} {|ω|}^{- 0.5} (w_{i} \times ω)$

where $ω$ is the particle velocity curl, $w_{i}$ is the slip velocity of the particle. SLC is Saffman constant with a default value of 1.615. This value can be modified while specifying the model inputs. $C_{l s}$ is the Saffman lift coefficient and is given by the expression:(31)
$C_{l s} = \{\begin{matrix} e^{- 0.1 {Re}_{i}} + 0.3314 \sqrt{γ_{i}} (1 - e^{- 0.1 {Re}_{i}}) {Re}_{i} \leq 40 \\ 0.0524 \sqrt{γ_{i} {Re}_{i}} {Re}_{i} > 40 \end{matrix}$
(32)
$γ_{i} = \frac{|ω| d_{i}}{2 |w_{i}|}$

The correlation for the Magnus force is given by:(33)
${\vec{f}}_{i}^{M a g n u s} = M L C \times C_{l m} w_{i}^{2} π d_{i}^{2} ρ_{f} \sqrt{μ_{f} ρ_{f}} \frac{(ω_{r} \times w_{i})}{|ω_{r}| |w_{i}|}$
(34)
$ω_{r} = \frac{1}{2} ω - ω_{i}$

where MLC is Magnus constant with a default value of 0.125. This value can be modified while specifying the model inputs. $C_{l m}$ is the Magnus lift coefficient and is given by the expression:(35)
$C_{l m} = d_{i} \frac{|ω_{r}|}{|w_{i}|} \{\begin{matrix} 1 {Re}_{i} \leq 1 \\ (0.178 + 0.822 {Re}_{i}^{- 0.522}) {Re}_{i} > 1 \end{matrix}$
Saffman-Magnus non-spherical lift
This is similar to the Saffman Magnus lift model for spherical lift, except the lift coefficients are replaced by the nonspherical lift coefficient which is explained in the nonspherical lift model section below.
Non-spherical lift
In this model, the lift coefficient is assumed to be proportional to the drag coefficient and the correlation for the lift force is given by:(36)
$F_{L} = \frac{1}{2} C_{L} ρ_{f} \frac{π}{4} d_{p}^{2} {|v_{f} - v_{i}|}^{2} \{\begin{matrix} \frac{C_{L}}{C_{D}} = \sin^{2} α \cdot \cos α (in Newton Law region) \\ \frac{C_{L}}{C_{D}} = \frac{\sin^{2} α \cdot \cos α}{0.65 + 40 {Re}^{0.72}} (30 < Re < 1500) \end{matrix}$

Here the drag coefficient is obtained from the drag force calculation. The direction of the lift force is given by:(37)
${\hat{e}}_{L_{o}} = \frac{u_{i} \cdot {v^{'}}_{f i}}{|u_{i} \cdot {v^{'}}_{f i}|} \frac{(u_{i} \times {v^{'}}_{f i}) \times {v^{'}}_{f i}}{(u_{i} \times {v^{'}}_{f i}) \times {v^{'}}_{f i}}$

where, $u_{i}$ is the particle principal axis and ${v^{'}}_{f i}$ is the relative velocity vector of the particle.

Torque Models

By using the torque models, the rotational drag force on the rotating particles due to the inertia of fluid can be considered. The three types of torque models available in AcuSolve are pitching torque, rotational torque and a combination of both.

Pitching torque
When the center of pressure x_cp acting on a non-spherical particle does not coincide with the center of mass of the particle x_cm, a hydrodynamic pitching torque (also known as offset torque) results and acts around the axis perpendicular to the plane of relative fluid velocity and particle orientation vector. The pitching torque can change the angle of incidence of the particle. The expression used for calculating the pitching torque is given by:(38)
$T = Δ x \times F$

where, $Δ x$ is the distance between the center of pressure and the center of mass of the particle and is given by the expression:(39)
$Δ x = \frac{L}{4} (1 - {(\sin α)}^{3})$

Where, $L$ is the length of the particle and $α$ is the angle of inclination of the particle. $F$ is the total aerodynamic force (sum of drag and lift force) acting on the particle.
Rotational torque
A particle experiences rotational torque, also known as rolling friction torque, when there is a difference between the local fluid rotation and the angular velocity of the particle. The rotational torque is applied at the center of mass of the particle and is given by the expression:(40)
$T = \frac{ρ_{f}}{2} {(\frac{d_{p}}{2})}^{5} c_{r} |\frac{1}{2} \nabla \times v_{f} - ω_{p}| (\frac{1}{2} \nabla \times v_{f} - ω_{p})$
(41)
$c_{r} = \{\begin{matrix} \frac{64 π}{{Re}_{r}} {Re}_{r} \leq 32 \\ \frac{12.9}{{Re}_{r}^{0.5}} + \frac{128.4}{{Re}_{r}} 32 ≪ {Re}_{r} ≪ 1000 \end{matrix}$
(42)
${Re}_{r} = \frac{ρ_{f} d_{p}^{2} |0.5 \nabla \times v_{f} - ω_{p}|}{μ}$

Here, $v_{f}$ is the velocity of the fluid and $ω_{p}$ is the angular velocity of the particle.
Pitching rotational torque
When the pitching rotational torque model is selected, both pitching torque and rotational torque are applied on the particle.

Heat Transfer Governing Equations

When heat transfer is active, in addition to the momentum equation, the energy conservation equations for the fluid and the particle are solved simultaneously to obtain the temperature of each phase. The governing equations for obtaining the temperatures of the particles and fluid are described below:(43)

\frac{\partial}{\partial t} (ρ ε_{f} C_{p f} T_{f}) + \nabla . (v_{f} ρ_{f} ε_{f} C_{p f} T_{f}) = \nabla . (ε_{f} k_{f} \nabla T_{f}) + Q_{p}

Where

Q_{p}

represents the heat source term resulting from the heat transfer from the particle.(44)

Q_{p} = \sum_{i = 0}^{N_{p}} h_{f p} A_{p} (T_{p} - T_{f})

A_{p}

is the surface area of the particle,

T_{p}

is the particle temperature and

T_{f}

is the fluid temperature. The heat transfer coefficient (

h_{f p}

) is calculated using the empirical correlation given by

{Nu}_{p} = \frac{h_{f p} d_{p}}{k_{f}}

where the Nusselt number is calculated using the following expression:(45)

{Nu}_{p} = (7 - 10 ε_{f} + 5 ε_{f}^{2}) [1 + 0.7 {Re}_{p}^{0.2} {Pr}^{0.33}] + (1.33 - 2.40 ε_{f} + 1.20 ε_{f}^{2}) {Re}_{p}^{0.7} {Pr}^{0.33}

Where $Pr$ is the Prandtl number given by $Pr = \frac{μ_{f} C_{p f}}{k_{f}}$ .