4.4. Force Field

The force field encompasses a broader set of operators and mechanisms responsible for computing forces acting on particles within the particle environment. It includes various types of forces such as gravitational forces, contact forces, and other external influences that affect particle dynamics.

4.4.1. Contact Force Laws

Contact's Law in the context of the Discrete Element Method (DEM) refers to the principle used to calculate forces between particles based on their relative displacements. In DEM simulations, Contact's Law is applied to model interactions between particles, enabling the simulation of elastic deformation and linear force behaviors within particle-based systems.

We distinguish two kinds of interaction :

  • pure contact interaction

  • cohesive interaction

Here after are given the main laws available in exaDEM :

Name

type

Description

hooke

pure contact

Default configuration : elastic linear normal force, normal viscosity force, Coulomb friction tangent force and rolling resistance

cohesive

Cohesive

Addition of a normal cohesive force depending on the distance between particules

dmt

Cohesive

Addition of an adhesive law at contact to model Van der Waals force for hard particles

The variables required to describe interactions are:

Variable

Description

\(cp\)

The contact position

\(r_i\)

The position of the particle i

\(r_j\)

The position of the particle j

\(v_i\)

The velocity of the particle i

\(v_j\)

The velocity of the particle j

\(vrot_i\)

The angular velocity of the particle i

\(vrot_j\)

The angular velocity of the particle j

\(m_i\)

The mass of the particle i

\(m_j\)

The mass of the particle j

\(\delta_n\)

The interpenetration / particle overlap

Each kind of interaction also requires additionnal constants, reflecting for most of them mechanical properties:

Constant

Description

\(\Delta_t\)

The timestep increment

\(\alpha_n\)

The damping rate

\(k_n\)

The normal stiffness coefficient

\(k_t\)

The tangential stiffness coefficient

\(v_t\)

The relative tangential velocity

\(k_r\)

The rotational stiffness coefficient

\(\mu\)

The coefficient of friction

4.4.1.1. Elastic linear force with friction tangent force (hooke)

In the Discrete Element Method (DEM), the equations of motion (translations and rotations) are discretized in time. Only the rigid-body displacements are considered. Small overlaps between the particles are allowed and used as strain variables. The total contact force between particle \(i\) and particle \(j\) is given by

\[\textbf{f}_{ij} = f_n \textbf{n} + \textbf{f}_t\]

where \(f_n\) the normal component of the contact force and \(\textbf{f}_t\) is the tangential force vector. These forces are expressed in the local contact frame \((\textbf{n},\textbf{t},\textbf{s})\) as a function of the overlaps and tangential displacements. They are calculated from force laws that generally describe frictional contact interactions. An important feature of DEM is to allow the particles to overlap. This overlap \(\delta_n\) represents a normal strain localized in the vicinity of the contact point. A simple linear relation is assumed between normal contact force and \(\delta_n\) . This is consistent with the fact that the overlaps allow for a penalty-based explicit formulation of particle motions, i.e., the elastic repulsion force is mobilized to prevent penalizing the overlap. The condition of particle undeformability implies that the overlaps must stay small compared to particle size. In this linear approximation, the normal component of the contact force is given by

\[f_n = - k_n \delta_n + \nu_n v_n\]

where \(k_n\) is the normal stiffness coefficient, \(v_n\) the normal component of the relative velocity (between particle \(i\) and particle \(j\)), and \(\nu_n\) is the viscous damping coefficient. \(\nu_n\) is related to the restitution coefficient \(e_n\) by

\[ \begin{align}\begin{aligned}\nu_n = \alpha_n \sqrt{2 m_{\text{eff}} v_n}\\\alpha_n = \frac{- \ln{e_n}}{\sqrt{\ln^2{e_n} + \pi^2}}\end{aligned}\end{align} \]

where:

  • \(\nu_n\) is the viscous damping rate during the collision

  • \(\alpha_n\) is the damping parameter, calculated from the restitution coefficient \(e_n\)

  • \(m_{\text{eff}}\) is the effective mass of the two colliding particles, defined by \(m_{\text{eff}} = \frac{m_i m_j}{m_i + m_j}\)

Note

