[CFD笔记1] 流体力学控制方程组

《计算流体力学基础及其应用》第二章笔记

最近开始学CFD了，先入个门，看的书是《计算流体力学基础及其应用》，通俗易懂，今天看完了第二章，做个记录；
包含了一些核心知识点和~~简陋的~~推导过程（自我理解），便于自己复习。

流体力学的控制方程组

两个概念

1. 物质导数

假设一运动的流体微团在笛卡尔坐标系下的某一个状态$(x{1},y{1},z{1},t{1})$，速度满足

$\overrightarrow{V}=u\overrightarrow{i}+v\overrightarrow{j}+w\overrightarrow{k}$

密度满足

$\rho_{1} = \rho\left(x_1,y_1,z_1,t_1\right)$

则在另一个状态$(x{2},y{2},z{2},t{2})$下，有：

$\rho_{2} = \rho\left(x_2,y_2,z_2,t_2\right)$

在$1$处使用泰勒展开，得到：

$\rho_2-\rho_1=\frac{\partial \rho}{\partial t}(t_2-t_1)+\frac{\partial \rho}{\partial x}(x_2-x_1)+\frac{\partial \rho}{\partial y}(y_2-y_1)+\frac{\partial \rho}{\partial z}(z_2-z_1)$

两边同时除以$(t_2-t_1)$，得到：

$\frac{\rho_2-\rho_1}{t_2-t_1}=\frac{\partial \rho}{\partial t}+\frac{\partial \rho}{\partial x}\frac{(x_2-x_1)}{t_2-t_1}+\frac{\partial \rho}{\partial y}\frac{(y_2-y_1)}{t_2-t_1}+\frac{\partial \rho}{\partial z}\frac{(z_2-z_1)}{t_2-t_1}$

当$t_2\rightarrow t_1$时，化简上式，得到：

$\frac{D\rho}{Dt}=\frac{\partial\rho}{\partial t}+u\frac{\partial\rho}{\partial x}+v\frac{\partial\rho}{\partial y}+w\frac{\partial\rho}{\partial z}$

其中$\frac{D\rho}{Dt}$为密度的物质导数，它等于密度的当地导数$\frac{\partial\rho}{\partial t}$和迁移导数$u\frac{\partial\rho}{\partial x}+v\frac{\partial\rho}{\partial y}+w\frac{\partial\rho}{\partial z}$（即$\overrightarrow{V}\cdot \overrightarrow{\nabla \rho}$）之和。

因此规定物质导数为$\frac{D}{Dt}$，当地导数为$\frac{\partial}{\partial t}$，迁移导数为$\overrightarrow{V}\cdot \overrightarrow{\nabla}$

有：

$\frac{D}{Dt}=\frac{\partial}{\partial t}+\overrightarrow{V}\cdot \overrightarrow{\nabla}$

当地导数的物理意义：流体微团在空间静止（不运动）时所对应的物理量的时间变化率，在数学上就表示为该物理量对时间$t$的偏导数；
迁移导数的物理意义：流体微团由于处于运动的状态（流畅空间的不均匀性）所导致的对应物理量的时间变化率。

2. 速度散度

场论里面，我们学过散度$div$或直接使用$\overrightarrow \nabla$算子点乘，而速度散度即

$\overrightarrow{\nabla}\cdot\overrightarrow{V}=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}$

考虑使用随流体运动的流体微团模型的性质（质量恒定），导出如下关系式：

$\overrightarrow{\nabla}\cdot\overrightarrow{V}=\frac{1}{\delta\mathscr{V}}\frac{D(\delta\mathscr{V})}{Dt}$

此式为速度散度赋予了物理意义：单位体积运动着的流体微团，体积变化的时间变化率。

质量守恒方程（连续性方程）

模型1：空间位置固定的有限控制体模型

核心原理：流体流入（出）控制体的净质量流量 = 流体质量的增加（减少）率
方程（积分形式/守恒形式）：

$\frac{\partial}{\partial t}\iiint_{\mathscr{V}}\rho d\mathscr{V}+\iint_S\rho \overrightarrow{V}\cdot d\overrightarrow{S}=0$

模型2：随流体流动的有限控制体模型

核心原理：控制体质量恒定
方程（积分形式/非守恒形式）：

