Skip to content

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{eqnarray} \p_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{eqnarray}

A possible algorithm, proposed by Chorin, is :

\begin{eqnarray} {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 \p_n \mathbf{u}|_{\Gamma_{out}} &=0 \end{eqnarray}
\begin{eqnarray} -\Delta p^{m+1} &= -\nabla\cdot \mathbf{u}^mo\mathbf{X}^m &\\ \p_n p^{m+1} &= 0 &\mbox{ on } \Gamma\\ p^{m+1} &= 0 &\mbox{ on } \Gamma_{out} \end{eqnarray}

where \mathbf{u}o\mathbf{X}(x) = \mathbf{u}(\mathbf{x}-\mathbf{u}(\mathbf{x})\delta t) since \p _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 \p_n u = 0), to compute a correction, q, to the pressure.

and define

\begin{eqnarray} \mathbf{u}^{m+1} &= \tilde{\mathbf{u}} + P \n q\delta t\\ p^{m+1} &= p^m-q \end{eqnarray}

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 \p_n u = 0.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
// Parameters
verbosity = 0;
int nn = 1;
real nu = 0.0025;
real dt = 0.2;
real epsv = 1e-6;
real epsu = 1e-6;
real epsp = 1e-6;

// Mesh
border a0(t=1, 0){x=-2; y=t; label=1;}
border a1(t=-2, 0){x=t; y=0; label=2;}
border a2(t=0, -0.5){x=0; y=t; label=2;}
border a3(t=0, 1){x=18*t^1.2; y=-0.5; label=2;}
border a4(t=-0.5, 1){x=18; y=t; label=3;}
border a5(t=1, 0){x=-2+20*t; y=1; label=4;}

mesh Th = buildmesh(a0(3*nn) + a1(20*nn) + a2(10*nn) + a3(150*nn) + a4(5*nn) + a5(100*nn));
plot(Th);

// Fespace
fespace Vh(Th, P1);
Vh w;
Vh u=0, v=0;
Vh p=0;
Vh q=0;

// Definition of Matrix dtMx and dtMy
matrix dtM1x, dtM1y;

// Macro
macro BuildMat()
{   /* for memory managenemt */
    varf vM(unused, v) = int2d(Th)(v);
    varf vdx(u, v) = int2d(Th)(v*dx(u)*dt);
    varf vdy(u, v) = int2d(Th)(v*dy(u)*dt);

    real[int] Mlump = vM(0, Vh);
    real[int] one(Vh.ndof); one = 1;
    real[int] M1 = one ./ Mlump;
    matrix dM1 = M1;
    matrix Mdx = vdx(Vh, Vh);
    matrix Mdy = vdy(Vh, Vh);
    dtM1x = dM1*Mdx;
    dtM1y = dM1*Mdy;
} //

// Build matrices
BuildMat

// Time iterations
real err = 1.;
real outflux = 1.;
for(int n = 0; n < 300; n++){
    // Update
    Vh uold=u, vold=v, pold=p;

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

    plot(u);

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

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

    //to have absolute epsilon in CG algorithm.
    epsv = -abs(epsv);
    epsu = -abs(epsu);
    epsp = -abs(epsp);

    p = pold-q;
    u[] += dtM1x*q[];
    v[] += dtM1y*q[];

    // Mesh adaptation
    if (n%50 == 49){
        Th = adaptmesh(Th, [u, v], q, err=0.04, nbvx=100000);
        plot(Th, wait=true);
        BuildMat // Rebuild mat.
    }

    // Error & Outflux
    err = sqrt(int2d(Th)(square(u-uold)+square(v-vold))/Th.area);
    outflux = int1d(Th)([u,v]'*[N.x,N.y]);
    cout << " iter " << n << " Err L2 = " << err << " outflux = " << outflux << endl;
    if(err < 1e-3) break;
}

// Verification
assert(abs(outflux)< 2e-3);

// Plot
plot(p, wait=1, ps="NSprojP.eps");
plot(u, wait=1, ps="NSprojU.eps");

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

Fig. 1: Rannacher's projection algorithm: result on an adapted mesh (top) showing the pressure (middle) and the horizontal velocity \mathbf{u} at Reynolds 400.
Navier Strokes projection Th
Navier Strokes projection P
Navier Strokes projection U