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
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