Discrete Element Modeling of ceramics

Micromechanics using discrete element modeling

Damien André*

* Univ. Limoges, IRCER, UMR 7315, F-87000 Limoges, damien.andre@unilim.fr

The basics of DEM

In a DEM model, each grain
⇢ has a position

$\vec{p_{i}}$ , a velocity

$\vec{\dot{p}_{i}}$ and an acceleration

$\vec{\ddot{p}_{i}}$
⇢ interacts with its neighbours through contacts

Contact forces $\vec{f_{ij}}$ depends on :
⇢ interpenetration

$\vec{f_{ij}^{c}} = F(\vec{\delta_{ij}})$
⇢ relative velocities

$\vec{f_{ij}^{v}} = G\left(\vec{v_{ij}} \right)$

Every grain obeys Newton's 2nd law

$\vec{a_{i}}=\frac{\sum\vec{f_{i}}}{m_{i}}$
This is a dynamic simulation that depends on time $t$ !

Focus on contact research

If the discrete domain contains $\text{N}$ discrete elements, this simple calculation should be repeated many times...

• for a single discrete element $(i)$ , the calculation should be repeated $(\text{N}-1)$ times

• for $\text{N}$ discrete elements, it gives $\frac{\text{N}}{2}(\text{N}-1)$ calculations ! $\rightarrow$ exponential complexity $\mathcal{O}(\text{N}^2)$ 😖

Computing contact reaction force

In DEM, a very popular contact model is the Hertz-Mindlin (see LIGGGHTS documentation ) one, where :

$\vec{f_{ij}} = \underbrace{\vec{n}\left(K_{n} \delta^{n}_{ij} - \gamma_{n}v^{n}_{ij}\right)}_{\text{normal force}} + \underbrace{\vec{t}\left(K_{t} \delta^{t}_{ij} - \gamma_{t}v^{t}_{ij}\right)}_{\text{tangential force}}$

⇢

$K_{n}$ and

$K_{t}$ depend on

$E$ ,

$\nu$ ,

$\delta^{n}$ and

$R$

⇢

$\gamma_{n}$ and

$\gamma_{t}$ depends on

$e$ (coefficient of restitution),

$E$ ,

$\nu$ ,

$\delta^{n}$ and

$M$

⇢ The Hertz-Mindlin model is an extension of the standard (non linear) Hertz model

⇢ It is based on material parameters !

Step #2 : Newton

At this step, forces $\sum\vec{f}$ and torques $\sum\vec{t}$ acting on DE are known

$\longrightarrow$ linear $\ddot{p}(t+\Delta t)$ and angular $\ddot{q}(t+\Delta t)$ accelerations with Newton's 2nd law:

$\boxed{\vec{\ddot{p}} = \frac{\sum\vec{f}}{m}} \text{ and } \boxed{\ddot{q} = \dot{q} \cdot \bar{q} \cdot \dot{q} + \frac{1}{2}q \cdot \left[\bar{\bar{I}}^{-1} \left( \sum\vec{\tau} - 4\left( \bar{q} \cdot \dot{q} \right) \wedge \bar{\bar{I}}\cdot\left(\bar{q} \cdot \dot{q} \right) \right)\right]}$

About quaternions

Quaternion is a quadruplet of numbers : $q = q_0 + q_1 \cdot i + q_2 \cdot j + q_3 \cdot k$

where $i, j,k$ are imaginary unit numbers. Another representation is : $q = (q_0, \vec{q})\quad \text{with} \quad \vec{q}(q_1,q_2,q_3)$

The norm of quaternion is defined as : $\left|\left|q\right|\right|= \sqrt{q_0^{2}+\left|\left|\vec{q}\right|\right|^{2}}$

Unit polar quaternions $\left|\left|q\right|\right|=1$ are used as : $q=(\cos\theta,\:\vec{q}\cdot\sin\theta) \quad \text{where} \quad q_0 = \cos\theta \quad \text{and} \quad \left|\left|\vec{q}\right|\right| = \sin\theta$

Rotation of rigid bodies can be described with quaternions where $\theta$ is the rotation angle and $\vec{q}$ is the rotational axis

Quaternions avoid Gimbal lock (that happens with Euler angles)

Step #3 : integration

