sandbox/M1EMN/Exemples/bingham_simple.c

    Free surface flow of a Bingham fluid with Navier Stokes solver

    An example of 2D complex flow over a plate with a free surface is presented here. The configuration is periodic what is injected to the left comes from the right. On the bottom solid wall there is a no slip boundary condition, on the top fixed free surface a slip condition is imposed. We use the centered Navier-Stokes solver with regularization for viscosity.

    Equations for a Bingham fluid

    The Bingham fluid is a non Newtonian fluid. The Bingham viscosity that we will put in the Navier-Stokes solver is such that if ||\tau|| \le \tau_y then there is no motion D=0 if the stress is high enough ||\tau|| > \tau_y note that ||\tau|| is the modulus defined as the Euclidian norm \sqrt{\frac{1}{2}{\tau_{ij} \tau_{ij}}}.

    It is not \sqrt{\tau_{11}^2 + \tau_{12}^2} as in Balmorth 06, which is the Frobenius norm.

    then stress tensor is linked to the shear rate tensor by \displaystyle \tau_{ij} = 2 \mu_0 D_{ij} + \tau_y \frac{D_{ij}}{||D||}

    where D_{ij} is the shear strain rate tensor (tenseur de taux de déformation) D_{ij}=(u_{i,j}+u_{j,i})/2, the components in 2D:

    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}

    in the Euclidian norm we have: \displaystyle ||D||=\sqrt{\frac{D_{ij}D_{ij}}{2}} The second invariant defined by D_2=\sqrt{D_{ij}D_{ij}} (this is the Frobenius norm) 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 and we have obviously ||D|| = D_2/\sqrt{2}

    We present here the formulation in Balmforth, he uses \dot \gamma which is by his definition \sqrt{\frac{1}{2}\dot \gamma_{ij} \dot \gamma_{ij}} and as \dot \gamma_{ij}=2 D_{ij} then \dot \gamma is \sqrt{2} D_2, that is why we have a \sqrt{2} in the equations. Factorising with 2 D_{ij} to obtain a equivalent viscosity \displaystyle \tau_{ij} = 2( \mu_0 + \frac{\tau_y}{2 ||D|| } ) D_{ij}=2( \mu_0 + \frac{\tau_y}{\sqrt{2} D_2 } ) D_{ij} as defined by Balmforth \displaystyle \tau_{ij} = 2 \mu_{eq} D_{ij} with \displaystyle \mu_{eq}= \mu_0 + \frac{\tau_y}{\sqrt{2} D_2 }

    Exact solution in the proposed case

    We look at an unidirectional flow, a pure shear flow u(y), v=0, so D_{11}=D_{22}=0 and D_{12}=D_{21}=(1/2)(\partial u/\partial y), this gives (mind square root of 2): D_2=\sqrt{D_{ij}D_{ij}} = \frac{1}{\sqrt{2}}\frac{\partial u}{\partial y}.

    \displaystyle \tau_{12} = 2 \mu D_{12} + 2 \tau_y \frac{D_{12}}{ \sqrt{2}D_2} = \mu \frac{\partial u}{\partial y} + \tau_y

    Equilibrium between pressure gradient and viscosity (writting \tau for a shorthand of \tau_{12}) \displaystyle 0=-\frac{\partial p}{\partial x} + \frac{\partial \tau}{\partial y} as there is no stress at the free surface y=h, the stress is \displaystyle \tau = (-\frac{\partial p}{\partial x})(h-y) the stress \tau increases from the free surface, as long as \tau<\tau_y, we are under the threshold, so shear is zero: \frac{\partial u}{\partial y} =0, hence velocity is constant, say it is U. Let us define Y=h-\tau_y/(-\frac{\partial p}{\partial x}), where \tau=\tau_y.

    So : Y<y<h, \tau<\tau_y, \frac{\partial u}{\partial y} =0, and u=U

    Then going down: 0<y<Y we have \tau =\tau_y +\mu \frac{\partial u}{\partial y}.

    This gives: \displaystyle \tau_y +\mu \frac{\partial u}{\partial y} = (-\frac{\partial p}{\partial x})(h-y) and this allows to solve for the velocity profile \displaystyle u = (-\frac{\partial p}{\partial x})(hy-\frac{1}{2}y^2) - \tau_y y which is indeed zero in y=0, and for y=Y, we have the plug flow u=U of value: \displaystyle U= \frac{(\tau_y - h(-\frac{\partial p}{\partial x}) )^2}{2 \mu (-\frac{\partial p}{\partial x})}

    as in y=0 u=0 and in y=Y, \tau=\tau_y, so \mu \frac{\partial u}{\partial y} =0 and is zero after, therefore \displaystyle u =U \left(1- \left( \frac{y-Y}{Y} \right)^2\right)\; with \;U= \frac{(\tau_y - h(-\frac{\partial p}{\partial x}) )^2}{2 \mu (-\frac{\partial p}{\partial x})} \; with \; Y = h-\tau_y/(-\frac{\partial p}{\partial x})

    flux \int_0^h u dy=\int_0^Y u dy + U(h-Y) is \displaystyle \int_0^h u dy= \frac{hU}{3}(1-Y/h)

    #include "navier-stokes/centered.h"
