sandbox/ghigo/src/test-navier-stokes/cylinder-steady.c
Fixed cylinder (moving) at the same speed as the surrounding inviscid flow
In this test case, both the fluid and the cylinder are moving at the same speed. The presence of the embedded boundary should not create any disturbance in the flow.
A similar test case we used in Gerris: hexagon.
We solve here the Euler equations and add the cylinder using an embedded boundary.
#include "../myembed.h"
#include "../mycentered.h"
#include "view.h"
Reference solution
#define d (0.753)
#define uref (0.912) // Reference velocity, uref
#define tref ((d)/(uref)) // Reference time, tref=d/u
We also define the shape of the domain.
#define cylinder(x,y) (sq (x) + sq (y) - sq ((d)/2.))
void p_shape (scalar c, face vector f)
{
vertex scalar phi[];
foreach_vertex()
phi[] = (cylinder (x, y));
boundary ({phi});
fractions (phi, c, f);
fractions_cleanup (c, f,
smin = 1.e-14, cmin = 1.e-14);
}
Setup
We define the mesh adaptation parameters.
#define lmin (7) // Min mesh refinement level (l=7 is 3pt/d)
#define lmax (10) // Max mesh refinement level (l=10 is 24pt/d)
#define cmax (1.e-2*(uref)) // Absolute refinement criteria for the velocity field
int main ()
{
The domain is 32\times 32.
L0 = 32.;
size (L0);
origin (-L0/2., -L0/2.);
We set the maximum timestep.
DT = 1.e-2*(tref);
We set the tolerance of the Poisson solver.
TOLERANCE = 1.e-4;
TOLERANCE_MU = 1.e-4*(uref);
We initialize the grid.
N = 1 << (lmax);
init_grid (N);
run();
}
Boundary conditions
We use inlet boundary conditions.
u.n[left] = dirichlet ((uref));
u.t[left] = dirichlet (0);
p[left] = neumann (0);
pf[left] = neumann (0);
u.n[right] = neumann (0);
u.t[right] = neumann (0);
p[right] = dirichlet (0);
pf[right] = dirichlet (0);
We give boundary conditions for the face velocity to “potentially” improve the convergence of the multigrid Poisson solver.
uf.n[left] = (uref);
uf.n[bottom] = 0;
uf.n[top] = 0;
Properties
Initial conditions
We use “third-order” face flux interpolation.
#if ORDER2
for (scalar s in {u, p, pf})
s.third = false;
#else
for (scalar s in {u, p, pf})
s.third = true;
#endif // ORDER2
We use a slope-limiter to reduce the errors made in small-cells.
#if SLOPELIMITER
for (scalar s in {u}) {
s.gradient = minmod2;
}
#endif // SLOPELIMITER
#if TREE
When using TREE and in the presence of embedded boundaries, we should also define the gradient of u at the cell center of cut-cells.
#endif // TREE
We initialize the embedded boundary.
#if TREE
When using TREE, we refine the mesh around the embedded boundary.
astats ss;
int ic = 0;
do {
ic++;
p_shape (cs, fs);
ss = adapt_wavelet ({cs}, (double[]) {1.e-30},
maxlevel = (lmax), minlevel = (1));
} while ((ss.nf || ss.nc) && ic < 100);
#endif // TREE
p_shape (cs, fs);
We also define the volume fraction at the previous timestep csm1=cs.
csm1 = cs;
We define the no-slip boundary conditions for the velocity.
u.n[embed] = dirichlet ((uref));
u.t[embed] = dirichlet (0);
p[embed] = neumann (0);
uf.n[embed] = dirichlet ((uref));
uf.t[embed] = dirichlet (0);
pf[embed] = neumann (0);
We initialize the velocity to speed-up convergence.
Embedded boundaries
Adaptive mesh refinement
#if TREE
event adapt (i++)
{
adapt_wavelet ({cs,u}, (double[]) {1.e-2,(cmax),(cmax)},
maxlevel = (lmax), minlevel = (1));
We also reset the embedded fractions to avoid interpolation errors on the geometry.
p_shape (cs, fs);
}
#endif // TREE
Outputs
event logfile (i++; t < 2.*(tref))
{
scalar e[], ef[], ep[];
foreach() {
if (cs[] <= 0.)
e[] = ef[] = ep[] = nodata;
else {
e[] = sqrt (sq (u.x[] - (uref)) + sq (u.y[]));
ep[] = cs[] < 1. ? e[] : nodata;
ef[] = cs[] >= 1. ? e[] : nodata;
}
}
boundary ((scalar *) {e, ef, ep});
fprintf (stderr, "%d %g %g %g %g %g %g %g %g\n",
i, t/(tref), dt/(tref),
normf(e).avg, normf(e).max,
normf(ep).avg, normf(ep).max,
normf(ef).avg, normf(ef).max
);
fflush (stderr);
}
Results
We plot the time evolution of the error. We observe a slow but steady increase of the error, compared to the test case cylinder-steady.c. This is, in my opinion, the propagation of errors made in small cells when using the variable divfnc in the function update_tracer. This errors could be reduced by using a slope limiter for the gradient of the velocity or by decreasing the timestep.
reset
set terminal svg font ",16"
set key top left spacing 1.1
set grid ytics
set xtics 0,1,10
set ytics format "%.0e" 1.e-18,1.e-4,1.e4
set xlabel 't/(d/u)'
set ylabel '||error||_{1}'
set xrange [0:2]
set yrange [1.e-18:1.e3]
set logscale y
plot 'log' u 2:($6) w l lw 2 lc rgb "black" t 'cut-cells', \
'' u 2:($8) w l lw 2 lc rgb "blue" t 'full cells', \
'' u 2:($4) w l lw 2 lc rgb "red" t 'all cells
set ylabel '||error||_{inf}'
plot 'log' u 2:($7) w l lw 2 lc rgb "black" t 'cut-cells', \
'' u 2:($9) w l lw 2 lc rgb "blue" t 'full cells', \
'' u 2:($5) w l lw 2 lc rgb "red" t 'all cells