At this step, we have :
$\rightarrow$ new linear $\vec{\ddot{p}}(t+\Delta t)$ and angular $\ddot{q}(t+\Delta t)$ accelerations
$\rightarrow$ the current linear $\vec{\dot{p}}(t)$ and angular $\dot{q}(t)$ velocities
$\rightarrow$ the current linear $\vec{p}(t)$ and angular $q(t)$ positions

Current statement :

	$(t)$			$(t+\Delta t)$
acceleration	$\vec{\ddot{p}}(t)$	$\ddot{q}(t)$	✔	$\vec{\ddot{p}}(t+\Delta t)$	$\ddot{q}(t+\Delta t)$	✔
velocity	$\vec{\dot{p}}(t)$	$\dot{q}(t)$	✔	$\boxed{\vec{\dot{p}}(t+\Delta t)}$	$\boxed{\dot{q}(t+\Delta t)}$	✖
position	$\vec{p}(t)$	$q(t)$	✔	$\boxed{\vec{p}(t+\Delta t)}$	$\boxed{q(t+\Delta t)}$	✖

The new positions and velocities at

$t+\Delta t$ must be deduced from accelerations thanks to the integration scheme

Integration with velocity Verlet

Linear and angular velocities : $\begin{align} \boxed{\vec{\dot{p}}(t+\Delta t)} &= \vec{\dot{p}}(t) + \frac{\Delta t}{2} \left[ \vec{\ddot{p}}(t) + \vec{\ddot{p}}(t+\Delta t) \right] \\ \boxed{\dot{q}(t+\Delta t)} &= \dot{q}(t) + \frac{\Delta t}{2} \left[ \ddot{q}(t) + \ddot{q}(t+\Delta t) \right] \end{align}$

Linear and angular displacements : $\begin{align} \boxed{\vec{p}(t+\Delta t)} &= \vec{p}(t) + \Delta t \: \vec{\dot{p}}(t)+\frac{\Delta t^{2}}{2}\vec{\ddot{p}}(t)\\ \boxed{q(t+\Delta t)} &= q(t) + \Delta t \: \dot{q}(t)+\frac{\Delta t^{2}}{2}\ddot{q}(t) \end{align}$

$\Delta t$ is the so-called time step.

$\Delta t$ should be small, critical time step is :

$\Delta t = \frac{1}{k} min(t_c) \qquad \text{with}\quad k \geq 2\pi \qquad \text{and} \quad t_c = \sqrt{\frac{M_i}{K_i}}$

About time step $\Delta t$

What happens if $\Delta t \geq t_c$ ?

The simulation is unstable !

⇢ Main issue of DEM methods
⇢ Inherent to explicit integration scheme
⇢ Some codes implement adaptive $\Delta t$
⇢ NSCD [1] is able to solve this problem
⇢ NSCD is not so popular (versatility problem)

[1] 1988, J.J. Moreau, P.D. Panagiotopoulos, and International Centre for Mechanical Sciences, “Nonsmooth mechanics and applications”

About computational performance

The latest discrete element codes use GPU-accelerated technologies. They are capable of simulating around 10,000,000 particles.

Rotating drum simulation with blaze dem

10,000,000 particles $\rightarrow 5 \times 5 \times 5\ mm^{3}$ (particle radii about $10\ \mu m$ )
⟶ We are far from engineering scales 😓
⟶ Active research to overcome this limitation (surrogate modelling)

Granular media vs "Continuum" media

Apparent mechanical behaviour of concrete shares similarities with granular media

Concrete material under confined compression tests (confined pressure $\sigma_L$ ) (*)

➡ For example, high sensibility to the confine pressure

(*) 2002, Sfer, Domingo, Carol, Ignacio, Gettu, Ravindra, and Etse, Guillermo, “Study of the behavior of concrete under triaxial compression”

Discrete media VS continuum media

How to consider such materials ? It depends on :
$\rightarrow$ the scale (macro-meso-micro) to consider
$\rightarrow$ the phenomenon to model

From andalusite brick (*) to its micro-structure

If the considered material exhibits an heterogeneous (micro)structure with large discontinuities such as cracks and damages... why not using DEM for modeling such (pseudo)continuum media ?

(*) photo from Zhengzhou Rongsheng Kiln Refractory Co., Ltd

Modeling brittle elastic media

In the literature, many bond models exist... Which one to choose ?