$\frac{D}{Dt}\iiint_{\mathscr{V}}\rho d\mathscr{V}=0$

模型3：空间位置固定的无穷小微团模型

核心原理：流体流入（出）流体微团的净质量流量 = 流体质量的增加（减少）率
方程（微分形式/守恒形式）：

$\frac{\partial \rho}{\partial t}+\overrightarrow{\nabla}\cdot (\rho \overrightarrow{V})=0$

模型4：随流体流动的无穷小微团模型

核心原理：控制体质量恒定
推导过程中使用了速度散度物理意义的公式
方程（微分形式/非守恒形式）：

$\frac{D\rho}{Dt}+\rho\overrightarrow{\nabla}\cdot\overrightarrow{V}=0$

四种模型的转换

首先只需要掌握两个公式（场论课都学过）：
1. 散度定理：
  $\begin{aligned}\iint_{S}\overrightarrow{V}\cdot d\overrightarrow{S}=\iiint_{\Omega}\left(\overrightarrow{\nabla}\cdot \overrightarrow{V}\right)d{\Omega}\end{aligned}$
2. 向量恒等式：
  $\overrightarrow{\nabla}\cdot\left(\rho \overrightarrow{V}\right)=\rho \overrightarrow{\nabla}\cdot \overrightarrow{V}+\overrightarrow{V}\cdot \overrightarrow{\nabla}\rho$

对于积分守恒形式方程，控制体空间位置固定，有

$\begin{aligned}&\frac{\partial}{\partial t}\iiint_{\mathscr{V}}\rho d\mathscr{V}+\iint_S\rho \overrightarrow{V}\cdot d\overrightarrow{S}\\&=\iiint_{\mathscr{V}}\frac{\partial\rho}{\partial t} d\mathscr{V}+\iiint_\mathscr{V}\overrightarrow{\nabla}\cdot(\rho\overrightarrow{V}) d{\mathscr{V}}\\&=\iiint_\mathscr{V}\left(\frac{\partial\rho}{\partial t}+\overrightarrow{\nabla}\cdot(\rho\overrightarrow{V})\right)d\mathscr{V}\\&=0\end{aligned}$

由于控制体是任取的，所以被积函数恒为零，导出微分守恒形式方程：

$\frac{\partial \rho}{\partial t}+\overrightarrow{\nabla}\cdot (\rho \overrightarrow{V})=0$

对于微分守恒形式方程，可通过如下方式导出微分非守恒形式方程：

$\begin{aligned} &\frac{\partial \rho}{\partial t}+\overrightarrow{\nabla}\cdot (\rho \overrightarrow{V})\\&=\frac{\partial \rho}{\partial t}+\rho\overrightarrow{\nabla}\cdot \overrightarrow{V}+\overrightarrow{V}\cdot\overrightarrow{\nabla}\rho\\ &=\left(\frac{\partial \rho}{\partial t}+\overrightarrow{V}\cdot\overrightarrow{\nabla}\rho\right)+\rho\overrightarrow{\nabla}\cdot \overrightarrow{V}\\ &=\frac{D\rho}{Dt}+\rho\overrightarrow{\nabla}\cdot \overrightarrow{V} \end{aligned}$

对于积分非守恒形式方程，可通过如下方法导出积分守恒形式方程

$\begin{aligned}\frac{D}{Dt}\iiint_{\mathscr{V}}\rho d\mathscr{V}&=\iiint_{\mathscr{V}}\frac{D\left(\rho d\mathscr{V}\right)}{Dt}\\&=\iiint_{\mathscr{V}}\left[\rho\frac{D(d\mathscr{V})}{Dt}+d\mathscr{V}\frac{D\rho}{Dt}\right]\\&=\iiint_{\mathscr{V}}\rho\frac{D(d\mathscr{V})}{Dt}\frac{1}{d\mathscr{V}}d\mathscr{V}+\iiint_{\mathscr{V}}\frac{D\rho}{Dt}d\mathscr{V}\\&=\iiint_{\mathscr{V}}\left(\rho\overrightarrow{\nabla}\cdot \overrightarrow{V}\right)d\mathscr{V}+\iiint_{\mathscr{V}}\left(\frac{\partial \rho}{\partial t}+\overrightarrow{V}\cdot \overrightarrow{\nabla}\rho\right)d\mathscr{V}\\&=\iiint_{\mathscr{V}}\left(\rho\overrightarrow{\nabla}\cdot \overrightarrow{V}+\overrightarrow{V}\cdot \overrightarrow{\nabla}\rho\right)d\mathscr{V}+\iiint_{\mathscr{V}}\frac{\partial \rho}{\partial t}d\mathscr{V}\\&=\iiint_{\mathscr{V}}\frac{\partial \rho}{\partial t}d\mathscr{V}+\iiint_\mathscr{V}\overrightarrow{\nabla}\cdot(\rho\overrightarrow{V}) d{\mathscr{V}}\\&=\frac{\partial}{\partial t}\iiint_{\mathscr{V}}d{\mathscr{V}}+\iint_{S}\rho \overrightarrow{V}\cdot d\overrightarrow{S}\end{aligned}$

