# Navier-Stokes equations

The Stokes equations are: for a given $\mathbf{f}\in L^2(\Omega)^2$,

$$\left.\begin{array}{cl} \label{eqn::Stokes} -\Delta \mathbf{u}+\nabla p & =\mathbf{f} \\ \nabla\cdot \mathbf{u} &=0 \end{array}\right\}\quad \hbox{ in }\Omega$$

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 \eqref{eqn::Stokes}:

Find $\mathbf{v}=(v_1,v_2)\in \mathbf{V}(\Omega)$

\mathbf{V}(\Omega)=\{\mathbf{w}\in H^1_0(\Omega)^2|\; \textrm{div}\mathbf{w}=0\}

which satisfy

\sum_{i=1}^2\int_{\Omega}\nabla u_i\cdot \nabla v_i=\int_{\Omega}\mathbf{f}\cdot \mathbf{w} \quad \textrm{for all }v\in V

Here it is used the existence $p\in H^1(\Omega)$ such that $\mathbf{u}=\nabla p$, if

\int_{\Omega}\mathbf{u}\cdot \mathbf{v}=0\quad \textrm{for all }\mathbf{v}\in V

Another weak form is derived as follows: We put

\begin{eqnarray*} \mathbf{V}=H^1_0(\Omega)^2;\quad W=\left\{q\in L^2(\Omega)\left|\; \int_{\Omega}q=0\right.\right\} \end{eqnarray*}

By multiplying the first equation in \eqref{eqn::Stokes} with $v\in V$ and the second with $q\in W$, subsequent integration over $\Omega$, and an application of Green's formula, we have

\begin{eqnarray*} \int_{\Omega}\nabla\mathbf{u}\cdot \nabla\mathbf{v}-\int_{\Omega}\textrm{div}\mathbf{v}\, p &=&\int_{\Omega}\mathbf{f}\cdot\mathbf{v}\\ \int_{\Omega}\textrm{div}\mathbf{u}\, q&=&0 \end{eqnarray*}

This yields the weak form of \eqref{eqn::Stokes}:

Find $(\mathbf{u},p)\in \mathbf{V}\times W$ such that

\begin{eqnarray} a(\mathbf{u},\mathbf{v})+b(\mathbf{v},p)&=&(\mathbf{f},\mathbf{v})\\ b(\mathbf{u},q)&=&0 \end{eqnarray}

for all $(\mathbf{v},q)\in V\times W$, where

\begin{eqnarray} a(\mathbf{u},\mathbf{v})&=&\int_{\Omega}\nabla \mathbf{u}\cdot \nabla\mathbf{v} =\sum_{i=1}^2\int_{\Omega}\nabla u_i\cdot \nabla v_i\\ b(\mathbf{u},q)&=&-\int_{\Omega}\textrm{div}\mathbf{u}\, q \end{eqnarray}

Now, we consider finite element spaces $\mathbf{V}_h\subset \mathbf{V}$ and $W_h\subset W$, and we assume the following basis functions

\begin{eqnarray*} &&\mathbf{V}_h=V_h\times V_h,\quad V_h=\{v_h|\; v_h=v_1\phi_1+\cdots +v_{M_V}\phi_{M_V}\},\\ &&W_h=\{q_h|\; q_h=q_1\varphi_1+\cdots +q_{M_W}\varphi_{M_W}\} \end{eqnarray*}

The discrete weak form is: Find $(\mathbf{u}_{h},p_{h}) \in \mathbf{V}_{h} \times W_{h}$ such that

$$\label{eqn::vfStokes} \begin{array}{cll} a(\mathbf{u}_h,\mathbf{v}_h)+b(\mathbf{v}_h,p) &= (\mathbf{f},\mathbf{v}_h) , &\forall \mathbf{v}_{h} \in \mathbf{V}_{h} \\ b(\mathbf{u}_h,q_h)&= 0, &\forall q_{h} \in W_{h} \end{array}$$

Note

Assume that:

1. 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\}
2. 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 \eqref{eqn::vfStokes} 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

\begin{eqnarray} A&=&(A_{ij}),\, A_{ij}=\int_{\Omega}\nabla \phi_j\cdot \nabla \phi_i\qquad i,j=1,\cdots,M_{\mathbf{V}}\\ \mathbf{B}&=&(Bx_{ij},By_{ij}),\, Bx_{ij}=-\int_{\Omega}\p \phi_j/\p x\, \varphi_i\qquad By_{ij}=-\int_{\Omega}\p \phi_j/\p y\, \varphi_i\nonumber\\ &&\qquad i=1,\cdots,M_W;j=1,\cdots,M_V\nonumber \end{eqnarray}

