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