Poiseuille flow of a yield-stress fluid
#include "grid/multigrid.h"
#include "yield.h"
double tau_y = 0.2;
int main() {
periodic (top);
size (0.5);
stokes = true;
TOLERANCE = 1e-5;
u.t[left] = dirichlet(0);
const face vector g[] = {0.,1.};
a = g;
const face vector muc[] = {1.,1.};
mu = muc;
const face vector tauc[] = {tau_y,tau_y};
tau_0 = tauc;
for (N = 8; N <= 64; N *= 2)
run();
}
double analytical (double x)
{
double h = 0.5, dpdy = 1., mu = 1.;
double U = sq(tau_y - h*dpdy)/(2.*mu*dpdy);
double X = h - tau_y/dpdy;
return (x < X ? dpdy*(h*x - sq(x)/2.) - tau_y*x : U);
}
scalar un[];
event logfile (t += 0.1; i <= 1000) {
double du = change (u.y, un);
if (i > 0 && du < 1e-7)
return 1; /* stop */
// fprintf (stderr, "%g %d %g %d %d\n", t, i, du, mgp.i, mgu.i);
}
event profile (t = end) {
if (N == 16) {
foreach()
fprintf (stdout, "%g %g %g %g %g %g %g\n", x, y, u.x[], u.y[], p[],
x - Delta/2., (u.y[] - u.y[-1])/Delta);
printf ("\n");
}
scalar e[];
foreach()
e[] = u.y[] - analytical(x);
norm n = normf (e);
fprintf (ferr, "%d %g %g %g\n", N, n.avg, n.rms, n.max);
}
set xlabel 'x'
set ylabel 'u.y'
h = 0.5
dpdy = 1.
tau_y = 0.2
mu = 1.
U = (tau_y - h*dpdy)**2/(2.*mu*dpdy)
X = h - tau_y/dpdy
u(x) = x < X ? dpdy*(h*x - x**2/2.) - tau_y*x : U
set xrange [0:0.5]
plot 'out' u 1:4 w p t 'basilisk', u(x) w l t 'analytical'
set xlabel 'x'
set ylabel 'du.y/dx'
plot 'out' u 6:7 w p t 'basilisk', (x<X?X-x:0) w l t 'analytical'
reset
set xlabel 'Resolution'
set ylabel 'Error norms'
set key bottom left
set logscale
set cbrange [1:2]
set xtics 16,2,128
set grid
ftitle(a,b) = sprintf("order %4.2f", -b)
f2(x)=a2+b2*x
fit f2(x) 'log' u (log($1)):(log($3)) via a2,b2
fm(x)=am+bm*x
fit fm(x) 'log' u (log($1)):(log($4)) via am,bm
plot exp (fm(log(x))) t ftitle(am,bm), \
exp (f2(log(x))) t ftitle(a2,b2), \
'log' u 1:4 t '|h|_{max}' ps 1.5, \
'log' u 1:3 t '|h|_2' ps 1.5
See also