sandbox/wmostert/poiseuille.c
We set up a 2D plane Poiseuille flow. This is in part an investigation of the mask function. This problem file is courtesy of J. Eggers.
#include "navier-stokes/centered.h"
Define Reynolds number, channel length, half-width, total time, and resolution.
#define Re (1.) // Reynolds number
#define L (2.) // length of channel
#define hw (L / 11.0) // half-width of channel
#define tfinal (5.) // length of run
#define LEVEL (9)
The domain is eight units long, centered vertically.
int main() {
L0 = L;
origin (-hw, -L0/2.);
init_grid(1 << LEVEL);
Set viscosity to unity.
Set the pressure drop over channel length, and the velocity amplitude.
double dp = 2.*L*Re/hw;
double U = Re / hw;
Set the inlet velocity at the left according to known solution. An outflow condition is used on the right boundary, with the exit pressure specified. No-slip is enforced.
u.n[left] = dirichlet(U*(sq(hw)-sq(y)));
p[left] = neumann(0.);
pf[left] = neumann(0.);
u.n[right] = neumann(0.);
p[right] = dirichlet(0.);
pf[right] = dirichlet(0.);
u.t[top] = dirichlet(0.);
u.t[bottom] = dirichlet(0.);
event init (t = 0) {
Define the masks at the channel halfwidths.
mask (y > hw ? top :
y < -hw ? bottom :
none);
Initialize the interior flow with a Poiseuille profile, although a quiescent initial flow should also work.
Output statistics.
event logfile (i++)
fprintf (stderr, "%d %g %d %d\n", i, t, mgp.i, mgu.i);
Output movies for pressure, horizontal and vertical velocity, and local refinement level.
event movies (i ++; t <= tfinal) {
scalar omega[];
vorticity (u, omega);
static FILE * fp = fopen ("p.ppm", "w");
static FILE * fp1 = fopen ("ux.ppm", "w");
static FILE * fp2 = fopen ("uy.ppm", "w");
static FILE * fp3 = fopen ("l.ppm", "w");
scalar l[];
foreach(){
l[] = level;
}
output_ppm (p, fp, box= {{-hw,-hw},{L-hw,hw}},
spread = 2, linear = true);
output_ppm (u.x, fp1, box = {{-hw,-hw},{L-hw,hw}},
linear = true, spread = 2);
output_ppm (u.y, fp2, n=1024, box= {{-hw,-hw},{L-hw,hw}},
linear = true, spread = 2);
output_ppm (l, fp3, n=1024, box= {{-hw,-hw},{L-hw,hw}} );
}
Output slices through the velocity field; time steps of 1 until final time tfinal.
event output_velocity_profiles (t += 1.; t<=tfinal) {
static FILE * fux = fopen ("poiseuille_profiles.dat", "w");
static FILE * fp = fopen ("poiseuille_pressure.dat", "w");
fprintf (fux, "\n");
fprintf (fux, "# time = %.4f\n",t);
fprintf (fp, "\n");
fprintf (fp, "# time = %.4f\n",t);
for (double yy = -hw; yy <= hw; yy += hw/100.0) {
double xslice = 1.;
fprintf (fux, "%.4f %.4f\n", yy, interpolate (u.x, xslice, yy));
}
for (double xx = -hw; xx <= L-hw; xx += 0.1) {
fprintf (fp, "%.4f %.4f\n", xx, interpolate (p, xx, 0.));
}
}
Adapt on the velocity field.