Mixed finite element method for stokes flow

Mixed finite element method for Stokes flow

Flow characteristics:

• Creeping flow
• Convective acceleration is negligible
• Incompressible

Strong form:

\[ \rho\frac{\partial\mathbf{v}}{\partial t}+\nabla p-\mu\nabla^{2}\mathbf{v}-\rho\mathbf{b}=\mathbf{0}, \]

\[ \nabla \cdot \mathbf{v} = 0 \]

Essential boundary condition:

\[ v_{i} = v_{ig}, \text{ on } \Gamma_{i}^{g} \]

Neumann boundary condition:

\[ -pn_{i}+\mu \frac{\partial v_i}{\partial x_j} n_{j}=h_i, \text{ on } \Gamma_{i}^{h} \]

here, note that \(h_{i}\) is not the real traction, it is called the pseudo traction.

Weak form:

Find \(\mathbf{v}\) and \(p\) such that for all \(\delta \mathbf{v}\) and \(\delta p\) following is true:

\[ B_{G}\left(\delta\mathbf{v},\delta p;\mathbf{v},p\right)=L_{G}(\delta\mathbf{v},\delta p) \]

where,

\[ B_{G}\left(\delta\boldsymbol{v},\delta p;\boldsymbol{v},p\right)=\int_{\Omega}\left(-\frac{\partial\delta v_{i}}{\partial x_{i}}p+\frac{\partial\delta v_{i}}{\partial x_{j}}\mu\frac{\partial v_{i}}{\partial x_{j}}+\delta p\frac{\partial v_{j}}{\partial x_{j}}\right)d\Omega \]

\[ L_{G}(\delta\mathbf{v},\delta p)=\int_{\Omega}\delta v_{i}\rho b_{i}d\Omega+\int_{\Gamma_{i}^{h}}\delta v_{i}h_{i}d\Omega \]

Proof of weak form:

The process of getting the above mentioned weak form is given below.

First, multiple the residual of momentum equation with the test function \(\delta v\) and the continuity equation with the \(\delta p\), and integrate over the entire domain.

\[ \int_{\Omega}\delta v_{i}\left\{ \frac{\partial p}{\partial x_{i}}-\mu\frac{\partial^{2}v_{i}}{\partial x_{j}\partial x_{j}}-\rho b_{i}\right\} d\Omega+\int_{\Omega}\delta p\frac{\partial v_{j}}{\partial x_{j}}d\Omega=0 \]

Now we use integration by parts and Green’s divergence theorem to pressure and shear stress gradient terms. Then we obtain the following

\[ \begin{split}-\int_{\Omega}\frac{\partial\delta v_{i}}{\partial x_{i}}pd\Omega+\int_{\Omega}\frac{\partial\delta v_{i}}{\partial x_{j}}\mu\frac{\partial v_{i}}{\partial x_{j}}d\Omega\\ -\int_{\Gamma}\delta v_{i}\left(-pn_{i}+\mu\frac{\partial v_{i}}{\partial x_{j}}\right)n_{j}d\Omega\\ -\int_{\Omega}\delta v_{i}\rho b_{i}d\Omega+\int_{\Omega}\delta p\frac{\partial v_{j}}{\partial x_{j}}d\Omega\\ =0 \end{split} \]

now noting that

\[ \delta v_{i}=0 \text{ on } \Gamma_{i}^{g} \]

and

\[ -pn_{i}+\mu \frac{\partial v_i}{\partial x_j} n_{j}=h_i, \text{ on } \Gamma_{i}^{h} \]

then,

\[ \begin{split}-\int_{\Omega}\frac{\partial\delta v_{i}}{\partial x_{i}}pd\Omega+\int_{\Omega}\frac{\partial\delta v_{i}}{\partial x_{j}}\mu\frac{\partial v_{i}}{\partial x_{j}}d\Omega\\ +\int_{\Omega}\delta p\frac{\partial v_{j}}{\partial x_{j}}d\Omega\\ -\int_{\Omega}\delta v_{i}\rho b_{i}d\Omega-\int_{\Gamma_{i}^{h}}\delta v_{i}h_{i}d\Omega\\ =0 \end{split} \]

or,

