couette_muw.c

Periodic couette flow with imposed shear stress at the lower wall

Periodic couette flow with imposed shear stress at the lower wall

This is a simple test of solution of Navier Stokes with mixed BC at the wall. No slip at upper wall were u=U_0, and we impose the shear stress at the lower wall \displaystyle \mu \frac{\partial u}{\partial y}|_0=\tau_w where \tau_w is given. In a classical plane Couette flow :u=U_0 at the upper wall and 0 at lower.

Analytical steady solution is just \tau(y)= \tau_w constant in the layer and \displaystyle u(y) = U_0 - \frac{\tau_w h}{\mu}(1-\frac{y}{h})

#include "navier-stokes/centered.h"
#define LEVEL 4 
double tauw;
scalar mu_eq[];
face vector muv[];

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

int main() {
  L0 = 1.0;
  origin (0., 0);
  tauw=.50;

Boundary conditions are periodic

    periodic (right);

classical couette OK no slip at the top wall which moves at unit velocity no slip at the bottom

#if 0
  u.t[top] = dirichlet(1);
  u.n[top] = dirichlet(0); 
  u.n[bottom] = dirichlet(0); 
  u.t[bottom] =  dirichlet(0);
#endif

impose slope at bottom: OK

#if 0
  u.t[top] = neumann(0);
  u.n[top] = dirichlet(0);
  u.n[bottom] = dirichlet(0); 
  u.t[bottom] = neumann(-tauw);
#endif

impose shear at the wall, note the minus sign (external normal), and the face value of \mu

#if 1
  u.t[top] =  dirichlet(1);
  u.n[top] = dirichlet(0);
  u.n[bottom] = dirichlet(0); 
  u.t[bottom] = neumann(mu.y[] ? -tauw/mu.y[] : 0.);
#endif

we tune the NS solver, as U=u(y)e_x. stokes is defined as true, so there is no CFL condition

 // DT = 0.1;
  NITERMAX = 100;
  TOLERANCE = 1e-3;
  stokes = true; 
  run(); 
}

event init (t = 0) {

prepare viscosity

  mu = muv;

pressure gradient, acceleration mdpdx if we want to do a Poiseuille flow \displaystyle -\frac{\partial p}{\partial x} = 0 \displaystyle -\frac{\partial p}{\partial y} = 0

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

Initialy at rest

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

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 += 1; i <= 10000) {
    double du = change (u.x, un);
    fprintf(stdout,"t= %g %g %g %g \n",t,interpolate (u.x, L0/2, .999),interpolate (u.x, L0/2, .999),du);
    if (i > 0 && du < 5.0e-5)
        return 1; /* stop */
}

Implementation of the Bagnold viscosity

event properties (i++) {

Compute viscosity, here constant and equal to one!

    foreach() {  
     mu_eq[] =1;
    }
    boundary ({mu_eq});
    foreach_face() {
        muv.x[] = (mu_eq[] + mu_eq[-1,0])/2.;
    }
    boundary ((scalar *){muv});
}

Save profiles

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

event profile (t = end)
{
  scalar shearS[];
  foreach()
    shearS[] = mu.y[]*(u.x[0,0] - u.x[0,-1])/(Delta);
  boundary ({shearS});
  for (double y = 0.; y < 1.0; y += 1./pow(2.,LEVEL))
    fprintf (stderr, "%g  %g  %g \n", y, interpolate (u.x, L0/2, y),interpolate (shearS, L0/2, y));
}

Run

To run the program

  qcc -g -O3 -o couette_muw  couette_muw.c -lm
 ./couette_muw

Results and plots

we compare here with the steady solution \displaystyle u(y) = U_0 + \frac{\tau_w h}{\mu}(\frac{y}{h}-1) where the stress is constant \displaystyle \tau(y)= \tau_w

Plots of the velocity and \tau with \tau_w=0.5, U_0=1, h=1, \mu=1

 set xlabel "y"
 set xlabel "U, tau"
 set key left
 p'xprof' u 1:2 t'U comp.' , 1-.5*(1-x) t'U anal.',''u 1:($3) t'tau comp.' , .5 t'imposed tau'

$Velocity, \tau= \mu du/dy profiles computed (script)$

Velocity, \tau= \mu du/dy profiles computed (script)