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$.
- 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.)
- 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}
$
- 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)
with
$ \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
|