then \eqref{eqn::vfStokes} is written by

where

Penalty method: This method consists of replacing \eqref{eqn::vfStokes} by a more regular problem: Find $(\mathbf{v}_h^{\epsilon},p_h^{\epsilon})\in \mathbf{V}_h\times \tilde{W}_{h}$ satisfying

$$\begin{array}{cll} a(\mathbf{u}_h^\epsilon,\mathbf{v}_h)+b(\mathbf{v}_h,p_h^{\epsilon}) &= (\mathbf{f},\mathbf{v}_h) , &\forall \mathbf{v}_{h} \in \mathbf{V}_{h} \\ b(\mathbf{u}_h^{\epsilon},q_h)-\epsilon(p_h^{\epsilon},q_h)&= 0, &\forall q_{h} \in \tilde{W}_{h} \end{array}$$

where $\tilde{W}_h\subset L^2(\Omega)$. Formally, we have

\textrm{div}\mathbf{u}_h^{\epsilon}=\epsilon p_h^{\epsilon}

and the corresponding algebraic problem

\begin{eqnarray*} \left( \begin{array}{cc} \mathbf{A}&B^*\\ B&-\epsilon I \end{array} \right) \left( \begin{array}{cc} \mathbf{U}_h^{\epsilon}\\ \{p_h^{\epsilon}\} \end{array} \right) = \left( \begin{array}{cc} \mathbf{F}_h\\ 0 \end{array} \right) \end{eqnarray*}

Note

We can eliminate $p_h^\epsilon=(1/\epsilon)BU_h^{\epsilon}$ to obtain

\begin{eqnarray} \label{eqn::StiffPvfStokes} (A+(1/\epsilon)B^*B)\mathbf{U}_h^{\epsilon}=\mathbf{F}_h^{\epsilon} \end{eqnarray}

Since the matrix $A+(1/\epsilon)B^*B$ is symmetric, positive-definite, and sparse, \eqref{eqn::StiffPvfStokes} can be solved by known technique. There is a constant $C>0$ independent of $\epsilon$ such that

\|\mathbf{u}_h-\mathbf{u}_h^\epsilon\|_{1,\Omega}+ \|p_h-p_h^{\epsilon}\|_{0,\Omega}\le C\epsilon

(see e.g. ROBERTS1993, 17.2)

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 \eqref{eqn::PvfStokes} 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

X_{h} = \left\{ \mathbf{v} \in H^{1}(]0,1[^2) \left|\; \forall K \in \mathcal{T}_{h}\quad v_{|K} \in P_{2}\right.\right\}

and

M_{h} = \left\{ v \in H^{1}(]0,1[^2) \left|\; \forall K \in \mathcal{T}_{h}\quad v_{|K} \in P_{1} \right.\right\}

The FE spaces and functions are constructed by

 1 2 3 4 5 fespace Xh(Th, P2); //definition of the velocity component space fespace Mh(Th, P1); //definition of the pressure space Xh u2, v2; Xh u1, v1; 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 \eqref{eqn::vfStokes} in freefem language is:

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 solve Stokes (u1, u2, p, v1, v2, q, solver=Crout) = int2d(Th)( ( dx(u1)*dx(v1) + dy(u1)*dy(v1) + dx(u2)*dx(v2) + dy(u2)*dy(v2) ) - p*q*(0.000001) - p*dx(v1) - p*dy(v2) - dx(u1)*q - dy(u2)*q ) + on(3, u1=1, u2=0) + on(1, 2, 4, u1=0, u2=0) ; 

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,

-\Delta\psi=\nabla\times u
  1 2 3 4 5 6 7 8 9 10 11 12 Xh psi, phi; solve streamlines (psi, phi) = int2d(Th)( dx(psi)*dx(phi) + dy(psi)*dy(phi) ) + int2d(Th)( - phi*(dy(u1) - dx(u2)) ) + on(1, 2, 3, 4, psi=0) ; 

Now the Navier-Stokes equations are solved

{\p {u}\over\p t} +u\cdot\nabla u-\nu \Delta u+\nabla p=0, \nabla\cdot u=0

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