At the IRCER, we work with :

$\rightarrow$ flat bond [1] available in the PFC software

$\rightarrow$ cohesive beam bond [2] available in the GranOO software

$\rightarrow$ distinct lattice spring [3] available in the GranOO software

[1] 2012, Potyondy, DO, and others, “A flat-jointed bonded-particle material for hard rock” [2] 2012, Damien André, Ivan Iordanoff, Jean-luc Charles, and Jérome Néauport, “Discrete element method to simulate continuous material by using the cohesive beam model” [3] 2019, André, Damien, Girardot, Jérémie, and Hubert, Cédric, “A novel DEM approach for modeling brittle elastic media based on distinct lattice spring model”

The cohesive beam model (CBM)

$\rightarrow$ Cohesive beams work in tension, compression, bending and torsion
$\rightarrow$ They are introduced between discrete elements

Geometrical		Material
Length	$L_{m}$	Young's modulus	$E_{m}$
Radius	$R_{m}$	Poisson's ratio	$\nu_{m}$

Analytical model

The response of the cohesive beam is defined as (force and torque):

$\mathbf{F}_1 = \begin{bmatrix} \frac{E S \Delta l}{l_{0}} & \text{tension}\newline -\frac{6 E I}{l_{0}^{2}} \left(\theta_{2z} + \theta_{1y}\right) &\text{bending} \newline -\frac{6 E I}{l_{0}^{2}} \left(\theta_{2y} + \theta_{1y}\right) &\text{bending} \end{bmatrix}$

$\mathbf{T}_1 = \begin{bmatrix} \frac{G I_0}{l_{0}} \left(\theta_{2x} - \theta_{1x}\right) & \text{twist} \newline -\frac{2 E I}{l_{0}} \left(\theta_{2y} + 2\theta_{1y}\right) & \text{bending}\newline +\frac{2 E I}{l_{0}} \left(\theta_{2z} + 2\theta_{1z}\right) & \text{bending} \end{bmatrix}$

Focus on non local criterion

A stress tensor is calculated for each discrete element (Love-Weber formula [1,2]) $\sigma_i = \frac{1}{4\Omega_i}\sum\limits_j{\mathbf{r}(ij) \otimes \mathbf{f}(ij)+\mathbf{f}(ij) \otimes {r}(ij)}$

Based on the Rankine criteria, a fracture occurs if : $max(\sigma_{\text{I}}, \sigma_{\text{II}}, \sigma_{\text{III}}) ≥ \sigma_f$
In such case, one (or more) cohesive beam that belong to the related DE is removed

Truong-Thi Nguyen validated this approach for complex loading [3]

[1] 2003, Zhou, Min, “A new look at the atomic level virial stress: on continuum-molecular system equivalence” [2] 1966, Weber, J, “Recherches concernant les contraintes intergranulaires dans les milieux pulvérulents” [3] 2019, Nguyen, Truong-Thi, André, Damien, and Huger, Marc, “Analytic laws for direct calibration of discrete element modeling of brittle elastic media using cohesive beam model”

The purpose of DLSM

DLSM was proposed in 2011 [1] → direct modeling of brittle elastic media $E, \nu, \sigma_f$

Original DLSM does not take into account rotations and works only for spherical elements and regular arrangements → not really compatible with DEM 😱

[1] 2011, Zhao, Gao-Feng, Fang, Jiannong, and Zhao, Jian, “A 3D distinct lattice spring model for elasticity and dynamic failure”

How it works ?

step #1 : computing displacements

Respecting the Voigt's hypothesis (linear disp.), local displacement fields give

$\vec{d}_i (x,y,z) = \begin{bmatrix} u(x,y,z)\newline v(x,y,z)\newline w(x,y,z) \end{bmatrix} = \begin{bmatrix} a_1 x + b_1 y + c_1 z\newline a_2 x + b_2 y + c_2 z\newline a_3 x + b_3 y + c_3 z \end{bmatrix} = \vec{a}\cdot x + \vec{b}\cdot y + \vec{c}\cdot z$

How to compute $\vec{a}, \vec{b}, \vec{c}$ from a point cloud ?

step #1 : computing displacements

Considering $O_i$ and considering the displacements of its $n$ neighbors, $\vec{a}, \vec{b}, \vec{c}$ can be determined…

