Stokes3D.edp

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

Problema considerado

$\Omega = \{(x,y,z)\in \mathbb{R}^3 : (x,y) \in\Omega_0, 0<z<z_{\max}\}$ $\Omega_0\subset \mathbb{R}^2$ $\partial\Omega_0$ $\partial\Omega = \gamma_{in}\cup\gamma_0\cup\gamma_{out}$ $\Gamma_{in}=\gamma_{in}\times (0,z_{\max})$ $\Gamma_{out} = \gamma_{out} \times (0,z_{\max})$ .

$\mathbf{u} = \mathbf{u}(x,y,z) = (u_1(x,y,z) ,u_2(x,y,z),u_3(x,y,z))$ $p = p(x,y,z)$ es la presión del fluido. Consideramos el sistema de Stokes con las siguientes condiciones de contorno:

\begin{matrix} (1) & {\begin{cases} - ν Δ u + \nabla p = f & en Ω, \\ \nabla \cdot u = 0 & en Ω, \\ u = u_{i n} & sobre Γ_{i n}, \\ (- p I d + ν D u) \cdot n = 0 & sobre Γ_{o u t}, \\ u = 0 & sobre \partial Ω ∖ (Γ_{i n} \cup Γ_{o u t}), \end{cases} \end{matrix}

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

Formulación débil

$\bar{\mathbf{u}}\in H^1(\Omega)^3$ $\bar{\mathbf{u}}=u_{in}$ $\Gamma_{in}$ $\bar{\mathbf{u}}=0$ $\partial\Omega\setminus(\Gamma_{in}\cup\Gamma_{out})$ $\nabla \cdot \bar{\mathbf{u}}=0$ $\Omega$ . Pongamos

\begin{matrix} (2) & W = {v \in H^{1} (Ω)^{3} : v = 0 sobre \partial Ω ∖ Γ_{o u t}} . \end{matrix}

El problema que queremos resolver es:

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

$\mathbb{R}^3$ . Como antes, no es es difícil ver que este es equivalente a:

\begin{matrix} (4) & {\begin{cases} Hallar u \in \bar{u} + W y 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_{Ω} f \cdot v d x d y \\ \forall v \in W, \forall q \in L^{2} (Ω) . \end{cases} \end{matrix}

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

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

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

Este problema regularizado posee solución única. De hecho, 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{―}{Ω})^{3}, v_{h} |_{T} \in P_{2} (T) \forall T \in T_{h}, v_{h} |_{\partial Ω ∖ Γ_{o u t}} = 0} \subset W \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}

$W$ $L^2(\Omega)$ .

$\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} - \int_{Ω} p_{h} (\nabla \cdot v_{h}) + \int_{Ω} q_{h} (\nabla \cdot u_{h}) + ε \int_{Ω} p_{h} q_{h} \\ = \int_{Ω} f \cdot v_{h} \forall v \in V_{h}, \forall q_{h} \in Q_{h} . \end{cases} \end{matrix}

Resolución con FreeFEM

$\varepsilon = 10^{-8}$ $\nu = 0.1$ ,

$\mathbf{u}_{in}(x,y) = \begin{pmatrix} y(1-y)z(1-z) \\ 0 \\ 0 \end{pmatrix}$ $\Gamma_{in}$ $\mathbf{f}(x,y) = \begin{pmatrix} 0\\ 0 \\ -1 \end{pmatrix}$ .

Código


//  Stokes3D.edp
//
//  3D stationnary Stokes equations
//
//  | -\nu \Delta u + \nabla p = f, \quad \div u = 0 in \Omega
//  |  u = u_{in} on \Gamma_{in}
//  |  u = 0 on \Gamma_w
//  |  (-p Id. + \nu Du)\cdot n = 0 on \Gamma_{out}
//
//  Mixed weak formulation
//
//  | \int_\Omega (\nabla u \nabla v + \nabla p v) = 
//  | \int_\Omega f v   \forall v \in ...
//  | \int_\Omega (\div u) q = 0  \forall q \in ...
//  |  u = u_{in} on \Gamma_{in}
//  |  u = 0 on \Gamma_w
//

verbosity = 0;
load "msh3"   

//
//  Data:
//

real cnu     =.1;         // eq. coeff. \nu
real coefin  = 1.;        // coeff. for u_{in}
real epsi = 1.e-8;        // varepsilon

func uin1 = coefin*y*(1.-y)*z*(1.-z);
func uin2 = 0.;
func uin3 = 0.;

func f1 =  0.;
func f2 =  0.;
func f3 = -1.;

//
//  Discretization data
//

int Lin  = 1000;
int Lout = 1001;
int Lw   = 1002;
int nu   = 3;      // \approx number of segments discretization per unit length (except elbow)

//
//  Borders
//

border C1(s=0., 3.){x=s; y=0.; label = Lw;}
border C2(t=-pi/2, 0.){x=3.+2.*cos(t); y=2.+2.*sin(t); label = Lw;}
border C3(s=2., 5.){x=5.; y=s; label = Lw;}
border Coutlet(s=5., 4.){x=s; y=5.; label = Lout;}
border C5(s=5., 2.){x=4.; y=s; label = Lw;}
border C6(t=0., -pi/2){x=3.+cos(t); y=2.+sin(t); label = Lw;}
border C7(s=3., 0.){x=s; y=1.; label = Lw;}
border Cinlet(s=1., 0.){x=0.; y=s; label = Lin;}

plot(C1(3*nu)+C2(3*nu)+C3(3*nu)+Coutlet(nu)+C5(3*nu)
    +C6(2*nu)+C7(3*nu)+Cinlet(nu),
    wait=true);

//
//  2D Mesh
//
mesh Th2d=buildmesh(C1(3*nu)+C2(8*nu)+C3(3*nu)+Coutlet(nu)
                   +C5(3*nu)+C6(5*nu)+C7(3*nu)+Cinlet(nu));
plot(Th2d, wait=true, dim = 2, ps = "Canal.eps");

//
//  3D Mesh
//

int Ltop =  Lw, Lbase = Lw;
int[int] rup = [0,Ltop];       // upper face --> Ltop
int[int] rdown = [0,Lbase];    // lower face --> Lbase

mesh3 Th3d = buildlayers(Th2d, nu, zbound=[0.,1.],
     labelup = rup, labeldown = rdown);
plot(Th3d,wait=true);

//
//  FE spaces
//
fespace Xh(Th3d,P2);
Xh u1, u2, u3, v1, v2, v3;
fespace Mh(Th3d,P1);
Mh p, q;

Xh f1h=f1, f2h=f2, f3h=f3;

cout << " " << endl;
cout << "Dimension of Xh x Mh = " << 3*Xh.ndof + Mh.ndof << endl;
cout << " " << endl;

//
//  Macros
//
macro Gradd(u1, u2, u3) [dx(u1), dy(u1), dz(u1),
                         dx(u2), dy(u2), dz(u2),
             dx(u3), dy(u3), dz(u3)]  // EOM
macro Grad(q) [dx(q),dy(q),dz(q)]                 // EOM
macro Div(u1,u2,u3) (dx(u1)+dy(u2)+dz(u3))        // EOM

//
//  Problem
//
problem Stokes([u1,u2,u3,p],[v1,v2,v3,q]) = 
          int3d(Th3d)( cnu*(Gradd(u1,u2,u3)'*Gradd(v1,v2,v3)) 
               - Div(v1,v2,v3)*p
           + Div(u1,u2,u3)*q
           + epsi*p*q )
    - int3d(Th3d)(f1h*v1+f2h*v2+f3h*v3)
        + on(Lin, u1=uin1,u2=uin2,u3=uin3)
        + on(Ltop,Lbase,Lw, u1=0.,u2=0.,u3=0.)
;
Stokes;

//
//  Plots
//

// 3D views
plot(Th3d, u1, wait=true, value=true, fill=false, cmm="x-component of the solution: u1");
plot(Th3d, u2, wait=true, value=true, fill=false, cmm="y-component of the solution: u2");
plot(Th3d, u3, wait=true, value=true, fill=false, cmm="z-component of the solution: u3");
plot(Th3d, p,  wait=true, value=true, fill=false, cmm="Pressure: p");

//
//  Additional 2D views: Horizontal cuts
//

fespace Vh2d(Th2d,P2); 
Vh2d uxy1, uxy2, uxy3, pxy;

real hz = .2, zz;
for(int i=1; i<5; i++)
{
  zz = i*hz;
  uxy1 = u1(x,y,zz);
  uxy2 = u2(x,y,zz);
  uxy3 = u3(x,y,zz);
  pxy  = p(x,y,zz);
    plot(uxy1, cmm="x-component, cut at z ="+ zz, value=true, fill=true, wait=true);
    plot(uxy2, cmm="y-component, cut at z ="+ zz, value=true, fill=true, wait=true);
    plot(uxy3, cmm="z-component, cut at z ="+ zz, value=true, fill=true, wait=true);
    plot([uxy1,uxy2], cmm="(x,y)-components, cut at z ="+ zz, value=true, fill=true, wait=true);
    plot(pxy, cmm="Pressure, cut at z ="+ zz, value=true, fill=true, wait=true);
}

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