至此，所有四种模型的方程（积分/微分/守恒/非守恒形式）可以互相转换，所以他们本质上就是一个方程：质量守恒方程

动量方程

模型

采用模型4：随流体运动的无穷小微团，利用牛顿第二定律：
$\overrightarrow{F}=m\overrightarrow{a}$
此处略去推导过程，只记录公式

动量方程式

非守恒形式 $x、y、z$方向：

$\begin{aligned} \rho\frac{Du}{Dt}&=-\frac{\partial p}{\partial x}+\frac{\partial \tau_{xx}}{\partial x}+\frac{\partial \tau_{yx}}{\partial y}+\frac{\partial \tau_{zx}}{\partial z}+\rho f_x\\ \rho\frac{Dv}{Dt}&=-\frac{\partial p}{\partial y}+\frac{\partial \tau_{xy}}{\partial x}+\frac{\partial \tau_{yy}}{\partial y}+\frac{\partial \tau_{zy}}{\partial z}+\rho f_y\\ \rho\frac{Dw}{Dt}&=-\frac{\partial p}{\partial z}+\frac{\partial \tau_{xz}}{\partial x}+\frac{\partial \tau_{yz}}{\partial y}+\frac{\partial \tau_{zz}}{\partial z}+\rho f_z \end{aligned}$

守恒形式 $x、y、z$方向：
$\begin{aligned} \frac{\partial(\rho u)}{\partial t}+\overrightarrow{\nabla}\cdot\left(\rho u \overrightarrow{V}\right)&=-\frac{\partial p}{\partial x}+\frac{\partial \tau_{xx}}{\partial x}+\frac{\partial \tau_{yx}}{\partial y}+\frac{\partial \tau_{zx}}{\partial z}+\rho f_x\\ \frac{\partial(\rho v)}{\partial t}+\overrightarrow{\nabla}\cdot\left(\rho v \overrightarrow{V}\right)&=-\frac{\partial p}{\partial y}+\frac{\partial \tau_{xy}}{\partial x}+\frac{\partial \tau_{yy}}{\partial y}+\frac{\partial \tau_{zy}}{\partial z}+\rho f_y\\ \frac{\partial(\rho w)}{\partial t}+\overrightarrow{\nabla}\cdot\left(\rho w \overrightarrow{V}\right)&=-\frac{\partial p}{\partial z}+\frac{\partial \tau_{xz}}{\partial x}+\frac{\partial \tau_{yz}}{\partial y}+\frac{\partial \tau_{zz}}{\partial z}+\rho f_z \end{aligned}$
注：关注到等号右边的几项，分别代表了流体微团所受的力，它包含：
- 表面力：压力 + 粘性力（正应力+切应力）
- 体积力

能量守恒方程

模型

采用模型4：随流体运动的无穷小微团，利用热力学第一定律：

微团总能量的变化率 = 流入的净热流率 + 体积力和表面力对微团做功的功率
注：此处总能量包括了分子和原子的所有能量（内能+动能），前者对于原子来说有平动能和电子能，对于分子（双原子及以上）来说又多出了转动能和振动能
此处略去推导过程，只记录公式

能量守恒方程式

