# Computational Fluid Dynamics

Please refresh the page if equations are not rendered correctly.
---------------------------------------------------------------

## 1.Intro

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

\vec{v} : velocity field, having different values in space

\frac{\partial \vec{v} }{\partial t} : unsteady term, indicating the change of velocity with time

(\vec{v}\cdot \nabla )\vec{v} : convective acceleration term, non-linear

\bigtriangledown p : gradient of thermodynamic pressure

\nu \nabla ^{2}\vec{v} : represent the viscous diffusion

Assumption:

• Newtonian fluid: stress and strain rate are linearly related

• The density \rho is constant：no compressibility effect

When we assume that \nu =0, the 2nd order in Navier-Stokes Equation turns to 0(inviscid fluid), which turns Navier-Stokes Equation into the Euler Equation.

In CFD, any assumption will result in errors, it's usual that there exist discrepancies between a CFD solution and an experiment(mathematical model and physical model)

Euler Equation:

\frac{\partial \vec{v} }{\partial t} +(\vec{v}\cdot \nabla )\vec{v}=-\nabla p

The acceration has two components, one component represents how velocity changes with respect to time, the other one represents how the velocity changes because a particle fluid moves from one place to another place(conviction).

### 1.1 The basic ingredients of CFD

#### 1.1.1 different mathematical models

The characteristics of these models/equations:
- have a set of partial differential equations or intergal differential equations

• associate with boundary conditions
• use target applications to simplify the model, like incompressible, inviscid, turbulent, 2d or 3d

#### 1.1.2 choose discretization method

To choose discretization method invovles two aspects:

• discretize geometry: the basic geometry is called grid/mesh(grid gives the vessel of the solution)
• discretize the model: transform the continum mathematical operators into arithmatic operators on the grid values

The discretization method is for approximation the PDEs by a system of algebraic equations, usually we have different discretization methods, as followed:

- finite difference(FD)
- finite volume(FV)
- finite elements(FE)
- spectral methods
- boundary dement method(BEM)
- particle methods

#### 1.1.3 analyze the numerical scheme

We use these schemes to obtain the grid point values of all the flow variables, usually they are
- time dependent situation

All the numercial schemes must satisfy certain conditions to be accepted:
- consistancy

• stability
• convegence

Also we must analyze accuracies of all the schemes.

## ## 2. The derivation of the differential form for the fluid equations

### 2.1.1 Conservation of Mass(a.k.a. Continuity)

• for a system:
\frac{\mathrm{d} M_{sys} }{\mathrm{d} t} =0
• for a control voulme:
\frac{\partial \int\limits_{c.v.}^{}\rho \mathrm{d}v }{\partial t} + \int\limits_{c.s.}^{}\rho \vec{v}\cdot \hat{n} \mathrm{d}tA=0

• $\frac{\partial \int\limits_{c.v.}^{}\rho \mathrm{d}v }{\partial t}$ rate of change of mass in control volume
• $\int\limits_{c.s.}^{}\rho \vec{v}\cdot \hat{n} \mathrm{d}tA$ net rate of mass flow across the control surface

Differential Form of volume- consider a small element \partial x\partial y\partial z
so we get new forms of seperate term above:
\frac{\partial \int\limits_{c.v.}^{}\rho \partial x \partial y \partial z }{\partial x}
Rate of mass flow:
mass flow is in the x direction, and in the center of this cube, the density is equal to \rho, so the x-component quantity is equal to \rho u, and we do a Talyor expansion, then the \rho u becomes (\rho u +\frac{\partial (\rho u)}{\partial x} \cdot \frac{\delta x}{2})\mathrm{d}y \mathrm{d}z

Net rate of mass outflow in x:\frac{\partial (\rho u)}{\partial x} \cdot \mathrm{d}x\mathrm{d}y\mathrm{d}z

similarly,
in y: \frac{\partial (\rho u)}{\partial y} \cdot \mathrm{d}x\mathrm{d}y\mathrm{d}z
in z: \frac{\partial (\rho u)}{\partial z} \cdot \mathrm{d}x\mathrm{d}y\mathrm{d}z

then we have differential form of conservation of mass:
\frac{\partial\rho}{\partial t}+\frac{\partial (\rho u)}{\partial x}+\frac{\partial (\rho u)}{\partial y} +\frac{\partial (\rho u)}{\partial z}=0
vector notation(continuity equation): \frac{\partial\rho}{\partial t}+\nabla \rho\cdot \vec{v} =0

