Discrete Element Method

The discrete element method (DEM) is designed for modeling the granular flow of materials such as sand, food particles, powders, capsules, and slurries. These flows are characterized by a high particle density where the inter-particle interaction is of importance.

Established by Cundall and Strack [720], the DEM model extends the Lagrangian formulation to account for inter-particle interaction in the particle equations of motion. These forces cannot be ignored for highly loaded flows, that is, many interacting particles.

The DEM particles can assume different shapes and volumes. Simcenter STAR-CCM+ offers several options for simulating non-spherical particles (see DEM Particle Types) including using arbitrary convex or concave shapes defined by CAD geometry. The collision between non-spherical particles and other particles or walls can have more than one force-transmitting contact point.

Simcenter STAR-CCM+ models DEM particles based on soft-particle formulation in which particles are allowed to develop an overlap. The calculated contact force is proportional to the overlap, as well as to the particle material and geometric properties. This contact force enters the Lagrangian momentum equation Eqn. (2958) as $F_{c}$ .

Discrete Element Method Particle Time Scale

To obtain the DEM particle trajectory, Eqn. (2958) is integrated over time using a characteristic DEM particle time scale.

The particle time scale, and therefore the time-step, depends on one of three times associated with the particle:

Rayleigh wave propagation time $τ_{1}$

The maximum time-step that is allowed for a DEM particle is constrained by the assumption that the force acting on a particle is only affected by the particle’s immediate neighbors during a single time-step duration. The time-step is therefore limited by the time it takes the Rayleigh wave to propagate across the surface of the particle to the opposite side [726] (from pole to pole in the case of a sphere):

Figure 1. EQUATION_DISPLAY

τ_{1} = π \frac{R}{V_{R a y l e i g h}}

(3244)

where $R$ is half the characteristic minimum dimension of the particle construed for its moment of inertia, (the radius in the case of a sphere).

The Rayleigh wave velocity depends on material properties, and the exact value is obtained as a solution to a secular equation [731]. Simcenter STAR-CCM+ uses an approximation of the solution that provides sufficient precision without incurring high computational cost [728], [729].

When the Linear Spring contact model is used,

Figure 2. EQUATION_DISPLAY

τ_{1} = 2 \sqrt{m / k}

(3245)

where $m$ is the particle mass and $k$ is the spring stiffness.

The maximum time-step must be no more than $t_{scale} τ_{1}$ , where $t_{scale}$ is a user-defined time-scale.

Impact duration $τ_{2}$

The second time-step limiting criterion applied to moving particles is the duration of impact of two perfectly elastic particles. For two identical spheres of radius

R

, this duration, assuming the Hertz contact theory, was derived by Timoshenko [735] to be:

Figure 3. EQUATION_DISPLAY

τ_{2} = 2.94 {(\frac{5 \sqrt{2} π ρ}{4} \frac{1 - ν^{2}}{E})}^{\frac{2}{5}} \frac{R}{\sqrt[5]{ν_{i m p a c t}}}

(3246)

The maximum time-step must be no more than $0.1 τ_{2}$ .

Particle transit time $τ_{3}$

The final restriction on a DEM particle time-step is geometric. It is based on the assumption that particles must not move too far within the time-step. This condition prevents missing contacts between DEM particles, as well as particles and walls. Therefore each particle is constrained such that it takes at least 10 time-steps for the particle to move the full length of the radius.

Figure 4. EQUATION_DISPLAY

τ_{3} = (\frac{R}{ν_{p a r t i c l e}})

(3247)

where $R$ is half the characteristic minimum dimension of the particle, (the radius in the case of a sphere).

The maximum time-step must be no more than $0.1 τ_{3}$ .

The final particle time-step is determined as a minimum of

t_{scale} τ_{1}, 0.1 τ_{2}, 0.1 τ_{3}

. In practice

τ_{1}

is typically the limiting factor while

τ_{2}

and

τ_{3}

only constrain particles that are moving fast, or the Young’s modulus of the material is set low in order to accelerate the simulation.

Artificial Viscosity

The artificial viscosity model is used to correct the results in cases where high impact velocities cause the numerical models to over-predict the overlap between the colliding objects. The energy that is dissipated by artificial viscosity represents unresolved processes of damage, deformation, and entropy.

The name "Artificial Viscosity" refers only to the similarity between the formulation of the Artificial Viscosity model and the formulations for viscosity. It does not refer to any physical effect resembling real viscosity.

This model introduces a force during compression (while the particles are overlapping and their centers are moving nearer) based on the von Neumann-Richtmyer artificial viscosity [738]:

Figure 5. EQUATION_DISPLAY

F_{v} = - f (Q_{2} | v_{n} | - Q_{1} c) \frac{m_{e q}}{R_{e q}} v_{n}

(3248)

where:

$F_{v}$ is the artificial viscosity force.
$v_{n}$ is the component of the velocity in the contact normal direction.
$m_{e q}$ is the equivalent (harmonic mean) mass of the particles.
$R_{e q}$ is the equivalent (harmonic mean) radius of the particles.
$c$ is the speed of sound in the material. It is estimated from the equivalent mass and radius, and from the contact stiffness for the specific contact model in use.
$Q_{1}$ is the linear coefficient.
$Q_{2}$ is the quadratic coefficient.

The factor $f$ is a linear ramp function that is scaled to the particle sizes as:

Figure 6. EQUATION_DISPLAY

f = {\begin{matrix} 0, & r / R_{e q} < p_{lower} \\ \frac{r / R_{e q} - p_{lower}}{p_{upper} - p_{lower}}, & p_{upper} \geq r / R_{e q} \geq p_{lower} \\ 1, & r / R_{e q} > p_{upper} \end{matrix}

(3249)

where:

$r$ is the penetration, the overlap distance between particles or between particle and wall.
$p_{upper}$ and $p_{lower}$ are the upper and lower ramp limits. They are defined in terms of the ratio of the penetration $r$ to the equivalent particle radius $R_{e q}$ .

This ramp function $f$ allows you to use artificial viscosity without affecting typical low penetration contacts.