if $n<3$ the system is under-determined, the solution can not be computed
if $n=3$ the system can be solved exactly by a linear system
if $n>3$ the system is over-determined, the solution is approximated thanks to least square method

step #2 : computing strains

The strain $\overline{\overline{\varepsilon}}$ tensor can be deduced from displacements $\vec{d}$ as:

$\overline{\overline{\varepsilon}} = \frac{1}{2}\left(\overline{\overline{grad}}(\vec{d}) + \overline{\overline{grad}}(\vec{d})^{T}\right) \text{ with } \vec{d}= \begin{bmatrix} a_1 x + b_1 y + c_1 z\newline a_2 x + b_2 y + c_2 z\newline a_3 x + b_3 y + c_3 z \end{bmatrix}$

So, the strain tensor at $O_i$ can be deduced from $\vec{a}, \vec{b}, \vec{c}$ with:

$\overline{\overline{\varepsilon}}_i(O_i) = \begin{bmatrix} \varepsilon_{xx} & \varepsilon_{xy} & \varepsilon_{xz}\newline \varepsilon_{xy} & \varepsilon_{yy} & \varepsilon_{yz}\newline \varepsilon_{xz} & \varepsilon_{yz} & \varepsilon_{zz} \end{bmatrix} = \frac{1}{2}\begin{bmatrix} 2 a_1&a_2+b_1&a_3+c_1\newline a_2+b_1&2 b_2&b_3 + c_2\newline a_3+c_1&b_3+c_2&2 c_3 \end{bmatrix}$

In most of case $n>3$
so

$(a_{1}, b_{1}, c_{1}, \cdots)$ should be approximated with (moving) least square method

step #3 : computing contact forces

In which interface… ? Voronoi tessellations are used.

Considering $(O_i, O_2)$ contact, reaction force is computed as $\boxed{\vec{F} = S\cdot\overline{\overline{\sigma}}(A)\cdot\vec{n}}$

$\overline{\overline{\sigma}}(A) = 2\mu\cdot\overline{\overline{\varepsilon}}(A) + \lambda\cdot\text{tr}(\overline{\overline{\varepsilon}}(A))\cdot\overline{\overline{\text{I}}} \quad \text{and} \quad \overline{\overline{\varepsilon}}(A) = \frac{\overline{\overline{\varepsilon}}(O_i) + \overline{\overline{\varepsilon}}(O_2)}{2}$

Forcing local strain information

⚠ the strain at the interface $(O_i, O_2)$ is non local. To force local information, the normal strain component is forced to be equal to the engineering strain :

$\overline{\overline{\varepsilon}}(A) = \frac{\overline{\overline{\varepsilon}}(O_{i}) + \overline{\overline{\varepsilon}}(O_{2})}{2}\quad \text{and} \quad \overline{\overline{\varepsilon}}(A)\cdot\vec{n}\cdot\vec{n} = \varepsilon_{nn} = \frac{\Delta L}{L}$

Displacements & local frame

When rigid body translations and/or rotations are expected, the displacements must be expressed in local frames. These frames are attached to discrete element.

$\vec{d}_{i|\mathcal{F}_{i}} (x,y,z) = \begin{bmatrix} u(x,y,z)\newline v(x,y,z)\newline w(x,y,z) \end{bmatrix} = \begin{bmatrix} a_1 x + b_1 y + c_1 z\newline a_2 x + b_2 y + c_2 z\newline a_3 x + b_3 y + c_3 z \end{bmatrix} = \vec{a}\cdot x + \vec{b}\cdot y + \vec{c}\cdot z$

⇒ Displacements of $O_0, O_1, \cdots, O_4$ must be expressed relatively to $\mathcal{F}_{i}$

About damping factor

Because the scheme is not strictly energy conservative, a small numerical damping is introduced into the Velocity Verlet integration

Linear Velocity Verlet

$\dot{\vec{\mathbf{p}}}(t+\Delta t) = \dot{\vec{\mathbf{p}}}(t) + \color{lime}{\varphi} \Delta t \left( \ddot{\vec{\mathbf{p}}}(t) + \ddot{\vec{\mathbf{p}}}(t+\Delta t) \right)$

Angular Velocity Verlet (with quaternion)

