Convection dominated flows

Cette partie est en anglais et nécessite d’être traduite.

Let $\Omega \subset \mathbb{R}^d, d=1,2,3$ , consider this equation defined in $\Omega$ , find $u$ a scalar field (e.g. temperature or concentration) such that

$\underbrace{\frac{\partial u}{\partial t}}_{time} \underbrace{- \nabla \cdot ( \epsilon \nabla u )}_{diffusion} + \underbrace{\mathbf{ \beta} \cdot \nabla u}_{advection/convection} + \underbrace{\mu u}_{reaction} = \underbrace{f}_{source},\quad \nabla \cdot \mathbf{ \beta} = 0$

with the given data

$\epsilon$ is the diffusion coefficient (matrix $d \times d$ ), $\epsilon(\mathbf{x}) > 0$
$\mathbf{\beta}$ velocity and $\nabla \cdot \mathbf{ \beta} \in L^2(\Omega)$
$\mu > 0$ reaction (or absorption) coefficient.

From now on we consider only steady problems ( $\frac{\partial u}{\partial t}=0$ )

1. Applications

Transport of temperature
Transport of pollutants
Transport of chemical species ( $O_2, CO_2,...$ )

These equations are often coupled with fluid flow equations such as incompressible Navier-Stokes equations (air, blood, water…) or darcy equations (porous media) to obtain the velocity field $\mathbf{ \beta}$ . In this lesson we consider $\mathbf{ \beta}$ given.

2. Boundary conditions

Boundary conditions can be of 2 types on parts of the boundary or the whole boundary of the domain $\Omega$ :

Dirichlet conditions : $\label{eq:2} u|_{{\partial \Omega}_D} = g$
Neumann conditions :
- Outflow $\label{eq:4} ( - \epsilon \nabla u + \mathbf{ \beta} u ) \cdot \mathbf{ n}|_{{\partial \Omega}_N} = (\mathbf{ \beta} u) \cdot \mathbf{ n}$
- (Heat) Flux $\label{eq:5} ( - \epsilon \nabla u ) \cdot \mathbf{ n}|_{{\partial \Omega}_N} = Q \quad (e.g. \mathbf{\beta} = \mathbf{0})$

3. Variational Formulation

The variational problem reads, Find $u \in V$ such that for all $v\ in V$

$a(u,v) \equiv \int_\Omega \epsilon \nabla u \cdot \nabla v + (\mathbf{ \beta} \cdot \nabla u ) v + \underbrace{\int_{\partial \Omega} (- \epsilon \nabla u ) \cdot \mathbf{ n}\ v}_{\mbox{apply boundary conditions}} = \int_\Omega f v \equiv \ell (v)$

where

$V$ is an Hilbert space
we used the identity below for the integration by part

$\nabla \cdot (\mathbf{ c} w ) = \mathbf{ c} \cdot \nabla w + \underbrace{w \nabla \cdot \mathbf{ c}}_{=0}$

Suppose that $-\frac{1}{2} \nabla \cdot ( \mathbf{ \beta} ) + \mu \geq 0$ almost everywhere in $\Omega$ then we can show that $a$ is coercive provide $\epsilon \geq \epsilon_0 > 0$

$a(v,v) \geq \alpha ||v||_{1,\Omega}, \quad \alpha = \frac{\epsilon_0}{1+C_\Omega^2}$

$a$ is continuous, there exists $M$ a constant such that

$|a(u,v)| \leq < M ||u||_{H^1{(\Omega)}} ||v||_{H^1{(\Omega)}}, \quad M = ||\mu||_{0,\Omega} + ||\epsilon||_{\infty,\Omega} + ||\mathbf{\beta}||_{\infty,\Omega}$

4. Discrete formulation

Let $V_h \subset V \equiv H^1(\Omega)$ a discrete space, consider the standard Galerkin approximation on $V_h$ for . The problem reads

Problem

Find $u_h \in V_h$ such that for all $v_h \in V_h$ we have:

$\int_\Omega \epsilon \nabla u_h \cdot \nabla v_h + (\mathbf{ \beta} \cdot \nabla u_h ) v_h + \underbrace{\int_{\partial \Omega} (- \epsilon \nabla u_h ) \cdot \mathbf{ n}\ v_h}_{\mbox{apply boundary conditions}} = \int_\Omega f v_h$

We can show that

$||u_h||_{1,\Omega} \leq \frac{1}{\alpha} ||f||_{0,\Omega}, \quad ||\nabla u_h||_{1,\Omega} \leq \frac{C_\Omega}{\epsilon_0} ||f||_{0,\Omega},$

which means that the Galerkin error inequality (best approximation) gives

$||u-u_h||_V \leq \frac{M}{\alpha} \mathrm{inf}_{v_h \in V_h} ||u-v_h||_V$