$$\begin{array}{cl} \frac{1}{\tau } (u^{n+1}-u^n\circ X^n) -\nu\Delta u^{n+1} + \nabla p^{n+1} &=0,\\ \nabla\cdot u^{n+1} &= 0 \end{array}$$

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 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 int i=0; real alpha=1/dt; problem NS (u1, u2, p, v1, v2, q, solver=Crout, init=i) = int2d(Th)( alpha*(u1*v1 + u2*v2) + nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1) + dx(u2)*dx(v2) + dy(u2)*dy(v2) ) - p*q*(0.000001) - p*dx(v1) - p*dy(v2) - dx(u1)*q - dy(u2)*q ) + int2d(Th)( - alpha*convect([up1,up2],-dt,up1)*v1 - alpha*convect([up1,up2],-dt,up2)*v2 ) + on(3, u1=1, u2=0) + on(1, 2, 4,u1=0, u2=0) ; // Time loop for (i = 0; i <= 10; i++){ // Update up1 = u1; up2 = u2; // Solve NS; // Plot if (!(i % 10)) plot(coef=0.2, cmm="[u1,u2] and p", p, [u1, u2]); } 

Notice that the stiffness matrices are reused (keyword init=i)

The complete script is available in Examples.

## 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 $\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 $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.

Stokes Uzawa

  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 // Mesh mesh Th = square(10, 10); // Fespace fespace Xh(Th, P2); Xh u1, u2; Xh bc1, bc2; Xh b; fespace Mh(Th, P1); Mh p; Mh ppp; //ppp is a working pressure // Problem varf bx (u1, q) = int2d(Th)(-(dx(u1)*q)); varf by (u1, q) = int2d(Th)(-(dy(u1)*q)); varf a (u1, u2) = int2d(Th)( dx(u1)*dx(u2) + dy(u1)*dy(u2) ) + on(3, u1=1) + on(1, 2, 4, u1=0) ; //remark: put the on(3,u1=1) before on(1,2,4,u1=0) //because we want zero on intersection matrix A = a(Xh, Xh, solver=CG); matrix Bx = bx(Xh, Mh); //B=(Bx, By)$matrix By = by(Xh, Mh); bc1[] = a(0,Xh); //boundary condition contribution on u1 bc2 = 0; //no boundary condition contribution on u2 //p_h^n -> B A^-1 - B^* p_h^n = -div u_h //is realized as the function divup func real[int] divup (real[int] & pp){ //compute u1(pp) b[] = Bx'*pp; b[] *= -1; b[] += bc1[]; u1[] = A^-1*b[]; //compute u2(pp) b[] = By'*pp; b[] *= -1; b[] += bc2[]; u2[] = A^-1*b[]; //u^n = (A^-1 Bx^T p^n, By^T p^n)^T$ ppp[] = Bx*u1[]; //ppp = Bx u_1 ppp[] += By*u2[]; //+ By u_2 return ppp[] ; } // Initialization p=0; //p_h^0 = 0 LinearCG(divup, p[], eps=1.e-6, nbiter=50); //p_h^{n+1} = p_h^n + B u_h^n // if n> 50 or |p_h^{n+1} - p_h^n| <= 10^-6, then the loop end divup(p[]); //compute the final solution 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