### 2.1.2 Conservation of the Momentum

#### 2.1.2.1 Conservation of Momentum

• for a system:
\vec{F} =\frac{\mathrm{d} \int\limits_{sys}^{} \vec{v}\mathrm{d}m }{\mathrm{d} x}
• for a control volume:
\sum \vec{F_{c.v.} } =\frac{\partial \int\limits_{c.v.}^{} \vec{v} \rho dv}{\partial t}+\int\limits_{c.s.}^{}\vec{v}\rho \vec{v} \cdot \hat{n} \mathrm{d}A

• ${\int\limits_{c.v.}^{} \vec{v} \rho dv}$: element of mass
• $\int\limits_{c.s.}^{}\vec{v}\rho \vec{v} \cdot \hat{n} \mathrm{d}A$: flux term, represents the flux of the momentum

#### 2.1.2.2 Infinitesimal fluid mass \mathrm{d}m

\mathrm{d}\vec{F} =\frac{\mathrm{d} (\vec{v}\mathrm{d}m) }{\mathrm{d}t} = \mathrm{d}m \frac{\mathrm{d} \vec{v}}{\mathrm{d} t}=\mathrm{d}m\cdot \vec{a}

#### 2.1.2.3 Types of forces and stresses:

• Forces:
• body forces: weight \mathrm{d}\vec{F} = \mathrm{d}m\cdot \vec{g}
• surface forces:

• normal: \mathrm{d}F_n
• tangential: \mathrm{d}F_1,\mathrm{d}F_2

有两种力，一种面内，一种面外，两种力两两正交

• Stresses:
• Normal Stress: $\sigma n=\lim{\mathrm{d}A \to 0} \frac{\mathrm{d} \vec{F_n}}{\mathrm{d} A}$
• Shearing Stresses:
• $\tau 1=\lim{\mathrm{d}A \to 0} \frac{\mathrm{d} \vec{F_1}}{\mathrm{d} A}$
• $\tau 2=\lim{\mathrm{d}A \to 0} \frac{\mathrm{d} \vec{F_2}}{\mathrm{d} A}$

With reference to a coordinate system:

Surface forces in terms of stresses:

In X direction：

A pair of normal stresses on the surface：
(\sigma _{xx}-\frac{\partial \sigma _{xx}}{\partial x} \frac{\mathrm{d}x}{2} )\mathrm{d}y \mathrm{d}z
(\sigma _{xx}+\frac{\partial \sigma _{xx}}{\partial x} \frac{\mathrm{d}x}{2} )\mathrm{d}y \mathrm{d}z

A pair of shear stresses on the surface：

(\tau _{zx}-\frac{\partial \tau _{zx}}{\partial z}\frac{\mathrm{d}z }{2} )\mathrm{d}x \mathrm{d}y
(\tau _{zx}+\frac{\partial \tau _{zx}}{\partial z}\frac{\mathrm{d}z }{2} )\mathrm{d}x \mathrm{d}y

A pair of shear stresses on the surface：

(\tau _{yx}-\frac{\partial \tau _{yx}}{\partial y}\frac{\mathrm{d}z }{2} )\mathrm{d}x \mathrm{d}z
(\tau _{yx}+\frac{\partial \tau _{yx}}{\partial y}\frac{\mathrm{d}z }{2} )\mathrm{d}x \mathrm{d}z

