sandbox/M1EMN/Exemples/couette_muw.c
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();
}
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
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'