Franz J. Vesely > CompPhys Tutorial > Selected Applications > Hydrodynamics  

< >
Part III: Ch. 8

8.2 Incompressible Flow with Viscosity

In (8.3) we put $d \rho / dt \equiv \partial \rho + \vec{v} \cdot \nabla \rho=0$ to find
$ \nabla \cdot \vec{v}=0 $     (8.52)
$\longrightarrow$ Flow is source-free!

A consequence of this is that

$ \nabla \cdot (\nabla \vec{v})+\nabla \cdot (\nabla \vec{v})^{T}= \nabla^{2} \vec{v} $     (8.53)

by which the Navier-Stokes equation reduces to the form

$ \frac{\textstyle \partial \vec{v}}{\textstyle \partial t} + (\vec{v} \cdot \nabla )\vec{v}= - \nabla \bar{p} + \nu \nabla^{2}\vec{v} $     (8.54)

with $\nu \equiv \mu/\rho$ and $\bar{p}\equiv p/\rho$.

Two techniques to solve (8.52), (8.54):

8.2.1 Vorticity Method

Take the rotation of

$ \frac{\textstyle \partial \vec{v}}{\textstyle \partial t} + (\vec{v} \cdot \nabla )\vec{v}= - \nabla \bar{p} + \nu \nabla^{2}\vec{v} $

to find

$ \frac{\textstyle \partial \vec{w}}{\textstyle \partial t} + (\vec{v} \cdot \nabla )\vec{w}= \nu \nabla^{2}\vec{w} $     (8.55)

with the vorticity $\vec{w}\equiv \nabla \times \vec{v}$.

$\longrightarrow$ Vorticity is transported both by an advective process ($\vec{v} \cdot \nabla \, \vec{w}$) and by viscous diffusion.

Write the (divergence-free) velocity as the rotation of a streaming function $\vec{u}$:

$ \vec{v}\equiv \nabla \times \vec{u} $     (8.56)

This does not determine $\vec{u}$ uniquely; $\longrightarrow$ require that $\nabla \cdot \vec{u}=0$.

The basic equations for the vorticity method are therefore:

$ \begin{eqnarray} \frac{\textstyle \partial \vec{w}}{\textstyle \partial t} + (\vec{v} \cdot \nabla )\vec{w}&=& \nu \nabla^{2}\vec{w} \\ \nabla^{2}\vec{u}&=&-\vec{w} \\ \vec{v}&=&\nabla \times \vec{u} \end{eqnarray} $          (8.57-8.59)

Two-dimensional case: $\vec{u}$ and $\vec{w}$ have only $z$-components:

$ \begin{eqnarray} \frac{\textstyle \partial w}{\textstyle \partial t} &=& \nu \nabla^{2}w -\left(v_{x}\partial_{y}-v_{y}\partial_{x}\right)w \\ \nabla^{2}u&=&-w \\ \vec{v}&=&u \nabla \times \vec{e}_{z} = \left( \begin{array}{r}\partial_{y}u \\ -\partial_{x}u \end{array} \right) \end{eqnarray} $

$\longrightarrow$ Apply Lax-Wendroff to these equations:

Vorticity method (2-dimensional):

