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:
A possible algorithm, proposed by Chorin, is:
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.
and define
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.