$\begin{eqnarray} -\nu\Delta\mathbf{u} + \nabla p &= \textbf{f} &\mbox{ in }\Omega\label{eqstokes:1}\\ % \nabla\cdot\mathbf{u} &= 0 &\mbox{ in }\Omega\label{eqstokes:2}\\ % \mathbf{u} &= \textbf{g}_1 &\mbox{ on }\Gamma_D\label{eqstokes:3}\\ % \nu\frac{\partial \mathbf{u}}{\partial \overrightarrow{n}} + p \cdot \overrightarrow{n} &= \textbf{g}_2 &\mbox{ on }\Gamma_N\label{eqstokes:4} \end{eqnarray}$

Let $\Omega\subset\mathbb{R}^3$ be a bounded domain with a Lipschitz-continuous boundary $\Gamma$. Suppose that $\Gamma$ consists of two measurable parts $\Gamma_D$ and $\Gamma_N$ where a Dirichlet and a Neumann condition are imposed respectively. Flows at low Reynold’s number satisfy the Stokes problem:

Where $\mathbf{u}=(u_1, u_2, u_3)$ is the fluid velocity, $p$ its pressure, $\nu$ its kinematic viscosity and $\mathbf{f}$ an applied force.

The differents operators are defined (in dimension $N$) by:

$\Delta\mathbf{u} = \sum_{i=1}^{N}{\frac{\partial^2\mathbf{u}}{\partial x_i^2}}$
$\nabla p = \left(\frac{\partial p}{\partial x_1}, \ldots, \frac{\partial p}{\partial x_N}\right)$

Variational form

We set $V=H_0^1(\Omega)^3=\left\{v\in H^1(\Omega)^3 | v_{|\Gamma_D}=0\right\}$ and $Q=L_0^2(\Omega)=\left\{q\in L^2(\Omega) | \int_{\Omega}{q}=0\right\}$. We then multiply the equation $\eqref{eqstokes:1}$ by $\mathbf{v}\in V$ and the equation $\eqref{eqstokes:2}$ by $q\in Q$, then integrate the equations over $\Omega$ , it holds:

$\begin{eqnarray*} - \mu\int_{\Omega}{\Delta\mathbf{u}\cdot\mathbf{v}} + \int_{\Omega}{\nabla p\cdot\mathbf{v}} &=& \int_{\Omega}{\textbf{f}\cdot\mathbf{v}}\\ \int_{\Omega}{\mathrm{div}(\mathbf{u})q} &=& 0 \end{eqnarray*}$

Applying the Green’s formula, we get:

$\begin{eqnarray*} \mu\int_{\Omega}{\nabla\mathbf{u}\cdot\nabla\mathbf{v}} - \int_{\Omega}{p\mathrm{div}(\mathbf{v})} - \int_{\Gamma}{\left(\mu\frac{\partial\mathbf{u}}{\partial\mathbf{n}}-p\mathbf{n}\right)\cdot\mathbf{v}} &=& \int_{\Omega}{\textbf{f}\cdot\mathbf{v}}\\ \int_{\Omega}{\mathrm{div}(\mathbf{u})q} &=& 0 \end{eqnarray*}$

Since $\mathbf{v}$ belongs to $V$, the variational form becomes: Find $\mathbf{u}\in V$, $p\in Q$ such as:

$\begin{eqnarray*} \mu\int_{\Omega}{\nabla\mathbf{u}\cdot\nabla\mathbf{v}} - \int_{\Omega}{p\mathrm{div}(\mathbf{v})} - \int_{\Gamma_N}{\textbf{g}_2\cdot\mathbf{v}} &=& \int_{\Omega}{\textbf{f}\cdot\mathbf{v}}\\ \int_{\Omega}{\mathrm{div}(\mathbf{u})q} &=& 0 \end{eqnarray*}$

Taylor-Hood spaces

Let $V_h\subset V$ and $Q_h\subset Q$ be subspaces of inite dimensions such as $V_h\rightarrow V$ and $Q_h\rightarrow Q$ when $h\rightarrow 0$. We consider the approximated Stokes problem:

Find $\mathbf{u}_h\in V_h$ and $p_h\in Q_h$ sush as:

$\begin{eqnarray*} \mu\int_{\Omega}{\nabla\mathbf{u}_h\cdot\nabla\mathbf{v}_h} - \int_{\Omega}{p_h\mathrm{div}(\mathbf{v}_h)} - \int_{\Gamma_N}{\textbf{g}_2\cdot\mathbf{v}_h} &=& \int_{\Omega}{\textbf{f}\cdot\mathbf{v}_h}\\ \int_{\Omega}{\mathrm{div}(\mathbf{u}_h)q_h} &=& 0 \end{eqnarray*}$

