# Flow past a rapidly pitching NACA0015 airfoil at Re=10000

This test case is inspired from the numerical work of Visbal et Shang., 1989. We reproduce here the results of Schneiders et al., 2013. The authors studied the laminar flow around a rapidly pitching NACA0015 airfoil using a adaptive grid where the smallest mesh size is equivalent to 2000pt/c (c is the NACA0015 airfoil corde length).

We solve here the 2D incompressible Navier-Stokes equations and add the NACA0015 airfoil using an embedded boundary for Re=10000 on an adaptive grid using double projection.

Vorticity \omega

Mesh level

#include "grid/quadtree.h"
#include "../myembed.h"
#include "../mycentered.h"
#include "view.h"

#define Re (10000.) // Reynolds number based on the cord length
#define p_w0 (0.6) // Pitch rate
#define p_t0 (1.) // Pitch characteristic time
#define p_ts (1.) // Pitch start time

#define uref (1.) // Reference velocity, uref=1
#define tref (1.) // Reference time, tref=c/U=1

#define lmin (2) // Min mesh refinement level
#define lmid (8) // Mesh refinement level (l=8 is 4pt/c)
#define lmax (14) // Max mesh refinement level (l=14 is 256pt/c)
#define cmax (1.e-2) // Mesh adaptation criterium for velocity

We define the NACA0015 airfoil’s pitch axis position p_pos, the pitch angle \theta, pitching rate and shape. The NACA0015 airfoil corde length is set to 1.

const coord p_pos = {0.25, 0.}; // Position of the pitch axis
double theta; // Instantaneous pitching angle of the NACA0015 airfoil
# define p_angle(t,ts,tau,w) ((t) < (ts) ? 0 : (w)*(((t) - (ts)) + (tau)/4.6*(exp(-4.6*((t) - (ts))/(tau)) - 1.))) // NACA0015 pitch angle
# define p_rotation(t,ts,tau,w) ((t) < (ts) ? 0 : (w)*(1. - exp(-4.6*((t) - (ts))/(tau)))) // NACA0015 pitch rotation rate

void particle (scalar c, face vector f, coord pos, double th)
{
double tt = 0.15;
vertex scalar phi[];
foreach_vertex() {
double xrot = pos.x + (x - pos.x)*cos (th) - (y - pos.y)*sin (th);
double yrot = pos.y + (x - pos.x)*sin (th) + (y - pos.y)*cos (th);
if (xrot >= 0. && xrot <= 1.)
phi[] = sq (yrot) - sq (5.*(tt)*(
0.2969*sqrt (xrot)
- 0.1260*(xrot)
- 0.3516*sq (xrot)
+ 0.2843*cube (xrot)
- 0.1015*pow (xrot, 4.)));
else
phi[] = 1.;
}
boundary ({phi});
fractions (phi, c, f);
boundary ((scalar *) {c, f});
restriction ((scalar *) {c, f});
}

We need a field for viscosity, so that the embedded boundary metric can be taken into account.

face vector muv[];

int main ()
{

The domain is 64\times 64.

  L0 = 64.;
size (L0);
origin (-L0/2., -L0/2.);

We set the viscosity field in the event properties to obtain a Reynolds number Re = 10000.

  mu = muv;

We set the maximum timestep.

  DT = 1.e-2*(tref);

We tune the Poisson solver. Note that the TOLERANCE of the multigrid solver is divided by \Delta t^2.

  TOLERANCE_MU = 1.e-3*(uref);
TOLERANCE = 1.e-3;

We use third order face_value and face_gradient reconstructions.

  for (scalar s in {u, p, pf})
s.third = true;

We initialize the grid at the minimum refinement level lmid that allows to describe the NACA0015 airfoil.

  N = 1 << (lmid);
init_grid (N);

run();
}

## Boundary conditions

We use an inlet velocity boundary condition of 1.

u.n[left] = dirichlet (1.);
u.t[left] = dirichlet (0);
uf.n[left] = dirichlet (1.);
uf.t[left] = dirichlet (0);
p[left] = neumann (0);
pf[left] = neumann (0);

u.n[right] = neumann (0.);
u.t[right] = neumann (0.);
uf.n[right] = neumann (0.);
uf.t[right] = neumann (0.);
p[right] = dirichlet (0.);
pf[right] = dirichlet (0.);

## Initial conditions