Let the flow field at time $t_{n}$ be given by $u_{i,j}^{n}$ and $w_{i,j}^{n}$. For simplicity, let $\Delta y=\Delta x \equiv \Delta l$.
  1. Auxiliary quantities:

    $ \begin{eqnarray} v_{x,i,j+1}^{n}&=&\frac{\textstyle 1}{\textstyle 2 \Delta l}\left( u_{i,j+2}^{n} -u_{i,j}^{n} \right) \\ v_{y,i,j+1}^{n}&=&-\frac{\textstyle 1}{\textstyle 2 \Delta l}\left( u_{i+1,j+1}^{n} -u_{i-1,j+1}^{n} \right) \\ && \\ w_{i,j+1}^{n+1/2}&=& \frac{\textstyle 1}{\textstyle 4} \left[ w_{i,j}^{n}+w_{i+1,j+1}^{n}+w_{i,j+2}^{n}+w_{i-1,j+1}^{n} \right] \\ && -\frac{\textstyle \Delta t}{\textstyle 4\Delta l}v_{x,i,j+1}^{n} \left( w_{i+1,j+1}^{n}-w_{i-1,j+1}^{n}\right) \\ && -\frac{\textstyle \Delta t}{\textstyle 4\Delta l}v_{y,i,j+1}^{n} \left( w_{i,j+2}^{n}-w_{i,j}^{n}\right) \end{eqnarray} $

    etc., for the $4$ lattice points nearest to $(i,j)$. Thus the viscous terms are being neglected for the time being (compare 8.57-8.59.)

  2. From the Poisson equation 8.58 the streaming function is also determined at half-step time, using diagonal differencing (see Section 1.3.4):

    $ u_{i,j+1}^{n+1/2}+u_{i,j-1}^{n+1/2}+u_{i+2,j+1}^{n+1/2}+u_{i+2,j-1}^{n} -4u_{i+1,j}^{n+1/2}=-w_{i+1,j}^{n+1/2} (\Delta l)^{2} $

  3. Now follows the integration step proper, the viscous term included:

    $ \begin{eqnarray} v_{x,i,j}^{n+1/2}&=&\frac{\textstyle 1}{\textstyle 2\Delta l} \left(u_{i,j+1}^{n+1/2}-u_{i,j-1}^{n+1/2}\right) \\ v_{y,i,j}^{n+1/2}&=&-\frac{\textstyle 1}{\textstyle 2\Delta l} \left(u_{i+1,j}^{n+1/2}-u_{i-1,j}^{n+1/2}\right) \\ && \\ w_{i,j}^{n+1}&=&w_{i,j}^{n} -\frac{\textstyle \Delta t}{\textstyle 2 \Delta l}v_{x,i,j}^{n+1/2} \left( w_{i+1,j}^{n+1/2}-w_{i-1,j}^{n+1/2}\right) -\frac{\textstyle \Delta t}{\textstyle 2 \Delta l}v_{y,i,j}^{n+1/2} \left( w_{i,j+1}^{n+1/2}-w_{i,j-1}^{n+1/2}\right) \\ && \;\;\;\; +\frac{\textstyle \nu\Delta t}{\textstyle 2(\Delta l)^{2}} \left( w_{i+1,j-1}^{n}+w_{i-1,j-1}+ w_{i-1,j+1}^{n}+w_{i+1,j+1}-4w_{i,j}^{n} \right) \end{eqnarray} $

Stability is again governed by the CFL condition
$ \Delta t \leq \frac{\textstyle 2 \Delta l}{\textstyle \sqrt{2}v_{max}} $

The presence of diffusive terms implies the additional restriction

$ \Delta t \leq \frac{\textstyle (\Delta l)^{2}}{\textstyle \nu} $

8.2.2 Pressure Method

Take the divergence (instead of rotation) of

$ \frac{\textstyle \partial \vec{v}}{\textstyle \partial t} + (\vec{v} \cdot \nabla )\vec{v}= - \nabla \bar{p} + \nu \nabla^{2}\vec{v} $

and use

$ \nabla \cdot \left( \vec{v} \cdot \nabla \right)\vec{v} = (\nabla \vec{v}):(\nabla \vec{v}) $

(with ${A}:{B} \equiv \sum_{i}\sum_{j}A_{ij}B_{ji}$) to obtain the basic equations

$ \begin{eqnarray} \frac{\textstyle \partial \vec{v}}{\textstyle \partial t} + (\vec{v} \cdot \nabla )\vec{v}&=& - \nabla \bar{p} + \nu \nabla^{2}\vec{v} \\ \nabla^{2} \bar{p}& =&-(\nabla \vec{v}):(\nabla \vec{v}) \end{eqnarray} $     (8.65-8.66)

Two-dimensional case, explicit:

$ \begin{eqnarray} \frac{\textstyle \partial v_{x}}{\textstyle \partial t}&=& -\frac{\textstyle \partial \bar{p}}{\textstyle \partial x}+\nu \left[\frac{\textstyle \partial^{2}v_{x}}{\textstyle \partial x^{2}} +\frac{\textstyle \partial^{2}v_{x}}{\textstyle \partial y^{2}}\right] -\frac{\textstyle \partial v_{x}^{2}}{\textstyle \partial x} -\frac{\textstyle \partial v_{x}v_{y}}{\textstyle \partial y} \\ \frac{\textstyle \partial v_{y}}{\textstyle \partial t}&=& -\frac{\textstyle \partial \bar{p}}{\textstyle \partial y}+\nu \left[\frac{\textstyle \partial^{2}v_{y}}{\textstyle \partial x^{2}} +\frac{\textstyle \partial^{2}v_{y}}{\textstyle \partial y^{2}}\right] -\frac{\textstyle \partial v_{y}^{2}}{\textstyle \partial y} -\frac{\textstyle \partial v_{x}v_{y}}{\textstyle \partial x} \\ \frac{\textstyle \partial^{2}\bar{p}}{\textstyle \partial x^{2}} +\frac{\textstyle \partial^{2}\bar{p}}{\textstyle \partial y^{2}}&=& -\left[\left(\frac{\textstyle \partial v_{x}}{\textstyle \partial x}\right)^{2} +2\left(\frac{\textstyle \partial v_{x}}{\textstyle \partial y}\right) \left(\frac{\textstyle \partial v_{y}}{\textstyle \partial x}\right) +\left(\frac{\textstyle \partial v_{y}}{\textstyle \partial y}\right)^{2} \right] \end{eqnarray} $     (8.67-8.69)