X direction:
$\mathrm{d}F_{sx} =(\frac{\partial{} \sigma {xx}}{\partial{x} }+\frac{\partial \tau{yx}}{\partial y}+\frac{\partial \tau_{zx}}{\partial z})\mathrm{d}x \mathrm{d}y \mathrm{d}z Y direction:\mathrm{d}F_{sy} =(\frac{\partial \tau_{xy}}{\partial x}+\frac{\partial{} \sigma {yy}}{\partial{y} }+\frac{\partial \tau{zy}}{\partial z})\mathrm{d}x \mathrm{d}y \mathrm{d}z Z direction:\mathrm{d}F_{sz} =(\frac{\partial \tau_{xz}}{\partial x}+\frac{\partial \tau_{yz}}{\partial y}+\frac{\partial{} \sigma _{zz}}{\partial{z} })\mathrm{d}x \mathrm{d}y \mathrm{d}z We know that\mathrm{d}\vec{F} =\mathrm{d}m\cdot \vec{a}and\mathrm{d}m=\rho\mathrm{d}x \mathrm{d}y\mathrm{d}zThen we have:\rho g_{x}+\frac{\partial{} \sigma {xx}}{\partial{x} }+\frac{\partial \tau{yx}}{\partial y}+\frac{\partial \tau_{zx}}{\partial z}=\rho (\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z} ) similarly in other directions:\rho g_{x}+\frac{\partial \tau_{xy}}{\partial x}+\frac{\partial{} \sigma {yy}}{\partial{y} }+\frac{\partial \tau{zy}}{\partial z}=\rho (\frac{\partial u}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z} ) \rho g_{x}+\frac{\partial \tau_{xz}}{\partial x}+\frac{\partial \tau_{yz}}{\partial y}+\frac{\partial{} \sigma _{zz}}{\partial{z} }=\rho (\frac{\partial u}{\partial t}+u\frac{\partial w}{\partial x}+v\frac{\partial w}{\partial y}+w\frac{\partial w}{\partial z} ) Now we have 3 equations above and continuity equation, there are 4 equations. In these equations, we don't know the three components of velocityu,v,wand all the stresses. To solve these equations ,we need to make some assumptions. Assumption: * For inviscid flow: no shearing stresses and all the normal stresses are equal to-p. Then the Euler Equation is equal to:\rho g_{x}-\frac{\partial{} p}{\partial{x} }=\rho (\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z} ) Vector Notation:\rho \vec{g}-\bigtriangledown p=\rho (\frac{\partial \vec{v}}{\partial t}+(\vec{v}\cdot \bigtriangledown)\vec{v})$

## 3. The Navier-Stokes Equations

### Newtonian Fluid

• Linear relationship between stresses
• Linear relatoinship between rates of deformation

Normal stresses:
\sigma_{xx}=-p+2\mu\frac{\partial u}{\partial x}
\sigma_{yy}=-p+2\mu\frac{\partial v}{\partial y}
\sigma_{zz}=-p+2\mu\frac{\partial w}{\partial z}
\sigma_{xx}+\sigma_{yy}+\sigma_{zz}=-3p

Shear stresses:
\tau_{xy}=\tau_{yx}=\mu(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x})
\tau_{yz}=\tau_{zy}=\mu(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y})
\tau_{zx}=\tau_{xz}=\mu(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z})

Stress terms, x-direction momentum equ:
\frac{\partial \sigma_{xx}}{\partial x}+\frac{\partial \tau_{yx}}{\partial y}+\frac{\partial \tau_{zx}}{\partial z}

## 4. Finite Difference Method(FDM)

### 3.2 FDM方法概述

#### 3.2.1 有限差分法的基本原理

u^{\prime}(x)=\lim _{\Delta x \rightarrow 0} \frac{u(x+\Delta x)-u(x)}{\Delta x} \approx \frac{u(x+\Delta x)-u(x)}{\Delta x} \quad (2)

\left{\begin{array}{l}\left(u_{i+1}-u_i\right) / h+c\left(x_i\right) u_i=f\left(x_i\right) \u_0=d

u_{i+1}=(1-hc_{i})u_{i}+hf_{i}
u_{0}=d

Au=F
\left[\begin{array}{cccc}1&&&\ C_1&1&& \&\ddots&\ddots&\&& C_{N-1}&1\end{array}\right]\left{\begin{array}{c}u_1 \ u_2 \ \vdots \ u_N\end{array}\right}=\left{\begin{array}{c}F_0 \ F_1 \ \vdots \ F_{N-1}\end{array}\right}-\left{\begin{array}{c}C_0 u_0 \ 0 \ \vdots \ 0\end{array}\right}

#### 3.2.2 微分的近似表示

u(x+h)=u(x)+h u^{\prime}(x)+\frac{h^2}{2} u^{\prime \prime}(x)+\cdots+\frac{h^k}{k !} u^{(k)}(x)+\dots \quad(9)

##### 3.2.2.1 一阶微分

u(x+h)=u(x)+h u^{\prime}(x)+\frac{h^2}{2} u^{\prime \prime}(x)+\frac{h^3}{6} u^{(3)}(x+\xi)+\dots \quad(10)