sandbox/ghigo/src/test-embed/pitching-naca0015.c

    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 "../mydoubleprojection.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 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 = 0.1/(p_w0);

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

      TOLERANCE = 1e-3;

    We use third order velocity and pressure face_value and face_gradient reconstructions.

      foreach_dimension()
        u.x.third = true;
      p.third = pf.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);
        s = adapt_wavelet ({cs}, (double[]){1e-2},
                           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);
    
      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 (t += 0.05/(p_w0))
    {
      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");
    }

    Adaptive mesh refinement

    event adapt (i++)
    {
      adapt_wavelet ({cs,u}, (double[]){1e-2,(cmax),(cmax)},
    		 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.

    fig. 18 from Schneiders al., 2013

    fig. 18 from Schneiders al., 2013

    Vorticity isolines for \theta = 44

    Vorticity isolines for \theta = 44

    fig. 19 from Schneiders al., 2013

    fig. 19 from Schneiders al., 2013

    Vorticity isolines for \theta = 55

    Vorticity isolines for \theta = 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,
    Drag and lift coefficients C_D and C_L as a function of the angle \theta (script)

    Drag and lift coefficients C_D and C_L as a function of the angle \theta (script)

    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'
    (i) Pressure coefficient C_p at \theta = 44 (script)

    (i) Pressure coefficient C_p at \theta = 44 (script)

    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'
    (ii) Pressure coefficient C_p at \theta = 55 (script)

    (ii) Pressure coefficient C_p at \theta = 55 (script)

    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.