测试----

李东岳

不只是N-S方程

谈到CFD，大部分人最开始想到的就是Navier-Stokes方程（N-S方程）。N-S方程从守恒定律推导而来：
\begin{equation}\tag{1}
\frac{\p \rho}{\p t}+\nabla\cdot(\rho\bfU)=0,
\label{C}
\end{equation}
\begin{equation}\tag{2}
\frac{\partial \rho\mathbf{U}}{\partial t}+\nabla \cdot (\rho\mathbf{U}\mathbf{U})=-\nabla p+\nabla \cdot\tau,
\label{M}
\end{equation}
其中$\rho$为密度，$\bfU$为速度，$p$为压力，$\tau$为剪切应力。N-S方程具有以下特点：

方程\eqref{M}中左边第二项是关于$\bfU$的偏导数，这种未知量和未知量乘积的问题构成非线性问题，CFD对非线性问题需要特殊处理。另一方面，非线性的双曲问题的解可能会存在间断（如激波）。激波通常存在于高超声速的欧拉问题求解中。同时，非线性项也是湍流在数学方程中的体现；

方程\eqref{M}的数学特征为抛物线。不同数学特征的问题需要调用不同的时间/空间离散格式，隐性时间格式更有利于求解抛物线问题。若方程\eqref{M}中省略若干项则会改变方程的数学特征，例如若将方程右侧置为$0$，则变为双曲特征的欧拉方程。欧拉方程得益于其双曲特性，可采用迎风类显性算法推进，各种基于有限体积法的高分辨率格式因此而生（交错网格中心格式、中心-迎风格式等）。同时，动量方程的对流项和扩散项的相对强弱，也会影响边界条件的设置（如DNS的出口无反射边界条件）；

在马赫数较大时（如大于0.3），方程\eqref{C}可用来求解密度，方程\eqref{M}可用来求解速度，同时附加能量方程求解温度以及状态方程求解压力，即密度基求解器。在马赫数较小时，并没有单独的压力方程，并且方程\eqref{C}缺少主要求解变量。这导致压力的求解需要特殊的策略。这也是CFD中压力基SIMPLE/PISO算法要处理的问题。在构建离散矩阵的情况下，对速度和压力整合处理还是单独处理，是分离求解与耦合算法的重要问题；

N-S方程之前传：N-S方程为宏观方程，调用了宏观假定，其可从玻尔兹曼方程附加Chapman-Enskog展开推导而来。N-S方程也可以看做从介尺度模型演化的宏观二阶矩模型。在无压力无粘性的条件下具备弱双曲特征。由于失去了高阶矩的统计学特征，因此N-S方程在某些情况下是不适用的，比如稀薄流、微流动、以及气固多相流中的颗粒轨迹交叉现象都不能用N-S方程来描述。其数学本源在于N-S方程假定流体为一种近似均衡状态，在某一个网格点存在一个单值的分布。然而，在分子/颗粒数量非常小的情况下，这些分子/颗粒基本无碰撞，会导致多值问题的产生。因此，在这种情况下，需要对统计学分布的高阶矩进行追踪，N-S方程作为二阶矩模型是往往不够的。

