Applied Mathematics AM 258

AM 258

- Computational Fluid Dynamics -

Instructor: Professor George Karniadakis

CONTENTS

(Chapter :1 - Chapter: 9)

Chapter 1:: Introduction (2 lectures)
$\bullet$ Navier-Stokes: (u,p), $(u,\omega)$ , Conservative Form, Approximations, B.C., Mathematical Nature, Computational Complexity.
Chapter 2:: The Finite Difference Method (4 lectures)
$\bullet$ Taylor Expansions, Difference Operators, Explicit Formulas, Fonberg's Construction, Implicit Formulas, Lele's Approach, Example, Thomas's Algorithm for Periodic Matrices, Two-dimensional Stencils.
Chapter 3:: Time-Stepping Algorithms (2 lectures)
$\bullet$ IVP, Consistency and Stability, Multi-step Methods (Derivation, Stability Regions), Polynomial Stability Methods.
Chapter 4:: Parabolic Equation (2 lectures)
$\bullet$ Explicit/Implicit Schemes, Discrete Perturbation Analysis, von Neumann Analysis, Matrix Methods, Spectra of Diffusion Operator, Iterative Methods of Elliptic Equations (Jacobi, Gauss Seidel, $\ldots$ ).
Chapter 5:: Hyperbolic Equations
$\bullet$ Linear Advection Equations, Review of Basic Properties of PDEs, Euler-Forward/Center-Difference Scheme, Euler-Forward Upwind-Differencing Scheme, Lax-Friedrichs Scheme, Second-Order Upwind Scheme, Lax-Wendroff Scheme, Crank-Nicolson/Center Differencing Scheme, Adam-Bashford/Center-Differencing Scheme, Effects of Boundary Conditions, Finite Element Discretizations.
Chapter 6:: Convection-Diffusion Equation
$\bullet$ The Peclet Number (Physical), Crank-Nicolson/Center Difference Scheme, Euler-Forward/Center Difference: Stability Condition, The (Numerical) Peclet Number Controversy; Mixed Explicit Schemes, Semi-Implicit Schemes, Boundary Conditions Normal Boundary Layers.

