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){
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.