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
  2verbosity = 0;
  3int nn = 1;
  4real nu = 0.0025;
  5real dt = 0.2;
  6real epsv = 1e-6;
  7real epsu = 1e-6;
  8real epsp = 1e-6;
  9
 10// Mesh
 11border a0(t=1, 0){x=-2; y=t; label=1;}
 12border a1(t=-2, 0){x=t; y=0; label=2;}
 13border a2(t=0, -0.5){x=0; y=t; label=2;}
 14border a3(t=0, 1){x=18*t^1.2; y=-0.5; label=2;}
 15border a4(t=-0.5, 1){x=18; y=t; label=3;}
 16border a5(t=1, 0){x=-2+20*t; y=1; label=4;}
 17
 18mesh Th = buildmesh(a0(3*nn) + a1(20*nn) + a2(10*nn) + a3(150*nn) + a4(5*nn) + a5(100*nn));
 19plot(Th);
 20
 21// Fespace
 22fespace Vh(Th, P1);
 23Vh w;
 24Vh u=0, v=0;
 25Vh p=0;
 26Vh q=0;
 27
 28// Definition of Matrix dtMx and dtMy
 29matrix dtM1x, dtM1y;
 30
 31// Macro
 32macro 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
 49BuildMat
 50
 51// Time iterations
 52real err = 1.;
 53real outflux = 1.;
 54for(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
120assert(abs(outflux)< 2e-3);
121
122// Plot
123plot(p, wait=1, ps="NSprojP.eps");
124plot(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