$\longrightarrow$ Apply finite difference scheme.

Note: Divergence condition $\nabla \cdot \vec{v}=0$ must stay intact in the course of the calculation.
$\longrightarrow$ Harlow and Welch procedure [Harlow 65], see also [Potter 80]:
-- Let the pressure $p_{i,j}$ be localized at the
centers of the Euler cells
-- Velocity components $v_{x,i,j}$ and $v_{y,i,j}$ are localized at the right and upper box sides, respectively

Figure 8.1: Grid structure in the pressure method

Now approximate the divergence of the velocity by

$ D_{i,j} \equiv \frac{\textstyle 1}{\textstyle \Delta l}\left[v_{x,i,j}-v_{x,i-1,j} \right] + \frac{\textstyle 1}{\textstyle \Delta l}\left[v_{y,i,j}-v_{y,i,j-1} \right] $     (8.70)

or, in "geographical" notation,

$ D_{C} \equiv \frac{\textstyle 1}{\textstyle \Delta l}\left[v_{x,C}-v_{x,W} \right] + \frac{\textstyle 1}{\textstyle \Delta l}\left[v_{y,C}-v_{y,S} \right] $     (8.71)

$\longrightarrow$ $D_{C}=0$ for vanishing divergence.

Apply the Lax scheme to (8.67-8.68) (all terms on the r.h.s. are taken at time $t_{n}$):

$ \begin{eqnarray} v_{x,C}^{n+1}&=&\frac{\textstyle 1}{\textstyle 4}\left[v_{x,N}+v_{x,E}+v_{x,S}+v_{x,W} \right] -\frac{\textstyle \Delta t}{\textstyle 2\Delta l}\left[v_{x,E}^{2}-v_{x,W}^{2} \right] \\ &&-\frac{\textstyle \Delta t}{\textstyle 2\Delta l}\left[\frac{\textstyle 1}{\textstyle 2}\left(v_{y,E}+v_{y,C}\right) \left(v_{x,N}+v_{x,C}\right)-\frac{\textstyle 1}{\textstyle 2}\left(v_{y,S}+v_{y,SE}\right) \left(v_{x,S}+v_{x,C}\right)\right] \\ &&-\frac{\textstyle \Delta t}{\textstyle \Delta l}\left(\bar{p}_{E}-\bar{p}_{C} \right) +\frac{\textstyle nu \Delta t}{\textstyle (\Delta l)^{2}} \left(v_{x,N}+v_{x,E}+v_{x,S}+v_{x,W}-4v_{x,C} \right) \\ && \\ v_{y,C}^{n+1}&=&\frac{\textstyle 1}{\textstyle 4}\left[v_{y,N}+v_{y,E}+v_{y,S}+v_{y,W} \right] -\frac{\textstyle \Delta t}{\textstyle 2\Delta l}\left[v_{y,N}^{2}-v_{y,S}^{2} \right] \\ &&-\frac{\textstyle \Delta t}{\textstyle 2\Delta l}\left[\frac{\textstyle 1}{\textstyle 2}\left(v_{x,N}+v_{x,C}\right) \left(v_{y,E}+v_{y,C}\right) -\frac{\textstyle 1}{\textstyle 2}\left(v_{x,NW}+v_{x,W}\right)\left(v_{y,W}+v_{y,C}\right)\right] \\ &&-\frac{\textstyle \Delta t}{\textstyle \Delta l}\left(\bar{p}_{N}-\bar{p}_{C} \right) +\frac{\textstyle nu \Delta t}{\textstyle (\Delta l)^{2}} \left(v_{y,N}+v_{y,E}+v_{y,S}+v_{y,W}-4v_{y,C} \right) \end{eqnarray} $     (8.72-8.73)

Insert the new velocity components in (8.71) to find