NS Uzawa Cahouet Chabart

  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 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 // Parameters verbosity = 0; real D = 0.1; real H = 0.41; real cx0 = 0.2; real cy0 = 0.2; //center of cylinder real xa = 0.15; real ya = 0.2; real xe = 0.25; real ye = 0.2; int nn = 15; //TODO real Um = 1.5; //max velocity (Rey 100) real nu = 1e-3; func U1 = 4.*Um*y*(H-y)/(H*H); //Boundary condition func U2 = 0.; real T=2; real dt = D/nn/Um; //CFL = 1 real epspq = 1e-10; real eps = 1e-6; // Variables func Ub = Um*2./3.; real alpha = 1/dt; real Rey = Ub*D/nu; real t = 0.; // Mesh border fr1(t=0, 2.2){x=t; y=0; label=1;} border fr2(t=0, H){x=2.2; y=t; label=2;} border fr3(t=2.2, 0){x=t; y=H; label=1;} border fr4(t=H, 0){x=0; y=t; label=1;} border fr5(t=2*pi, 0){x=cx0+D*sin(t)/2; y=cy0+D*cos(t)/2; label=3;} mesh Th = buildmesh(fr1(5*nn) + fr2(nn) + fr3(5*nn) + fr4(nn) + fr5(-nn*3)); // Fespace fespace Mh(Th, [P1]); Mh p; fespace Xh(Th, [P2]); Xh u1, u2; fespace Wh(Th, [P1dc]); Wh w; //vorticity // Macro macro grad(u) [dx(u), dy(u)] // macro div(u1, u2) (dx(u1) + dy(u2)) // // Problem varf von1 ([u1, u2, p], [v1, v2, q]) = on(3, u1=0, u2=0) + on(1, u1=U1, u2=U2) ; //remark : the value 100 in next varf is manualy fitted, because free outlet. varf vA (p, q) = int2d(Th)( grad(p)' * grad(q) ) + int1d(Th, 2)( 100*p*q ) ; varf vM (p, q) = int2d(Th, qft=qf2pT)( p*q ) + on(2, p=0) ; varf vu ([u1], [v1]) = int2d(Th)( alpha*(u1*v1) + nu*(grad(u1)' * grad(v1)) ) + on(1, 3, u1=0) ; varf vu1 ([p], [v1]) = int2d(Th)(p*dx(v1)); varf vu2 ([p], [v1]) = int2d(Th)(p*dy(v1)); varf vonu1 ([u1], [v1]) = on(1, u1=U1) + on(3, u1=0); varf vonu2 ([u1], [v1]) = on(1, u1=U2) + on(3, u1=0); matrix pAM = vM(Mh, Mh, solver=UMFPACK); matrix pAA = vA(Mh, Mh, solver=UMFPACK); matrix AU = vu(Xh, Xh, solver=UMFPACK); matrix B1 = vu1(Mh, Xh); matrix B2 = vu2(Mh, Xh); real[int] brhs1 = vonu1(0, Xh); real[int] brhs2 = vonu2(0, Xh); varf vrhs1(uu, vv) = int2d(Th)(convect([u1, u2], -dt, u1)*vv*alpha) + vonu1; varf vrhs2(v2, v1) = int2d(Th)(convect([u1, u2], -dt, u2)*v1*alpha) + vonu2; // Uzawa function func real[int] JUzawa (real[int] & pp){ real[int] b1 = brhs1; b1 += B1*pp; real[int] b2 = brhs2; b2 += B2*pp; u1[] = AU^-1 * b1; u2[] = AU^-1 * b2; pp = B1'*u1[]; pp += B2'*u2[]; pp = -pp; return pp; } // Preconditioner function func real[int] Precon (real[int] & p){ real[int] pa = pAA^-1*p; real[int] pm = pAM^-1*p; real[int] pp = alpha*pa + nu*pm; return pp; } // Initialization p = 0; // Time loop int ndt = T/dt; for(int i = 0; i < ndt; ++i){ // Update brhs1 = vrhs1(0, Xh); brhs2 = vrhs2(0, Xh); // Solve int res = LinearCG(JUzawa, p[], precon=Precon, nbiter=100, verbosity=10, veps=eps); assert(res==1); eps = -abs(eps); // Vorticity w = -dy(u1) + dx(u2); plot(w, fill=true, wait=0, nbiso=40); // Update dt = min(dt, T-t); t += dt; if(dt < 1e-10*T) break; } // Plot plot(w, fill=true, nbiso=40); // Display cout << "u1 max = " << u1[].linfty << ", u2 max = " << u2[].linfty << ", p max = " << p[].max << endl; 

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. 24: The vorticity at Reynolds number 100 a time 2s with the Cahouet-Chabart method.

## References#

[TEMAM1977] TEMAM, Roger. Navier-Stokes equations: theory and numerical analysis. 1977.

[ROBERTS1993] ROBERTS, J. E. et THOMAS, J. M. Mixed and Hybrid Methods, Handbook of Numerical Anaysis, Vol. II. North-Holland, 1993, vol. 183, p. 184.

[GLOWINSKI1979] GLOWINSKI, R. et PIRONNEAU, O. On numerical methods for the Stokes problem. In: Energy methods in finite element analysis.(A79-53076 24-39) Chichester, Sussex, England, Wiley-Interscience, 1979, p. 243-264., 1979, p. 243-264.

[GLOWINSKI1985] GLOWINSKI, Roland et ODEN, J. Tinsley. Numerical methods for nonlinear variational problems. Journal of Applied Mechanics, 1985, vol. 52, p. 739.

[GLOWINSKI2003] GLOWINSKI, Roland. Finite element methods for incompressible viscous flow. Handbook of numerical analysis, 2003, vol. 9, p. 3-1176.