正如上面讨论的，N-S方程中存在大量的数学问题，初学者也看的不明不白。这无关紧要。本章仅仅抛出一块砖头，让大家了解隐藏在CFD中的方程的美妙含义。这些方程中充满火星符号（如$\nabla$，$\p$），也让人摸不清头脑。\href{https://www.cfd-china.com/topic/2397}{东岳流体CFD课程}因此分为两步，第一步是基本的N-S方程入门，第二步是N-S方程求解。第一步需要同学们通过本笔记进行预习。第二步将在课程上讲授。

学习CFD理论首先遇到的就是各种各样的偏微分算符。能看懂、拆分CFD方程是研究算法的最基本步骤。给出一个CFD方程（如动量方程），可以不知道是怎么推出来的，但要能看懂。本章介绍CFD方程的各种写法，不涉及到任何的CFD算法。

矢量标识法

本笔记中主要采用矢量标识法讨论CFD方程，在矢量标识法中，标量全部采用斜体，如压力$p$。矢量采用正体加粗，如速度矢量$\bfU$，其具有三个分量$u_1$，$u_2$，$u_3$或$u$，$v$，$w$。二阶张量也采用正体\footnote{如果难以理解二阶张量的含义，可以这样尝试：矢量是一阶张量，具有三个分量，二阶张量则具有9个分量。}，比如应力张量$\boldsymbol{\tau}$，其具备9个分量，其可表示为：
\begin{equation}
\boldsymbol{\tau}=
\begin{bmatrix}
\tau_{11} & \tau_{12} & \tau_{13} \
\tau_{21} & \tau_{22} & \tau_{23} \
\tau_{31} & \tau_{32} & \tau_{33}
\end{bmatrix}
\end{equation}
下面粘度$\nu$为$1$的不可压缩流体动量方程为例对其进行拆分，这个方程若采用矢量标识法可以写为：
\begin{equation}
\frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot (\mathbf{U}\mathbf{U})=-\nabla \frac{p}{\rho}+\nabla \cdot(\nabla \mathbf{U})
\label{mom}
\end{equation}
其中的$\mathbf{U}$为速度矢量，$p$为压力，$\rho$为密度。一般来讲，CFD文献中通常采用方程\eqref{mom}的形式，而并不进行展开进而更加紧凑。下面介绍如何将其展开为3个方程。
\begin{itemize}
\item 方程\eqref{mom}第一项表示$\bfU$对时间的偏导数，因为$\mathbf{U}$为矢量，故其导数的分量形式为：
\begin{equation}\label{time}
\frac{\partial \mathbf{U}}{\partial t}=
\left[
\begin{matrix}
\frac{\partial u_1}{\partial t} \
\frac{\partial u_2}{\partial t} \
\frac{\partial u_3}{\partial t} \
\end{matrix}
\right]
\end{equation}
其中$u_1$表示$x$方向速度，$u_2$表示$y$方向速度，$u_3$表示$z$方向速度。这样拆分之后的方程，即为各个方向的速度针对时间的偏导数。如果理解方程\eqref{time}有困难，那么有必要预习一下《高等数学》（同济大学版）第二章。
\item 方程\eqref{mom}第二项$\nabla \cdot (\mathbf{U}\mathbf{U})$中的$\mathbf{U}\mathbf{U}$是一种简写，完整形式为$\mathbf{U}\otimes \mathbf{U}$，$\otimes$是一个张量运算符。依据$\otimes$的定义，$\mathbf{U}\mathbf{U}$可以写为：
\begin{equation}
\mathbf{U} \otimes \mathbf{U}=\mathbf{U}\mathbf{U}=\mathbf{U}\cdot\mathbf{U}^T
=\left[\begin{matrix}
u_1\
u_2\
u_3
\end{matrix}\right][u_1, u_2, u_3]=\left[
\begin{matrix}
u_1 u_1 & u_1 u_2 & u_1 u_3\
u_2 u_1 & u_2 u_2 & u_2 u_3\
u_3 u_1 & u_3 u_2 & u_3 u_3
\end{matrix}
\right]
\end{equation}
接下来看 $\nabla \cdot$，其为散度算符，有时用$\mathrm{div}$来表示。对一个矢量（$1$阶张量）做散度的结果为一个标量（$0$阶张量），对一个$2$阶张量做散度的结果为矢量（$1$阶张量）。因此，对任意$n$阶张量做散度操作之后，结果为$n-1$阶张量。举例，对一个矢量$\bfU$做散度有：$\nabla \cdot \mathbf{U} = \frac{\partial u_1}{\partial x}+\frac{\partial u_2}{\partial y}+\frac{\partial u_3}{\partial z}$。因此，方程\eqref{mom}中的第二项$\nabla \cdot (\mathbf{U}\mathbf{U})$即为对一个$2$阶张量做散度：
\begin{equation}\label{conv}
\nabla \cdot \left(\mathbf{U}\mathbf{U}\right) = \nabla \cdot \left[
\begin{matrix}
u_1 u_1 & u_1 u_2 & u_1 u_3\
u_2 u_1 & u_2 u_2 & u_2 u_3\
u_3 u_1 & u_3 u_2 & u_3 u_3
\end{matrix}
\right]=\left[
\begin{matrix}
\frac{\partial u_1 u_1}{\partial x}+\frac{\partial u_2 u_1}{\partial y}+\frac{\partial u_3 u_1}{\partial z} \
\frac{\partial u_1 u_2}{\partial x}+\frac{\partial u_2 u_2}{\partial y}+\frac{\partial u_3 u_2}{\partial z} \
\frac{\partial u_1 u_3}{\partial x}+\frac{\partial u_2 u_3}{\partial y}+\frac{\partial u_3 u_3}{\partial z}
\end{matrix}
\right]
\end{equation}

warnerchang

@李东岳东岳老师，方程双曲性质具体是什么？不是特别懂，有时是双曲，有时又是抛物线，椭圆的，他们之间是什么关系呢？

其次，为啥式（2）中的运动粘度为1呢？式（2）是openfoam中求解方程的形式吗？

李东岳

@warnerchang 方程2是简写，省略粘度符号，要不看起来太麻烦

warnerchang

@李东岳噢噢，好的