event init (t = 0)
{

We decrease the value of the CFL (similar value then with the VOF algorithm).

  CFL = 0.5;

We define the initial pitch angle (in rad).

  theta = (p_angle(0., (p_ts), (p_t0), (p_w0)));

We refine the grid down to level lmax around the airfoil.

  astats s;
int ic = 0;
do {
ic++;
particle (cs, fs, p_pos, theta);
fractions_cleanup (cs, fs);
maxlevel = (lmax), minlevel = (lmin));
} while ((s.nf || s.nc) && ic < 100);

We initialize the embedded boundaries.

  particle (cs, fs, p_pos, theta);
fractions_cleanup (cs, fs);
foreach_face()
if (uf.x[] && !fs.x[])
uf.x[] = 0.;
boundary ((scalar *) {uf});

The initial boundary condition is zero velocity on the embedded boundary.

  u.n[embed] = dirichlet (0);
u.t[embed] = dirichlet (0);
p[embed] = neumann (0);

uf.n[embed] = dirichlet (0);
uf.t[embed] = dirichlet (0);
pf[embed] = neumann (0);
}

The viscosity is set to 1/Re and needs to take the embedded boundary metric into account.

event properties (i++)
{
foreach_face()
muv.x[] = fm.x[]/(Re);
boundary ((scalar *) {muv});
}

## Embedded boundaries

The airfoil’s position is advanced to time t + \Delta t and the appropriate embedded boundary conditions (for both u and uf) are set.

NOTE: the rotation rate is defined in the opposite trigonometric direction, the velocity is given by v = - \Omega \times \bm{x - x_{pitch\, axis}}.

event embed_advection (i++)
{
theta = (p_angle(t + dt, (p_ts), (p_t0), (p_w0)));

particle (cs, fs, p_pos, theta);
fractions_cleanup (cs, fs);
}

event embed_boundary (i++)
{
u.n[embed] = dirichlet ((y - p_pos.y)*(p_rotation(t + dt, (p_ts), (p_t0), (p_w0))));
u.t[embed] = dirichlet (-(x - p_pos.x)*(p_rotation(t + dt, (p_ts), (p_t0), (p_w0))));
p[embed] = neumann (sqrt(sq (x - (p_pos.x)) + sq (y - (p_pos.y)) + SEPS)*sq ((p_rotation(t + dt, (p_ts), (p_t0), (p_w0)))));

uf.n[embed] = dirichlet ((y - p_pos.y)*(p_rotation(t + dt, (p_ts), (p_t0), (p_w0))));
uf.t[embed] = dirichlet (-(x - p_pos.x)*(p_rotation(t + dt, (p_ts), (p_t0), (p_w0))));
pf[embed] = neumann (sqrt(sq (x - (p_pos.x)) + sq (y - (p_pos.y)) + SEPS)*sq ((p_rotation(t + dt, (p_ts), (p_t0), (p_w0)))));
}

## Log files

The log file contains the time, timestep, the pitch angle, the pressure and viscous force components exerted by the fluid on the cylinder.

event logfile (i++)
{
coord Fp, Fmu;
embed_force (p, u, mu, &Fp, &Fmu);
fprintf (ferr, "%d %g %g %g %g %g %g %g\n",
i, t, dt,
theta,
Fp.x, Fp.y, Fmu.x, Fmu.y);
}

## Surface pressure coefficient and vorticity

We compute here the distribution of the pressure coefficient C_p and vorticity \omega at the surface of the NACA0015 airfoil when it reaches the pitch angles \theta = 44,\, 55 approximately.

void cpout (FILE * fcp)
{
foreach_leaf() {
if (cs[] > 0. && cs[] < 1.) {
coord b, n;
double area = embed_geometry (point, &b, &n);
double xe = x + b.x*Delta, ye = y + b.y*Delta;
double xcord = p_pos.x + (xe - p_pos.x)*cos (theta)
- (ye - p_pos.y)*sin (theta);

fprintf (fcp, "%g %g %g %g %g %g %g\n",
xe - p_pos.x, // 1
ye - p_pos.y, // 2
xcord, // 3
M_PI + atan2(ye - p_pos.y, xe - p_pos.x), // 4
embed_interpolate (point, p, b), // 5
area*Delta, // 6
embed_vorticity (point, u, b, n) // 7
);
fflush (fcp);
}
}
}

event surface_profile (t = {(p_ts) + 1.4971*(p_t0),
(p_ts) + 1.8172*(p_t0)})
{
char name[80];
sprintf (name, "cp-%d-pid-%d",
(int) ((p_angle(t, (p_ts), (p_t0), (p_w0)))/M_PI*180.), pid());
FILE * fcp = fopen (name, "w");
cpout (fcp);
fclose (fcp);
}