非守恒形式：
$\begin{aligned} \rho \frac{D}{Dt}\left(e+\frac{V^2}{2}\right)&=\rho\dot{q}+\frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial x}\right)+\frac{\partial}{\partial y}\left(k\frac{\partial T}{\partial y}\right)+\frac{\partial}{\partial z}\left(k\frac{\partial T}{\partial z}\right)\\ &-\frac{\partial (up)}{\partial x}-\frac{\partial (vp)}{\partial y}-\frac{\partial (wp)}{\partial z}\\ &+\frac{\partial (u\tau_{xx})}{\partial x}+\frac{\partial (u\tau_{yx})}{\partial y}+\frac{\partial (u\tau_{zx})}{\partial z}\\ &+\frac{\partial (v\tau_{xy})}{\partial x}+\frac{\partial (v\tau_{yy})}{\partial y}+\frac{\partial (v\tau_{zy})}{\partial z}\\ &+\frac{\partial (w\tau_{xz})}{\partial x}+\frac{\partial (w\tau_{yz})}{\partial y}+\frac{\partial (w\tau_{zz})}{\partial z}\\ &+\rho\overrightarrow{f}\cdot\overrightarrow{V} \end{aligned}$
守恒形式：
$\begin{aligned} &\frac{\partial}{\partial t}\left[\rho\left(e+\frac{V^2}{2}\right)\right]+\overrightarrow{\nabla}\cdot\left[\rho\left(e+\frac{V^2}{2}\right)\overrightarrow{V}\right] \\ &=\rho\dot{q}+\frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial x}\right)+\frac{\partial}{\partial y}\left(k\frac{\partial T}{\partial y}\right)+\frac{\partial}{\partial z}\left(k\frac{\partial T}{\partial z}\right)\\ &-\frac{\partial (up)}{\partial x}-\frac{\partial (vp)}{\partial y}-\frac{\partial (wp)}{\partial z}\\ &+\frac{\partial (u\tau_{xx})}{\partial x}+\frac{\partial (u\tau_{yx})}{\partial y}+\frac{\partial (u\tau_{zx})}{\partial z}\\ &+\frac{\partial (v\tau_{xy})}{\partial x}+\frac{\partial (v\tau_{yy})}{\partial y}+\frac{\partial (v\tau_{zy})}{\partial z}\\ &+\frac{\partial (w\tau_{xz})}{\partial x}+\frac{\partial (w\tau_{yz})}{\partial y}+\frac{\partial (w\tau_{zz})}{\partial z}\\ &+\rho\overrightarrow{f}\cdot\overrightarrow{V} \end{aligned}$
可以发现，能量方程和动量方程（守恒和非守恒形式）的等号左边各式的形式是一致的，而从非守恒形式转化到守恒形式时之后，仅仅是等号左边有变化而已。其转化的核心步骤如下：
$\begin{aligned} \rho \frac{Du}{Dt}&=\rho \frac{\partial u}{\partial t}+\rho\overrightarrow{V}\cdot \overrightarrow{\nabla}u\\ &=\frac{\partial (up)}{\partial t}-u\frac{\partial \rho}{\partial t}+\overrightarrow{\nabla}(\rho u\overrightarrow{V})-u\overrightarrow{\nabla}\cdot(\rho\overrightarrow{V})\\ &=\frac{\partial(\rho u)}{\partial t}+\overrightarrow{\nabla}(\rho u \overrightarrow{V})-\left[u\frac{\partial \rho}{\partial t}+u\overrightarrow{\nabla}\cdot(\rho\overrightarrow{V})\right] \end{aligned}$
注意到中括号内为空间位置固定的流体微团模型的质量守恒方程，故有：
$\rho \frac{Du}{Dt}=\frac{\partial(\rho u)}{\partial t}+\overrightarrow{\nabla}(\rho u \overrightarrow{V})$
同样的，能量方程的守恒和非守恒式之间的关系也可由此导出。

无粘流的欧拉方程（Euler Equation）

要点

忽略了耗散、粘性输运、质量扩散以及热传导的流动

方程式

连续性方程

非守恒式：
$\frac{D\rho}{Dt}+\rho\overrightarrow{\nabla}\cdot\overrightarrow{V}=0$
守恒式：
$\frac{\partial \rho}{\partial t}+\overrightarrow{\nabla}\cdot (\rho \overrightarrow{V})=0$
动量方程

