粘性通量离散方法

粘性通量及其离散方法

Viscous Flux

对于一般三维 NS 方程，粘性通量表示为

$$\begin{equation} \bold{F_v} = \begin{bmatrix} 0 \\ n_x \tau_{xx} + n_y \tau_{xy} + n_z \tau_{xz} \\ n_x \tau_{yx} + n_y \tau_{yy} + n_z \tau_{yz} \\ n_x \tau_{zx} + n_y \tau_{zy} + n_z \tau_{zz} \\ n_x \Theta_{x} + n_y \Theta_{y} + n_z \Theta_{z} \\ \end{bmatrix} \end{equation}$$

式中

$$\begin{equation} \begin{gather*} \Theta_{x} = u \tau_{xx} + v \tau_{xy} + w \tau_{xz} + k \frac{\partial T}{\partial x} \\ \Theta_{y} = u \tau_{yx} + v \tau_{yy} + w \tau_{yz} + k \frac{\partial T}{\partial y} \\ \Theta_{z} = u \tau_{zx} + v \tau_{zy} + w \tau_{zz} + k \frac{\partial T}{\partial z} \\ \end{gather*} \end{equation}$$

对于组分方程

$$\begin{equation} \bold{F_v^{(m)}} = - \rho Y_m V_j^{(m)} \end{equation}$$

式中 $V_j^{(m)}$ 为扩散速度

$$\begin{equation} V_j^{(m)} = - \frac{D}{Y_m} \frac{Y_m}{x_j} \end{equation}$$

Diffusion	Transport Poperties	Sign	Law
Momutum	Dynamic Viscosity	$\mu$	Stokes
Heat	Thermal Conductivity	$k$	Fourier
Species	Diffusion Coefficient	$D$	Fick

Viscous Stress

粘性应力通过应力张量 $\bar{\bar{\tau}}$，笛卡尔坐标系中一般形式为

$$\begin{equation} \bar{\bar{\tau}} = \begin{bmatrix} \tau_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & \tau_{yy} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \tau_{zz} \\ \end{bmatrix} \qquad \tau_{ij} \Rightarrow \begin{cases} i = j, & \text{normal stress} \\ i \ne j, & \text{shear stress} \\ \end{cases} \end{equation}$$

对于牛顿流体介质，粘性应力张量 (Stokes) 包含两部分

• 体积膨胀率 volumetric dilatation, $\lambda \text{div}(\vec{v})$
• 线膨胀率 linear dilatation, $\mu(\partial u/\partial x)$

$$\begin{equation} \begin{gather*} \tau_{xx} = \lambda \left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} \right) + 2 \mu \frac{\partial u}{\partial x} \\ \tau_{yy} = \lambda \left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} \right) + 2 \mu \frac{\partial v}{\partial y} \\ \tau_{zz} = \lambda \left( \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} \right) + 2 \mu \frac{\partial w}{\partial z} \\ \tau_{xy} = \tau_{yx} = \mu \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) \\ \tau_{xz} = \tau_{zx} = \mu \left( \frac{\partial u}{\partial z} + \frac{\partial w}{\partial x} \right) \\ \tau_{yz} = \tau_{zy} = \mu \left( \frac{\partial v}{\partial z} + \frac{\partial w}{\partial y} \right) \\ \end{gather*} \end{equation}$$

式中 $\lambda$ 为 second viscosity, $\mu$ 为 dynamic viscosity, 定义 kinematic viscosity

$$\begin{equation} \nu = \mu / \rho \end{equation}$$

为封闭正应力表达式，Stokes 引入 bulk viscosity 假设

$$\begin{equation} \lambda + \frac{2}{3} \mu = 0 \end{equation}$$

上式表征了有限速率体积变化时均匀温度流体中的能量耗散特性，目前除了极端高温或高压情况，尚未有实验证据表明上述关系不成立，因此应力张量可表示为

$$\begin{equation} \tau_{ij} = \mu \left( 2 S_{ij} - \frac{2}{3} \delta_{ij} S_{kk} \right), \quad S_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \end{equation}$$

Heat Conduction

热导率可以通过普朗特数 Prandtl 计算

$$\begin{equation} Pr = \frac{\nu}{\alpha} = \frac{\mu / \rho}{k / \rho c_p} = c_p \frac{\mu}{k} \end{equation}$$

因此

$$\begin{equation} k \frac{\partial T}{\partial x} = \frac{\mu}{Pr} \frac{\partial(c_p T)}{\partial x} = \frac{\mu}{(\gamma - 1) Pr} \frac{\partial c^2}{\partial x} \end{equation}$$

常温下空气 Prandtl 数 $Pr = 0.71$

Discretization

有限差分法离散

$$\begin{equation} \frac{\partial\bold{F_v}}{\partial\xi} = \frac{(\bold{F_v})_{I+1/2} - (\bold{F_v})_{I-1/2}}{\Delta{\xi}} \end{equation}$$

有限体积法

$$\begin{equation} \begin{align*} \frac{(\bold{F_v} S)_{I+1/2} - (\bold{F_v} S)_{I-1/2}}{\Omega_I} = \frac{1}{\Omega_I}(\bold{F_v}(\bold{Q}_{I+1/2}) S_{I+1/2} - \bold{F_v}(\bold{Q}_{I-1/2}) S_{I-1/2}) \\ \end{align*} \end{equation}$$

粘性通量一般采用中心差分格式实现

• 粘性通量描述守恒量扩散过程，物理上各向同性
• 粘性会产生物理耗散，而迎风格式引入数值耗散，会污染边界层发展和湍流耗散
• 扩散方程为抛物型方程，特征速度无穷大，不存在特征方向

Gradient

计算变量梯度一般有两种方法

• 有限差分法 (Finite Difference)
• 格林公式 (Green’s Theorem)

Finite Difference

将笛卡尔坐标系 $(x, y, z)$ 转换为曲线坐标系 $(\xi, \eta, \zeta)$, 然后在曲线坐标系上应用有限差分法

$$\begin{equation} \frac{\partial U}{\partial x} = \frac{\partial U}{\partial \xi} \frac{\partial \xi} {\partial x} + \frac{\partial U}{\partial \eta} \frac{\partial \eta} {\partial x} + \frac{\partial U}{\partial \zeta} \frac{\partial \zeta}{\partial x} \end{equation}$$

Green’s Theorem

基于控制体面构建辅助控制体 $\Omega'$ 再应用格林公式求控制面处的流动变量 $U$ 梯度

$$\begin{equation} \frac{\partial U}{\partial x} = \frac{1}{\Omega'} \int_{\partial \Omega'} U dS_x' = \frac{1}{\Omega'} \sum_{m=1}^{N_F} U_m S_{x,m}' \end{equation}$$

对于 Cell-Centered 网格，辅助控制体如图所示

|           |           |           |
 ----------- ----- ..... ..... -----
|      I-1  |     :I    |     :I+1  |
|     *     |     *     o     *     |
|           |     :     |I+1/2:     |
 ----------- ----- ..... ......-----
|           |           |           |

对于 Cell-Vertex 网格，辅助控制体如图所示

* --------- * --------- * --------- *
|           |           |           |
|            ...........            |
|I-1        :I          :I+1        |
* --------- * --- o --- * --------- *
|           :     I+1/2 :           |
|            ...........            |
|           |           |           |
* --------- * --------- * --------- *

Reference

注意方程采用了两种表示法，后续需要进行统一，个人习惯采用张量表示法