Chapter 7:: Vorticity-Based Formulations
Vorticity Streamfunction Formulation $(\Psi - \omega)$ : Governing Equations, Boundary Condition (Thom's Pearson's, Woods's) 1D Model $\Psi - \omega$ Formulation: Explicit Euler, Semi-Implicit, Full-Implicit; Green's Functions; Vorticity, Velocity Formulations: Green's Functions for 1D Model Problem, Penalty Method; Application in 2D: Flow Over a Backward Facing Step, Streamfunction-Only Formulation.
Chapter 8:: Navier-Stokes in Primitive Variables
Governing Equations-Pressure Interpretation, Stokes Model Problem, Staggered Mesh/Implicit Discretization: Discrete Poisson Equation, Green's Functions Implementation, MAC Grid; Non-Staggered Mesh (Implicit Discretization): Consistent 2D Poisson Equation, Discretization of Stoke Model Problem, Pressure Boundary Conditions, Classical Splitting Scheme, High-Order Splitting Scheme, Other Splitting Schemes; Staggered Mesh (Explicit Discretization), The Artificial Compressibility Methods: Basic Idea, MAC Discretization, The Choice of Parameter C, The Turkel Preconditioner, Stability of Artificial Compressibility Method, The Original Chorin Scheme.
Chapter 9:: Examples/Applications

Chapter 1: Introduction

CFD is often confused with Numerical Analysis and sometimes with Computer Science or Fluid Dynamics but it is perhaps the intersection of all three:

$\begin{figure} \centerline{\psfig{file=lec1.1.eps,height=1in}} \end{figure}$
There is no official birth date of CFD but the 1965 article of Harlow and Fromm in Scientific American (vol. 212, No. 3, pp 104-140) set the stage for computer experiments in Fluid Dynamics and drew the analogy with the modern ``thought experiments'' of Fermi and Einstein.
You are supposed to know, for example, the eigenvalues of

$\begin{displaymath}\left [ \begin{array}{llll} 1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1\\ 1 & 1 & 1 & 1\end{array}\right ]\end{displaymath}$

and properties of linear systems but you should also be able to distinguish between convective and diffusive time scales in order to develop appropriate non-dimensionalizations.
Also, simultaneous advances in computer hardware often suggest approriate algorithms. Case in point is the Jacobi algorithm, obviously obsolete on serial computers in favor of the Gauss-Seidel algorithm in the seventies. However, it was re-discovered in the eighties since straight Gauss-Seidel (G-S) is not easily parallelizable (although the black-red G-S is!)
In this course we will learn how to numerically formulate and solve the incompressible Navier-Stokes equations in two- and three-dimensions, covering fundamentals of different discretizations algorithms as well as formulations.

Incompressible Navier-Stokes Equations

I. Velocity-Pressure Formulation (primitive variables):

$\begin{figure} \centerline{\psfig{file=lec1.2.eps,width=3in}} \end{figure}$

Eulerian description

$\begin{displaymath}\hspace*{1in} \left \{ \begin{array}{llll} \displaystyle{\fr... ...{\rm constrained set of PDEs/mixed IVP/BVP}}\end{array}\right .\end{displaymath}$

where

$\begin{displaymath}\tau_{ij} = \mu (\frac{\partial u_j}{\partial x_i} + \frac{\partial u_i}{\partial x_j})\end{displaymath}$

is the stress tensor for incompressible fluids.

*

Pressure is not a thermodynamic quantity, i.e. $p=C\rho^{\gamma}$ is not valid. It is a constraint, projecting the solution $\vec{u}(x,y,z,t)$ onto a divergence-free space. Otherwise it would overdetermine the system. Note that $\sqrt{\frac{dp}{d\rho}} \rightarrow \infty$ (sound speed). Therefore, disturbances propagate infinitely fast in incompressible media.

*

Viscous terms (assuming linear Newtonian shear stress law): It may be that $\nu =\mu/\rho = \nu (x,y,z,t)$ , then the full stress tensor should be used instead of $\nu \nabla^2 \vec{u}$ , i.e. $\nabla \cdot [\nu(\nabla \vec{u} + \nabla^{T} \vec{u})]$ . For example, in LES formulations, $(\nu = \nu_{0} + \nu_{e})$ , or temperature-dependent viscosity in stability studies.

*

Convective/Advective Terms: $\frac{D\vec{u}}{Dt}$ : =

Non-conservative: $\frac{\partial \vec{u}}{\partial t} + (\vec{u} \cdot \nabla)\vec{u}$ , convective-form
Consevative (Lax '54): $\frac{\partial \vec{u}}{\partial t} + \nabla \cdot (\vec{u}\vec{u})$ , flux-form, (basis of finite volume schemes)
Rotational form: $\frac{\partial \vec{u}}{\partial t} - \vec{u} \times \vec{\omega} + \nabla (1/2 \vert u\vert^{2})$
Skew-symmetric: $\frac{\partial \vec{u}}{\partial t} + \frac{1}{2} [ (\vec{u}\cdot \nabla )\vec{u} + \nabla \cdot (\vec{u}\vec{u})]$ , used for transport of species, e.g. temperature.

Note that the continuum forms are all equivalent but semi-discrete forms may differ.

* Semi-discrete forms for Non-Linear Terms in incompressible Navier-Stokes:

Ref: ``On the rotation and skew-symmetric forms for incompressible flow simulation,'' T.A. Zang, Appl. Num. Math.

For $\nu = 0$ (inviscid, frictionless) flows and no source term $(\vec{F} = 0)$ , the N-S conserve linear momentum $\int u\; d\; \Omega$ and kinetic energy $\int\frac{1}{2} \vert u\vert^{2}d \Omega$ . However, for semi-discrete system (i.e. continuous in time but discrete in space) we have:

Integral Forms (Galerkin) conserve both.
Pointwise forms (collocation):

$\rightarrow$
rotational/skew-symmetric: conserve both
$\rightarrow$
divergence (flux): conserves only linear momentum
$\rightarrow$
convection (non-cons): unstable, conserves neither
$\rightarrow$
alternate (divergence & convective): stable

*Computational Complexity:

-: Rotation: 6 derivatives
-: Skew- symm: 18 derivatives
-: Convective/divergence/alternate: 9 derivatives

*Conservation Property of Skew-Symmetric Form:

$\begin{displaymath}\frac{\partial \phi}{\partial t} + {\cal L}_{N} \phi = 0 \Rightarrow (\phi,\phi_{t})+ (\phi,{\cal L}_{N} \phi) = 0\end{displaymath}$

But

$\begin{figure} \centerline{\psfig{file=lec1.3.eps,width=3in}} \end{figure}$

$\begin{eqnarray*}(\phi, {\cal L}_{N}\phi) = - ({\cal L}_{N}\phi, \phi) = 0 & \Ri... ...}) = 0\\ & \Rightarrow & \int \phi^{2} d\Omega = {\rm const.} \end{eqnarray*}$

Note: Skew symmetric matirx

$\begin{displaymath}A = - A^T \;(n\times n)\end{displaymath}$

Note: If n = odd $\Rightarrow$

$\begin{displaymath}\Rightarrow \vert A\vert = 0\end{displaymath}$

Review: The general Navier-Stokes equations are:

$\begin{displaymath}\rho \frac{D\vec{u}}{Dt} = -\nabla p + \nabla \cdot \stackrel{=}{\tau} + \rho \vec{F}\end{displaymath}$

Newtonian:

$\begin{displaymath}\stackrel{=}{\tau}_{ij} = \mu \left (\frac{\partial u_{j}}{\p... ...x_{j}}\right ) + \lambda (\nabla \cdot \vec{u}) \delta _{ij}\end{displaymath}$

Stokes hypothesis: $2\mu + 3\lambda = 0$ , from experimental evidence (equilibrium). Thus,

$\begin{displaymath}\stackrel{=}{\tau}_{ij} = \mu \left[ \left (\frac{\partial u_... ...ight ) - \frac{2}{3} (\nabla \cdot \vec{u}) \delta_{ij}\right ]\end{displaymath}$

constant $\mu$ :

$\begin{displaymath}\nabla \cdot \stackrel{=}{\tau} = \mu \left [\nabla^{2} \vec{... ...{u})\right ] = \mu \nabla^2 \vec{u}, \; {\rm (incompressible)}\end{displaymath}$

* Viscoelastic fluids: $\stackrel{=}{\tau} = \tau^{p} {\rm (polymeric) + \tau^{s}}\;{\rm (Newtonian)}$ ; $\tau^{s} = \mu (\nabla u + \nabla^{T} u)$

In Non-Newtonian fluids (e.g. blood flow), we then need constitutive equations for $\tau^{p}$ , the polymeric contributions, which are typically of hyperbolic nature, e.g.,

$\begin{displaymath}\lambda \left \{ \frac{\partial \tau^p}{\partial t} + \vec{u}... ...}\right \} = \mu_{p} [ \nabla \vec{u} + \nabla^{T} \vec{u}],\end{displaymath}$

where $\mu_p$ is the polymeric shear viscosity.

II. Velocity-Vorticity Formulation: $\vec{\omega} = \nabla \times \vec{u}$

$\begin{displaymath}\frac{D\vec{\omega}}{Dt} = (\vec{\omega}\cdot \nabla)\vec{u} + \nu \nabla^{2} \vec{\omega}\end{displaymath}$

$\begin{displaymath}(*) \left \{ \begin{array}{lll} &\displaystyle{\frac{\partia... ...l \Omega} \vec{u} \cdot \hat{n} \; ds = 0} \end{array}\right .\end{displaymath}$

System (*) is equivalent to the $(\vec{u},p)$ formulation. The proof, for simply-connected domains, is done in 3 steps.

1.: $\nabla \cdot \vec{\omega} = 0\; {\rm in}\; \Omega$ ?
2.: $\vec{\omega} = \nabla \times \vec{u} \; {\rm on}\; \partial \Omega \Rightarrow \vec{\omega} = \nabla \times \vec{u} \; {\rm in}\; \Omega$
3.: $\nabla\cdot \vec{u} = 0\; {\rm in}\; \Omega$ ?

$\begin{displaymath}\mbox{{\rm Let us define}}\; Q \equiv \nabla \cdot \omega \ri... ...rrow \frac{\partial Q}{\partial t} = \nu \nabla^{2} Q\leqno{1.}\end{displaymath}$

$\begin{displaymath}\partial \Omega \colon \; Q = \nabla \cdot \omega =\nabla \cd... ...\times \vec{u}) = 0 \Rightarrow Q \equiv 0 \; {\rm in}\; \Omega\end{displaymath}$

$\begin{displaymath}\nabla^{2} \vec{u} = - \nabla \times \vec{\omega} = \nabla (\... ...cdot \vec{u}) - \nabla \times (\nabla \times \vec{u})\leqno{2.}\end{displaymath}$

$\begin{displaymath}\rightarrow \nabla \times [\vec{\omega} - \nabla \times \vec{... ...imes \vec{u}] = - \nabla \times \nabla (\nabla \cdot \vec{u})=0\end{displaymath}$

$\begin{displaymath}\rightarrow \nabla \{\nabla \cdot [\vec{\omega} - \nabla \tim... ...vec{u}]\} - \nabla^{2} [\vec{\omega} - \nabla \times \vec{u}]=0\end{displaymath}$

$\begin{displaymath}\left . \begin{array}{ll} \Rightarrow \nabla^{2} (\vec{\omeg... ...arrow \vec{\omega} = \nabla \times \vec{u}\; {\rm in}\; \Omega \end{displaymath}$

$\begin{displaymath}\left .\begin{array}{lll} \nabla^{2} \vec{u} & = &\nabla (\n... ...tarrow \nabla \cdot \vec{u} = 0\; {\rm in} \; \Omega.\leqno{3.}\end{displaymath}$

III. Non-Dimensionalization:

$\begin{figure} \centerline{\psfig{file=lec1.4.eps,width=3in}} \end{figure}$

We define $\bar{u} = \frac{u}{u_{0}}; \; \bar{x} = \frac{x}{L_{0}}; \; \bar{t} = \frac{t}{(L/u_{0})}$

$\bullet$ Convective time scale

$\begin{displaymath}\bar{\omega} = \frac{\omega}{u_{0}/L_{0}}\end{displaymath}$

then $\frac{\partial \omega}{\partial \bar{t}} +\nabla \cdot (\vec{u} \omega) = Re^{-1} \nabla^{2} \omega$ in 2D, which is appropriate for modest to large Re, where

$\begin{displaymath}Re = \frac{u_{0}L_{0}}{\nu} = \frac{{\rm inertial}}{{\rm viscous}}.\end{displaymath}$

For $Re \rightarrow 0$ , we introduce the diffusive time scale $t^{'} = \frac{t}{L^{2}/\nu}$ , therefore the 2D-vorticity equation is: $\frac{\partial \omega }{\partial t} + Re \nabla \cdot (\vec{u}\omega) = \nabla^{2} \omega$ , which is more appropriate for low Re (creeping) flows.

*Velocity-pressure: $\frac{\partial u}{\partial t} + u\cdot \nabla u = - \frac{\nabla p}{\rho} + \nu \nabla^{2} u$ (dimensional)

* $\frac{L_{0}}{u^{2}_{0}} \;\; \left (\frac{\partial\bar{u}}{\partial \bar{t}} + ... ...{0}L_{0}} \nabla ^{2} \bar{u} = - \nabla \bar{p} + Re^{-1} \nabla^{2} \bar{u}$ , where $\bar{p} = \frac{p}{\rho u^{2}_{0}}$

In summary: convective time units (fast): $\left (\frac{t}{L/u_{0}}\right )$ , and difffusive time scale (slow): $\left ( \frac{t}{L^{2}/\nu}\right )$

IV. Boundary Conditions:

$\begin{figure} \centerline{\psfig{file=lec1.5.eps,width=3in}} \end{figure}$

$\bullet$ No-slip (Dirichlet) U_s = U_w implies perfect momentum exhange of molecules with wall, molecules stick on the wall. Valid for continuum media. ( $\lambda \ll L_{0}$ ; the mean free path (in gases) much less than characteristic length scale)

Define $K_{n} \equiv \frac{\lambda}{L}$ , Knudsen number

If K_n < 10^-3, continuum $\Rightarrow$ Direct Methods, (e.g., finite differences, finite elements, spectal methods). However, use Monte Carlo, Boltzmann for K_n > 10^-3; rarefied allows slip and thus, model b.c. as

$\begin{displaymath}U_{s} -U_{w} = K_{n} \frac{\partial U_{s}}{\partial \hat{n}}.\end{displaymath}$

Applications: high-altitude, high-speed flows $(\lambda$ Large), micro-mechanis (L small) $\sim 1\mu m.$ ( $\lambda =$ 65 nm at STP). $\bullet$ No-penetration (Dirichlet) For viscous or inviscid flows, in the norml direction we always have

$\begin{displaymath}U_s \cdot \hat{u} = U_w \cdot \hat{n}\end{displaymath}$

$\bullet$ Flux (Neumann b.c.)

$\begin{displaymath}\underbrace{ \frac{\partial \theta}{\partial n} }_{\nabla \theta\cdot \hat{n}} = \underbrace{g(x,y,z,t)}_{{\rm known}}\end{displaymath}$

$\begin{figure} \centerline{\psfig{file=lec1.6.eps,height=1in}} \end{figure}$

$\bullet$ Robin b.c. - (heat transfer; mass transfer) $k \frac{\partial\theta}{\partial n} = h (\theta-\theta_{\infty})$ , where h = heat transfer coefficient.

$\bullet$ Outflow boundary conditions: Flow past an airfoil:

$\begin{figure} \centerline{\psfig{file=lec1.7.eps,height=2in}} \end{figure}$

outflow:

mathematically required, not physically
unknown state
extrapolate, e.g. $\frac{\partial u}{\partial n} = 0$ (with $p=0 \Rightarrow$ zero normal stress)
$\frac{\partial u}{\partial t} + (u\cdot \nabla) u=0$ in a sponge or a buffer layer.
parabolize the N.S. making singular the Laplacian operator, e.g. $\nabla \cdot [\nu(x)(\nabla u+\nabla^{\tau}u)]$

$\begin{figure} \centerline{\psfig{file=lec1.8.eps,height=1in}} \end{figure}$

NOTE: More difficult for incompressible than compressible where use of characteristics may work (not rigorous...)

Special Cases

2D Flow (does it exist?)

$\begin{figure} \centerline{\psfig{file=lec1.9.eps,height=2in}} \end{figure}$

$\bullet$ 2D geometry does not guarantee 2D flow. 3D for Re > Re^c, bifurcation. Instabilities (bifurcations) may develop. Vorticity $\vec{\omega}_{z} = \hat{k} \nabla \times \vec{u}$ , perpendicular to 2D plane in 2D flow. Stretching term $(\vec{\omega} \cdot \nabla) \vec{\omega} \equiv 0$ in 2D thus,

$\begin{displaymath}\frac{\partial \omega_{z}}{\partial t} + (\vec{u} \cdot \nabla ) \omega_{z} = \nu \nabla^{2} \omega_{z}\end{displaymath}$

Define $\psi$ : (stream function) 4th order for $\psi$ (3rd for BL): $\vec{u} = -\nabla \times (\psi \hat{k})$ , and thus $u=- \frac{\partial \psi}{\partial y} ; v = \frac{\partial \psi}{\partial x}$ Continuity: $\nabla \cdot \vec{u} \equiv 0$ is satisfied automatically and $\nabla^{2} \psi = - \omega$ , expresses a kinematic constraint (it is not a time evolution equation).

For $\nu = 0\colon \; \frac{\partial \omega}{\partial t} + \left ( - \frac{\partial ... ...ght ) = 0 \Rightarrow \frac{\partial \omega}{\partial t} + J(\psi, \omega) = 0$

$\bullet$ steady state $\Rightarrow J(\psi, \omega) = 0 \Rightarrow \omega = f_{n} (\psi)$ that is, in inviscid 2D steady flow, vorticity is constant along steamlines. Streamlines $\equiv$ isovorticity lines

* Boundary Conditions: Too many for $\psi$ /Thom's b.c. None for $\omega$ - create b.c. through local expansions and use the excessive $\psi$ b.c.

$\bullet$ Advantage: Incompressibility, the main headache in incompressible formulations, is honored automatically within discretization error. (Very popular method in the seventies.)

Stokes (creeping) flow:

$\begin{displaymath}\left ( \vec{u}\cdot \nabla \right ) \vec{u} \ll \nu \nabla^{2} \vec{u},\end{displaymath}$

i.e. very viscous flows $\rightarrow$ linear equations.

Applications: dispersion, combustion, particle flows small size $\rightarrow Re \ll 1\rightarrow$ approximate as Stokes flow.

parabolic/elliptic problems.

*: Stokes drag of a sphere : $D=6 \pi R \mu U_{\infty}$

$\begin{figure} \centerline{\psfig{file=lec1.10.eps,height=1in}} \end{figure}$

This is the oldest known analytical solution for a creeping flow (Stokes, 1851)
*: Here, we obtain symmetric solutions if geometry is symmetric. Contrast with symmetry-breaking bifurcations for finite Re.

Other complex shapes (particles): Boundary element methods

Inviscid flow:

$\begin{displaymath}\nu \nabla^{2} \vec{u} \ll (\vec{u}\cdot \nabla)\vec{u}, \; {\rm for}\; Re \gg 1\end{displaymath}$

e.g. high-speed flows

$\bullet$: hyperbolic equations
$\bullet$: There is some debate in the current literataure as to whether there exist singular solutions in incompressible flows, i.e., $\omega \sim (t-t_{0})^{m}, \; m < 0$ for t₀ = finite.
$\bullet$: Standard vortex methods (Chorin/Leonard) work well for simulations of inviscid flows.

V. Potential Flow: $\vec{\omega} = 0$ , $\vec{u} = \nabla \phi$ /incompressible, inviscid, irrotational flow, purely elliptic

$\begin{eqnarray*}\nabla^{2} \psi &= &0\\ \nabla^{2} \phi &= & 0\; \mbox{{\rm with B.C.}} \end{eqnarray*}$

does it exist?

$\begin{figure} \centerline{\psfig{file=lec1.11.eps,height=1in}} \end{figure}$

at sufficiently large distances $u\rightarrow u_{\infty}$ in freestream and $\omega \simeq 0$ with $Re \gg 1$

$\bullet$ Compressible:

$\begin{displaymath}\left ( 1 - \frac{u^{2}}{c_{0}} \right ) \frac{\partial \phi}... ...^{2}}{c_0}\right ) \frac{\partial^{2} \phi}{\partial y^{2}} = 0\end{displaymath}$

+ energy equation

$\begin{displaymath}(\vec{u} = \nabla \phi)\end{displaymath}$

unknowns, $(\phi, c_0).$ Boundary Layer Flow

$\begin{figure} \centerline{\psfig{file=lec1.12.eps,height=1in}} \end{figure}$

neglect streamwise diffusion compared to normal (cross) diffusion $\delta /x \ll 1$ , e.g.

$\begin{displaymath}\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial ... ...tial p}{\partial x} + \nu \frac{\partial ^{2}u}{\partial y^{2}}\end{displaymath}$

great simplification, computation complexity in y-only; smaller domain.

*sometimes called the thin shear layer approximation

$\begin{figure} \centerline{\psfig{file=lec1.13.eps,height=2in}} \end{figure}$

In steady state we also call them: Parapolized N-S if the predominant main flow direction equation is parabolic in $x \leftarrow$ plays the role of time $\rightarrow$ advance explicitly. VI. Mathematical Nature of the Flow Equations

The balance of various terms in Navier-Stokes describes dominant diffusion or convection fluxes and thus the mathematical character of N-S will be varying from parabolic to hyperbolic depending on the flow problem.

$\bullet$ unsteady N-S: parabolic-hyperbolic; steady: elliptic-hyperbolic. The divergence-free constraint imposes the hyperbolicity.

*Examples: Stokes equation/creeping flow

$\begin{displaymath}\frac{-u^{2}_{0}T_{0}}{\nu} \frac{\partial u}{\partial t} + \nabla^{2} u=Re \left ( \frac{\partial p}{\partial x}\right )\end{displaymath}$

$\bullet$ steady state $\longrightarrow$ purely elliptic (for fixed $\frac{\partial p}{\partial x}$ ). $\bullet$ unsteady $\longrightarrow$ parabolic.

Outside the B.L. for $\nu = 0$ (high-speed flows) the governing equations are:

$\begin{displaymath}\left \{ \begin{array}{lll} \displaystyle{\frac{\partial u}{... ...frac{1}{\rho} \frac{\partial p}{\partial y}}\end{array}\right .\end{displaymath}$

essentially a hyperbolic equation in space and time describing a propagation phenomenon. Inside the B.L. the governing equations are parabolic/hyperbolic.

PDEs of 2 $^{{\rm nd}}$ -order (Review)

$\begin{displaymath}a \frac{\partial ^{2}\phi}{\partial x ^{2}} + 2b \frac{\parti... ...ial ^{2} \phi}{\partial y^{2}} = 0; \; {\rm {\bf quasi-linear}}\end{displaymath}$

Where $a,b,c = {\cal F} \left (x,y,\phi, \frac{\partial \phi}{\partial x}, \frac{\partial \phi}{\partial y}\right )$ (not 2nd derivatives). For example, $u = \frac{\partial \phi}{\partial x}; v = \frac{\partial \phi}{\partial y}$ , where $\phi$ is the potential.

$\begin{displaymath}\left \{ \begin{array}{ll} \displaystyle{a \frac{\partial u}... ...rtial y} \left ( \begin{array}{cc}u\\ v\end{array} \right ) = 0\end{displaymath}$

$\begin{displaymath}\Rightarrow A_{1} \frac{\partial U}{\partial x} + A_{2} \frac{\partial U}{\partial y} = 0\end{displaymath}$

$\bullet$ Solution: wave of the form $U = \hat{U} e^{i\vec{\eta}\cdot \vec{x}}, \vec{\eta} = (\eta_{x}, \eta_{y})$ , $\vec{\eta} =$ propagation vector. Therefore, by substitution

$\begin{displaymath}(A_{1} \eta_{x} + A_{2} \eta_{y} )\hat{U} = 0 \Rightarrow \det \vert A_{1} \eta_{x} + A_{2} \eta_{y}\vert=0\end{displaymath}$

$\begin{displaymath}\left \vert \begin{array}{ll} a\eta_{x} + 2b\eta_{y} & c\eta... ...ht )^{2} + 2b \left ( \frac{\eta_{x}}{\eta_{y}}\right ) + c = 0\end{displaymath}$

$\begin{displaymath}\begin{array}{ll} \mbox{{\rm Wave-like solution}}\\ \mbox... ...\\ \mbox{{\rm parabolic}} & b^{2}-ac = 0 \end{array}\right .\end{displaymath}$

Example: Transonic-Inviscid Flow/Steady

$\begin{figure} \centerline{\psfig{file=lec1.14.eps,width=3in}} \end{figure}$

Local Mach $\left (\frac{u}{c_{0}}\right )$ variation; c₀ = speed of sound

$\begin{displaymath}\left (\underbrace{1-\frac{u^{2}}{c_{0}^{2}}}_{a}\right ) \fr... ...}^{2}}}_{c}\right ) \frac{\partial^{2}\phi}{\partial y^{2}} = 0\end{displaymath}$

Determinant:

$\begin{displaymath}b^{2}-ac = \frac{u^{2}+v^{2}}{c^{2}_{0}}-1 = M^{2}-1\end{displaymath}$

$\begin{displaymath}\begin{array}{ll} \bullet M < 1, \;\; \mbox{{\rm subsonic fl... ...rm supersonic flow}} \rightarrow \; {\rm hyperbolic}\end{array}\end{displaymath}$

This mixed nature of the transonic potential equation has been a great challenge in CFD since the transition line between the subsonic and the supersonic regions is part of the solution. Note that for $M = 1 \rightarrow$ parabolic.

$\bullet$ Non-linear wave propagation (general model) $U=\hat{U} e^{iS(x,y)}; S (x,y)$ is the wave characteristic surface, if real $\Rightarrow$ hyperbolic system. We thus consider, $\det (A_{1}S_{x} + A_{2} S_{y})=0$ . If we interpret $\vec{\eta} = \nabla S$ , normal to surface = direction of propagation, then $\Rightarrow U = \hat{U} e^{i(\vec{x}\cdot \nabla S)} = \hat{U} e^{i(xS_{x}+yS_{y})}$

Then the tangent vectors of the surface are defined as:

$\begin{displaymath}dS = \vec{\nabla} S\cdot d\vec{x} = \frac{\partial S}{\partial x} dx + \frac{\partial S}{\partial y} dy = 0\end{displaymath}$

Thus, the direction of the chacteristic surface is given by: $\frac{dy}{dx}= - \frac{S_{x}}{S_{y}} = - \frac{\eta_{x}}{\eta_{y}}$

$\begin{figure} \centerline{\psfig{file=lec1.15.eps,height=1in}} \end{figure}$

Note: The surface S(x,y) = const is a wavefront surface that separates the points already influenced by the propagating disturbance from the points not yet reached by the wave. Example: Potential flow for a thin airfoil:

$\begin{displaymath}(1-M^{2}_{\infty} ) \frac{\partial^{2}\phi}{\partial x^{2}} + \frac{\partial ^{2} \phi}{\partial y^{2}} = 0\end{displaymath}$

$\begin{displaymath}\det = 0\Rightarrow a\left ( \frac{\eta_{x}}{\eta_{y}}\right ... ...arrow \frac{\eta_{y}}{\eta_{x}} = \pm \sqrt{M^{2}_{\infty} - 1}\end{displaymath}$

Directon of Characteristic

$\begin{figure} \centerline{\psfig{file=lec1ab.eps,width=3in}} \end{figure}$

$\begin{displaymath}\frac{\eta_{x}}{\eta_{y}} = \mp \frac{1}{\sqrt{M^{2}_{\infty} -1}} = \tan \mu\end{displaymath}$

$\begin{displaymath}\Rightarrow \sin \mu = \frac{1}{M_{\infty}}, \;\; \mbox{{\rm Mach lines}}\end{displaymath}$

*System of $1^{{\rm st}}$ -order equations/Generalization e.g. inviscid flows conservation (flux) form:

$\begin{displaymath}\frac{\partial}{\partial x^{k}} F^{k}_{i} = Q_{i}, \begin{array}{ll}i=1,\ldots n\\ k=1, \ldots m\end{array}\end{displaymath}$

n-equations
m-dimensions (including time)

$\rightarrow$ Quasi-linear form via Jacobian:

$\begin{displaymath}A^{k}_{ij} = \frac{\partial F^{k}_{i}}{\partial u_{j}}\end{displaymath}$

$\begin{displaymath}\Rightarrow A^{k}_{ij} \frac{\partial u_{j}}{\partial x^{k}} ... ...} i,j = 1, \ldots, n\\ k = 1, \ldots , m \end{array}\right.\end{displaymath}$

matrix form: (summation)

$\begin{displaymath}A^{k} \frac{\partial U}{\partial x^{k}}= 0, \; k=1,\ldots, m\end{displaymath}$

A^k, Q = F(x^k,U) do not depend on derivatives of U!! A plane wave solution will exist with propagating direction $\vec{\eta}$ if $[A^{k}\eta_{k}]\hat{U} = 0 \Rightarrow \det (A^{k}\eta_{k})=0$ , where $\eta_{k} =$ components of $\vec{\eta}$ . The solutions $\vec{\eta} = [\eta_{k}], k=1,=\ldots, m$ are (n) altogether (equal to # of 1st order eqn) $\Rightarrow$ define characteristic surfaces from $\eta_{k} = \partial S/\partial x^{k} \Rightarrow$ at most (n) characteristic surfaces.

Definition:

hyperbolic: for all n characteristic surfaces are real and linearly independent
elliptic: all complex
hybrid: some elliptic, some hyperbolic
parabolic: if matrix $(A^{k}\eta_{k})$ is not of rank (n), e.g. u₁ has derivatives w.r.t. x¹ missing $\Rightarrow A^1_{i1} = 0 \;\forall\; i$ .

Example: 2 equations:

$\begin{displaymath}\left \{ \begin{array}{lll} \displaystyle{a\frac{\partial u}... ...} + d \frac{\partial u}{\partial y} = F_{2}}\end{array}\right .\end{displaymath}$

$\begin{displaymath}\Rightarrow \left \vert \underbrace{\begin{array}{ll} a & 0\\... ... ) = \left ( \begin{array}{ll} F_{1}\\ F_{2}\end{array}\right )\end{displaymath}$

$\begin{displaymath}\vert A^{1}\underbrace{\eta_{1}}_{\eta_{x}} + A^{2} \eta_{2} ... ...\\ d & \frac{b\eta_{x}}{\eta_{y}}\end{array} \right \vert = 0\end{displaymath}$

$\begin{displaymath}\Rightarrow \left ( \frac{\eta_{x}}{\eta_{y}}\right )^{2} = \frac{cd}{ab}\end{displaymath}$

if $cd/ab > 0 \Rightarrow$ hyperbolic (e.g. a = b=1, c=d=1),

$\begin{displaymath}% latex2html id marker 1808 F_{1}=F_{2} = 0 \rightarrow \frac... ...\partial ^{2}u}{\partial y^{2}},\; \mbox{{\rm (wave equation)}}\end{displaymath}$
if $cd/ab<0 \Rightarrow$ elliptic $\rightarrow$ (e.g a=b=1 c=-d=-1)

$\begin{displaymath}% latex2html id marker 1818 \rightarrow \frac{\partial ^{2}u}... ...^{2}u}{\partial y^{2}} = 0,\; \mbox{{\rm (Laplace's equation)}}\end{displaymath}$
if $b=0 \Rightarrow$ parabolic, a=1, b=0, c=-d=-1, F₁ =0, F₂ = v

$\begin{displaymath}% latex2html id marker 1824 \Rightarrow \frac{\partial u}{\pa... ...al ^{2}u}{\partial y^{2}}\;\; \mbox{{\rm (diffusion equation)}}\end{displaymath}$

Exercise: Show that the Cauchy-Riemann equations are of elliptic type

$\begin{displaymath}C-R\left ( \begin{array}{ll} \displaystyle{\frac{\partial u}... ...tial x} - \frac{\partial u}{\partial y} = 0}\end{array}\right )\end{displaymath}$

Example: 2D Stationary shallow-water equations

$\begin{displaymath}\left \{ \begin{array}{llll}\displaystyle{u \frac{\partial h}... ...{ccc}h\\ u\\ v\end{array}\right \} \; \mbox{{\rm state vector}}\end{displaymath}$

$\begin{displaymath}\left \vert\underbrace{\begin{array}{ccc}u & h & 0\\ g & u & ... ...= 0;\;\; {\rm where}\; \lambda \equiv \frac{\eta_{x}}{\eta_{y}}\end{displaymath}$

$\begin{displaymath}\left \vert \begin{array}{ccc} u \lambda + v & h\lambda & h\... ...da + v & 0\\ g & 0 & u\lambda + v\end{array} \right \vert = 0\end{displaymath}$

where

$\begin{displaymath}\lambda^{(1)} = - \frac{v}{u}; \;\;\lambda^{(2),(3)} = \frac{-uv \pm \sqrt{u^{2}+v^{2} - gh}}{u^{2} - gh}\end{displaymath}$

$\sqrt{gh} \approx$ ``sonic'' speed, $u^{2} + v^{2} > gh \Rightarrow$ hyperbolic

$\begin{displaymath}\Rightarrow \; \mbox{{\rm else hybrid}}\end{displaymath}$

$\bullet$ Characteristic surface $S^{1} \colon \vec{\eta}^{(y)} \cdot \vec{v} \Rightarrow \eta_{x} \cdot u + \eta_{y} v = (-v)u+(u)v=0$ $\Rightarrow$ characteristic surface $\equiv$ streamline surface.

Example: Inviscid, Incompressible Equations (Euler)/steady Define $U= \left ( \begin{array}{cc} u\\ v\\ p\end{array}\right )$ then

$\begin{displaymath}\left [ \begin{array}{ccc} u & 0 & 1/\rho\\ 0 & u & 0\\ 1 &0... ... 0 & 1 & 0\end{array}\right ] \frac{\partial U}{\partial y} = 0\end{displaymath}$

Thus,

$\begin{displaymath}\det \left [ \begin{array}{ccc} u + v \lambda & 0 & 1/\rho\\ ... ...ambda & \lambda/\rho\\ 1 & \lambda & 0\end{array}\right ] = 0\end{displaymath}$

$\begin{displaymath}\Rightarrow (u+v\lambda)(\lambda^{2} + 1) = 0 \Rightarrow \be... ...m real}\\ \lambda_{2,3} = \pm i,\; {\rm imaginary}\end{array}\end{displaymath}$

Threfore, this system is hybrid since it has mixed eigenvalues.

*Time-dependent flows We single out time from the system, e.g. x^m = t,

$\begin{displaymath}\Rightarrow \frac{\partial U}{\partial t} + A^{k} \frac{\part... ...; \underbrace{k=1,\ldots, m-1}_{\mbox{{\rm spatial variables}}}\end{displaymath}$

$\begin{displaymath}\det [\eta_{t} + \underbrace{A^{k}\eta_{k}}_{K_{ij}}]=0\end{displaymath}$

equivalent to eigenvalue problem

$\begin{displaymath}\det [K - \lambda I] = 0\end{displaymath}$

where

$\begin{displaymath}\eta_{t} = - \lambda\end{displaymath}$

(compare with $U = \hat{U} e^{i(\vec{k}x-\omega t)} \rightarrow$ where, $\omega =$ frequency either real or complex, and $\vec{\eta}_{k} \rightarrow \vec{k}$

Example: 1D Time-dependent shallow water equations

$\begin{displaymath}\left \{ \begin{array}{lll} \displaystyle{\frac{\partial h}{... ...ght . ; U = \left \{ \begin{array}{c} h\\ u\end{array}\right \}\end{displaymath}$

$\begin{displaymath}\Rightarrow \frac{\partial U}{\partial t} + \left \vert \begi... ... u\\ \end{array} \right \vert \frac{\partial U}{\partial x} = 0\end{displaymath}$

$\begin{displaymath}\Rightarrow \left \vert \begin{array}{ccc} u -\lambda & h\\ ... ... 0 \Rightarrow \lambda _{1,2,} = u \pm \sqrt{gh}\; \epsilon \Re\end{displaymath}$

$\begin{displaymath}\Rightarrow \; \mbox{{\rm hyperbolic system}}\end{displaymath}$

Question: What is the 2D time-dependent shallow water? Domain of Dependence/Influence

$\begin{figure} \centerline{\psfig{file=lec4.9.eps,width=4in}} \end{figure}$

$\Theta_{x} + \Theta_{y} = 0$	$\Theta_{x} = \Theta_{yy}$	$\Theta_{xx} + \Theta_{yy} = q$
$\bullet$ hyperbolic	$\bullet$ parabolic	$\bullet$ elliptic
$\bullet$ IVP	$\bullet$ IVP/BVD	$\bullet$ BVP

Computational Complexity of Navier-Stokes

*At high Re Kolmogorov's scaling theory; also Foias's argument

$\begin{figure} \centerline{\psfig{file=lec1.17.eps,width=3in}} \end{figure}$

Assume that L = size of a large eddy and that (u_rms) represents $\{rms\}$ value. Also $\epsilon$ is the dissipation $(\epsilon \sim u^{3}/L)$ and $\eta =$ Kolmogorov's scale, which is the smallest scale of turbulence.

1D: disparity of scales Number of Degrees of freedom $\displaystyle{\sim \frac{L}{\eta} Re^{3/4}}$ $Re_{\eta} = \displaystyle{\frac{\eta u_{rms}}{\nu}} = 1$
3D: scales carry energy $\sim Re^{9/4}$ $\bullet$ Time scales: fast: $\displaystyle{\frac{\eta}{v_{rms}}}$ ; slow: $\displaystyle{\frac{L}{v_{rms}}}$

$\bullet$ Explicit integration in time dictates that we honor the CFL condition, i.e.,

$\begin{eqnarray*}\Delta t & < & CFL< {\cal O} (1)\\ & \Rightarrow & \Delta t < \frac{\eta}{u_{rms}} \end{eqnarray*}$

$\bullet$ Need to integrate for long time $T = \frac{L}{u_{rms}}$ (one large eddy-turnover time). Thus, number of time steps:

$\begin{displaymath}\frac{T}{\Delta t} \sim \frac{L}{\eta} = Re^{3/4}\end{displaymath}$

Thus total work $\sim (Re^{9/4}) (Re^{3/4}) = Re^{3}$ , so increasing Re by a factor of 2 leads to an order of magnitude increase in the CPU requirement.

Recommended Reading:

(1): ``Computer Experiments in Fluid Dynamics,'' Harlow, F.H. and Fromm, J.E., Scientific American, vol. 212, No. 3, pp 104-110.
(2): ``Nodes, Modes and Flow Codes'', Karniadakis, G.E. and Orszag, S.A., Physics Today, March 1993, pp 34-42.
(3): G. Batchelor's book: Derivation of the Navier-Stokes Equations and Conservation Laws
(4): P-L. Lions, Mathematical Topics in Fluid Mechanics, Oxford Science Publications, 1996.
(5): P. Constantin and C. Foias, Navier-Stokes Equations, University of Chicago Press, 1988.
(6): V.A. Solonnikov and A.V. Kazhmikov, ``Existence theorems for the equations of a compressible viscous fluid,'' Ann. Rev. Fluid Mech., vol. 13, p. 79, 1982.