Stokes2D.edp

Resumen: En este ejemplo se considera el problema estacionario de Stokes bidimensional.

Problema considerado

$\Omega \subset \mathbb{R}^2$ $\mathbf{u} = \mathbf{u}(x,y)$ $p = p(x,y)$ $f = f(x,y)$ $\mathbf g(x,y)$ $(x,y)\in \partial\Omega$ .

$\mathbf{u} = \mathbf{u}(x,y) = (u_1(x,y) ,u_2(x,y))$ $\mathbf{u}$ $p$ deben verificar el sistema de Stokes:

\begin{matrix} (1) & {\begin{cases} - ν Δ u + \nabla p = f (x, y) & en Ω, \\ \nabla \cdot u = 0 & en Ω, \\ u = g (x, y) & sobre \partial Ω, \end{cases} \end{matrix}

$\nu$ es una constante positiva dada (la viscosidad cinemática).

Formulación débil

$L^2_0(\Omega):=\{q\in L^2(\Omega): \displaystyle\int_\Omega q = 0\}$ $\eqref{eq.Stokes2}$ $\mathbf v\in H^1_0(\Omega)^2$ $q\in L^2_0(\Omega)$ $\Omega$ $\eqref{eq.Stokes2}$ :

\begin{matrix} (2) & {\begin{cases} Hallar u \in \bar{u} + H_{0}^{1} (Ω)^{2} y p \in L_{0}^{2} (Ω) tales que \\ ν \int_{Ω} \nabla u : \nabla v d x d y - \int_{Ω} p (\nabla \cdot v) d x d y = \int_{Ω} f \cdot v d x d y \forall v \in H_{0}^{1} (Ω)^{2} \\ \int_{Ω} q (\nabla \cdot u) d x d y = 0 \forall q \in L_{0}^{2} (Ω), \end{cases} \end{matrix}

$\bar{\mathbf{u}}$ $\bar{\mathbf{u}}\in H^1(\Omega)^2$ $\nabla \cdot \bar{\mathbf{u}}=0$ $\Omega$ .

$\eqref{eq.StokesFV}$ que se conoce como formulación débil mixta. No es difícil ver que es equivalente a:

\begin{matrix} (3) & {\begin{cases} Hallar u \in \bar{u} + H_{0}^{1} (Ω)^{2} y p \in L_{0}^{2} (Ω) tales que \\ ν \int_{Ω} \nabla u : \nabla v d x d y - \int_{Ω} p (\nabla \cdot v) d x d y + \int_{Ω} q (\nabla \cdot u) d x d y = \int_{Ω} f \cdot v d x d y \\ \forall v \in H_{0}^{1} (Ω)^{2}, \forall q \in L_{0}^{2} (Ω) . \end{cases} \end{matrix}

$\eqref{eq.StokesFVe}$ $H_0^1(\Omega)^2\times L_0^2(\Omega)$ . Por tanto, dicha formulación variacional no está en las condiciones del Teorema de Lax-Milgram.

$\eqref{eq.StokesFVe}$ $\varepsilon>0$ :

\begin{matrix} (4) & {\begin{cases} Hallar u \in \bar{u} + H_{0}^{1} (Ω)^{2}, p \in L^{2} (Ω) tales que \\ ν \int_{Ω} \nabla u : \nabla v d x d y - \int_{Ω} p (\nabla \cdot v) d x d y + \int_{Ω} q (\nabla \cdot u) d x d y + ε \int_{Ω} p q d x d y = \int_{Ω} f \cdot v d x d y \\ \forall v \in H_{0}^{1} (Ω)^{2}, \forall q \in L^{2} (Ω) . \end{cases} \end{matrix}

Este problema regularizado posee solución única, puesto que ahora sí se verifican las hipótesis del Teorema de Lax-Milgram.

Problema aproximado

$P_2$ $P_1$ -Lagrange para la presión. Consideramos los espacios

\begin{matrix} V_{h} = {v_{h} | v_{h} \in C^{0} (\overset{―}{Ω})^{2}, v_{h} |_{T} \in P_{2} (T) \forall T \in T_{h}, v_{h} |_{\partial Ω} = 0} \subset H_{0}^{1} (Ω)^{2} \end{matrix}

\begin{matrix} Q_{h} = {q_{h} \in C^{0} (Ω) | q_{h} |_{T} \in P_{1} (T) \forall T \in T_{h}} \subset L^{2} (Ω) . \end{matrix}

$\eqref{eq.StokesFVep}$ es el siguiente:

\begin{matrix} (PV_{h}) & {\begin{cases} Hallar u_{h} \in \bar{u} + V_{h} y p_{h} \in Q_{h} tales que \\ ν \int_{Ω} \nabla u_{h} : \nabla v_{h} d x d y - \int_{Ω} p_{h} (\nabla \cdot v_{h}) d x d y + \int_{Ω} q_{h} (\nabla \cdot u_{h}) d x d y + ε \int_{Ω} p_{h} q_{h} d x d y \\ = \int_{Ω} f \cdot v_{h} d x d y \forall v \in V_{h}, \forall q_{h} \in Q_{h} . \end{cases} \end{matrix}

Resolución con FreeFEM

$\Omega = (0,1)\times(0,1)$ $\varepsilon = 10^{-8}$ $\nu = 1$ ,

$\mathbf{u}(x,y) = \begin{pmatrix} 10^5 \\ 0 \end{pmatrix}$ $\Gamma_3$ $\mathbf{u}(x,y) = \begin{pmatrix} 0 \\ 0 \end{pmatrix}$ $\Gamma_1\cup\Gamma_2\cup\Gamma_4$ $\mathbf{f}(x,y) = \begin{pmatrix} 1\\ 0 \end{pmatrix}$ .

Código


x
//
//  Data:
//
real xnu = 1., f1 = 1.,f2 = 0.;
real Uh = 1.e5;
real epsr = 1.e-8;

//
//  Mesh:
//
mesh Th = square(10,10,flags=3); 
plot(Th,wait = 1);

//
//  Finite element space
//  [u1, u2, p]  : unknown
//  [v1, v2, p]  : function test
//
fespace Vh(Th,P2); 
fespace Ph(Th,P1); 

Vh u1,u2,v1,v2;
Ph p,q;

//
//  Macro definitions:
//  Macros can have arguments (like those ones) and MUST finish by  //
macro Grad(u,v) [dx(u), dy(v), dx(v), dy(u)] // 
macro div(u1,u2) (dx(u1)+dy(u2)) //

//
//  Problem + Solution:
//
 solve Stokes( [u1,u2,p],[v1,v2,q] ) = int2d(Th) ( xnu*Grad(u1,u2)'*Grad(v1,v2) 
          - p*div(v1,v2)
          + div(u1,u2)*q + epsr*p*q )
          + on(1,2,4, u1 = 0.,u2 = 0.) + on(3, u1 = Uh,u2 = 0.); 
//
//  Plots
//
plot([u1,u2],wait = true,value = true);
plot(p,wait = true,value = true ,fill = true ,nbiso = 40);

Anna Doubova - Rosa Echevarría - Dpto. EDAN - Universidad de Sevilla