## Vorticity isolines

We plot here the vorticity isolines when the NACA0015 airfoil reaches the pitch angles \theta = 44,\, 55 approximately to compare them with fig. 18 and 19 of Schneiders et al., 2013.

event images_Schneiders2013 (t = {(p_ts) + 1.4971*(p_t0),
(p_ts) + 1.8172*(p_t0)})
{
scalar omega[];
vorticity (u, omega);

clear();
view (fov = 0.4,
tx = -(p_pos.x + 0.2)/L0, ty = -(p_pos.y - 0.1)/L0,
bg = {1,1,1},
width = 400, height = 400);
draw_vof ("cs", "fs", filled = -1, lw = 5);
squares ("omega", map = cool_warm, min = -100, max = 100);

char name[80];
sprintf (name, "vorticity-%d.png",
(int) ((p_angle(t, (p_ts), (p_t0), (p_w0)))/M_PI*180.));
save (name);
}

## Animation

We display here the time evolution of the vorticity \omega.

event movie_vorticity (i += 25)
{
scalar omega[];
vorticity (u, omega);

view (fov = 0.3,
tx = -(p_pos.x + 0.2)/L0, ty = -(p_pos.y - 0.1)/L0,
bg = {1,1,1},
width = 400, height = 200);
char legend[80];
sprintf (legend, "theta=%g", theta);

draw_vof ("cs", "fs", filled = -1, lw = 5);
squares ("omega", map = cool_warm, min = -100, max = 100);
draw_string (legend, pos = 1, lw = 2);
save ("vorticity.mp4", opt = "-r 4");

draw_vof ("cs", "fs", lw = 5);
squares ("level", min = (lmid), max = (lmax), map = cool_warm);
cells ();
draw_string (legend, pos = 1, lw = 2);
save ("mesh.mp4", opt = "-r 4");
}

event adapt (i++)
{
maxlevel = (lmax), minlevel = (lmin));
}

## End of simulation

event stop_run (t = (p_ts) + 1.85*(p_t0))
{
return 1.;
}

## Results for Re=10000 and level=14 using double projection

#### Vorticity

We compare here the vorticity isolines with those of fig. 18 and 19 from Schneiders et al., 2013, obtained when the NACA0015 airfoil reaches the pitch angles \theta = 44,\, 55.

#### Drag coefficient C_D

We next plot the drag and lift coefficents C_D and C_L as a function of the pitch angle \theta. We compare the results with those of fig. 17 from Schneiders et al., 2013.

reset
set key left top
set xlabel 'theta'
set ylabel 'C_{D,L}'
set xrange[0:55]
set yrange[0:5]
plot 'Schneiders2013/Schneiders2013-fig17-CD.csv' u 1:2 w p ps 2 pt 6 lc rgb "black" t "Schneiders et al., 2013, fig. 17, C_D", 'Schneiders2013/Schneiders2013-fig17-CL.csv' u 1:2 w p ps 2 pt 7 lc rgb "black" t "Schneiders et al., 2013, fig. 17, C_L", 'log' u ($4/pi*180):(2.*($5+$7)) w l lw 3 lc rgb "blue" t "l=14 (256pt/D), 2p", 'log' u ($4/pi*180):(2.*($6+$8)) w l lw 3 lc rgb "blue" notitle,

#### Surface pressure coefficient C_p around the NACA0015 airfoil

We finally plot the distribution of the pressure coefficient C_p at the surface of the NACA0015 airfoil when it reaches the pitch angles \theta = 44,\, 55 approximately.

reset
set xlabel 'x/c'
set ylabel 'C_p'
set xrange[0:1]
set yrange[-20:5]
plot '< cat cp-44-pid-* | sort -k3,3' u ($3):(2.*$5) w p pt 6 ps 1.5 lc rgb "blue" t 'l=14 (256pt/D), 2p'
reset
set xlabel 'x/c'
set ylabel 'C_p'
set xrange[0:1]
set yrange[-20:5]
plot '< cat cp-54-pid-* | sort -k3,3' u ($3):(2.*$5) w p pt 6 ps 1.5 lc rgb "blue" t 'l=14 (256pt/D), 2p'

## References

 [Schneiders2013] L. Schneiders, D. Hartmann, M. Meinke, and W. Schroder. An accurate moving boundary formulation in cut-cell methods. Journal of Computational Physics, 235:786–809, 2013. [Visbal1989] M.R. Visbal and J.S. Shang. Investigation of the flow structure around a rapidly pitching airfoil. AIAA Journal, 27:1044–1051, 1989.