Discrete Element Method

The Discrete Element Method (DEM) is a computational technique widely used in physics, engineering, and applied mathematics for simulating the mechanical behavior of granular materials, such as powders, sand, soil, or rock, as well as other discontinua. Unlike continuum mechanics that treats materials as continuous, DEM considers individual particles or elements and their interactions. This approach provides detailed insights into the micro-mechanical behavior of materials, making it particularly valuable in fields such as geomechanics, material science, and mechanical engineering.

Fundamental Principles

The core idea behind DEM is the discretization of a material system into a finite set of distinct, interacting mass elements (particles). These elements (particles) can vary in shape, size, and properties, and they interact with each other and possibly with their boundaries through contact forces and potential fields. The motion and behavior of each mass element are governed by Newton's laws of motion, accounting for the forces and moments acting upon them.

API

TrixiParticles.DEMSystem — Type

DEMSystem(initial_condition, contact_model; damping_coefficient=0.0001,
          acceleration=ntuple(_ -> 0.0, ndims(initial_condition)), source_terms=nothing,
          radius=nothing)

Constructs a Discrete Element Method (DEM) system for numerically simulating the dynamics of granular and particulate matter. DEM is employed to simulate and analyze the motion, interactions, and collective behavior of assemblies of discrete, solid particles, typically under mechanical loading. The model accounts for individual particle characteristics and implements interaction laws that govern contact forces (normal and tangential), based on specified material properties and contact mechanics.

Arguments

initial_condition: Initial condition of the system, encapsulating the initial positions, velocities, masses, and radii of particles.
contact_model: Contact model used for particle interactions.

Keywords

acceleration: Global acceleration vector applied to the system, such as gravity. Specified as an SVector of length NDIMS, with a default of zero in each dimension.
source_terms: Optional; additional forces or modifications to particle dynamics not captured by standard DEM interactions, such as electromagnetic forces or user-defined perturbations.
damping_coefficient=0.0001: Set a damping coefficient for the collision interactions.
radius=nothing: Specifies the radius of the particles, defaults to initial_condition.particle_spacing / 2.

Experimental Implementation

This is an experimental feature and may change in a future releases.

References

[33], [34], [35]

source

Contact Models

TrixiParticles.HertzContactModel — Type

HertzContactModel(; elastic_modulus, poissons_ratio)

Non-linear contact model based on Hertzian contact theory ([35]).

This model calculates the normal contact force between two spherical particles (or a particle and a boundary represented by an equivalent sphere) based on their material properties and the overlap $\delta$. The elastic part of the force is given by:

\[F_{\text{elastic}} = \frac{4}{3} E^* \sqrt{R^*} \delta^{3/2},\]

where $E^*$ is the effective Young's modulus and $R^*$ is the effective radius.

The effective Young's modulus $E^*$ is calculated from the Young's moduli $E_a, E_b$ and Poisson's ratios $\nu_a, \nu_b$ of the two contacting bodies:

\[E^* = \left( \frac{1 - \nu_a^2}{E_a} + \frac{1 - \nu_b^2}{E_b} \right)^{-1}.\]

The effective radius $R^*$ is calculated from the radii of the two particles $R_a$ and $R_b$:

\[R^* = \left( \frac{1}{R_a} + \frac{1}{R_b} \right)^{-1} = \frac{R_a R_b}{R_a + R_b}.\]

For particle-wall interactions, $R_b \to \infty$, so $R^* = R_a$.

The implementation also includes a damping force based on the approach described in [35], proportional to the normal component of the relative velocity $v_{\text{rel,n}}$:

\[F_{\text{damping}} = C_{\text{damp}} \gamma_c v_{\text{rel,n}},\]

where $C_{\text{damp}}$ is the user-provided damping coefficient (damping ratio), and $\gamma_c$ is a non-linear critical damping coefficient:

\[\gamma_c = 2 \sqrt{m^* K_{\text{nonlin}}}\]

with $m^*$ being the effective mass and $K_{\text{nonlin}}$ being a non-linear stiffness term related to the current state:

\[K_{\text{nonlin}} = \frac{F_{\text{elastic}}}{\delta} = \frac{4}{3} E^* \sqrt{R^* \delta}.\]

The total normal force is $F_n = F_{\text{elastic}} + F_{\text{damping}}$.

Fields

elastic_modulus::Float64: Material Young's modulus $E$.
poissons_ratio::Float64: Material Poisson's ratio $\nu$.

source

TrixiParticles.LinearContactModel — Type

LinearContactModel(; normal_stiffness)

Linear spring-dashpot contact model ([34]).

This model calculates the normal contact force between two objects based on a linear spring law for the elastic component and a linear viscous damping law for the dissipative component. The total normal force $F_n$ is given by

\[F_n = k_n \delta + \gamma_d v_{\text{rel,n}},\]

where $k_n$ is the normal stiffness, $\delta$ is the overlap between the objects, $v_{\text{rel,n}}$ is the normal component of the relative velocity, and $\gamma_d$ is the damping coefficient.

The damping coefficient $\gamma_d$ is calculated based on the critical damping coefficient $\gamma_c$ and a user-provided damping ratio $C_{\text{damp}}$:

\[\gamma_d = C_{\text{damp}} \gamma_c,\]

where the critical damping for this linear system is

\[\gamma_c = 2 \sqrt{m^* k_n}\]

and $m^*$ is the effective mass of the colliding pair.

The total force is applied along the normal direction connecting the centers of the contacting objects.

Fields

normal_stiffness::Real: Constant spring stiffness $k_n$ for the normal direction.

source