This problem is well-defined if the spaces Vh and Qh verify the inf-sup condition:

$\exists\beta_h>0, \inf_{q_h\in Q_h}{\sup_{v_h\in Vh}{\frac{\int_{\Omega}{q_h\mathrm{div}(\mathbf{v}_h)}}{\Vert q_h\Vert_Q\Vert\mathbf{v}_h\Vert_V}}}\geq\beta_h$

We can note that Stokes equations only use secoind derivatives for the velocity $u$ and first derivative for the pressure $p$, so it makes sense to use $k\geq 1$ polynomial degree for $V_h$ and $k−1$ degree for $Q_h$. Here, we put $k=2$.

A regular triangulation $Th$ of $\Omega$, such as $\bar{\Omega}=\bigcup_{K\in T_h}{K}$, verifies:

1. Each $K$ is a polyhedron such as $int(K)\ne 0$

2. $\text{int}(K_1)\cap \text{int}(K_2)=\emptyset$ for each $K_1$ distincts to $K_2$ .

3. if $E=K_1\cap K_2$, then $E$ is either a shared face, a shared edge or a share vertex.

Figure 1. Example of triangulation

Finally, we define the Taylor-Hood space in the following way:

$V_h = \left\{X_h^2\cup H_0^1(\Omega)\right\}^3$
$Q_h = \left\{X_h^1\cup L_0^2(\Omega)\right\}$

where $X_h^r = \left\{g\in \mathcal{C}(\bar{\Omega}), \forall K\in T_h, g_{|K}\in\mathbb{P}^r\right\}$

The Taylor-Hood spaces verify the inf-sup conditions and according to [GR86] we have the convergence estimation:

$\Vert\mathbf{u}-\mathbf{u}_h\Vert_V+\Vert p-p_h\Vert_Q\leq Ch^s(\Vert\mathbf{u}\Vert_{s+1}+\Vert p\Vert_s), s=1, 2$

FreeFem++ code

In the Dirichlet condition only case, it is necessary to fix the constant part of the pressure. We add a stabilization term to the weak form:

$\begin{eqnarray} \mu\int_{\Omega}{\nabla\mathbf{u}\cdot\nabla\mathbf{v}} - \int_{\Omega}{p\mathrm{div}(\mathbf{v})} &=& \int_{\Omega}{\textbf{f}\cdot\mathbf{v}}\\ \int_{\Omega}{\mathrm{div}(\mathbf{u})q} - \varepsilon\int_{\Omega}{pq} &=& 0 \end{eqnarray}$

Where $\epsilon$ is small, typically between $10−6$ and $10−10$. The following program solves the Stokes equations on a cube with a given velocity:

``````load "msh3"

//Mesh
int nn = 8;
mesh T2d = square(nn, nn);
mesh3 Th = buildlayers(T2d, nn, zbound=[0., 1.]);

//Parameters
real mu = 1.;

//Easy ktest
func f1 = 1.;
func f2 = 1.;
func f3 = 1.;
func uEx = x;
func vEx = y;
func wEx = -2*z;
func pEx = x+y+z-3/2.;

//Taylor Hood spaces
fespace Vh(Th,P2);
Vh u1, u2, u3;
Vh v1, v2, v3;
fespace Qh(Th, P1);
Qh p, q;

//Weak formulation
solve Stokes ([u1,u2,u3,p], [v1,v2,v3,q], solver=UMFPACK)
= int3d(Th)(
mu * (
dx(u1)*dx(v1) + dy(u1)*dy(v1) + dz(u1)*dz(v1)
+ dx(u2)*dx(v2) + dy(u2)*dy(v2) + dz(u2)*dz(v2)
+ dx(u3)*dx(v3) + dy(u3)*dy(v3) + dz(u3)*dz(v3)
)
- p * q * 1e-10
- p*dx(v1) - p*dy(v2) - p*dz(v3)
- dx(u1)*q - dy(u2)*q - dz(u3)*q
)
- int3d(Th)(
f1*v1 + f2*v2 + f3*v3
)
+ on(1, 2, 3, 4, 5, 6, u1=uEx, u2=vEx, u3=wEx)
;

plot(p, fill=true);``````

Uzawa method

We have to solve the following problem:

$\begin{eqnarray} \mu\int_{\Omega}{\nabla\mathbf{u}_h\cdot\nabla\mathbf{v}_h} - \int_{\Omega}{p_h\mathrm{div}(\mathbf{v}_h)} - \int_{\Gamma_N}{\textbf{g}_2\cdot\mathbf{v}_h} &=& \int_{\Omega}{\textbf{f}\cdot\mathbf{v}_h}\\ \int_{\Omega}{\mathrm{div}(\mathbf{u}_h)q_h} &=& 0 \end{eqnarray}$

This problem is equivalent to the matricial form:

$\left( \begin{array}{cc} A & B^t\\ B & 0\\ \end{array} \right) \left( \begin{array}{c} \mathbf{u}_h\\ p_h \end{array} \right) = \left( \begin{array}{c} F\\ 0 \end{array} \right)$

Where A, B and F are the respectives representatives matrices of the bilinear forms:

$a(\mathbf{u}, \mathbf{v}) = \mu\int_{\Omega}{\nabla\mathbf{u}\cdot\nabla\mathbf{v}}$
$b(p, \mathbf{v}) = \int_{\Omega}{p\mathrm{div}(\mathbf{v})}$
$f(\mathbf{u}, \mathbf{v}) = \int_{\Omega}{\textbf{f}\cdot\mathbf{v}} - \int_{\Gamma_N}{\textbf{g}_2\cdot\mathbf{v}}$

The Uzawa [RT93] algorithm consists of minimizing the functionnal:

$J(\mathbf{u}) = \frac{1}{2}(A,\mathbf{u})-(F,\mathbf{u})$

Under the constraint $B=0$

The Uzawa algorithm reads:

$\begin{array}{ll} \textbf{Initialisation} & \\ & \mbox{Let }p_h^1\mbox{ be known}\\ & k=1.\\ \textbf{Step 1} & \\ & \mbox{Find }\mathbf{u}_h^k\mbox{ such as:}\\ & A\mathbf{u}_h^k + B^tp_h^k=F\\ \textbf{Step 2} & \\ & \mbox{Find }p_h^{k+1}\mbox{ such as:}\\ & p_h^{k+1} = B\mathbf{u}_h^k\\ \textbf{Step 3} & \\ & k=k+1.\\ & \mbox{Return to Step 1 until convergence}\\ \end{array}$

The following program solves the Stokes equations using the Uzawa method on a cube with a given velocity:

``````load "msh3"

//Mesh
int nn = 8;
mesh T2d = square(nn, nn);
mesh3 Th = buildlayers(T2d, nn, zbound=[0., 1.]);

//Parameters
real mu = 1.;

//Easy ktest
func f1 = 1.;
func f2 = 1.;
func f3 = 1.;
func uEx = x;
func vEx = y;
func wEx = -2*z;
func pEx = x+y+z-3/2.;

//Taylor Hood spaces
fespace Vh(Th,P2);
Vh u1, u2, u3;
fespace Qh(Th, P1);
Qh p;

//Weak formulation
macro grad(a) [dx(a), dy(a), dz(a)] //

varf vA(u, v)
= int3d(Th)(
)
+ on(1, 2, 3, 4, 5, 6, u=0.)
;

varf vBX(q, v)
= int3d(Th)(
q * dx(v)
)
;

varf vBY(q, v)
= int3d(Th)(
q * dy(v)
)
;

varf vBZ(q, v)
= int3d(Th)(
q * dz(v)
)
;

varf vOnX(u, v)
= on(1, 2, 3, 4, 5, 6, u=uEx)
+ int3d(Th)(
f1 * v
)
;

varf vOnY(u, v)
= on(1, 2, 3, 4, 5, 6, u=vEx)
+ int3d(Th)(
f2 * v
)
;

varf vOnZ(u, v)
= on(1, 2, 3, 4, 5, 6, u=wEx)
+ int3d(Th)(
f3 * v
)
;

matrix A = vA(Vh, Vh, solver=sparsesolver);
matrix BX = vBX(Qh, Vh);
matrix BY = vBY(Qh, Vh);
matrix BZ = vBZ(Qh, Vh);

real[int] OnX = vOnX(0, Vh);
real[int] OnY = vOnY(0, Vh);
real[int] OnZ = vOnZ(0, Vh);

func real[int] Uzawa (real[int] &pp){
real[int] bX = OnX; bX += BX * pp;
real[int] bY = OnY; bY += BY * pp;
real[int] bZ = OnZ; bZ += BZ * pp;

u1[] = A^-1 * bX;
u2[] = A^-1 * bY;
u3[] = A^-1 * bZ;

pp  = BX' * u1[];
pp += BY' * u2[];
pp += BZ' * u3[];
pp  = -pp;

return pp;
}

p=0.;
LinearCG(Uzawa, p[], eps=1.e-6);

plot(p, fill=true);``````