非守恒式：
$\rho\frac{Du/v/w}{Dt}=-\frac{\partial p}{\partial x/y/z}+\rho f_{x/y/z}$
守恒式：
$\frac{\partial(\rho u/v/w)}{\partial t}+\overrightarrow{\nabla}\cdot\left(\rho u/v/w \overrightarrow{V}\right)=-\frac{\partial p}{\partial x/y/z}+\rho f_{x/y/z}\\$
能量方程

非守恒式：
$\begin{aligned} \rho \frac{D}{Dt}\left(e+\frac{V^2}{2}\right)=\rho\dot{q} -\frac{\partial (up)}{\partial x}-\frac{\partial (vp)}{\partial y}-\frac{\partial (wp)}{\partial z}+\rho\overrightarrow{f}\cdot\overrightarrow{V} \end{aligned}$
守恒形式：
$\frac{\partial}{\partial t}\left[\rho\left(e+\frac{V^2}{2}\right)\right]+\overrightarrow{\nabla}\cdot\left[\rho\left(e+\frac{V^2}{2}\right)\overrightarrow{V}\right] \\=\rho\dot{q} -\frac{\partial (up)}{\partial x}-\frac{\partial (vp)}{\partial y}-\frac{\partial (wp)}{\partial z}+\rho\overrightarrow{f}\cdot\overrightarrow{V}$

附加方程

我们发现，无论是$NS$方程还是$Euler$方程，都只有5个方程式（质量守恒1个、动量守恒3个、能量守恒1个），而其中却有6个未知数（$p、u、v、w、\rho、e$）
因此，我们需要补充方程，在空气动力学中，假设气体是完全气体（理想气体），则有状态方程：
$p=\rho RT$
但此方程又引入了新变量$T$，因此还需要补充状态参量间的热力学关系（量热状态方程），比如对定比热容的完全气体，有：
$e=c_vT$
其中$c_v$是比定容热容。

边界条件

对于任意特定的流动情况，其控制方程是相同的，而决定他们之间不同的因素便是边界条件（包括初始条件）

适合粘性流动的边界条件

无滑移条件

在物面处，有：
$u=v=w=0$
若壁面温度已知，则有：
$T=T_w$
若未知，由傅里叶导热定律，有：
$\dot{q_w}=-\left(k\frac{\partial T}{\partial n}\right)_w$
若为绝热壁面，有：
$\left(\frac{\partial T}{\partial n}\right)_w=0$
等等……

适合无粘流动的边界条件

无滑移条件

在物面处，有：
$\overrightarrow{V}\cdot\overrightarrow{n}=0$
等等……

适合CFD的控制方程

那必须是守恒式，其实理论上控制方程的守恒式与非守恒式是无差别的，但由于CFD的出现，导致这两者出现了区分，非守恒式所求出的解往往不尽如人意（会出现震荡or跑偏），而守恒式得出的解则较为合理。
观察$NS$方程，发现五个方程有共通的地方，可以归结为如下形式：
$\frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}+\frac{\partial G}{\partial y}+\frac{\partial H}{\partial z}=J$
而其中的$U、F、G、H、J$可以从$NS$方程中得出并以列向量的形式给出
通常情况采用时间推进算法解方程（后面会学到）
通常我们把$U$命名为解向量，因为该项的分量通常就是每一时间步中直接被求解的未知函数，而$F G H$为通量项，当然我们也可以自己指定解向量

弱守恒&强守恒

强守恒：所有东西都写进了导数里面，没有任何流动变量单独留在导数之外（如上式）
弱守恒：出现了有流动变量单独留在导数之外的情形

激波装配法&激波捕捉法

激波捕捉法：将激波作为CFD方程的计算结果，直接计算出激波的流场

优点：不需要事先知道激波位置，以及激波的边界条件

缺点：可能会与真实情况相差甚远，误差起伏大
激波装配法：利用兰金-许贡钮关系式，将激波人为引入到流场解中，而流动控制方程只用来解除激波之外的解

优点：与真实情况接近，激波位置确定

缺点：需要事先知道激波位置，并利用斜激波关系式引入

心得：orz这一章的东西真TM硬核啊，看了一整天，算是把大二下的工流没弄明白的东西全给整明白了，~~我吐了~~。

物理

流体力学 CFD

本博客所有文章除特别声明外，均采用 CC BY-SA 4.0 协议，转载请注明出处！

快速幂&矩阵快速幂上一篇

链式前向星解释下一篇