FreeFEM Documentation on GitHub

stars - forks

A projection algorithm for the Navier-Stokes equations

Summary : Fluid 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:

\[\begin{split}\begin{array}{rcl} \partial_t \mathbf{u} + \mathbf{u}\cdot\nabla\mathbf{u} + \nabla p - \nu\Delta\mathbf{u} &= 0 &\hbox{ in } \Omega\times ]0,T[\\ \nabla\cdot\mathbf{u} &= 0 &\hbox{ in } \Omega\times ]0,T[\\ \mathbf{u}|_{t=0} &= \mathbf{u}^0\\ \mathbf{u}|_\Gamma &= \mathbf{u}_\Gamma \end{array}\end{split}\]

A possible algorithm, proposed by Chorin, is:

\[\begin{split}\begin{array}{rcl} {1\over \delta t}[\mathbf{u}^{m+1} - \mathbf{u}^mo\mathbf{X}^m] + \nabla p^m -\nu\Delta \mathbf{u}^m &=& 0\\ \mathbf{u}|_\Gamma &=& \mathbf{u}_\Gamma\\ \nu \partial_n \mathbf{u}|_{\Gamma_{out}} &=&0 \end{array}\end{split}\]
\[\begin{split}\begin{array}{rcl} -\Delta p^{m+1} &= -\nabla\cdot \mathbf{u}^mo\mathbf{X}^m &\\ \partial_n p^{m+1} &= 0 &\mbox{ on } \Gamma\\ p^{m+1} &= 0 &\mbox{ on } \Gamma_{out} \end{array}\end{split}\]

where \(\mathbf{u}o\mathbf{X}(x) = \mathbf{u}(\mathbf{x}-\mathbf{u}(\mathbf{x})\delta t)\) since \(\partial_t \mathbf{u} + \mathbf{u}\cdot\nabla \mathbf{u}\) is approximated by the method of characteristics, as in the previous section.

We use the Chorin’s algorithm with free boundary condition at outlet (i.e. \(p=0,\nu \partial_n u = 0\)), to compute a correction, q, to the pressure.

\[\begin{split}\begin{array}{rcl} -\Delta q &= \nabla\cdot\mathbf{u}&\\ q &= 0 \mbox{ on } &\Gamma_{out} \end{array}\end{split}\]

and define

\[\begin{split}\begin{array}{rcl} \mathbf{u}^{m+1} &=& \tilde{\mathbf{u}} + P \nabla q\delta t\\ p^{m+1} &=& p^m-q \end{array}\end{split}\]

where \(\tilde{\mathbf{u}}\) is the \((\mathbf{u}^{m+1}, v^{m+1})\) of Chorin’s algorithm, and where \(P\) is the \(L^2\) projection with mass lumping ( a sparse matrix).

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.

Note

