A Projection Algorithm for the Navier-Stokes equations :

SummaryFluid flows require good algorithms and good triangultions. We show here an example of a complex algorithm and or first example of mesh adaptation.
An incompressible viscous fluid satisfies: $$ \partial _t u + u\cdot\nabla u + \nabla p - \nu\Delta u = 0,\quad \nabla\cdot u=0 \quad \hbox{ in } \Omega\times ]0,T[, $$ $$ u|_{t=0} = u^0,\quad u|_\Gamma = u_\Gamma. $$ A possible algorithm, proposed by Chorin, is $$ {1\over \delta t}[u^{m+1} - u^moX^m] + \nabla p^m -\nu\Delta u^m= 0,\quad u|_\Gamma = u_\Gamma, $$ $$ -\Delta p^{m+1} = -\nabla\cdot u^moX^m, \quad \partial _n p^{m+1} = 0, $$ where \(uoX(x) = u(x-u(x)\delta t)\) since \(\partial _t u + u\cdot\nabla u \) is approximated by the method of characteristics, as in the previous section.
An improvement over Chorin's algorithm, given by Rannacher, is to compute a correction, q, to the pressure (the overline denotes the mean over \(\Omega\)) $$ -\Delta q= \nabla\cdot\vec u - \overline{\nabla\cdot\vec u} $$ and define $$ u^{m+1}=\tilde u + \nabla q\delta t,   p^{m+1}=p^m-q-\overline{p^m-q} $$ where \(\tilde u\) is the \((u^{m+1},v^{m+1})\) of Chorin's algorithm.

The backward facing step

The geometry is that of a channel with a backward facing step so that the inflow section is smaller than the outflow section. This geometry produces a fluid recirculation zone that must be captured correctly.

This can only be done if the triangulation is sufficiently fine, or well adapted to the flow.

Remark (FH), The are a technical difficulty is the example, the discret flow flux must be \(0\) and in the previous version this is not the case, the correction is not so simple.

Example[NSprojection.edp]

 // file NSprojection.edp

 border a0(t=1,0){ x=0;      y=t;      label=1;}
 border a1(t=0,1){ x=2*t;    y=0;        label=2;}
 border a2(t=0,1){ x=2;      y=-t/2;       label=2;}
 border a3(t=0,1){ x=2+18*t^1.2;  y=-0.5;       label=2;}
 border a4(t=0,1){ x=20;     y=-0.5+1.5*t;   label=3;}
 border a5(t=1,0){ x=20*t; y=1;        label=4;}
 int n=1;
 mesh Th= buildmesh(a0(3*n)+a1(20*n)+a2(10*n)+a3(150*n)+a4(5*n)+a5(100*n));
 plot(Th);
 fespace Vh(Th,P1);
 real nu = 0.0025, dt = 0.2; // Reynolds=200
 func uBCin =  4*y*(1-y)*(y>0)*(x<2) ; 
 func uBCout =  4./1.5*(y+0.5)*(1-y) *(x>19);
 Vh w,u = uBCin, v =0, p = 0, q=0;
 real area= int2d(Th)(1.);
 Vh ubc  = uBCin + uBCout; 
 real influx0  = int1d(Th,1) (ubc*N.x), // FH add
       outflux0 = int1d(Th,3) (ubc*N.x); // FH add
 verbosity=1;
 for(int n=0;n<300;n++){

   Vh uold = u,  vold = v, pold=p;
   Vh f=convect([uold,vold],-dt,uold); 
   real outflux = int1d(Th,3) (f*N.x); // FH add
   f = f - (influx0+outflux)/outflux0 * uBCout; // FH add
   outflux = int1d(Th,3) (f*N.x); // FH add
   assert( abs(influx0+outflux) < 1e-10); // WARNING the flux must be 0 ..

   solve pb4u(u,w,init=n,solver=LU)
         =int2d(Th)(u*w/dt +nu*(dx(u)*dx(w)+dy(u)*dy(w)))
         -int2d(Th)((convect([uold,vold],-dt,uold)/dt-dx(p))*w)
         + on(1,u = 4*y*(1-y)) + on(2,4,u = 0) + on(3,u=f);
   plot(u);

   solve pb4v(v,w,init=n,solver=LU)
         = int2d(Th)(v*w/dt +nu*(dx(v)*dx(w)+dy(v)*dy(w)))
         -int2d(Th)((convect([uold,vold],-dt,vold)/dt-dy(p))*w)
         +on(1,2,3,4,v = 0);

  real meandiv = int2d(Th)(dx(u)+dy(v))/area;

  solve pb4p(q,w,init=n,solver=LU)= int2d(Th)(dx(q)*dx(w)+dy(q)*dy(w))
     - int2d(Th)((dx(u)+ dy(v)-meandiv)*w/dt)+ on(3,q=0);

  real meanpq = int2d(Th)(pold - q)/area;
  if(n%50==49){
    Th = adaptmesh(Th,[u,v],q,err=0.04,nbvx=100000);
     plot(Th, wait=true);
     ubc = uBCin + uBCout; // reinterpolate B.C. 
     influx0 = int1d(Th,1) (ubc*N.x); // FH add
     outflux0 = int1d(Th,3) (ubc*N.x); // FH add
  }
  p = pold-q-meanpq;
  u = u + dx(q)*dt;
  v = v + dy(q)*dt;
  real err = sqrt(int2d(Th)(square(u-uold)+square(v-vold))/Th.area) ;
  cout << " iter " << n << " Err L2 = " << err << endl;
  if(err < 1e-3) break;
 }
 plot(p,wait=1,ps="NSprojP.eps");
 plot(u,wait=1,ps="NSprojU.eps");

[htbp]

[figNSproj] Rannacher's projection algorithm: result on an adapted mesh (top) showing the pressure (middle) and the horizontal velocity \(u\) at Reynolds 400.

We show in figure [figNSproj] the numerical results obtained for a Reynolds number of 400 where mesh adaptation is done after 50 iterations on the first mesh.