Logo
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ on Travis-CI Feel++ on Twitter Feel++ on YouTube Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Pages
Stokes Tutorial

Table of Contents

Author
Feel++ Consortium
Date
2013-02-19




Theory

Let solve the stokes equation considering a Poiseuille flow profile.
We have the following system of equations,

\( \left\{ \begin{aligned} -\mu\Delta \bf u + \nabla p & = \bf f & \text{on}\; \Omega \;, \\ \nabla\cdot\bf u & = 0 & \text{on}\; \Omega \;, \\ \bf u & = g & \text{on}\; \Gamma \;, \\ \end{aligned} \right \)


where \(u\in [H_g^1(\Omega)]^d\) denotes the flow speed, \(p\in [L_0^2(\Omega)]\) the fluid pressure, \(\mu\) the fluid viscosity.
The last boundary condition expresses a null pressure fixed on the outlet.
The Poiseuille profile on the boundary is,

\( g(x,y)=\left( \begin{aligned} y(1-y) \\ 0 \\ \end{aligned} \right) \)


The method used to obtain the strong formulation is closed to the one used for the laplacian (see section Laplacian Examples ).
We multiply the first equation by a test function \(v\in H^1(\Omega)\) and we integrate on the domain \(\Omega\),

\( \begin{aligned} \left( \int_\Omega \mu \nabla \mathbf u : \nabla \mathbf v -\int_{\partial\Omega} \frac{\partial \mathbf u}{\partial \mathbf n} \cdot \mathbf v \right) +\int_\Omega ( \nabla\cdot(p \mathbf v) - \mathbf p \nabla\cdot v ) =\int_\Omega \mathbf f \cdot \mathbf v \;. \end{aligned} \right) \)


where \(n\) denotes a normal vector on the boundary.
The divergence theorem (or Gauss's theorem) gives,

\( \begin{aligned} \int_\Omega \nabla\cdot(p \mathbf v) = \int_{\partial\Omega} p \mathbf v\cdot \mathbf n \;. \end{aligned} \right) \)


We have to add a consistency terms to the equation to guaranty the symmetry of the bilinear form.
This term is provided by the second equation. We multiply this equation by a test function \(q\in L_2(\Omega)\) and we integrate on the domain \(\Omega\),

\( \begin{aligned} \int_{\Omega} \nabla\cdot\mathbf u q = 0 \;, \end{aligned} \right) \)


Finally, we deduce from the equations and after rearranging the integrals the variationnal formulation,

\( \begin{aligned} \int_\Omega \mu \nabla \mathbf u :\nabla \mathbf v +\int_\Omega \left( \nabla\cdot\mathbf u q - p \nabla\cdot\mathbf v \right) + \int_{\partial\Omega} \left( p \mathbf n - \frac{\partial \mathbf u}{\partial \mathbf n} \right) \cdot \mathbf v =\int_\Omega \mathbf f \cdot \mathbf v \end{aligned} \right) \)


Let us assume now that \((\mathbf v,q) \in [H_0^1(\Omega)]^d \times L_0^2(\Omega)\), the variationnal formulation leads to: Find \((\mathbf u,p)\in [H_g^1(\Omega)]^d\times L_0^2(\Omega) \) such that for all \((\mathbf v,q) \in [H_0^1(\Omega)]^d \times L_0^2(\Omega)\)

\( \begin{aligned} \int_\Omega \mu \nabla \mathbf u :\nabla \mathbf v +\int_\Omega \left( \nabla\cdot\mathbf u q - p \nabla\cdot\mathbf v \right) =\int_\Omega \mathbf f \cdot \mathbf v \end{aligned} \right) \)


Or equivalently:

\( \begin{aligned} a((\mathbf u,p),(\mathbf v,q)) = l((\mathbf v,q)) \end{aligned} \right) \)


where \(a\) is a bilinear form, continuous, coercive and where \(l\) is a linear form.

Implementation

Let's see the Feel++ code corresponding to this mathematical statement (source "doc/manual/tutorial/mystokes.cpp").
We suppose for this example the viscosity \(\mu=1\) and \(\mathbf f = 0\).

The procedure to create the mesh is very simple.
You have to provide to the command line (or via the cfg file) the gmsh.filename option.
You can provide a .geo or a .msh file (created via gmsh).

As for the laplacian problem, the code is very closed to the mathematical formulation.
We define the product of function spaces for the flow speed and the flow pressur using THch<order>()(TH stands for Taylor-Hoods) function which is Pch<N+1> \(\times\) Pch<N> for respectively flow speed and pressure spaces.
We take an element \(U=\left( \begin{array}{c} u \\ p \\ \end{array} \right) \) in this space. Then we define the integrals of the variationnal formulation for the left and the right hand side. Finally, we apply the Poiseuille profile on the boundary.
We call the solver to resolve the problem (Algebraic solutions).

int main(int argc, char**argv )
{
Environment env( _argc=argc, _argv=argv,
_about=about(_name="mystokes",
_author="Feel++ Consortium",
_email="feelpp-devel@feelpp.org"));
// create the mesh
auto mesh = loadMesh(_mesh=new Mesh<Simplex< 2 > > );
// function space
auto Vh = THch<2>( mesh );
// element U=(u,p) in Vh
auto U = Vh->element();
auto u = U.element<0>();
auto p = U.element<1>();
// left hand side
auto a = form2( _trial=Vh, _test=Vh );
a = integrate(_range=elements(mesh),
_expr=trace(gradt(u)*trans(grad(u))) );
a+= integrate(_range=elements(mesh),
_expr=-div(u)*idt(p)-divt(u)*id(p));
auto syms = symbols<2>();
auto u1 = parse( option(_name="functions.alpha").as<std::string>(), syms );
auto u2 = parse( option(_name="functions.beta").as<std::string>(), syms );
matrix u_exact = matrix(2,1);
u_exact = u1,u2;
auto p_exact = parse( option(_name="functions.gamma").as<std::string>(), syms );
auto f = -laplacian( u_exact, syms ) + grad( p_exact, syms ).transpose();
LOG(INFO) << "rhs : " << f;
// right hand side
auto l = form1( _test=Vh );
l = integrate(_range=elements(mesh),
_expr=trans(expr<2,1,5>( f, syms ))*id(u));
a+=on(_range=boundaryfaces(mesh), _rhs=l, _element=u,
_expr=expr<2,1,5>(u_exact,syms));
// solve a(u,v)=l(v)
a.solve(_rhs=l,_solution=U);
//# endmarker_main #
double mean_p = mean(_range=elements(mesh),_expr=idv(p))(0,0);
double mean_p_exact = mean(_range=elements(mesh),_expr=expr(p_exact,syms))(0,0);
double l2error_u = normL2( _range=elements(mesh), _expr=idv(u)-expr<2,1,5>( u_exact, syms ) );
double l2error_p = normL2( _range=elements(mesh), _expr=idv(p)-mean_p-(expr( p_exact, syms )-mean_p_exact) );
LOG(INFO) << "L2 error norm u: " << l2error_u;
LOG(INFO) << "L2 error norm p: " << l2error_p;
// save results
auto e = exporter( _mesh=mesh );
e->add( "u", u );
e->add( "p", p );
e->save();
}