$\dot{\mathbf{q}}(t+\Delta t) = \dot{\mathbf{q}}(t) + \color{lime}{\varphi} \Delta t \left( \ddot{\mathbf{q}}(t) + \ddot{\mathbf{q}}(t+\Delta t) \right)$

Values of $\varphi$ higher than $1$ give stable simulations. Commonly, we use $\varphi \in [1; 1.3]$

Quantitative checking (static)

Apparent values of Young's modulus $\langle E \rangle$ and Poisson's ratio $\langle \nu \rangle$

Input		Output
$E$ (GPa)	$\nu$ (-)	$\langle E \rangle$ (GPa)	$\langle \nu \rangle$ (-)
20	0.49	20.973	0.483
20	0.4	20.164	0.394
20	0.3	20.141	0.292
20	0.2	20.117	0.192
20	0.1	20.095	0.093
20	0.0	20.078	-0.006

$\left. \begin{align} \langle \sigma_{xx} \rangle = \frac{f}{S} &= \frac{1}{2S}\left( \sum_{i\ \in\ \text{x}^{+}}^{n^{\text{x}^{+}}} \vec{f_{i}} - \sum_{i\ \in\ \text{x}^{-}}^{n^{\text{x}^{-}}} \vec{f_{i}} \right)\bullet \vec{x}\\ \langle \varepsilon_{xx} \rangle = \frac{\Delta l}{l_{0}} &= \frac{1}{l_{0}} \left( \frac{1}{n^{\text{x}^{+}}}\sum_{i\ \in\ \text{x}^{+}}^{n^{\text{x}^{+}}} \vec{d_{i}} - \frac{1}{n^{\text{x}^{-}}}\sum_{i\ \in\ \text{x}^{-}}^{n^{\text{x}^{-}}} \vec{d_{i}} \right)\bullet \vec{x} \end{align} \right. \Rightarrow \begin{array}{ll} \langle E \rangle = \frac{\langle\sigma_{xx}\rangle}{\langle\varepsilon_{xx}\rangle} \\ \langle \nu \rangle = -\frac{\langle\varepsilon_{yy}\rangle}{\langle\varepsilon_{xx}\rangle} \end{array}$

The basics

Let's consider a discrete domain

We can associate to each discrete element $(i)$ a temperature $T_{i}$

The temperature of a $(ij)$ bond is simply $T_{ij} = \frac{1}{2}\left( T_{i}+T_{j}\right)$

$\rightarrow$ Thermal expansion is managed by the bonds

Thermal expansion - uniform T°

lattice structure in a relaxing state at $\color{cyan}{T_{0}} \rightarrow l_{0}(\color{cyan}{T_{0}})$ relaxed lengths of bonds

the structure is then heated to $\color{red}{T_{1}}$

the structure expands due to the heating and keep its relaxing state

$l_{0}(\color{red}{T_{1}})$ are the new relaxed lengths of bonds

thermal expansion is taken into account by varying relaxed length of bonds $l_{0}(T)$

The formula is quite simple $\rightarrow l_{0}(T) = l_{0}^{T_{0}} + \color{lime}{\alpha}\left(T-T_{0}\right)$ (no calibration)

$\rightarrow$ Let's see anisotropic thermal expansion

Anisotropic thermal expansion

now, let's consider non-isotropic thermal expansion $\vec{\alpha}$

first, let's consider an elongation along $\vec{x}$

and secondly a contraction along $\vec{y}$

for anisotropic thermal expansion, $l_{0}(T)$ depends on the bond direction

$\rightarrow$ the spherical shape is not very well adapted

Thermal conduction

Applying the heat equation, temperature variation of the element

$i$ due to its interaction with

$j$ is

$\Delta T^{ij} = \frac{\left(T^{j}-T^{i}\right) S_{t} \color{yellow}{\lambda} } { d_{ij} \color{yellow}{c_{p} \rho_{i}} V_{i}} \Delta t$
Total temperature variation of

$i$ due its interaction with its

$N$ neighbors

$\Delta T^{i} = \sum^{N} \Delta T^{ij}$

Special attention must be paid on surface transmission $S_{t}$ calculation [1] [1] 2013, I. Terreros, I. Iordanoff, and J.L. Charles, “Simulation of continuum heat conduction using DEM domains”

$S_{t}$ is transmission surface between $i$ and $j$
$\Delta t$ is time step
$V_{i}$ is the volume of $i$

