Contact Force

Hertz-Mindlin No-Slip Contact Model

The Hertz-Mindlin contact model is a variant of the non-linear spring-dashpot contact model that is based on the Hertz-Mindlin contact theory [726], [722]. The forces between two spheres, A and B, are described by the following set of equations:

where $F_n$ and $F_t$ are the magnitudes of the normal and tangential components, respectively.

The normal direction is defined by the following:

• Normal force:

  $$F_n={-K}_n{d}_n-N_n{v}_n$$

  (3251)
• Normal spring stiffness:

  $$K_n=\frac{4}{3}{E}_{eq}\sqrt{d_n{R}_{eq}}$$

  (3252)
• Normal damping:

  $$N_n=\sqrt{\left({5K}_n{M}_{eq}\right)}{N}_{n\text{ damp}}$$

  (3253)

  where $N_{n\text{ damp}}$ is the normal damping coefficient as given in Eqn. (3257).

The tangential direction is defined by the following:

• The tangential force is defined by $-K_t{d}_t-N_t{v}_t$ if $\left|K_t{d}_t\right|<\left|K_n{d}_n\right|C_{fs}$ , where $C_{fs}$ is a static friction coefficient. Otherwise:

  $$F_t=\frac{\left|K_n{d}_n\right|C_{fs}{d}_t}{\left|d_t\right|}$$

  (3254)
• Tangential spring stiffness:

  $$K_t=8G_{eq}\sqrt{d_n{R}_{eq}}$$

  (3255)
• Tangential damping:

  $$N_t=\sqrt{\left({5K}_t{M}_{eq}\right)}{N}_{t\text{ damp}}$$

  (3256)

  where $N_{t\text{ damp}}$ is the tangential damping coefficient as given in Eqn. (3258).

Here, $C_{n\text{ rest}}$ and $C_{t\text{ rest}}$ are the normal and tangential coefficients of restitution, model parameters set by the user. If $C_{n\text{ rest}}=0$, then $N_{n\text{ damp}}=1$, and if $C_{t\text{ rest}}=0$, then $N_{t\text{ damp}}=1$.

The equivalent radius is defined as:

The equivalent particle mass is:

The equivalent Young’s modulus is expressed as:

The equivalent shear modulus is:

where:

• $M_A$ and $M_B$ are masses of sphere A and sphere B.
• $d_n$ and $d_t$ are overlaps in the normal and tangential directions at the contact point.
• $R_A$ and $R_B$ are the radii of the spheres.
• $E_A$ and $E_B$ are the Young’s modulus of the spheres.
• ${\nu }_A$ and ${\nu }_B$ are the Poisson's ratios.
• ${\nu }_n$ and ${\nu }_t$ are the normal and tangential velocity components of the relative sphere surface velocity at the contact point.

For particle-wall collisions, the formulas stay the same, but the wall radius and mass are assumed to be $R_{wall}=\infty$ and $M_{wall}=\infty$ , so the equivalent radius is reduced to $R_{eq}=R_{particle}$ and equivalent mass to $M_{wall}=M_{particle}$ .

Di Renzo and Di Maio [722] proposed several formulations for detailed treatment of tangential micro-slip and tangential force computation. Simcenter STAR-CCM+ uses the formulation that Tsuji proposed [737], where the tangential force is assumed to be non-linear but the detail micro-slip tracking is replaced by analytical expression. The resulting code is computationally efficient while still matching experimental data.

Critical Velocity

When the particle-wall impact velocity exceeds a critical velocity, the particle has enough energy to pass through a wall boundary. Simcenter STAR-CCM+ calculates the critical velocity as:

$$v_{\text{crit}}=\sqrt{\frac{4}{5}\frac{E}{\pi \rho (1-{\nu }^2)}}$$

(3263)

where $E$ is the Young's modulus, $\rho$ is the density, and $\nu$ is the Poisson ratio.

Linear Spring Contact Model

This model is based on the work of Cundall and Strack [720]. The normal and tangential forces ${\mathbf{F}}_n$ and ${\mathbf{F}}_t$ are defined as they are in Eqn. (3251) and Eqn. (3254), but parameters $K_n$ , $K_t$ , $N_n$ , and $N_t$ have the following definitions:

• $K_n$ is the normal spring constant.
• $K_t$ is the tangential spring constant.
• $N_n=2N_{n\mathrm{d}\mathrm{a}\mathrm{m}\mathrm{p}}\sqrt{K_n{M}_{eq}}$ where:
  • $N_{n\mathrm{d}\mathrm{a}\mathrm{m}\mathrm{p}}$ is the normal damping coefficient as in Eqn. (3257).
  • $M_{eq}$ is the equivalent particle mass as in Eqn. (3260).
• $N_t=2N_{t\mathrm{d}\mathrm{a}\mathrm{m}\mathrm{p}}\sqrt{K_t{M}_{eq}}$ where $N_{t\mathrm{d}\mathrm{a}\mathrm{m}\mathrm{p}}$ is the tangential damping coefficient as in Eqn. (3258).

When the Linear Spring model uses the particle material based method, the spring stiffness is estimated based on the assumption that the linear model is a linearized Hertz model, in which the maximum normal contact force is:

where:

• ${\delta }_{\mathit{max}\mathit{}}$ is the maximum overlap.

• $E_{eq}$ is the equivalent Young’s modulus defined in Eqn. (3261).

• $R_{eq}$ is the equivalent radius defined in Eqn. (3259).

To calculate the spring constants, let:

This gives the normal spring constant $K_n$ as:

The tangential spring constant $K_t$ is:

where $G_{eq}$ is the equivalent shear modulus defined in Eqn. (3262).

Walton Braun Hysteretic Contact Model

The Walton Braun hysteretic plastic-elastic contact system is characterized by the following relations:

when the contact is loading, and:

during unloading, where:

$E$ is the material Young modulus, $\delta$ is the contact normal overlap, $\mathbf{\text{n}}$ is the contact normal vector, and Yield Stress Fraction is a user-controllable model property defining the onset of plastic deformation. $E_f$ is a user-controllable energy fraction that defines the amount of energy that is recovered during the unloading and ${\delta }_{\text{max}}$ is the maximum overlap that is reached during the contact loading.

Contact Cycle

When particle contact is formed, the state of the contact is tracked at each time-step as follows:

• Contact is considered loading when ${\mathbf{\text{v}}}_r\cdot \mathbf{\text{n}}\ge 0$ and $\delta >{\delta }_{\text{max}}$ , where ${\mathbf{\text{v}}}_r$ is relative surface velocity.

  When contact is loading, the force is calculated using $K_1$ stiffness and ${\delta }_{\text{max}}$ is updated to a value of $\delta$ .

• If the contact is not loading, it is either unloading ( ${\mathbf{\text{v}}}_r\cdot \mathbf{\text{n}}<0$ ) or reloading ( $\delta \le {\delta }_{\text{max}}$ ). During unloading, the $K_2$ stiffness is used to compute the normal force and the value of ${\delta }_{\text{max}}$ is not updated.
• If the particles separate completely, the contact tracking information is removed and subsequent contact is treated as new.

Geometric Interpretation

The loading/unloading cycle can be geometrically represented as follows: