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

Table of Contents

Author
Christophe Prud'homme
Date
2013-05-01

Theory

We are interested in solving \( -\nu \Delta \bvec{u} + \nabla p = \bvec{f},\quad \nabla \cdot \bvec{u} = 0 \) in curl-curl formulation in 2D.

First we introduce the following notations:

\[ u \times n = u_1 n_2 - u_2 n_!, \quad u = (u_1, u_2),\ n = (n_1, n_2) \]

The curl of a vectorial field \(\phi\)

\[ \mathrm{curl} \phi = (-\frac{\partial \phi_2}{\partial x_1}, -\frac{\partial \phi_1}{\partial x_2} ) \]

The curl of a scalar field $ \(\psi\)

\[ \nabla \times \psi = \frac{\partial \psi}{\partial x_1} -\frac{\partial \psi}{\partial x_2} \]

Using the above notations we easily verify that:

\[ \mathrm{curl}( \nabla \times \psi) = \Delta \psi + \nabla (div \psi) \]

We also need a partial integration formulae:

\[ \int_\Omega \bvec{v}\cdot \mathrm{curl}w = \int_\Omega w ( \nabla \times \bvec{v}) -\int_{\partial\Omega} w (\bvec{v}\times n) \]

Using the formulae above, the system can now be written as:

\( \left\{ \begin{aligned} \mu \mathrm{curl}( \nabla \times \bvec{u}) - \nabla (div \bvec{u}) + \nabla p & = \bvec{f} & \text{on}\; \Omega \;, \ \ \nabla\cdot \bvec{u} & = 0 & \text{on}\; \Omega \; \end{aligned} \right \)


Then we recall the same method used to obtain the strong formulation for the laplacian problem (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\)

\[ \mu \int_\Omega \mathrm{curl}( \nabla \times \bvec{u}) \cdot \bf \bvec{v} - \mu \int_\Omega \nabla (div \bvec{u}) \cdot \bvec{v} + \nabla p \cdot \bvec{v} & = \int_\Omega \bvec{f}\cdot \bvec{v} \]


Then we apply the above partial integration formulae on the first term, and the green formulae on the rest of the terms, we obtain:

\[ \mu \int_\Omega ( \nabla \times \bvec{u})( \nabla \times \bvec{v}) - \mu \int_{\partial\Omega} (\nabla \times \bvec{u}) \cdot ( \bvec{v} \times \bvec{n})+ \mu \int_\Omega ( \nabla \cdot \bvec{u})( \nabla \times \bvec{v}) - \mu \int_{\partial\Omega} (\nabla \cdot \bvec{u})n \cdot \bvec{v} + \int_\Omega p \nabla \cdot \bvec{v} -\int_{\partial\Omega} p \mathbf n \cdot \bvec{v} = \int_\Omega \bvec{f} \cdot \bvec{v} \]


Implementation

int main(int argc, char**argv )
{
Environment env( _argc=argc, _argv=argv,
_desc=feel_options(),
_about=about(_name="stokes_curl",
_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=curlxt(u)*curlx(u) + divt(u)*div(u) );
a+= integrate(_range=elements(mesh),
_expr=-div(u)*idt(p)-divt(u)*id(p));
a+= integrate(_range=markedfaces(mesh,"inlet"),
_expr=( curlxt(u)*cross(id(u),N()) +
curlx(u)*cross(idt(u),N()) +
-trans(divt(u)*N())*id(u) +
30*cross(id(u),N())*cross(idt(u),N())/hFace() ) );
a+= integrate(_range=markedfaces(mesh,"outlet"),
_expr=( curlxt(u)*cross(id(u),N()) +
curlx(u)*cross(idt(u),N()) +
-trans(divt(u)*N())*id(u) +
30*cross(id(u),N())*cross(idt(u),N())/hFace() ) );
// right hand side
auto l = form1( _test=Vh );
l += integrate(_range=markedfaces(mesh,"inlet"),
_expr=2*trans(id(u))*N() );
l += integrate(_range=markedfaces(mesh,"outlet"),
_expr=1*trans(id(u))*N() );
// Dirichlet condition
a+=on(_range=markedfaces(mesh,"wall"), _rhs=l, _element=u,
_expr=vec(cst(0.),cst(0.)));
// solve a(u,v)=l(v)
a.solve(_rhs=l,_solution=U);
// save results
auto e = exporter( _mesh=mesh );
e->add( "u", u );
e->add( "p", p );
e->save();
}