$\color{yellow}{\lambda}$ is thermal conductivity
$\color{yellow}{\rho_{i}}$ is density of $i$
$\color{yellow}{c_{p}}$ is thermal capacity

$d_{ij}$ is distance between $i$ and $j$

1st case - andalousite castable

Stress/strain curves at room temperature (from [1])

Mechanical behavior strongly depends on thermal history

→ Thermal history affects the microstructure

→ High non linear behavior $\approx$ high fracture toughness $\approx$ high thermal shock strength

Discrete model of andalousite microstructure

[1] 2007, Mahdi Ghamessi Kakroudi, “Comportement thermomécanique en traction de bétons réfractaires : influence de la nature des agrégats et de l’histoire thermique”

2nd case : alumina incl. / glass matrix

To understand these micro-macro relationships, experimental studies on model materials were conducted [1].

Only the cooling part, starting from $T^{\circ}$ , is considered by the DEM simulations [2]

[1] 2003, Nicolas Tessier-Doyen, “Etude expérimentale et numérique du comportement thermomécanique de matériaux réfractaires modèles” [2] 2017, André, Damien, Levraut, Bertrand, Tessier-Doyen, Nicolas, and Huger, Marc, “A discrete element thermo-mechanical modelling of diffuse damage induced by thermal expansion mismatch of two-phase materials”

	$E$ (GPa)	$\nu$	$\sigma_f$ (MPa)	$\alpha$ (K $^{-1}$ )
Glass	$72.0$	$0.23$	$50.0$	$11.6 \times 10^{-6}$
Alumina	$340.0$	$0.24$	$380.0$	$7.6 \times 10^{-6}$

3rd case - aluminum titanate

The aluminum titanate exhibits very high thermal expansion anisotropy

$\begin{align} \color{red}{\alpha_{a}} &= \color{red}{-2.38} & \times 10^{-6} \text{K}^{-1} \\ \color{green}{\alpha_{b}} &= \color{green}{11.97} & \times 10^{-6} \text{K}^{-1} \\ \color{blue}{\alpha_{c}} &= \color{blue}{20.80} & \times 10^{-6} \text{K}^{-1} \end{align}$
The distribution of crystallographic grain orientation is considered as random

Again, it induces thermal internal stresses leading to microcracks

Quantitative 3D simulation - Overview

Al2TiO5 DEM microstructure given by Voronoi tesselation

Mechanical input paramters
→ $E = 170$ GPa
→ $\nu = 0.28$
→ $\sigma_f = 375$ MPa

Thermal expansion parameters
→ $\alpha_{a} = -2.38 \times 0.8 \times 10^{-6}$ K $^{-1}$
→ $\alpha_{b} = 11.97 \times 0.8 \times 10^{-6}$ K $^{-1}$
→ $\alpha_{c} = 20.80 \times 0.8 \times 10^{-6}$ K $^{-1}$

Remarks
→ 3D periodic conditions
→ $0.8$ factor comes from glassy phase
→ No mechanical anisotropy
→ Grain random orientation
→ ~40 DEs per grain
→ Total number of grains: ~550

About PlugIn

Users are free to assemble these bricks to build their own simulations

<GRANOO TotIteration="20000" TimeStep="2e-5" OutDir="Results"/>

<STEP Label="PreProcessing">
 <NEW-FRAME Center="(0, 0, 0)" Quat="(0, 0, 1)(90)" ID="F1"/>
 <NEW-GROUND Type="Rectangle" DimY="2." DimZ="0.2" FrameID="F1" ID="G1"/>
...
</STEP>

<STEP Label="Processing">
 <AddElement/>
 <CLEAR-LOAD/>
 <APPLY-GRAVITY X="0." Y="-10." Z="0."/>
 <MANAGE-COLLISION Between="Body/Ground" .../>
 <MANAGE-COLLISION Between="Body/Body".../>
 <INTEGRATE-ACCELERATION Linear="Yes" Angular="Yes"/>
</STEP>

</GRANOO>

A lot of macro-commands are provided by GranOO
Users can develop their own macro-commands thanks to plugin

1/155

Discrete Element Modeling of ceramics Micromechanics using discrete element modeling Damien André* * Univ. Limoges, IRCER, UMR 7315, F-87000 Limoges, damien.andre@unilim.fr