#define LEVEL 6
double tau_y,mu_0,mumax;

    The domain is one unit long. 0<x<1 0<y<1

    int main() {
  L0 = 1.;
  origin (0., 0);

    Values of yeld stress and viscosity

      tau_y=0.25;
  mu_0 = 1;

    the regularisation value of viscosity

      mumax=1000;

    Boundary conditions are periodic

        periodic (right);

    slip at the top

        u.t[top] = neumann(0);
    u.n[top] = neumann(0);
    uf.n[top] = neumann(0);

    no slip at the bottom

        u.n[bottom] = dirichlet(0);
    uf.n[bottom] = dirichlet(0);
    u.t[bottom] = dirichlet(0);

    note the pressure

        p[top] = neumann(0);
    pf[top] = neumann(0);
    p[bottom] = neumann(0);
    pf[bottom] = neumann(0);

    the \Delta t_{max} should be enough small

      DT = 0.05;
  run(); 
}
    face vector muv[];
scalar D2c[];

event init (t = 0) {

    prepare viscosity

      mu = muv;

    presure gradient mdpdx \displaystyle -\frac{dp}{dx} = 1 

      const face vector mdpdx[] = {1.0,0.};
  a = mdpdx;

    Initialy at rest

      foreach() {
    u.x[] = 0;
    u.y[] = 0;
  }
}

    We check the number of iterations of the Poisson and viscous problems.

    event logfile (i++)
  fprintf (stderr, "%d %g %d %d\n", i, t, mgp.i, mgu.i);

    old value of the velocity is saved

    scalar un[];
event init_un (i = 0) {
    foreach()
    un[] = u.x[];
}

    so that when it does not more change we are converged

    event conv (t += 0.1; i <= 10000) {
    double du = change (u.x, un);
    fprintf(stdout,"t= %g \n",t);
    if (i > 0 && du < 1e-5)
        return 1; /* stop */
}

    Implementation of the Bingham viscosity

    event bingham(i++) {
    scalar muB[];

    Compute 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 where the second invariant is D_2=\sqrt{D_{ij}D_{ij}} hence \displaystyle D_2^2= D_{ij}D_{ij}= D_{11}D_{11} + D_{12}D_{21} + D_{21}D_{12} + D_{22}D_{22} the equivalent viscosity is \displaystyle \mu_{eq}= \mu_0 + \frac{\tau_y}{\sqrt{2} D_2 } the final one is the min of of \mu_{eq} and a large \mu_{max}, then the fluid flows always, it is not a solid, but a very viscous fluid.

        foreach() {
        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(D2)/(2.*Delta);
                D2c[] = D2;
                muB[] =  min(tau_y/(sqrt(2.)*D2) + mu_0 , mumax );    }
            }
    boundary ({muB});
    foreach_face() {
        muv.x[] = (muB[] + muB[-1,0])/2.;
    }
    boundary ((scalar *){muv});
}

    Save profiles

    event profiles (t += .1)
{
    FILE * fp = fopen("xprof", "w");
    scalar shear[];
    foreach()
    shear[] = (u.x[0,1] - u.x[0,-1])/(2.*Delta);
    boundary ({shear});
    for (double y = 0.; y < 1.0; y += 1./pow(2.,LEVEL))
        fprintf (fp, "%g %g %g %g\n", y, interpolate (u.x, L0/2, y), interpolate (shear, L0/2, y),
                  interpolate (D2c, L0/2, y));
    fclose (fp);
}

    We adapt according to the error on the velocity field.

    event adapt (i++) {
//  adapt_wavelet ({u}, (double[]){3e-2,3e-2}, 7, 4);
}

    Results and plots

    To run the program

     qcc -g -O3 -o bingham_simple bingham_simple.c -lm
 ./bingham_simple
 
 
 lldb bingham_simple

    Plot of the velocity and shear

     set xlabel "y"
 set ylabel "u, shear"
 Y = 1-.25
 U = Y*Y/2/1
 p[:][:]'xprof' u 1:2 t'U computed' w lp ,''u 1:3 t'tau comp.',(x<Y ?U*(1-((x-Y)/Y)**2):U) t 'U exact',(x<Y?Y-x:0) t'shear exact'
    Velocity and stress profiles for Bingham flow (script)

    Velocity and stress profiles for Bingham flow (script)

    Paris Avril 2015

    Bibliography