Tutorial 3 - gnomeCreative/HYBIRD GitHub Wiki

logo

Granular flow down an incline with LBM

This tutorial reproduces a steady, free-surface granular flow on a simple infinite incline of constant slope $\theta$, as shown in the figure below. The flow thickness $H$ is assumed constant, and the inclination is set by tilting the principal axis with respect to an imposed pressure gradient $g$, i.e. $g_x = g \sin \theta$, $g_z = g \cos \theta$. The boundary conditions are no-slip at the bottom wall and zero pressure at the top ($z=H$).

We assume that the granular flow behaviour can be well represented by the $\mu(I)$ rheology. Following Jop et al. (2006), the rheology is implemented by locally defining the viscosity as:

$\nu = \frac{p \mu(I)}{\dot{\gamma}}$,

where $\dot{\gamma}$ is the second invariant of the shear rate tensor, and $I$ is the inertial number, defined as:

$I = \frac{\dot{\gamma}d}{P / \rho_\mathrm{s}}$,

where $d$ and $\rho_\mathbf{s}$ are the diameter and density of the grains, respectively. The dependence of the friction coefficient on the inertial number is phenomenologically set as:

$\mu(I) = \mu_\mathrm{s} + \frac{\mu_2 - \mu_\mathrm{s}}{I_0 / I + 1}$.

The $\mu(I)$ rheology needs a set of five material-dependent parameters: $d$, $\rho_\mathrm{s}$, the critical friction coefficient $\mu_\mathrm{s}$, its asymptotic value $\mu_2$, and the constant $I_0$. The three last parameters can be obtained using an imposed-pressure rheometer.

In LBM, in analogy to other non-Newtonian rheologies, the $\mu(I)$ rheology can be reproduced, but a few restrictions apply. Whenever the shear rate is very low, the viscosity diverges. This is the way the fluid reproduces the yield threshold implicit to the Coulomb law. Areas where this threshold is not reached are static, hence the divergence of the viscosity. A second problem arises when pressure is close to zero, e.g., at the free surface. In this case, viscosity reduces to zero.

show_inclineFlow Geometry of the granular flow on an incline.

On an infinite incline, and at steady state, the conservation of momentum in the two directions leads to:

$\tau_{xz} = \rho g \sin \theta H (1 - z / H)$,

$p = \rho g \cos \theta H (1 - z / H)$,

which implies that the friction coefficient, defined as $\mu(I) = \tau / p$, is constant for a given slope $\theta$ and equal to $\mu(I) = \tan \theta$, and that the inertial number is itself constant, and equal to

$I(\theta) = I_0 \frac{\tan \theta - \mu_0}{\mu_2 - \tan \theta}$.

For more details, please see Lagree et al. (2011). From the definition of the inertial number, one obtains the shear rate as

$\frac{\partial u}{\partial z} = \frac{I(\theta)}{d} \left(g H \cos \theta\right)^{\frac{1}{2}} \left(1 - z / H\right)^{\frac{1}{2}}$,

which can be integrated to yield the velocity profile:

$u(z) = \frac{2}{3} \frac{I(\theta)}{d} \left(g H^3 \cos \theta\right)^{\frac{1}{2}} \left(1 - \left(1 - z / H\right)^{\frac{3}{2}}\right)$,

and the maximum velocity (at the surface) is

$u_\textrm{max} = u(H) = \frac{2}{3} \frac{I(\theta)}{d} \left(g H^3 \cos \theta\right)^{\frac{1}{2}}$.

In this tutorial, we retrieve this solution by applying the same geometrical setup as in Tutorial 2. The simulation is halted after stationary conditions are reached. The velocity profile $u(z)$ that should be obtained is illustrated, for a slope with $\theta = 0.43$.

mui_profile Profiles for velocity, pressure, and equivalent viscosity.

If $\tan \theta \lt \mu_0$, no flow is expected to occur, and a vanishing inertial number is predicted for $\tan \theta \rightarrow \mu_0$. The regularization of viscosity and pressure implies that a slow creeping flow will be reproduced when $\tan \theta < \mu_0$. However, the magnitude of this flow is small enough to be considered negligible. Fig. \ref{fig } shows the velocity at the free surface (the maximum) for variable $\theta$.

umax_profile Variation of the free-surface velocity as a function of the slope incline.

Setting the environment

All configuration files should be set inside the tutorials/granular_flow folder. It is up to the user whether to visualize the flow using Eulerian or Lagrangian files. The configuration file is already ready for use (granular_flow.cfg), and is very similar to the one used for the previous tutorial. Only the differences are described in the following: Domain Setup

  • forceX = 3.99: The imposed pressure gradient parallel to the incline $g \sin(\theta)$, with $\theta = 0.42$ (4^\circ$).
  • forceZ = -8.95: The imposed pressure gradient orthogonal to the incline $g \cos(\theta)$, with $\theta = 0.42$ (4^\circ$).

Physical Parameters

  • rheologyModel = MUI: We use a simple Newtonian model.
  • fluidDensity = 1.5e3: The bulk granular density $\rho_\textrm{b}$ (a typical value).
  • initVisc = 0.2: The initial dynamic viscosity $\nu$ (will be used to initialize viscosity).
  • frictionCoefFluid = 0.382: The base friction coefficient, $\mu_0 = \tan(20.9^\circ)$.
  • deltaFriction = 0.42: The difference $\mu_2 - \mu_0 = \tan(32.76^\circ) - \tan(20.9^\circ)$.
  • baseInertial = 0.279: The base inertial number $I_0$.
  • particleDiameterFluid = 0.00054: The diameter of the particles constituting the granular continuum.
  • turbulenceSolver = 0: The turbulence model should not be used for non-Newtonian models.

Code Configuration

  • fluidTimeStep = 1e-5: The time step $\Delta t$. In this case, it cannot be chosen automatically, since it is a key choice with respect to stability and accuracy.

Based on the choices of time step $\Delta t$ and lattice spacing $\Delta S$, the code will work with the same boundaries for the dynamic viscosity as in the previous tutorial. The tutorial folder contains an Excel sheet scaling_granular.xlsx, which can be used to explore the effect of a change in the parameter setup on the viscosity range.