PFEM–DEM for particle-laden flows with free surface

Abstract

This work proposes a fully Lagrangian formulation for the numerical modeling of free-surface particle-laden flows. The fluid phase is solved using the particle finite element method (PFEM), while the solid particles embedded in the fluid are modeled with the discrete element method (DEM). The coupling between the implicit PFEM and the explicit DEM is performed through a sub-stepping staggered scheme. This work only considers suspended spherical particles that are assumed not to affect the fluid motion. Several tests are presented to validate the formulation. The PFEM–DEM results show very good agreement with analytical solutions, laboratory tests and numerical results from alternative numerical methods.

Appendix

All vectors and matrices introduced in Sect. 3 are listed below. The variables used in the following refer to the fluid domain, hence the subindex f is here omitted (e.g. $\rho ={\rho }_f$).

The fully-discretized form of the linear momentum equations (Eq. (1)) reads

$$\begin{aligned} \varvec{K} \varvec{\varDelta } {\bar{{\varvec{u}}}}=\varvec{R} \end{aligned}$$

(23)

where $\varvec{K}$ is computed as the sum of stiffness and mass matrices as

$$\begin{aligned} \varvec{K}= \varvec{K}^{m} + \varvec{K}^{\rho } \end{aligned}$$

(24)

with

$$\begin{aligned} \displaystyle \varvec{K}^{\rho }_{IJ}= & {} \varvec{I}\int _{\varOmega } N_I \frac{2\rho }{\varDelta t} N_J d\varOmega \end{aligned}$$

(25)

$$\begin{aligned} \displaystyle \varvec{K}^{m}_{IJ}= & {} \int _{\varOmega } \varvec{B}^{T}_I \varvec{c} \varvec{B}_J d\varOmega \end{aligned}$$

(26)

and $\varvec{R}$ is the residual of the discretized linear momentum equations computed at each iteration i as

$$\begin{aligned} \displaystyle {R^i_{Ii}}= & {} \int _\varOmega N_I \rho N_J d\varOmega \ {\bar{{\dot{u}}}^i_{Ji}}+ \int _\varOmega \frac{\partial N_I}{\partial x_j} {\sigma }^i_{ij} d\varOmega \nonumber \\&- \int _\varOmega {N}_I \rho {g}_i d\varOmega \end{aligned}$$

(27)

where $N_I$ are the linear shape functions and, for a 3D problem, $\varvec{B}$ and $\varvec{c}$ are defined as follows

$$\begin{aligned} \displaystyle \mathbf{B}_I= & {} \left[ \begin{array}{ccc} \displaystyle {\partial N_I \over \partial x} &{}0 &{}0 \\ \displaystyle 0 &{} \displaystyle {\partial N_I \over \partial y} &{}0 \\ \displaystyle 0 &{} 0 &{} \displaystyle {\partial N_I \over \partial z} \\ \displaystyle \displaystyle {\partial N_I \over \partial y} &{}\displaystyle {\partial N_I \over \partial x} &{} 0 \\ \displaystyle \displaystyle {\partial N_I \over \partial z} &{} 0 &{} \displaystyle {\partial N_I \over \partial x} \\ 0 &{}\displaystyle {\partial N_I \over \partial z} &{} \displaystyle \displaystyle {\partial N_I \over \partial y} \end{array} \right] ,~ \\ \displaystyle {\varvec{c}}= & {} \left[ \begin{array}{cccccc} \kappa \varDelta t + \frac{4 \mu }{3 } &{} \kappa \varDelta t - \frac{2 \mu }{3} &{} \kappa \varDelta t - \frac{2 \mu }{3} &{} 0 &{} 0 &{} 0 \\ \kappa \varDelta t - \frac{2 \mu }{3}&{} \kappa \varDelta t+ \frac{4 \mu }{3 } &{} \kappa \varDelta t - \frac{2 \mu }{3} &{} 0 &{} 0 &{} 0 \\ \kappa \varDelta t - \frac{2 \mu }{3} &{} \kappa \varDelta t - \frac{2 \mu }{3} &{} \kappa \varDelta t+ \frac{4 \mu }{3 } &{} 0 &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} \mu &{} 0 &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} \mu &{} 0\\ 0 &{} 0 &{} 0 &{} 0 &{} 0 &{} \mu \\ \end{array} \right] \end{aligned}$$

The fully-discretized form of the continuity equation (Eq. (4)) reads

$$\begin{aligned} \varvec{H}\bar{\varvec{p}}=\varvec{F}_{p} \end{aligned}$$

(28)

where

$$\begin{aligned} \varvec{H}= & {} \left( \frac{1}{ \varDelta t} \mathbf{M }_{1} + \frac{1}{ {\varDelta t}^2} \mathbf{M }_{2} + \mathbf{L}+\mathbf{M}_{b} \right) \end{aligned}$$

(29)

$$\begin{aligned} \varvec{F}_{p}= & {} \frac{\mathbf{M }_{1}}{\varDelta t}{\bar{\varvec{p}}^n } + \frac{ \mathbf{M }_{2}}{{\varDelta t}^2} \left( {{\bar{{\varvec{p}}}}^n }+ { \bar{\dot{\varvec{p}}}^n} \varDelta t \right) +{\varvec{f}}_{p} - \mathbf{Q}^T\ {\bar{{{\varvec{v}}}}^{n+1}} \end{aligned}$$

(30)

with

$$\begin{aligned} \displaystyle {M}_{1_{IJ}}= & {} \int _{\varOmega } \frac{1}{\kappa }{N}_I {N}_J d\varOmega \end{aligned}$$

(31)

$$\begin{aligned} \displaystyle {M}_{2_{IJ}}= & {} \int _{\varOmega } \tau \frac{\rho }{\kappa }{N}_I {N}_J d\varOmega \end{aligned}$$

(32)

$$\begin{aligned} {L}_{{IJ}}= & {} \int _{\varOmega } \tau ({\pmb \nabla }^T {N}_I) {\pmb \nabla }{N}_J d\varOmega \end{aligned}$$

(33)

$$\begin{aligned} \displaystyle {M}_{b_{IJ}}= & {} \int _{\varGamma _t} \frac{2\tau }{h_n} {N}_I {N}_J d\varGamma \end{aligned}$$

(34)

where $\varGamma _t$ is the free-surface contour and the stabilization parameter $\tau$ is defined as

$$\begin{aligned} \displaystyle \tau = \left( \frac{8 \mu }{h^2}+ \frac{2\rho }{\delta } \right) ^{-1} \end{aligned}$$

(35)

where h and $\delta$ are characteristic distances in space and time [52].

Finally,

$$\begin{aligned} \displaystyle {f}_{p_I}= & {} \int _{\varGamma _t}\tau {N}_I \left[ \rho \frac{Du_n}{Dt}-\frac{2}{h_n} (2 \mu d_n - {\hat{t}})\right] d\varGamma \nonumber \\&- \int _{\varOmega } \tau {\pmb \nabla }^T {N}_I {\rho {\varvec{g}}} d\varOmega \end{aligned}$$

(36)

$$\begin{aligned} \mathbf{Q }_{IJ}= & {} \int _{\varOmega } \mathbf{B }_I^T \mathbf{m } N_J d\varOmega \end{aligned}$$

(37)

where in 3D $\mathbf{m } = [1,1,1,0,0,0]^T$.

Details on the derivation of above matrices and vectors can be found in [52].