which means that given $M$ and $\alpha$ , the estimate $\epsilon_0 \rightarrow 0$ In such case, the standard Galerkin method can yield to inaccurate solutions unless:

we use a very small $h$ (cost!!)
we use a stabilisation method (loss of $\frac{1}{2}$ in convergence rate)

5. Stabilisation methods for convection dominated flows

We now introduce

$\mathrm{Pe}=\frac{|\mathbf{ \beta}|L}{2 \epsilon}$ the global number. $L$ is the characteristic size of the domain.
the local Péclet number $\mathrm{Pe}_h=\frac{|\mathbf{\beta}|h}{2 \epsilon}$

The dominating convection and inacurrate behavior occurs when, on at least some cells, $\mathrm{Pe} > 1$ or $\mathrm{Pe}_h > 1$ .

A remedy is to use a sufficiently small

$h$ same to ensure

$\mathrm{Pe}_h < 1$ . For example if

$|\mathbf{beta}| = 1$ and

$\epsilon=5e-4$ , we should take

$h=1e-4$ .

Another remedy is to use a different approximation space for the unknown

$u_h$ and the test functions

$v_h$ . We talk about Petrov-Galerkin approximation.

Problem

Find $u_h \in V_h$ such that

$a_h(u_h,v_h) = \ell_h(v_h),\quad \forall v_h \in V_h$

with

$a_h(u_h,v_h) = a(u_h,v_h) + b_h(u_h,v_h),\quad \ell_h(v_h) = \ell(v_h) + g_h(v_h),\quad$

The purpose of

$b_h$ and

$g_h$ is to eliminate(or reduce) the numerical oscillations produced by the standard Galerkin method: they are called stabilisation methods and depend on

$h$ . The term stabilisation method is not exact. Indeed the Galerkin method is already stable (i.e. continuity). Here it is to be understood as the aim of reducing (or elimination) numerical oscillations when

$\mathrm{Pe} > 1$ .

Without doing anything wiggles occur.

There are remedies so called stabilisation techniques, here some some examples:

Artificial diffusion (streamline diffusion) (SDFEM)
Galerkin Least Squares method (GaLS)
Streamline Upwind Petrov Galerkin (SUPG)
Continuous Interior Penalty methods (CIP)

5.1. Artificial diffusion (or streamline diffusion) (SDFEM)

Method The method consists in adding an

$\epsilon_h =\epsilon(1+\phi(\mathrm{Pe}))$

with $\phi(\mathrm{Pe}) \rightarrow 0$ as $h \rightarrow 0$ , e.g. $\phi(\mathrm{Pe}) = \mathrm{Pe}-1+B(2*\mathrm{Pe})$ where $B$ is the so-called Bernoulli function $B(t) = \frac{t}{e^t-1}$ if $t > 0$ and $B(0) = 1$ (also exponential fitting scheme)

$b_h(u_h,v_h) = \int_\Omega \epsilon \Phi(\mathrm{Pe}) \nabla u_h \cdot \nabla v_h, \quad g_h(v_h) = 0$

Theorem

for a given $\epsilon$ and for $h$ tending to $0$ , we have for $u \in H^{r+1}(\Omega)$

$||u-u_h||_{1,\Omega} \leq C_1 \Big[ h^r||u||_{r+1,\Omega} + \phi(\mathrm{Pe})||u||_{1,\Omega}\Big]$

and for a given $h$ and $\epsilon$ tending to 0,

$||u-u_h||_{1,\Omega} \leq C_1 \Big[ h^{r-1}||u||_{r+1,\Omega} + ||u||_{1,\Omega}\Big]$

If $\phi(\mathrm{Pe})=\frac{|\mathbf{ \beta}|h}{2 \epsilon}$ , the convergence is linear, with the exponential fitting scheme it is quadratic if $r \geq 2$ .

5.2. GaLS and SUPG

