$$ \newcommand{\uex}{u_{\small\mbox{e}}} \newcommand{\Aex}{A_{\small\mbox{e}}} \newcommand{\half}{\frac{1}{2}} \newcommand{\halfi}{{1/2}} \newcommand{\Ddt}[1]{\frac{D #1}{dt}} \newcommand{\xpoint}{\pmb{x}} \newcommand{\normalvec}{\pmb{n}} \newcommand{\x}{\pmb{x}} \newcommand{\X}{\pmb{X}} \renewcommand{\u}{\pmb{u}} \renewcommand{\v}{\pmb{v}} \newcommand{\w}{\pmb{w}} \newcommand{\V}{\pmb{V}} \newcommand{\e}{\pmb{e}} \newcommand{\f}{\pmb{f}} \newcommand{\F}{\pmb{F}} \newcommand{\stress}{\pmb{\sigma}} \newcommand{\strain}{\pmb{\varepsilon}} \newcommand{\stressc}{{\sigma}} \newcommand{\strainc}{{\varepsilon}} \newcommand{\I}{\pmb{I}} \newcommand{\T}{\pmb{T}} % Unit vectors \newcommand{\ii}{\pmb{i}} \newcommand{\jj}{\pmb{j}} \newcommand{\kk}{\pmb{k}} \newcommand{\ir}{\pmb{i}_r} \newcommand{\ith}{\pmb{i}_{\theta}} \newcommand{\iz}{\pmb{i}_z} % Finite elements \newcommand{\basphi}{\varphi} \newcommand{\refphi}{\tilde\basphi} \newcommand{\phib}{\pmb{\varphi}} \newcommand{\sinL}[1]{\sin\left((#1+1)\pi\frac{x}{L}\right)} \newcommand{\xno}[1]{x_{#1}} %\newcommand{\xno}[1]{x^{(#1)}} \newcommand{\Xno}[1]{X_{(#1)}} \newcommand{\yno}[1]{y_{#1}} \newcommand{\Yno}[1]{Y_{(#1)}} \newcommand{\xdno}[1]{\pmb{x}_{#1}} % FEniCS commands \newcommand{\dX}{\, \mathrm{d}X} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\Integerp}{\mathbb{N}} \newcommand{\Integer}{\mathbb{Z}} $$ INF5620 Lecture: Numerical solution of the Navier-Stokes equations

INF5620 Lecture: Numerical solution of the Navier-Stokes equations

Dec 6, 2012

The physical and mathematical problem
The Navier-Stokes equations
Boundary conditions
The classical splitting method
A simple, naive approach
A working scheme
Summary
Boundary conditions
Spatial discretization by the finite element method
Increasing the implicitness
Methods based on slight compressibility
Applications

The physical and mathematical problem

Applications involving fluid flow:

Weather (flow in the atmosphere)

Ocean currents

Flight

Drag on cars

Blood circulation

Breathing

The Navier-Stokes equations

Assumptions:

Incompressible flow (velocity $ < 1/3 $ of the speed of sound)

Laminar flow

Simple fluids (constant viscosity $ \nu $)

Primary unknowns:

velocity $ \u(\x,t) $

pressure $ p(\x,t) $

Figure 1: Flow around a cylinder.

Momentum balance (Newton's 2nd law):

$$ \begin{equation} \u_t + (\u\cdot\nabla) \u = -\frac{1}{\varrho}\nabla p + \nu\nabla^2\u + \f, \label{ns:NS:mom} \end{equation} $$

Mass balance (eq. of continuity):

$$ \begin{equation} \nabla\cdot\u = 0 \label{ns:NS:mass} \thinspace . \end{equation} $$

Boundary conditions

Dirichlet conditions: components of $ \u $ are known

Neumann conditions:

Stress condition: components of the stress vector $ \stress\cdot\normalvec $ are prescribed

Outflow or symmetry condition: $ \partial\u/\partial n=0 $ (or components of this vector are zero)

Pressure known at a single point

The classical splitting method

Idea: split the N-S equations into simpler problems (operator splitting).

A simple, naive approach

The equation for $ \u $ looks like a diffusion equation...why not a Forward Euler scheme?

$$ \u_t + (\u\cdot\nabla) \u = -\frac{1}{\varrho}\nabla p + \nu\nabla^2\u + \f $$

$$ \begin{equation} \frac{\u^{n+1}-\u^n}{\Delta t} + (\u^n\cdot\nabla)\u^n = -\frac{1}{\varrho} \nabla p^n + \nu\nabla^2\u^n + \f^n, \end{equation} $$

$$ \begin{equation} \u^{n+1} = \u^n - \Delta t(\u^n\cdot\nabla)\u^n -\frac{\Delta t}{\varrho} \nabla p^n + \Delta t\,\nu\nabla^2\u^n + \Delta t\f^n \label{ns:split:FE0} \thinspace . \end{equation} $$

Two fundamental problems:

$ \nabla\cdot\u^{n+1}\neq 0 $ (that equation is not used!)

no computation of $ p^{n+1} $

A working scheme

Idea: Forward Euler in time, but evaluate $ \nabla p $ at $ t_{n+1} $ and enforce $ \nabla\cdot\u^{n+1} = 0 $.

$$ \begin{align*} \u^{n+1} &= \u^n - \Delta t(\u^n\cdot\nabla)\u^n -\frac{\Delta t}{\varrho} \nabla p^{n+1} + \Delta t\,\nu\nabla^2\u^n + \Delta t\f^n,\\ \nabla\cdot\u^{n+1} &= 0 \thinspace \end{align*} $$

Implicit system for $ \u^{n+1} $ and $ p^{n+1} $

Splitting technique:

Use old $ \beta\nabla p^n $ for $ \nabla p^{n+1} $ and advance to intermediate velocity $ \u^* $

Correct the $ \u^* $ velocity by $ \nabla\cdot\u^{n+1} = 0 $

Intermediate velocity (Forward Euler):

$$ \begin{equation} \u^{*} = \u^n - \Delta t(\u^n\cdot\nabla)\u^n - \beta\frac{\Delta t}{\varrho} \nabla p^{n} + \Delta t\,\nu\nabla^2\u^n + \Delta t\f^n \label{ns:split:FE2} \thinspace \end{equation} $$

Seek correction $ \delta\u $ such that

$$ \begin{equation} \u^{n+1} = \u^* + \delta\u, \label{ns:split:corr:du} \end{equation} $$ fulfills

$$ \nabla\cdot \u^{n+1} = 0 \thinspace .$$

Subtract $ \u^* $ equation from original $ \u^{n+1} $ equation to find $ \delta\u $:

$$ \begin{equation} \delta\u = \u^{n+1} - \u^* = -\frac{\Delta t}{\varrho}\nabla\Phi, \label{ns:split:du} \end{equation} $$ where

$$ \begin{equation} \Phi = p^{n+1} - \beta p^n \label{ns:split:Phi} \thinspace . \end{equation} $$

The oldest methods had $ \beta=0 $, but $ \beta\neq 0 $ gives in general better speed and accuracy.

$ \nabla\cdot \u^{n+1} = 0 $ implies

$$ \nabla\cdot\delta\u = - \nabla\cdot \u^* ,$$ which gives

$$ \begin{equation} \nabla^2 \Phi = \frac{\varrho}{\Delta t}\nabla\cdot\u^* \label{ns:split:Poisson} \thinspace . \end{equation} $$

When $ \Phi $ is computed,

$$ \begin{equation} \u^{n+1} = \u^* -\frac{\Delta t}{\varrho}\nabla\Phi, \label{ns:split:u:np1} \end{equation} $$ and

$$ \begin{equation} p^{n+1} = \Phi + \beta p^n \label{ns:split:p:np1} \thinspace . \end{equation} $$

Summary

Compute the intermediate velocity $ \u^* $ from \eqref{ns:split:FE2}

Solve the Poisson equation \eqref{ns:split:Poisson} for $ \Phi $

Update the velocity: $ \u^{n+1} = \u^* -\frac{\Delta t}{\varrho}\nabla\Phi $

Update the pressure: $ p^{n+1} = \Phi + \beta p^n $

Basically, we have $ u=f $ approximation problems (1, 3, 4) and a Poisson equation to solve.

Boundary conditions

Problem: $ p $ condition at one point only in the original N-S equations. Now we need boundary conditions for $ \Phi $ along the whole boundary (Poisson equation).

Use conditions for $ \u $ also for $ \u^* $

Known pressure: known $ \Phi $

Known pressure gradient: known $ \partial\Phi/\partial n $

Otherwise $ \partial\Phi/\partial n = 0 $

Spatial discretization by the finite element method

$ \u^*,\u^{n+1}\in V^{(u)} $ (modulo nonzero Dirichlet cond.)

$ p^{n+1}\in V^{(\Phi)} $ (modulo nonzero Dirichlet cond.)

Test function $ \v^{(u)}\in V^{(u)} $ for vector equations (velocity)

Test function $ v^{(\Phi)}\in V^{(\Phi)} $ for scalar equations (pressure)

Take inner product of vector equation and $ \v^{(u)} $

Integrate $ \nabla^2\u\cdot \v^{(u)} $ by parts

Integrate $ \nabla p\cdot \v^{(u)} $ by parts (optional)

Notation: $ \u $ is $ \u^{n+1} $, $ \u_1 $ is $ \u^n $, $ p $ is $ p^{n+1} $, $ p_1 $ is $ p^n $ (as in code)

$$ \begin{align} \int_\Omega \bigl( \u^*\cdot\v^{(u)} &+ \Delta t((\u_1\cdot\nabla)\nabla\u_1)\cdot\v^{(u)} - \frac{\Delta t}{\varrho}p\nabla\cdot\v^{(u)} + \nonumber\\ & \Delta t\,\nu\nabla\u_1\cdot\nabla\v^{(u)} - \Delta t f_1\bigr)\dx + \int_{\partial\Omega_{N,u}}\left( \nu\frac{\partial\u}{\partial n} - p\normalvec\right)\cdot\v^{(u)}\ds, \label{ns:split:FE2:vf} \end{align} $$ $ \forall\v^{(u)}\in V^{(u)} $.

Natural boundary condition:

$$ \nu\frac{\partial\u}{\partial n} - p\normalvec \quad (=0)$$ Usually $ {\partial\u /\partial n}=0 $ and $ p=0 $ at outlets.

Pressure Poisson equation:

$$ \begin{equation} \int_\Omega\nabla\Phi\cdot\nabla v^{(\Phi)}\dx = \frac{\varrho}{\Delta t}\int_\Omega \nabla\cdot\u^*\, v^{(\Phi)}\dx + \int_{\partial\Omega_{N,p}} \frac{\partial\Phi}{\partial n}v^{(\Phi)}\ds, \quad\forall v^{(\Phi)}\in V^{(\Phi)} \thinspace . \label{ns:split:Poisson:vf} \end{equation} $$

Velocity update:

$$ \begin{equation} \int_\Omega \u\cdot\v^{(u)}\dx = \int_{\Omega} (\u^* -\frac{\Delta t}{\varrho} \nabla\Phi)\cdot\v^{(u)}\dx,\quad\forall\v^{(u)}\in V^{(u)} \thinspace . \label{ns:split:u:np1:vf} \end{equation} $$

Pressure update:

$$ \begin{equation} \int_\Omega p v^{(\Phi)}\dx = \int_{\Omega} (\Phi + \beta p_1)v^{(\Phi)}\dx, \quad\forall v^{(\Phi)}\in V^{(\Phi)} \thinspace . \label{ns:split:p:np1:vf} \end{equation} $$

Increasing the implicitness

Stability (due to Forward Euler-style scheme):

$$ \begin{equation} \Delta t \leq \frac{h^2}{2\nu + Uh} \thinspace . \end{equation} $$ $ h $: minimum element size, $ U $: typical velocity.

Better stability by a Backward Euler scheme:

$$ \begin{align} \u^{n+1} &= \u^n - \Delta t(\u^{n+1}\cdot\nabla)\u^{n+1} -\frac{\Delta t}{\varrho} \nabla p^{n+1} + \Delta t\,\nu\nabla^2\u^{n+1} + \Delta t\f^{n+1}, \label{ns:split:BE}\\ \nabla\cdot\u^{n+1} &= 0 \thinspace . \end{align} $$

Intermediate velocity ($ \nabla p^{n+1}\rightarrow \beta p^n $):

$$ \u^{*} = \u^n - \Delta t(\u^{*}\cdot\nabla)\u^{*} - \beta\frac{\Delta t}{\varrho} p^{n+1} + \Delta t\,\nu\nabla^2\u^{*} + \Delta t\f^{n+1} \label{ns:split:BE1} \thinspace $$

Deal with nonlinearity in a simple way (1 Pickard it.):

$$ \begin{equation} (\u^{*}\cdot\nabla)\u^{*} \approx (\u^{n}\cdot\nabla)\u^{*} \thinspace . \end{equation} $$

Then we have a linear problem for $ \u^* $:

$$ \u^{*} = \u^n - \Delta t(\u^{n}\cdot\nabla)\u^{*} - \beta\frac{\Delta t}{\varrho} \nabla p^{n} + \Delta t\,\nu\nabla^2\u^{*} + \Delta t\f^{n+1} \label{ns:split:BE2} \thinspace $$

Correction (assume $ \u^{n+1}-\u^* $ small):

$$ \delta\u = \Delta t((\u^{n+1}\cdot\nabla)\u^{n+1} - (\u^{n}\cdot\nabla)\u^{*}) -\frac{\Delta t}{\varrho} \nabla \Phi + \Delta t\,\nu (\nabla^2 (\u^{n+1} - \u^*) \approx -\frac{\Delta t}{\varrho} \nabla \Phi \thinspace .$$

So, as before,

$$ \nabla^2 \Phi = \frac{\varrho}{\Delta t}\nabla\cdot\u^* $$

Need to solve an implicit equation (linear system) for $ \u^* $

The other steps are the same

A Crank-Nicolson method is more accurate (also implicit scheme for $ \u^* $)

Methods based on slight compressibility

$ \nabla\cdot\u=0 $ is problematic. Allow slight compressibility in the fluid:

$$ p_t + c^2\nabla\cdot\u =0 \thinspace . $$ $ c $: speed of sound.

Now we have evolution equations for $ \u $ and $ p $:

$$ \begin{align} \u_t &= -(\u\cdot\nabla) \u -\frac{1}{\varrho}\nabla p + \nu\nabla^2\u + \f,\\ p_t &= - c^2\nabla\cdot\u \thinspace . \end{align} $$

Forward Euler:

$$ \begin{align} \u^{n+1} &= \u^n - \Delta t (\u^n\cdot\nabla) \u^n -\frac{\Delta t}{\varrho}\nabla p^n + \Delta t\,\nu\nabla^2\u^n + \Delta t\f^n,\\ p^{n+1} &= p^n - \Delta t c^2\nabla\cdot\u^n \thinspace . \end{align} $$

Stability requires $ \Delta t\sim h^2/c $

$ c $ is large

Remedy: choose smaller $ c $ ("pseudo sound speed")

Can use Runge-Kutta methods

Can use implicit methods (but then the other scheme with $ \nabla\cdot\u=0 $ is not more complicated)

Applications

Figure 2: Flow in a channel.

Figure 3: Flow in a half a channel with a symmetry line.

Figure 4: Flow over a backward facing step.

Figure 5: Flow around a cylinder.

INF5620 Lecture: Numerical solution of the Navier-Stokes equations

Dec 6, 2012

Table of contents

The physical and mathematical problem

The Navier-Stokes equations

Boundary conditions

The classical splitting method

A simple, naive approach

A working scheme

Summary

Boundary conditions

Spatial discretization by the finite element method

Increasing the implicitness

Methods based on slight compressibility

Applications