An unresolved FEM-DEM model
for granular suspension
M.Henry, F.Dubois, J.Lambrechts, V.Legat
www.migflow.be
Two families of suspension
-
Coilloidal suspension : 1nm - 5$\mu$m
-
Paints, Inks, ...
-
Granular suspension : > 5$\mu$m
-
Submarine avalanches, mudflow, ...
A brief history of granular rheology
Direct approaches to deal with two phases
-
Body-fitted mesh
-
Fluid-Structure interaction at the boundary
-
Immersed boundary method
-
Rigidity constraint at the solid position
Unresolved DEM-FEM model
- Solid phase
Discrete element
Lagrangian representation
Insight in contacts physics
- Fluid phase
Continuous medium
Eulerian representation
Computational convenience
- Fluid-grains interaction
Empirical momentum transfer
Non-smooth contact dynamics for dense regime
- Newton's second law \[ \frac{d (m_k \mathbf{u}_k)}{d t} = m_k\mathbf{g} - \mathbf{f}_k \]
- Dynamics equation \[ \mathbf{M}\Delta \mathbf{u} = \mathbf{p} \]
- Local Frame \[ \mathbf{U} = \mathbf{H}\mathbf{u} \hspace{5mm} \mathbf{P} = \mathbf{H} \mathbf{p} \]
- Solve for $\mathbf{P}$ \[ \begin{pmatrix} \Delta U_n \\ \Delta U_t \end{pmatrix} = \underbrace{\begin{pmatrix} w_{nn} & 0 \\ 0 & w_{tt} \end{pmatrix}}_{ \mathbf{H}\mathbf{M}^{-1}\mathbf{H}^T} \begin{pmatrix} P_n \\ P_t \end{pmatrix} \]
- Normal and tangential components \[ \begin{aligned} P_n &= \frac{1}{w_{nn}}\text{max}(0, U_n - \frac{\delta}{\Delta t})\\ P_t &= \frac{1}{w_{tt}}\text{min}(|U_t|, \frac{w_{tt}}{w_{nn}}\mu |U_n|)\frac{U_t}{|U_t|} \end{aligned} \]
- Update velocities \[ \mathit{u} \rightarrow \mathit{u} + \mathbf{M}^{-1}\mathbf{H}^T \mathbf{P} \]
Parametrizations are required
- Averaged Navier-Stokes equations, $\mathbf{u} = \epsilon \mathbf{v}$, $\mathbf{u}^s = \phi \mathbf{v}^s$ \[\begin{aligned} &\color{red}{\nabla \cdot \mathbf{u}^s} + \nabla \cdot \mathbf{u} = 0\\ &\rho \left[\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \left(\frac{\mathbf{u}\mathbf{u}}{\color{red}{\epsilon}}\right)\right] = -\color{red}{\epsilon} \nabla p + \nabla \cdot (2\eta^*\color{red}{\epsilon}\mathbf{d}) - \color{red}{\mathbf{f}} + \color{red}{\epsilon} \rho \mathbf{g} \end{aligned} \]
- Drag parametrization \[ \mathbf{f} = g(\epsilon) C_d \pi r^2 \frac{\rho}{2} \left\lVert \mathbf{v^s} - \mathbf{v} \right\rVert \left( \mathbf{v}^s - \mathbf{v} \right) \]
- Subgrid scale viscous model \[ \eta^* = \eta\ \underbrace{(1+2.5\phi + 5\phi^2)}_{\text{Batchelor model}} \]
Comparison
to experimental observations
\[
\eta_r = \left(1 - \frac{\phi}{\phi_c}\right)^\alpha
\]
A model applicable to sensible experiments
Instabilities may arise once immerged
-
Fluid-particle interaction force
\[ \mathbf{f}_{fp} = -V_p \nabla p - \gamma \left(\mathbf{v} - \mathbf{u}\right) \] From Newton's second law :
\[ \frac{\text{d}}{\text{dt}}(m \mathbf{v}) = m \mathbf{g} + \mathbf{f}_{fp} + \mathbf{f}_{pp} \] From semi implicit scheme, a particle velocity predication is made :
\[ \frac{m}{\Delta t} (\mathbf{v}^* - \mathbf{v}^n) = m \mathbf{g} - V_p \nabla p^{n+1} - \gamma^n \left(\mathbf{v}^* - \mathbf{u}^{n+1}\right) + \mathbf{f}_{pp}^n \] The fluid-particle force can lead to instabilities : \[ \mathbf{f}_{fp}^{n+1} = F(\mathbf{u}^{n+1}, \mathbf{p}^{n+1}, \mathbf{v}^{n}, \mathbf{f}_{pp}^n) \]
How to avoid pressure instability
-
Decompose the pressure :
\[
\begin{aligned}
\nabla p &= \nabla P_{h} + \nabla P_{d} \\
&= -\rho \mathbf{g} + \epsilon \gamma \left(\mathbf{v} - \mathbf{u}\right)
\end{aligned}
\]
-
Add an extra-diffusivity as a stabilization term (PSPG-SUPG)
NSCD with cohesive particles
-
Forces can be added into a particle-particle contact to study agglomeration phenomenon.
Force derivated from Morse's potential : \[U_{morse}(r) = D_e \left[1 - e^{-\beta d \left(\frac{r}{d}-1\right)}\right]^2\]
Constitutive law proposed by A.Yazdani[2017]
Predefined bond for particle chains
-
Each chain is pre-generated contacts with a cohesive spring force based on a FENE model :
\[
\begin{aligned}
F_{s} = \frac{k \Delta \mathbf{x}}{1 - \left(\frac{\Delta\mathbf{x}}{Q}\right)^2}
\end{aligned}
\]
Small timestep avoid wide reactions
-
Monodisperse chains submitted to a fixed-volume simple shear
-
An over-simplified model for platelet interaction.
A droplet can dig into a granular media
Granular bed is fluidized by Leidenfrost effect
The droplet digs up to 10 times its radius
New developments towards thermal particulate flows
-
Main difficulties :
- Couple the temperature with the fluid solver
- Solve a FEM-DEM energy equation
- A model for vapor creation
- Track the bubble
A weak coupling by the Boussinesq's approximation
-
Navier-Stokes equation :
\[
\begin{aligned}
&\nabla \cdot \mathbf{u} = 0\\
&\rho \left[\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \left(\mathbf{u}\mathbf{u}\right)\right] = -\nabla p + \nabla \cdot (2\eta\mathbf{d}) + \color{red}{\rho (1-\alpha T) \mathbf{g}}
\end{aligned}
\]
Energy Equation : \[ c_p \rho \left(\frac{\partial T}{\partial t} + \nabla \cdot \left(\mathbf{u} T\right)\right) = \nabla \cdot (k \nabla T) \]
An evaporation model for the vapor
-
Gas mixture composed of two ideal gas mixture,
\[ p_g = p_v + p_a \]
At the interface, from the Clausius-Clapeyron law, the saturated vapor pressure leads to \[ p_v\vert_{\Gamma} = p_g \exp\left[ -\frac{L_{vap}M_1}{R} \left(\frac{1}{T\vert_{\Gamma}} - \frac{1}{T_{sat}}\right) \right] \] This is linked to the energy equation by an interface condition \[ \left[ k \nabla T \cdot \mathbf{n} \right]_{\Gamma} = \dot{m} \left(L_{vap} + (C^p_{f} - C^p_{v})\ (T_{sat} - T\vert_{\Gamma})\right) \]
The local mass flow rate is deduced
from evolution vapor
-
The conservation of vapor is given by
\[
\rho \frac{\partial Y}{\partial t} + \nabla \cdot \mathbf{u}Y = \nabla \cdot \left(\rho D \nabla Y \right)
\]
with the interface condition
\[
\begin{aligned}
&\left[\rho D \nabla Y \cdot \mathbf{n}\right]_{\Gamma} = -\dot{m} \left[Y\right]_{\Gamma}\\
&Y_{\Gamma} = \frac{p_v\vert{\Gamma}M_v}{p_v\vert{\Gamma}M_v + (p_g - p_v)M_a}
\end{aligned}
\]
From the evolution of vapor, the local mass flow rate can be deduced
The evaporation is tracked by a mesh deformation and adaptation
The arbitrary Lagrangian Euler representation is used
Energy Fluid-grains coupling
-
Averaged Navier-Stokes equation :
\[
\begin{aligned}
&\color{red}{\nabla \cdot \mathbf{u}^s} + \nabla \cdot \mathbf{u} = 0\\
&\rho \left[\frac{\partial \mathbf{u}}{\partial t} + \nabla \cdot \left(\frac{\mathbf{u}\mathbf{u}}{\color{red}{\epsilon}}\right)\right] = -\color{red}{\epsilon} \nabla p + \nabla \cdot (2\eta^*\color{red}{\epsilon}\mathbf{d}) - \color{red}{\mathbf{f}} + \color{red}{\epsilon} \rho (1-\alpha T) \mathbf{g}
\end{aligned}
\]
Energy Equation : \[ c_p \rho \left(\frac{\partial (\color{red}{\epsilon} T)}{\partial t} + \nabla \cdot \left(\mathbf{u} T\right)\right) = \nabla \cdot (k \color{red}{\epsilon} \nabla T) + \color{red}{Q_p} \] with \[ \color{red}{Q_p} = \sum_{i \in \mathcal{P}} S_i h_i (T_i - T_i^s) \]