sandbox/M1EMN/Exemples/granular.h

    Granular rheology, \mu(I)

    Purpose

    This is the implementation of \mu(I) rheology (GDR MiDi, Jop Pouliquen Forterre).

    Implementation

    We use

    #include "navier-stokes/centered.h"

    with two phases,

    #include "vof.h"

    One phase is a passive fluid of small density to preserve 0 pressure at the surface of the other phase, the granular one. f is the VOF field.

    possible values define RHOF 1e-4 define mug 1e-5, here we take:

    #define RHOF 1e-3
#define mug  1e-4

face vector alphav[];
face vector muv[];
face vector eta[];
scalar rhov[];
// note the v
scalar f[];
scalar * interfaces = {f};

    parameters of the rheology

    by definition (By GDR MiDi) from a scalar analysis of a unidirectional shear flow \overrightarrow{u}=u(y)\overrightarrow{e}_x, the tangential and normal stress are linked through a Coulomb relation : \tau = \mu(I ) p, the friction law is: \displaystyle \mu = \mu_s+ \frac{\Delta \mu}{1+I/I_0}\text{ and where } I=D \frac{\partial u/\partial y}{\sqrt(p)} with in the definition of I, the diameter of grains D=1./30 (30 diameters of grains in the unit), and \partial u/\partial y the shear of velocity and p pressure supposed to be zero at the free surface. The friction coefficients are \mu_s=0.32, \Delta \mu=.28, and I_0=.4, the cohesion is C_o=0, it can be changed.

    In this file, we code the tensorial version \tau_{ij}, the shear will be replaced by \sqrt{2} \sqrt{D_{ij}D_{ji}} where D_{ij} is the shear strain rate tensor (see for analytical solutions the basic Bagnold flow and Bagnold with cohesion

    All those numerical values are to be changed in the main.

    double Co=0,D=1./30,mus=0.32,dmu=.28,I0=.4;
double etamax=10000;

    Do not forget Boundary conditions for granular flow, pressure must be zero at the surface: p[top] = dirichlet(-RHOF*LDOMAIN);

    Total density

    #define rho(f) ((f) + RHOF*(1. - (f)))

    beware that \alpha = 1/\rho and the \mu is the dynamic viscosity or absolute viscosity. mu is a reserved variable of basilisk. To avoid confusion we use symbol \eta in papers and in the redefinition of dynamic viscosity.

    event init_granul (t = 0) {
    alpha = alphav;
    mu = muv;
    rho = rhov;
}

    computing the viscosity

    The total stress is written in a pressure part and a dissipative part \displaystyle \sigma_{ij}= - p \delta_{ij} + \tau_{ij}

    To write this in tensorial form, \tau_{ij} is linked with D_{ij} the shear strain rate tensor (tenseur des taux de déformation) D_{ij}=(u_{i,j}+u_{j,i})/2,

    \displaystyle \tau_{ij} =2 \eta D_{ij}

    where \eta is the equivalent dynamic viscosity

    The components of D_{ij} in 2D are: D_{11}=\frac{\partial u}{\partial x}, D_{12} =\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}), D_{21} =D_{12} =\frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x}), D_{22}=\frac{\partial v}{\partial y}

    And the second invariant is D_2=\sqrt{D_{ij}D_{ij}} hence \displaystyle D_2^2= D_{ij}D_{ij}= ( \frac{\partial u}{\partial x})^2 + (\frac{\partial v}{\partial y})^2 + \frac{1}{2}( \frac{\partial u}{\partial y}+ \frac{\partial v}{\partial x})^2

    The \mu(I) is a non newtonain viscosity, it looks like the Bingham one, but with a threshold depending on the pressure.

    The viscosity is computed as a fonction of D_2=\sqrt{D_{ij}D_{ji}} and p;

    The inertial number I=D \sqrt{2} D_2/\sqrt(p) and \mu(I) = \mu_s+ \frac{\Delta \mu}{1+I/I_0}, the equivalent dynamic viscosity is \displaystyle \eta = \frac{\mu(I)p}{\sqrt{2} D_2} by Jop et al hypothesis, tangential stress is linked to the shear rate tensor by \displaystyle \tau_{ij} = 2 \mu(I) p \frac{D_{ij}}{ \sqrt{2} D_2}

    This equivalent viscosity will be coded next. If the shear is too small, we have to do a regularisation to avoid infinite viscosity at rest. so \eta=min(\eta,\eta_{max}).

    Note that if \eta is too small an artificial small viscosity \rho D \sqrt{gD} is taken see Lagrée et al. 11 § 2.3.1

    Note

    In the pure shear flow D_{11}=D_{22}=0 et D_{12}=D_{21}=(1/2)(\partial u/\partial y), so that D_2=\sqrt{D_{ij}D_{ij}} =\sqrt{ 2 D_{12}^2} = \frac{\partial u}{ \sqrt{2}\partial y}. In a pure shear flow, \partial u/\partial y= \sqrt{2} D_2. The inertial number I is indeed D \partial u/\partial y/\sqrt(p)

    event properties (i++) {
    
    trash ({alphav});
    scalar eta[];
    
    scalar fa[];
    foreach()
    fa[] = (4.*f[] +
            2.*(f[-1,0] + f[1,0] + f[0,-1] + f[0,1]) +
            f[1,1] + f[-1,1] + f[1,-1] + f[-1,-1])/16.;
    boundary ({fa});
    
    foreach() {
        eta[] = mug;
        if (p[] > 0.) {
            double D2 = 0.;
            foreach_dimension() {
                double dxx = u.x[1,0] - u.x[-1,0];
                double dxy = (u.x[0,1] - u.x[0,-1] + u.y[1,0] - u.y[-1,0])/2.;
                D2 += sq(dxx) + sq(dxy);
            }
            if (D2 > 0.) {
                D2 = sqrt(2.*D2)/(2.*Delta); // this D2 is sqrt(2) D2
                double In = D2*D/sqrt(p[]);
                double muI = mus + dmu*In/(I0 + In);
                double etamin = sqrt(D*D*D);
                eta[] = max((Co + muI*p[])/D2, etamin);// this D2 is sqrt(2) D2
                eta[] = min(eta[],etamax);      }
        }
    }
    boundary ({eta});

    foreach_face() {
        double fm = (fa[] + fa[-1,0])/2.;
        muv.x[] = (fm*(eta[] + eta[-1,0])/2. + (1. - fm)*mug);
        // mu.x[] = 1./(2.*fm/(eta[] + eta[-1,0]) + (1. - fm)/mug); was in Gerris
        alphav.x[] = 1./rho(fm);
    }
    foreach()
    rhov[] = rho(fa[]);
    boundary ({muv,alphav,rhov});
}

    alternate possibility

    http://basilisk.fr/sandbox/vatsal/GenaralizedNewtonian/LidDrivenBingham.c

    event properties (i++) {

    scalar fa[]; foreach() fa[] = (4.f[] + 2.(f[-1,0] + f[1,0] + f[0,-1] + f[0,1]) + f[1,1] + f[-1,1] + f[1,-1] + f[-1,-1])/16.; foreach() fa[]= (2*f[] + (f[-1,0] + f[1,0] + f[0,-1] + f[0,1]))/6; boundary ({fa});

    trash ({alphav}); foreach_face() { eta.x[] = mug; if (p[] > 0.) { double D2 = 0.; foreach_dimension() { double dxx = u.x[1,0] - u.x[-1,0]; double dxy = (u.x[0,1] - u.x[0,-1] + u.y[1,0] - u.y[-1,0])/2.; D2 += sq(dxx) + sq(dxy); } if (D2 > 0.) { D2 = sqrt(2.D2)/(2.Delta); // this D2 is sqrt(2) D2 double In = D2D/sqrt(p[]); double muI = mus + dmuIn/(I0 + In); double etamin = sqrt(DDD); eta.x[] = max((Co + muI*p[])/D2, etamin);// this D2 is sqrt(2) D2 eta.x[] = min(eta.x[],100000); } }

    double fm = (fa[] + fa[-1,0])/2.; muv.x[] = (fmeta.x[] + (1. - fm)mug); alphav.x[] = 1./rho(fm); } boundary ({eta});

    foreach() rhov[] = rho(fa[]); boundary ({muv,alphav,rhov}); }

    Biblio

    • Jop

    • GDR MiDi

    Exemples of use

    Exemples of implemetations not yet in granular.h, some en axi

    Confiné 04/2020