The formulas are identical to those used in Rockable (see Rockable Force Laws) but are implemented differently to align with the exaDEM data structure.

The tangential force represents the frictional resistance between particles when they slide against each other. This force is calculated based on the relative tangential velocity (\(v_t\)) and a tangential stiffness parameter (\(k_t\), keyword ktContact).

The Coulomb friction model is used to limit the tangential force, ensuring that it does not exceed the product of the friction coefficient \(\mu\) and the normal force (\(f_n\)):

\[f_t \leq \mu \cdot f_n\]

The tangential force \(f_t\) is incrementally updated at each time step according to the following relation, starting from the onset of contact:

\[f_t = f_t + k_t \cdot v_t \cdot \Delta t\]

where: - \(f_t\) is the tangential force, - \(k_t\) is the tangential stiffness, - \(v_t\) is the relative tangential velocity, and - \(\delta t\) is the time step.

The friction force is reset to zero as soon as contact is lost.

4.4.1.2. cohesive normal law

This interaction requires two more parameters that are :

Constant

Description

\(fc\)

Cohesive force threshold

\(dncut\)

Distance cutoff for cohesive interaction

Three possibilities depending on the value of the interpenetration between \(\delta_n\) two particles:

  • \( \delta_n < -dncut \)

  • \( -dncut < \delta_n < 0 \)

  • \( 0 < \delta_n < dncut \)

Warning

Cohesive forces (dncut) are only applied if you use the operators contact_[*]_[*]_[*]_cohesive. Otherwise we only consider the case \(\\d_n < 0.0\).

Formula between particle i and particle j if \( \delta_n < -dncut \) :

\[\textbf{f}_{ij} = -k_n . \delta_n + \alpha_n \sqrt{2.m_{eff}} v_n + k_t . v_t . \Delta_t\]

with the effective mass:

\[m_{eff} = \frac{m_i.m_j}{m_i+m_j}\]

and the relative velocity norm:

\[v_n = (v_i - (cp - r_i) \wedge vrot_i) - (v_j - (cp - r_j) \wedge vrot_j)\]

Formula between particle i and particle j if \( -dncut < \delta_n < 0 \) :

\[\textbf{f}_{ij} = f_n + k_t . v_t . \Delta_t\]

with:

