Navier-Stokes equations
The Stokes equations are: for a given \(\mathbf{f}\in L^2(\Omega)^2\):
where \(\mathbf{u}=(u_1,u_2)\) is the velocity vector and \(p\) the pressure. For simplicity, let us choose Dirichlet boundary conditions on the velocity, \(\mathbf{u}=\mathbf{u}_{\Gamma}\) on \(\Gamma\).
In [TEMAM1977], Theorem 2.2, there is a weak form of (61):
Find \(\mathbf{v}=(v_1,v_2)\in \mathbf{V}(\Omega)\):
which satisfy:
Here it is used the existence \(p\in H^1(\Omega)\) such that \(\mathbf{u}=\nabla p\), if:
Another weak form is derived as follows: We put:
By multiplying the first equation in (61) with \(v\in V\) and the second with \(q\in W\), subsequent integration over \(\Omega\), and an application of Green’s formula, we have:
This yields the weak form of (61):
Find \((\mathbf{u},p)\in \mathbf{V}\times W\) such that:
for all \((\mathbf{v},q)\in V\times W\), where:
Now, we consider finite element spaces \(\mathbf{V}_h\subset \mathbf{V}\) and \(W_h\subset W\), and we assume the following basis functions:
The discrete weak form is: Find \((\mathbf{u}_{h},p_{h}) \in \mathbf{V}_{h} \times W_{h}\) such that:
Note
Assume that:
There is a constant \(\alpha_h>0\) such that:
\[a(\mathbf{v}_h,\mathbf{v}_h)\ge \alpha\| \mathbf{v}_h\|_{1,\Omega}^2\quad \textrm{for all }\mathbf{v}_h\in Z_h\]where:
\[Z_h=\{\mathbf{v}_h\in \mathbf{V}_h|\; b(\mathbf{w}_h,q_h)=0\quad \textrm{for all }q_h\in W_h\}\]There is a constant \(\beta_h>0\) such that:
\[\sup_{\mathbf{v}_h\in \mathbf{V}_h}\frac{b(\mathbf{v}_h,q_h)}{\| \mathbf{v}_h\|_{1,\Omega}} \ge \beta_h\| q_h\|_{0,\Omega}\quad \textrm{for all }q_h\in W_h\]Then we have an unique solution \((\mathbf{u}_h,p_h)\) of (62) satisfying:
\[\| \mathbf{u}-\mathbf{u}_h\|_{1,\Omega}+\| p-p_h\|_{0,\Omega} \le C\left( \inf_{\mathbf{v}_h\in \mathbf{V}_h}\| u-v_h\|_{1,\Omega} +\inf_{q_h\in W_h}\| p-q_h\|_{0,\Omega}\right)\]with a constant \(C>0\) (see e.g. [ROBERTS1993], Theorem 10.4).
Let us denote that:
then (62) is written by:
where:
Penalty method: This method consists of replacing (62) by a more regular problem:
Find \((\mathbf{v}_h^{\epsilon},p_h^{\epsilon})\in \mathbf{V}_h\times \tilde{W}_{h}\) satisfying:
where \(\tilde{W}_h\subset L^2(\Omega)\). Formally, we have:
and the corresponding algebraic problem:
Note
We can eliminate \(p_h^\epsilon=(1/\epsilon)BU_h^{\epsilon}\) to obtain:
Since the matrix \(A+(1/\epsilon)B^*B\) is symmetric, positive-definite, and sparse, (64) can be solved by known technique. There is a constant \(C>0\) independent of \(\epsilon\) such that:
(see e.g. [ROBERTS1993], 17.2)
Tip
Cavity
The driven cavity flow problem is solved first at zero Reynolds number (Stokes flow) and then at Reynolds 100. The velocity pressure formulation is used first and then the calculation is repeated with the stream function vorticity formulation.
We solve the driven cavity problem by the penalty method (63) where \(\mathbf{u}_{\Gamma}\cdot \mathbf{n}=0\) and \(\mathbf{u}_{\Gamma}\cdot \mathbf{s}=1\) on the top boundary and zero elsewhere (\(\mathbf{n}\) is the unit normal to \(\Gamma\), and \(\mathbf{s}\) the unit tangent to \(\Gamma\)).
The mesh is constructed by:
1 mesh Th = square(8, 8);
We use a classical Taylor-Hood element technique to solve the problem:
The velocity is approximated with the \(P_{2}\) FE (\(X_{h}\) space), and the pressure is approximated with the \(P_{1}\) FE (\(M_{h}\) space), where:
and:
The FE spaces and functions are constructed by:
1 fespace Xh(Th, P2); //definition of the velocity component space
2 fespace Mh(Th, P1); //definition of the pressure space
3 Xh u2, v2;
4 Xh u1, v1;
5 Mh p, q;
The Stokes operator is implemented as a system-solve for the velocity \((u1,u2)\) and the pressure \(p\). The test function for the velocity is \((v1,v2)\) and \(q\) for the pressure, so the variational form (62) in freefem language is:
1 solve Stokes (u1, u2, p, v1, v2, q, solver=Crout)
2 = int2d(Th)(
3 (
4 dx(u1)*dx(v1)
5 + dy(u1)*dy(v1)
6 + dx(u2)*dx(v2)
7 + dy(u2)*dy(v2)
8 )
9 - p*q*(0.000001)
10 - p*dx(v1) - p*dy(v2)
11 - dx(u1)*q - dy(u2)*q
12 )
13 + on(3, u1=1, u2=0)
14 + on(1, 2, 4, u1=0, u2=0)
15 ;
Each unknown has its own boundary conditions.
If the streamlines are required, they can be computed by finding \(\psi\) such that \(rot\psi=u\) or better:
1 Xh psi, phi;
2
3 solve streamlines (psi, phi)
4 = int2d(Th)(
5 dx(psi)*dx(phi)
6 + dy(psi)*dy(phi)
7 )
8 + int2d(Th)(
9 - phi*(dy(u1) - dx(u2))
10 )
11 + on(1, 2, 3, 4, psi=0)
12 ;
Now the Navier-Stokes equations are solved:
with the same boundary conditions and with initial conditions \(u=0\).
This is implemented by using the convection operator convect
for the term \({\p u\over\p t} +u\cdot\nabla u\), giving a discretization in time
The term \(u^n\circ X^n(x)\approx u^n(x-u^n(x)\tau )\) will be computed by the operator convect
, so we obtain:
1 int i=0;
2 real alpha=1/dt;
3 problem NS (u1, u2, p, v1, v2, q, solver=Crout, init=i)
4 = int2d(Th)(
5 alpha*(u1*v1 + u2*v2)
6 + nu * (
7 dx(u1)*dx(v1) + dy(u1)*dy(v1)
8 + dx(u2)*dx(v2) + dy(u2)*dy(v2)
9 )
10 - p*q*(0.000001)
11 - p*dx(v1) - p*dy(v2)
12 - dx(u1)*q - dy(u2)*q
13 )
14 + int2d(Th)(
15 - alpha*convect([up1,up2],-dt,up1)*v1
16 - alpha*convect([up1,up2],-dt,up2)*v2
17 )
18 + on(3, u1=1, u2=0)
19 + on(1, 2, 4,u1=0, u2=0)
20 ;
21
22 // Time loop
23 for (i = 0; i <= 10; i++){
24 // Update
25 up1 = u1;
26 up2 = u2;
27
28 // Solve
29 NS;
30
31 // Plot
32 if (!(i % 10))
33 plot(coef=0.2, cmm="[u1,u2] and p", p, [u1, u2]);
34 }
Notice that the stiffness matrices are reused (keyword init=i
)
The complete script is available in cavity example.
Uzawa Algorithm and Conjugate Gradients
We solve Stokes problem without penalty. The classical iterative method of Uzawa is described by the algorithm (see e.g. [ROBERTS1993], 17.3, [GLOWINSKI1979], 13 or [GLOWINSKI1985], 13):
Initialize: Let \(p_h^0\) be an arbitrary chosen element of \(L^2(\Omega)\).
Calculate :math:`mathbf{u}_h`: Once \(p_h^n\) is known, \(\mathbf{v}_h^n\) is the solution of:
\[\mathbf{u}_h^n = A^{-1}(\mathbf{f}_h-\mathbf{B}^*p_h^n)\]Advance :math:`p_h`: Let \(p_h^{n+1}\) be defined by;
\[p_h^{n+1}=p_h^n+\rho_n\mathbf{B}\mathbf{u}_h^n\]
There is a constant \(\alpha>0\) such that \(\alpha\le \rho_n\le 2\) for each \(n\), then \(\mathbf{u}_h^n\) converges to the solution \(\mathbf{u}_h\), and then \(B\mathbf{v}_h^n\to 0\) as \(n\to \infty\) from the Advance \(p_h\). This method in general converges quite slowly.
First we define mesh, and the Taylor-Hood approximation. So \(X_{h}\) is the velocity space, and \(M_{h}\) is the pressure space.
Tip
Stokes Uzawa
1 // Mesh
2 mesh Th = square(10, 10);
3
4 // Fespace
5 fespace Xh(Th, P2);
6 Xh u1, u2;
7 Xh bc1, bc2;
8 Xh b;
9
10 fespace Mh(Th, P1);
11 Mh p;
12 Mh ppp; //ppp is a working pressure
13
14 // Problem
15 varf bx (u1, q) = int2d(Th)(-(dx(u1)*q));
16 varf by (u1, q) = int2d(Th)(-(dy(u1)*q));
17 varf a (u1, u2)
18 = int2d(Th)(
19 dx(u1)*dx(u2)
20 + dy(u1)*dy(u2)
21 )
22 + on(3, u1=1)
23 + on(1, 2, 4, u1=0) ;
24 //remark: put the on(3,u1=1) before on(1,2,4,u1=0)
25 //because we want zero on intersection
26
27 matrix A = a(Xh, Xh, solver=CG);
28 matrix Bx = bx(Xh, Mh); //B=(Bx, By)
29 matrix By = by(Xh, Mh);
30
31 bc1[] = a(0,Xh); //boundary condition contribution on u1
32 bc2 = 0; //no boundary condition contribution on u2
33
34 //p_h^n -> B A^-1 - B^* p_h^n = -div u_h
35 //is realized as the function divup
36 func real[int] divup (real[int] & pp){
37 //compute u1(pp)
38 b[] = Bx'*pp;
39 b[] *= -1;
40 b[] += bc1[];
41 u1[] = A^-1*b[];
42 //compute u2(pp)
43 b[] = By'*pp;
44 b[] *= -1;
45 b[] += bc2[];
46 u2[] = A^-1*b[];
47 //u^n = (A^-1 Bx^T p^n, By^T p^n)^T
48 ppp[] = Bx*u1[]; //ppp = Bx u_1
49 ppp[] += By*u2[]; //+ By u_2
50
51 return ppp[] ;
52 }
53
54 // Initialization
55 p=0; //p_h^0 = 0
56 LinearCG(divup, p[], eps=1.e-6, nbiter=50); //p_h^{n+1} = p_h^n + B u_h^n
57 // if n> 50 or |p_h^{n+1} - p_h^n| <= 10^-6, then the loop end
58 divup(p[]); //compute the final solution
59
60 plot([u1, u2], p, wait=1, value=true, coef=0.1);
NSUzawaCahouetChabart.edp
In this example we solve the Navier-Stokes equation past a cylinder with the Uzawa algorithm preconditioned by the Cahouet-Chabart method (see [GLOWINSKI2003] for all the details).
The idea of the preconditioner is that in a periodic domain, all differential operators commute and the Uzawa algorithm comes to solving the linear operator \(\nabla. ( (\alpha Id + \nu \Delta)^{-1} \nabla\), where \(Id\) is the identity operator. So the preconditioner suggested is \(\alpha \Delta^{-1} + \nu Id\).
To implement this, we do:
Tip
NS Uzawa Cahouet Chabart
1 // Parameters
2 verbosity = 0;
3 real D = 0.1;
4 real H = 0.41;
5 real cx0 = 0.2;
6 real cy0 = 0.2; //center of cylinder
7 real xa = 0.15;
8 real ya = 0.2;
9 real xe = 0.25;
10 real ye = 0.2;
11 int nn = 15;
12
13 //TODO
14 real Um = 1.5; //max velocity (Rey 100)
15 real nu = 1e-3;
16
17 func U1 = 4.*Um*y*(H-y)/(H*H); //Boundary condition
18 func U2 = 0.;
19 real T=2;
20 real dt = D/nn/Um; //CFL = 1
21 real epspq = 1e-10;
22 real eps = 1e-6;
23
24 // Variables
25 func Ub = Um*2./3.;
26 real alpha = 1/dt;
27 real Rey = Ub*D/nu;
28 real t = 0.;
29
30 // Mesh
31 border fr1(t=0, 2.2){x=t; y=0; label=1;}
32 border fr2(t=0, H){x=2.2; y=t; label=2;}
33 border fr3(t=2.2, 0){x=t; y=H; label=1;}
34 border fr4(t=H, 0){x=0; y=t; label=1;}
35 border fr5(t=2*pi, 0){x=cx0+D*sin(t)/2; y=cy0+D*cos(t)/2; label=3;}
36 mesh Th = buildmesh(fr1(5*nn) + fr2(nn) + fr3(5*nn) + fr4(nn) + fr5(-nn*3));
37
38 // Fespace
39 fespace Mh(Th, [P1]);
40 Mh p;
41
42 fespace Xh(Th, [P2]);
43 Xh u1, u2;
44
45 fespace Wh(Th, [P1dc]);
46 Wh w; //vorticity
47
48 // Macro
49 macro grad(u) [dx(u), dy(u)] //
50 macro div(u1, u2) (dx(u1) + dy(u2)) //
51
52 // Problem
53 varf von1 ([u1, u2, p], [v1, v2, q])
54 = on(3, u1=0, u2=0)
55 + on(1, u1=U1, u2=U2)
56 ;
57
58 //remark : the value 100 in next varf is manualy fitted, because free outlet.
59 varf vA (p, q) =
60 int2d(Th)(
61 grad(p)' * grad(q)
62 )
63 + int1d(Th, 2)(
64 100*p*q
65 )
66 ;
67
68 varf vM (p, q)
69 = int2d(Th, qft=qf2pT)(
70 p*q
71 )
72 + on(2, p=0)
73 ;
74
75 varf vu ([u1], [v1])
76 = int2d(Th)(
77 alpha*(u1*v1)
78 + nu*(grad(u1)' * grad(v1))
79 )
80 + on(1, 3, u1=0)
81 ;
82
83 varf vu1 ([p], [v1]) = int2d(Th)(p*dx(v1));
84 varf vu2 ([p], [v1]) = int2d(Th)(p*dy(v1));
85
86 varf vonu1 ([u1], [v1]) = on(1, u1=U1) + on(3, u1=0);
87 varf vonu2 ([u1], [v1]) = on(1, u1=U2) + on(3, u1=0);
88
89 matrix pAM = vM(Mh, Mh, solver=UMFPACK);
90 matrix pAA = vA(Mh, Mh, solver=UMFPACK);
91 matrix AU = vu(Xh, Xh, solver=UMFPACK);
92 matrix B1 = vu1(Mh, Xh);
93 matrix B2 = vu2(Mh, Xh);
94
95 real[int] brhs1 = vonu1(0, Xh);
96 real[int] brhs2 = vonu2(0, Xh);
97
98 varf vrhs1(uu, vv) = int2d(Th)(convect([u1, u2], -dt, u1)*vv*alpha) + vonu1;
99 varf vrhs2(v2, v1) = int2d(Th)(convect([u1, u2], -dt, u2)*v1*alpha) + vonu2;
100
101 // Uzawa function
102 func real[int] JUzawa (real[int] & pp){
103 real[int] b1 = brhs1; b1 += B1*pp;
104 real[int] b2 = brhs2; b2 += B2*pp;
105 u1[] = AU^-1 * b1;
106 u2[] = AU^-1 * b2;
107 pp = B1'*u1[];
108 pp += B2'*u2[];
109 pp = -pp;
110 return pp;
111 }
112
113 // Preconditioner function
114 func real[int] Precon (real[int] & p){
115 real[int] pa = pAA^-1*p;
116 real[int] pm = pAM^-1*p;
117 real[int] pp = alpha*pa + nu*pm;
118 return pp;
119 }
120
121 // Initialization
122 p = 0;
123
124 // Time loop
125 int ndt = T/dt;
126 for(int i = 0; i < ndt; ++i){
127 // Update
128 brhs1 = vrhs1(0, Xh);
129 brhs2 = vrhs2(0, Xh);
130
131 // Solve
132 int res = LinearCG(JUzawa, p[], precon=Precon, nbiter=100, verbosity=10, veps=eps);
133 assert(res==1);
134 eps = -abs(eps);
135
136 // Vorticity
137 w = -dy(u1) + dx(u2);
138 plot(w, fill=true, wait=0, nbiso=40);
139
140 // Update
141 dt = min(dt, T-t);
142 t += dt;
143 if(dt < 1e-10*T) break;
144 }
145
146 // Plot
147 plot(w, fill=true, nbiso=40);
148
149 // Display
150 cout << "u1 max = " << u1[].linfty
151 << ", u2 max = " << u2[].linfty
152 << ", p max = " << p[].max << endl;
Warning
Stop test of the conjugate gradient
Because we start from the previous solution and the end the previous solution is close to the final solution, don’t take a relative stop test to the first residual, take an absolute stop test (negative here).

Fig. 164 The vorticity at Reynolds number 100 a time 2s with the Cahouet-Chabart method.