$ \begin{eqnarray} D_{C}^{n+1}&=& \frac{\textstyle 1}{\textstyle 4}\left(D_{N}^{n}+D_{E}^{n}+D_{S}^{n}+D_{W}^{n} \right) -\frac{\textstyle \Delta t}{\textstyle 2(\Delta l)^{2}}S_{C}^{n} \\ &&-\frac{\textstyle \Delta t}{\textstyle (\Delta l)^{2}} \left(\bar{p}_{N}^{n} +\bar{p}_{E}^{n}+\bar{p}_{S}^{n}+\bar{p}_{W}^{n} -4\bar{p}_{C}^{n} \right) \\ &&+\frac{\textstyle nu \Delta t}{\textstyle (\Delta l)^{2}} \left(D_{N}^{n}+D_{E}^{n}+D_{S}^{n}+D_{W}^{n}-4D_{C}^{n} \right) \end{eqnarray} $     (8.74)


$ \begin{eqnarray} S_{C} &\equiv&\left(v_{x,E}^{2}-v_{x,C}^{2}-v_{x,W}^{2}+v_{x,WW}^{2} \right) +\left(v_{y,N}^{2}-v_{y,C}^{2}-v_{y,S}^{2}+v_{y,SS}^{2} \right) \\ &&+\frac{\textstyle 1}{\textstyle 2}\left(v_{y,E}+v_{y,C}\right)\left(v_{x,N}+v_{x,C}\right) -\frac{\textstyle 1}{\textstyle 2}\left(v_{y,S}+v_{y,SE}\right)\left(v_{x,S}+v_{x,C}\right) \\ &&-\frac{\textstyle 1}{\textstyle 2}\left(v_{x,NW}+v_{x,W}\right)\left(v_{y,C}+v_{y,W}\right) +\frac{\textstyle 1}{\textstyle 2}\left(v_{x,W}+v_{x,SW}\right)\left(v_{y,ES}+v_{y,SW}\right) \end{eqnarray} $     (8.75)

Now solve the Poisson equation (8.69) to obtain the pressures.

Problem: After applying the Poisson solver the pressures contain small errors which cause the central values $D_{i,j}^{n}$ and $D_{i,j}^{n+1}$ deviate from zero. The simple ansatz

$ \bar{p}_{N}+\bar{p}_{E}+\bar{p}_{S}+\bar{p}_{W}-4\bar{p}=-S_{C} $     (8.76)

for the pressures is therefore not usable. Rather, we write

$ \begin{eqnarray} \bar{p}_{N}+\bar{p}_{E}+\bar{p}_{S}+\bar{p}_{W}-4\bar{p}_{C}&=&-S_{C} +\frac{\textstyle (\Delta l)^{2}}{\textstyle 4\Delta t}\left(D_{N}+D_{E}+D_{S}+D_{W} \right) \\ &&+\nu \left(D_{N}+D_{E}+D_{S}+D_{W}-4D_{C}\right) \end{eqnarray} $     (8.77)

Stability: Again, the conditions are

$ \Delta t \leq \frac{\textstyle \Delta l}{\textstyle \sqrt{2}|v|_{max}}\,\;\;\;\;
{\rm and} \;\;\; \Delta t \leq \frac{\textstyle 1}{\textstyle 2} \, \frac{\textstyle (\Delta l)^{\textstyle 2}}
{nu} $     (8.78)

8.2.3 Free Surfaces: Marker-and-Cell Method (MAC)

Introduce appropriate boundary conditions to handle such an open surface: "marker" particles distinguish between liquid-filled and empty Euler cells.
The hydrodynamic equations with gravity,

$ \begin{eqnarray} \frac{\textstyle \partial \vec{v}}{\textstyle \partial t} &=& -\nabla \bar{p}-(\vec{v}\cdot \nabla) \vec{v} + \nu \nabla^{2}\vec{v}+\vec{g} \\ \nabla \cdot \vec{v}&=&0 \end{eqnarray} $     (8.79-8.80)

are integrated by any of the foregoing techniques (usually the pressure method). The marker particles in each cell move along according to $\vec{r}^{n+1}=\vec{r}^{n}+\vec{v}^{n}\Delta t$, where $\vec{v}^{n}$ is a particle velocity calculated by interpolation (with suitable weights) from the velocities $v_{x},v_{y}$ in the adjacent Euler cells [Harlow 65].

Treatment of boundary cells: Four possible types -- see Figure 8.2. for the respective velocity boundary conditions. The pressure boundary conditions are the same in all cases: $p=p_{vac}$, where $p_{vac}$ is the "external" pressure in the empty Euler cells.

Figure 8.2: MAC method: the 4 types of surface cells and the appropriate boundary conditions for $v_{x}, v_{y}$ (see [Potter 80].)

vesely 2006

< >