First we decompose our operators into a symmetric ( $<Lu,v> = <u,Lv>$ and skew symmetric ( $<L u, v> = -<u,L v>$ ) contributions, we start with

$L u = -\epsilon \Delta u + \nabla \cdot (\mathbf{ \beta} u ) + \mu u$

$L u = \underbrace{-\epsilon \Delta u + \Big[ \mu + \frac{1}{2} \nabla \cdot \mathbf{ \beta} \Big] u}_{L_S u} + \underbrace{\frac{1}{2}\Big[ \nabla \cdot ( \mathbf{ \beta} u) + \mathbf{ \beta} \cdot \nabla u \Big]}_{L_{SS} u}$

Consistent schemes

We say that a method is consistent when adding a term to a problem such as: Find $u_h \in V_h$ such that

$a(u_h,v_h) + \mathcal{L}_h(u_h,f;v_h) = (f,v_h), \quad \forall v_h \in V_h$

the term added statisfies

$\mathcal{L}_h(u,f;v_h) = 0, \forall v_h \in V_h$

5.2.1. Choice for consistent methods

A possible choice for $\mathcal{L}_h$ is the following

$\mathcal{L}_h(u_h,f;v_h) = \mathcal{L}^{(\rho)}_h(u_h,f;v_h) = \sum_{K \in \mathcal{T}_h} \delta (L u_h - f, \mathcal{S}^{(\rho)}_K(v_h))_{0,\Omega}$

where

$(\cdot,\cdot)_{0,\Omega}$ is the $L^2$ scalar product
$\rho$ and $\delta$ are parameters

and we have set

$\mathcal{S}^{(\rho)}_K(v_h) = \frac{h_K}{|\mathbf{\beta}|}\Big[ L_{SS} v_h + \rho L_S v_h\Big]$

Galerkin Least-Square: if $\rho = 1$ we have the Galerkin Least Square method (GaLS)

$\mathcal{S}^{(\rho)}_K(v_h) = \frac{h_K}{|\mathbf{ \beta}|}\Big[ L v_h\Big]$

Streamline Upwind Petrov-Galerkin: if $\rho = 0$ we have the Streamline Upwind Petrov-Galerkin (SUPG)

$\mathcal{S}^{(0)}_K(v_h) = \frac{h_K}{|\mathbf{ \beta}|}\Big[ L_{SS} v_h\Big]$

Douglas and Wang: if $\rho = -1$ we have the Douglas and Wang (DW)

$\mathcal{S}^{(-1)}_K(v_h) = \frac{h_K}{|\mathbf{ \beta}|}\Big[ (L_{SS} -L_S )v_h\Big]$

We define the $\rho$ Norm

$||v||_{(\rho)} = \Big\{\epsilon ||\nabla u||^2_{0,\Omega} + ||\sqrt{\gamma} v||^2_{0,\Omega} + \sum_{K \in \mathcal{T_h}} \delta \Big( (L_{SS}+\rho L_S )v, \mathcal{S}^{(\rho)}_K(v) \Big)_{0,\Omega} \Big\}^{1/2}$

where $\gamma$ is a positive constant such that $-\frac{1}{2} \nabla \cdot \mathbf{\beta} + \mu \geq \gamma > 0$

We have the following result

Theorem

if $u \in H^{k+1}(\Omega)$ , then the following error estimates hold:

${\|u-u_h\|_{(\rho)}} \leq C {h^{k+1/2}} |u|_{k+1,\Omega}$

GaLS

In practice for GaLS ( $\rho = 1$ ) we take $\delta$ such that

$\delta(h_K,\epsilon) \frac{h_K}{|\mathbf{ \beta}|} = \Big( \frac{1}{h_K} + \frac{\epsilon}{h^2_K} \Big)^{-1}$

and we can prove the following estimates if $u\in H^{k+1}(\Omega)$ ,

$\forall \epsilon \quad {\|u-u_h\|_{0,\Omega}} \leq c {h^{k+1/2}} \|u\|_{k+1,\Omega}$

$\forall \epsilon \geq c h \quad {\|u-u_h\|_{1,\Omega}} \leq c {h^{k}} \|u\|_{k+1,\Omega}$

and finally if the family $\{\mathcal{T}_h\}_{h > 0}$ is quasi-uniform and $\epsilon \leq c h$ , then

$\| \beta \cdot \nabla (u -u_h) \|_{0,\Omega} \leq c h^k \|u \|_{k+1,\Omega}$

5.3. Continuous Interior Penalty

In the continuous interior penalty we add the following term

$\sum_{F \in \Gamma_\mathrm{int} } \int_{F} \eta\ h_F^2\ |\mathbf{ \beta} \cdot \mathbf{n}|\ \jump{\nabla u} \jump{\nabla v}$

where

$\Gamma_\mathrm{int}$ is the set of internal faces
the $\mathrm{Pe}>>1$ (typically it is applied to all internal faces)
we have

\jump{\nabla u} = \nabla u \cdot \mathbf{n}|_1 + \nabla u \cdot \mathbf{n}|_2

is the so called jump of $\nabla u$ (scalar valued) across the face.

In the case of scalar valued functions

$\jump{u} = u \mathbf{n}|_1 + u \mathbf{n}|_2$

Choice for $\eta$ $\eta$ can be taken in the range $[1e-2;1e-1$ ]. A typical value is $\eta=2.5e-2$ . A similar error estimate $O(h^{r+1/2})$ holds for CIP.

Example CIP

// define the stabilisation coefficient expression
auto stab_coeff = (eta*abs(idv(beta))*abs(trans(N())*idv(beta)))*vf::pow(hFace(),2.0));

// assemble the stabilisation operator
form2( Xh, Xh, M ) +=
 integrate( internalfaces(Xh->mesh()), // faces of the mesh
            stab_coeff*(trans(jumpt(gradt(u)))*jump(grad(v))));

cpp