There is a technical difficulty in the example: the output B.C. Here we put \(p=0\) and \(\nu \partial_n u = 0\).

  1 // Parameters
  2 verbosity = 0;
  3 int nn = 1;
  4 real nu = 0.0025;
  5 real dt = 0.2;
  6 real epsv = 1e-6;
  7 real epsu = 1e-6;
  8 real epsp = 1e-6;
  9 
 10 // Mesh
 11 border a0(t=1, 0){x=-2; y=t; label=1;}
 12 border a1(t=-2, 0){x=t; y=0; label=2;}
 13 border a2(t=0, -0.5){x=0; y=t; label=2;}
 14 border a3(t=0, 1){x=18*t^1.2; y=-0.5; label=2;}
 15 border a4(t=-0.5, 1){x=18; y=t; label=3;}
 16 border a5(t=1, 0){x=-2+20*t; y=1; label=4;}
 17 
 18 mesh Th = buildmesh(a0(3*nn) + a1(20*nn) + a2(10*nn) + a3(150*nn) + a4(5*nn) + a5(100*nn));
 19 plot(Th);
 20 
 21 // Fespace
 22 fespace Vh(Th, P1);
 23 Vh w;
 24 Vh u=0, v=0;
 25 Vh p=0;
 26 Vh q=0;
 27 
 28 // Definition of Matrix dtMx and dtMy
 29 matrix dtM1x, dtM1y;
 30 
 31 // Macro
 32 macro BuildMat()
 33 {   /* for memory managenemt */
 34     varf vM(unused, v) = int2d(Th)(v);
 35     varf vdx(u, v) = int2d(Th)(v*dx(u)*dt);
 36     varf vdy(u, v) = int2d(Th)(v*dy(u)*dt);
 37 
 38     real[int] Mlump = vM(0, Vh);
 39     real[int] one(Vh.ndof); one = 1;
 40     real[int] M1 = one ./ Mlump;
 41     matrix dM1 = M1;
 42     matrix Mdx = vdx(Vh, Vh);
 43     matrix Mdy = vdy(Vh, Vh);
 44     dtM1x = dM1*Mdx;
 45     dtM1y = dM1*Mdy;
 46 } //
 47 
 48 // Build matrices
 49 BuildMat
 50 
 51 // Time iterations
 52 real err = 1.;
 53 real outflux = 1.;
 54 for(int n = 0; n < 300; n++){
 55     // Update
 56     Vh uold=u, vold=v, pold=p;
 57 
 58     // Solve
 59     solve pb4u (u, w, init=n, solver=CG, eps=epsu)
 60         = int2d(Th)(
 61               u*w/dt
 62             + nu*(dx(u)*dx(w) + dy(u)*dy(w))
 63         )
 64         -int2d(Th)(
 65                 convect([uold, vold], -dt, uold)/dt*w
 66             - dx(p)*w
 67         )
 68         + on(1, u=4*y*(1-y))
 69         + on(2, 4, u=0)
 70         ;
 71 
 72     plot(u);
 73 
 74     solve pb4v (v, w, init=n, solver=CG, eps=epsv)
 75         = int2d(Th)(
 76               v*w/dt
 77             + nu*(dx(v)*dx(w) + dy(v)*dy(w))
 78         )
 79         -int2d(Th)(
 80                 convect([uold,vold],-dt,vold)/dt*w
 81             - dy(p)*w
 82         )
 83         +on(1, 2, 3, 4, v=0)
 84         ;
 85 
 86     solve pb4p (q, w, solver=CG, init=n, eps=epsp)
 87         = int2d(Th)(
 88               dx(q)*dx(w)+dy(q)*dy(w)
 89         )
 90         - int2d(Th)(
 91               (dx(u)+ dy(v))*w/dt
 92         )
 93         + on(3, q=0)
 94         ;
 95 
 96     //to have absolute epsilon in CG algorithm.
 97     epsv = -abs(epsv);
 98     epsu = -abs(epsu);
 99     epsp = -abs(epsp);
100 
101     p = pold-q;
102     u[] += dtM1x*q[];
103     v[] += dtM1y*q[];
104 
105     // Mesh adaptation
106     if (n%50 == 49){
107         Th = adaptmesh(Th, [u, v], q, err=0.04, nbvx=100000);
108         plot(Th, wait=true);
109         BuildMat // Rebuild mat.
110     }
111 
112     // Error & Outflux
113     err = sqrt(int2d(Th)(square(u-uold)+square(v-vold))/Th.area);
114     outflux = int1d(Th)([u,v]'*[N.x,N.y]);
115     cout << " iter " << n << " Err L2 = " << err << " outflux = " << outflux << endl;
116     if(err < 1e-3) break;
117 }
118 
119 // Verification
120 assert(abs(outflux)< 2e-3);
121 
122 // Plot
123 plot(p, wait=1, ps="NSprojP.eps");
124 plot(u, wait=1, ps="NSprojU.eps");

Rannacher’s projection algorithm: result on an adapted mesh, Fig. 22, showing the pressure, Fig. 23, and the horizontal velocity Fig. 24 for a Reynolds number of 400 where mesh adaptation is done after 50 iterations on the first mesh.

Table of content