\[\begin{split}f_n =\left \{ \begin{array}{lcl} -k_n . \delta_n + \alpha_n \sqrt{2.m_{eff}} v_n & if & -k_n . \delta_n + \alpha_n \sqrt{2.m_{eff}} v_n >= -fc \\ & & \\ -fc & if & -k_n . \delta_n + \alpha_n \sqrt{2.m_{eff}} v_n < -fc \end{array} \right.\end{split}\]

Formula between particle i and particle j if \( 0 < \delta_n < dncut \) :

\[\textbf{f}_{ij} = (\frac{fc}{dncut} . \delta_n - fc) . \textbf{n}\]

with n normalized vector from particle i to particle j

4.4.1.3. dmt adhesive normal law

DMT (Derjaguin–Muller–Toporov) is usually used to compute pull-off forces (force needed to split two rigid objects in contact). It is computed with global energy calculation. To define it, it needs an additional parameter representing an energy :

Constant

Description

\(\gamma\)

Adhesion energy per unit of surface

In case of contact, an additional force \(f_{DMT}\) is added to the normal force such that :

\[f_n = f_n - f_{DMT}\]

with:

\[f_{DMT} = 2 \pi R_{eff} \gamma\]

where the effective radius \(R_{eff}\) is obtained from the respective radii \(R_{i}\) and \(R_{j}\) of the particles i and j:

\[R_{eff} = \frac{R_i.R_j}{R_i+R_j}\]

and \(\gamma\) is the surface energy (linked to Van der Waals forces).

Note

Since DMT force law is added to the normal force, it requires the definition of a pure contact law, such as hooke one.

4.4.2. Contact Law Operators

Contact Law Operators define the type of contact law and the parameters associated with.

The operator naming convention follows the rule below:

contact_[material_mode]_[particle_type]_[pure_contact_law]_[cohesion_law]

Where each field is defined as follows:

Operator name fields

Field

Possible values

Description

[material_mode]

singlemat, multimat

Material interaction mode

[particle_type]

sphere, polyhedron

Type of particles

[pure_contact_law]

hooke

Pure contact law

[cohesion_law]

cohesive, dmt, none

Cohesion / adhesion law

Note

For convenience, simplified operator names are available. They are strict equivalents of the full naming convention.

Sphere operators

  • contact_sphere is equivalent to contact_singlemat_sphere_hooke or contact_singlemat_sphere_hooke_none

  • contact_multimat_sphere is equivalent to contact_multimat_sphere_hooke_none

Polyhedron operators

  • contact_polyhedron is equivalent to contact_singlemat_polyhedron_hooke or contact_singlemat_polyhedron_hooke_none

  • contact_multimat_polyhedron is equivalent to contact_multimat_polyhedron_hooke_none

In a contact operator, the following parameters can be defined:

symetric

Activate or disable symmetric updates (do not disable it with polyhedron).

config

Data structure that contains contact force parameters (dncut, kn, kt, kr, fc, mu, damp_rate). Type = exaDEM::ContactParams. No default parameter.

config_driver

Data structure that contains contact force parameters (dncut, kn, kt, kr, fc, mu, damp_rate). Type = exaDEM::ContactParams. This parameter is optional.

save_interactions

Store interactions into the classifier data structure. Default is false.

Here are 4 examples with YAML:

- contact_sphere:
   symetric: true
   config: { kn: 100000, kt: 100000, kr: 0.1, mu: 0.9, damp_rate: 0.9}

- contact_polyhedron:
   config: { kn: 10000, kt: 10000, kr: 0.1, mu: 0.1, damp_rate: 0.9}
   config_driver: { kn: 10000, kt: 10000, kr: 0.1, mu: 0.3, damp_rate: 0.9}

- contact_singlemat_sphere_hooke_cohesive:
   symetric: true
   config: { dncut: 0.1 m, kn: 100000, kt: 100000, kr: 0.1, fc: 0.05, mu: 0.9, damp_rate: 0.9}

- contact_singlemat_polyhedron_hooke_cohesive:
   config: { dncut: 0.1 m, kn: 10000, kt: 10000, kr: 0.1, fc: 0.05, mu: 0.1, damp_rate: 0.9}
   config_driver: { dncut: 0.1 m, kn: 10000, kt: 10000, kr: 0.1, fc: 0.05, mu: 0.3, damp_rate: 0.9}

Note

It is important to check that interaction lists have been built with this option enabled. By default, exaDEM always builds interaction lists using the symmetry option to limit the number of calculations.

Note

Contact Force With Cohesion ([cohesion_law] = cohesive) operator includes a cohesion force from rcut to rcut+dncut with the cohesion force parameter fc.

Note

  • The contact_[*]_sphere_[*]_[*] operators are designed to process interactions built in nbh_sphere (please, include the config_spheres.msp file).

  • The contact_[*]_polyhedron_[*]_[*] operators are designed to process interactions built in nbh_polyhedron (please, include the config_polyhedra.msp file).

4.4.3. Multi-Material

In the previous section, the contact law used the same parameters for all interactions ([material_mode] = singlemat). It is also possible to specify the contact law depending on the type of interaction ([material_mode] = multimat) . In this section, we introduce how to define the values of the contact law between particles, as well as between particles and drivers.

Note

Unlike in some other DEM codes, the coefficients are not derived from the material properties (such as Poisson’s ratio and Young’s modulus).

To handle multiple particle types, you must use either the contact_multimat_[*]_[*]_[*] operators as illustrated below.

  • YAML example for polyhedra:

compute_force:
  - gravity_force
  - contact_multimat_polyhedron

  • YAML example for spheres:

compute_force:
  - gravity_force
  - contact_multimat_sphere:
     symetric: true

The following examples illustrate the definition of contact parameters for two particle types (Type1, Type2) and a driver identified by 0.

  • Operator Name: multimat_contact_params

  • Description: This operator defines the contact law parameters between different particle types.

Parameter

Description

mat1

List of the first particle type(s).

mat2

List of the second particle type(s).

kn

Normal force coefficient for the specified interaction type.

kt

Tangential force coefficient for the specified interaction type.

kr

Rolling resistance coefficient for the specified interaction type.

mu

Friction coefficient for the specified interaction type.

damprate

Damping rate coefficient for the specified interaction type.

default_config

Applies the same parameter set to all undefined interaction configurations.

YAML example:

- multimat_contact_params:
   mat1:      [  Type1, Type1, Type2 ]
   mat2:      [  Type1, Type2, Type2 ]
   kn:        [   5000, 10000, 15000 ]
   kt:        [   4000,  8000, 12000 ]
   kr:        [    0.0,   0.0,   0.0 ]
   mu:        [    0.1,   0.2,   0.3 ]
   damprate:  [  0.999, 0.999, 0.999 ]

With default_config:

- multimat_contact_params:
   mat1:      [  Type1 ]
   mat2:      [  Type1 ]
   kn:        [   5000 ]
   kt:        [   4000 ]
   kr:        [    0.0 ]
   mu:        [    0.1 ]
   damprate:  [  0.999 ]
   default_config: { kn: 15000, kt: 12000, kr: 0.0, mu: 0.3, damp_rate: 0.999 }

A complete example is available (please report if the link does not work): rotating-multimat.msp

  • Operator Name: drivers_contact_params

  • Description: This operator defines the contact law parameters between particles and drivers.

  • Parameters:

Parameter

Description

mat

List of particle type(s) concerned by the interaction.

driver_id

Identifier(s) of the driver(s) interacting with the particles.

kn

Normal force coefficient for the specified interaction type.

kt

Tangential force coefficient for the specified interaction type.

kr

Rolling resistance coefficient for the specified interaction type.

mu

Friction coefficient for the specified interaction type.

damprate

Damping rate coefficient for the specified interaction type.

default_config

Applies the same parameter set to all undefined interaction configurations.

 
- drivers_contact_params:
   mat:       [  Type1, Type2 ]
   driver_id: [      0,     0 ]
   kn:        [  10000, 15000 ]
   kt:        [   8000, 12000 ]
   kr:        [    0.0,   0.1 ]
   mu:        [    0.5,   0.5 ]
   damprate:  [  0.999, 0.999 ]

A complete example is available (please report if the link does not work): rotating-multimat.msp

4.4.4. External Forces

External forces are additional influences acting on particles within the simulation environment, originating from sources outside the particle system itself. These forces can include environmental factors like wind, fluid flow, or magnetic fields, as well as user-defined forces applied to specific particles or regions.

Formula:

(1)\[\textbf{f} = m.\textbf{g}\]

With f the forces, m the particle mass, and g the gravity constant.

  • Operator Name: gravity_force

  • Description: This operator computes forces related to the gravity.

  • Parameter:

gravity

Define the gravity constant in function of the gravity axis, default value are x axis = 0, y axis = 0 and z axis = -9.807

YAML example:

- gravity_force:
   gravity: [0,0,-0.009807]

Formula:

\[\textbf{f} = -\mu.cx.\|v\|.\textbf{v}\]

With f the particle forces, cx the aerodynamic coefficient, and \(\mu\) the drag coefficient, ||v|| the norm of the particle velocity, and v the particle velocity.

  • Operator Name: quadratic_force

  • Description: External forces that model air or fluid, f = - mu * cx * norm(v) * vector(v).

  • Parameter:

cx

aerodynamic coefficient, default value is for air = 0.38.

mu

drag coefficient. default value is for air = 0.000015.

YAML example: see example quadratic-force-test/QuadraticForceInput.msp

- quadratic_force:
   cx: 0.38
   mu: 0.0000015

  • Operator Name: sphere_fluid_friction

  • Description: External forces that model a fluid computed from a grid such as: f = ||fv - pv||

\[dv = fv - pv\]
\[f = cx . dv . ||dv|| . \pi . r . r.\]

With fv the fluid velocity, pv the particle velocity, r the particle radius, and cx a coefficient set to 1 by default.

Note

The fluid velocity fv for each point of the grid has been defined by the operator set_cell_values (pure exaNBody operator).

4.4.5. Inner Bond Forces

Some words:

  • Interaction type = 13

  • Parameters: kn, kt, damprate, g