\[ \begin{aligned}\int_{\Omega}\left(-\frac{\partial\delta v_{i}}{\partial x_{i}}p+\frac{\partial\delta v_{i}}{\partial x_{j}}\mu\frac{\partial v_{i}}{\partial x_{j}}+\delta p\frac{\partial v_{j}}{\partial x_{j}}\right)d\Omega\\ =\int_{\Omega}\delta v_{i}\rho b_{i}d\Omega+\int_{\Gamma_{i}^{h}}\delta v_{i}h_{i}d\Omega \end{aligned} \]

This completes the proof.

Finite element approximation

For test and trial function of pressure field use the following finite element approximation:

\[ \delta p=\sum_{I=1}^{n_{n}}\delta p_{I}N_{p}^{I} \]

\[ p=\sum_{I=1}^{n_{n}} p_{I}N_{p}^{I} \]

where, \(\delta p_I\) and \(p_I\) are the nodal values of \(\delta p\) and \(p\), and \(N^I_p\) is the shape function.

For test and trial function of velocity field we use the following finite element approximation:

\[ \delta v_{iI}=\sum_{I=1}^{n_{n}}\delta v_{iI}N_{v}^{I} \]

\[ v_{iI}=\sum_{I=1}^{n_{n}} v_{iI}N_{v}^{I} \]

where, \(\delta v_{iI}\) and \(v_{iI}\) are the nodal value of \(i\)th component of velocity field, and \(N_{v}^I\) is the shape function of velocity field.

Finite element discretization

\[ \begin{aligned}\int_{\Omega}-\frac{\partial\delta v_{i}}{\partial x_{i}}pd\Omega & =-\delta v_{iI}\left[\int_{\Omega}\frac{\partial N_{v}^{I}}{\partial x_{i}}N_{p}^{J}d\Omega\right]p_{J}\\ & =\delta v_{iI}\left[G_{i1}^{IJ}\right]p_{J} \end{aligned} \]

\[ \begin{aligned}\int_{\Omega}\frac{\partial\delta v_{i}}{\partial x_{j}}\mu\frac{\partial v_{i}}{\partial x_{j}}d\Omega & =\delta v_{iI}\int_{\Omega}\frac{\partial N_{v}^{I}}{\partial x_{j}}\mu\frac{\partial v_{i}}{\partial x_{j}}d\Omega\\ & =\delta v_{iI}\left[\int_{\Omega}\frac{\partial N_{v}^{I}}{\partial x_{j}}\mu\frac{\partial N_{v}^{J}}{\partial x_{j}}d\Omega\right]v_{iJ}\\ & =\delta v_{iI}\left[K_{ij}^{IJ}\right]v_{jJ} \end{aligned} \]

\[ \begin{aligned}\int_{\Omega}\delta p\frac{\partial v_{i}}{\partial x_{i}}d\Omega & =\delta p_{I}\left[\int_{\Omega}N_{p}^{I}\frac{\partial N_{v}^{J}}{\partial x_{i}}d\Omega\right]v_{iJ}\\ & =\delta p_{I}\left[-G_{1i}^{JI}\right]v_{iJ} \end{aligned} \]

\[ \begin{aligned}\int_{\Omega}\delta v_{i}\rho b_{i}d\Omega+\int_{\Gamma_{i}^{h}}\delta v_{i}h_{i}d\Omega & =\delta v_{iI}\left\{ \int_{\Omega}N_{v}^{I}\rho b_{i}d\Omega+\int_{\Gamma_{i}^{h}}N_{v}^{I}h_{i}ds\right\} \\ & =\delta v_{iI}\left\{ F_{iI}\right\} \end{aligned} \]

where,

\[ G_{i1}^{IJ}=-\int_{\Omega}\frac{\partial N_{v}^{I}}{\partial x_{i}}N_{p}^{J}d\Omega \]

\[ K_{ij}^{IJ}=\delta_{ij}\int_{\Omega}\frac{\partial N_{v}^{I}}{\partial x_{k}}\mu\frac{\partial N_{v}^{J}}{\partial x_{k}}d\Omega \]

\[ F_{iI}=\int_{\Omega}N_{v}^{I}\rho b_{i}d\Omega+\int_{\Gamma_{i}^{h}}N_{v}^{I}h_{i}ds \]

System of linear equation

The system of linear equation is given by following. \[ \left[\begin{array}{cc} \mathbf{K} & \mathbf{G}\\ \mathbf{G}^{T} & \mathbf{0} \end{array}\right]\left\{ \begin{array}{c} \mathbf{V}\\ \mathbf{P} \end{array}\right\} =\left\{ \begin{array}{c} \mathbf{F}\\ \mathbf{0} \end{array}\right\} \]