在对CFD即计算流体力学的学习过程中，Navier-Stokes方程是绕不过去的。在较常见的物理场中，Navier-Stokes 方程还是适用的。下面主要整理记录一下Navier-Stokes方程的基本推导过程。
\[\frac{\partial \rho}{\partial t}+\nabla \cdot(\rho \mathbf{U})=0\]
\begin{equation}\frac{\partial \rho \mathbf{U}}{\partial t}+\nabla \cdot(\rho \mathbf{U U})=-\nabla p+\nabla \cdot \tau\end{equation}
基础知识
在对Navier-Stokes方程进行推导之前，还是有必要明确相关数学概念，毕竟对于一个八百年没用过高数的人来说，单靠回忆怕是不怎么靠谱。
泰勒展开
泰勒展开式应该在大一高数有讲过，简单来说就是一个用函数在某点的信息描述其附近取值的公式。既可用于一元函数也可用于多元函数，下列公式摘自维基百科，是带有拉格朗日型余项的泰勒公式：
\begin{equation}
f(x)=f(a)+\frac{f^{\prime}(a)}{1 !}(x-a)+\frac{f^{(2)}(a)}{2 !}(x-a)^{2}+\cdots+\frac{f^{(n)}(a)}{n !}(x-a)^{n}+\frac{f^{(n+1)}(\theta)}{(n+1) !}(x-a)^{(n+1)}
\end{equation}
一般来说，在下述Navier-Stokes方程的推导过程中，只用到了一阶表达形式。至于为什么不用到高阶的形式，个人怀疑是因为针对有限体积法来说，一阶形式的精度可能就足够了？下面分别是一元及多元函数的泰勒展开式：
一元泰勒展开式：
\begin{equation}f\left(x_{1}\right)=f\left(x_{0}\right)+\left.\frac{\partial f(x)}{\partial x}\right|_{x=x_{0}}\left(x_{1}-x_{0}\right)+\ldots\end{equation}
多元泰勒展开式：
\begin{equation}\begin{array}{r}f(x, y, z, t)=f\left(x_{1}, y_{1}, z_{1}, t_{1}\right)+\left.\frac{\partial f(x, y, z, t)}{\partial x}\right|_{x=x_{1}}\left(x-x_{1}\right)+\left.\frac{\partial f(x, y, z, t)}{\partial y}\right|_{y=y_{1}}\left(y-y_{1}\right)+ \\\left.\frac{\partial f(x, y, z, t)}{\partial z}\right|_{z=z_{1}}\left(z-z_{1}\right)+\left.\frac{\partial f(x, y, z, t)}{\partial t}\right|_{t=t_{1}}\left(t-t_{1}\right)+\ldots\end{array}\end{equation}
数学符号
显然，在文章开头所展示的Navier-Stokes方程有一点晦涩，主要是因为它的表达形式上引入了散度、梯度等内容，如果将其替换为常见的微分形式，相信起码能让人有一点亲近感。
Nabla算子
首先引入Nabla算子，也称为倒三角算子，这可能是因为它的写法。其形式为：
\begin{equation}\nabla=\frac{d}{d r}\end{equation}
梯度
可直接使用Nabla算子表达梯度，当其直接作用于函数时，无论其是否为矢量或是标量，都表示其梯度。标量函数的梯度为向量，向量的梯度为二阶张量。
\begin{equation}\nabla f=\left(\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y}, \frac{\partial f}{\partial z}\right)=\frac{\partial f}{\partial x} \mathbf{i}+\frac{\partial f}{\partial y} \mathbf{j}+\frac{\partial f}{\partial z} \mathbf{k}\end{equation}
散度
散度是向量场的一种强度性质，就如同密度、浓度、温度一样，它对应的广延性质是一个封闭区域表面的通量，所以说散度是通量的体密度。流体力学中，散度为零的流体称为不可压缩流体，也就是说此流体中不会有一部分凭空消失或突然产生，每个微小时间间隔中流入一个微小体元的流体总量都等于在此时间间隔内流出此体元的流体总量。
$$\nabla \cdot \mathbf{A}=\frac{\partial A_{x}}{\partial x}+\frac{\partial A_{y}}{\partial y}+\frac{\partial A_{z}}{\partial z}$$
实际上，散度是一个标量。
旋度
旋度是向量场的一种强度性质，就如同密度、浓度、温度一样，它对应的广延性质是向量场沿一个闭合曲线的环量，所以说旋度是环量的面密度。如果一个向量场中处处的旋度都是零，则称这个场为无旋场或保守场。
不同于梯度和散度，旋度不能简单的推广到其他维度；某些推广是可能的，但是只有在三维中，在几何上定义的向量场旋度还是向量场。
\begin{equation}\boldsymbol{\nabla} \times \mathbf{A}=\left(\frac{\partial A_{z}}{\partial y}-\frac{\partial A_{y}}{\partial z}\right) \mathbf{i}+\left(\frac{\partial A_{x}}{\partial z}-\frac{\partial A_{z}}{\partial x}\right) \mathbf{j}+\left(\frac{\partial A_{y}}{\partial x}-\frac{\partial A_{x}}{\partial y}\right) \mathbf{k}\end{equation}
当然，也可表示为行列式。
\begin{equation}\left|\begin{array}{ccc}\mathbf{i} & \mathbf{j} & \mathbf{k} \\\frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\A_{x} & A_{y} & A_{z}\end{array}\right|\end{equation}
不过今天的推导过程中，倒没有涉及到旋度。
N-S方程
前言
流体与固体有一个明显的区别，就是在承受剪应力时将会发生连续变形。当流体受到剪应力变形或拉伸应力时，其会产生一定阻力，这样的阻力被称为粘滞力或者粘性力。我们用黏度（Viscosity）来表征流体的黏性，黏度也被称为运动粘度、粘滞系数等。当一流体黏度恒定时，称其为牛顿流体，当一流体黏度不恒定时，称其为非牛顿流体。
流体主要包括液体和气体，至于等离子体暂不做讨论。液体与气体相比，其可压缩性大大降低。通常来看，将液体视为不可压缩流体，其密度为常数。至于是否将气体视为可压缩流体，这个需要结合实际分析情况。
同时，在使用N-S方程的一个重要前提是连续介质假定（the concept of continuum）。此时，一个单元体内有足够的流体分子，在宏观上充分小，而在微观上充分大。
总体而言，为描述流体的运动，我们将使用以下三种方程。其一是连续性方程，满足质量守恒定律；其二是动量方程，满足牛顿第二定律；在考虑传热的情况下，还应附加能量方程，满足能量守恒定律。
多数情况下，CFD求解方法为有限体积法（FVM）。
流动模型
在讨论N-S方程前，我们要明确流动模型。
流动模型分为有限控制体模型和无穷小微团模型，同时，也可进一步分为空间位置固定和随流线运动，对应的描述为欧拉描述和拉格朗日描述。在本文中，为方便推导，使用的是空间位置固定的无穷小微团模型。此时，质量、体积的描述可使用dV、dm来表示。
并引入物质导数的概念。描述t1时刻笛卡尔坐标系下的一微团速度U1(x1,y1,z1,t1)，随微团运动， t2时刻笛卡尔坐标系下的一微团温度U2(x2,y2,z2,t2) 。基于多元函数的泰勒展开，显然有如下方程。
\begin{equation}\begin{aligned}U_{2} &=U_{1}+\left.\frac{\partial T}{\partial x}\right|_{x=x_{1}}\left(x_{2}-x_{1}\right)+\left.\frac{\partial T}{\partial y}\right|_{y=y_{1}}\left(y_{2}-y_{1}\right)+\left.\frac{\partial T}{\partial z}\right|_{z=z_{1}}\left(z_{2}-z_{1}\right)+\left.\frac{\partial T}{\partial t}\right|_{t=t_{1}}\left(t_{2}-t_{1}\right)\end{aligned}\end{equation}
两侧分别除以(t2– t1)：
\begin{equation}\frac{u_{2}-u_{1}}{t_{2}-t_{1}}=\left.\frac{\partial U}{\partial x}\right|_{x=x_{1}} \frac{x_{2}-x_{1}}{t_{2}-t_{1}}+\left.\frac{\partial U}{\partial y}\right|_{y=y_{1}} \frac{y_{2}-y_{1}}{t_{2}-t_{1}}+\left.\frac{\partial U}{\partial z}\right|_{z=z_{1}} \frac{z_{1}-z_{1}}{t_{x}-t_{1}}+\left.\frac{\partial U}{\partial t}\right|_{t=t_{1}} \end{equation}
当 t2无限接近 t1时，根据速度的定义即可得到速度U 在笛卡尔坐标系下的物质导数的定义 ，
\begin{equation}\frac{d U}{d t}=\frac{\partial U}{\partial x} \cdot u+\frac{\partial U}{\partial y} \cdot v+\frac{\partial U}{\partial t} \cdot w+\frac{\partial U}{\partial t}\end{equation}
其余性质，例如特征线等，留待下次记录。
\begin{equation}\underbrace{\frac{\mathrm{D}}{\mathrm{D} t}}_{\text {Lagrangian }}=\underbrace{\frac{\partial}{\partial t}+\mathrm{U} \cdot \nabla}_{\text {Euler }}\end{equation}
连续性方程
以较为简单的欧拉法对连续性方程进行推导，欧拉观点下的质量守恒意味着位置固定无穷小微团质量的变化= 流入无穷小微团质量−流出无穷小微团质量。也就是这两者变化率相同。
首先表达无穷小微团质量的变化率，显然流体微团体积为\(\mathrm{d} V=\mathrm{d} x \mathrm{~d} y \mathrm{~d} z\)，基于密度\(\rho\)，可得到其对应的无穷小微团质量：
\[\mathrm{d} m=\rho \mathrm{d} x \mathrm{~d} y \mathrm{~d} z\]
其变化率为：
\[\frac{\partial \mathrm{d} m}{\partial t}=\frac{\partial \rho \mathrm{d} x \mathrm{~d} y \mathrm{~d} z}{\partial t}=\frac{\partial \rho}{\partial t} \mathrm{~d} x \mathrm{~d} y \mathrm{~d} z\]
随后表达（流入无穷小微团质量−流出无穷小微团质量）的变化率，其中以\(x\)方向为例，此时单位时间内的流入质量为：
\[\rho u \mathrm{~d} y \mathrm{~d} z\]
根据泰勒公式的一阶展开，能得到单位时间内的流出质量为：
\[\left(\rho u+\frac{\partial \rho u}{\partial x} \mathrm{~d} x\right) \mathrm{d} y \mathrm{~d} z\]
因此，可得到（流入无穷小微团质量−流出无穷小微团质量）的变化率为：
\[\rho u \mathrm{~d} y \mathrm{~d} z-\left(\rho u+\frac{\partial \rho u}{\partial x} \mathrm{~d} x\right) \mathrm{d} y \mathrm{~d} z=-\frac{\partial \rho u}{\partial x} \mathrm{~d} x \mathrm{~d} y \mathrm{~d} z\]
同理， \(y\)方向上的质量变化率为：
\[ -\frac{\partial \rho v}{\partial y} \mathrm{~d} x \mathrm{~d} y \mathrm{~d} z\]
\(z\)方向上的质量变化率为：
\[ -\frac{\partial \rho w}{\partial z} \mathrm{~d} x \mathrm{~d} y \mathrm{~d} z\]
此时，联立上述方程，有
\[\frac{\partial \rho}{\partial t} \mathrm{~d} x \mathrm{~d} y \mathrm{~d} z=-\left(\frac{\partial \rho u}{\partial x}+\frac{\partial \rho v}{\partial y}+\frac{\partial \rho w}{\partial z}\right) \mathrm{d} x \mathrm{~d} y \mathrm{~d} z\]
也即为连续方程：
\[\frac{\partial \rho}{\partial t}+\frac{\partial \rho u}{\partial x}+\frac{\partial \rho v}{\partial y}+\frac{\partial \rho w}{\partial z}=0\]
基于散度，得到下式：
\begin{equation}\frac{\partial \rho}{\partial t}+\nabla \cdot(\rho \mathbf{U})=0\end{equation}
动量方程
动量方程的推导基于空间位置移动的无穷小微团，这是因为其比较简单。
首先对无穷小微团进行受力分析。其受力可分为体积力和表面力。体积力作用在无穷小微团的全部体积上，例如重力、电磁力；表面力作用在无穷小微团的面上，包括压力、表面张力。在当前分析推导中，最重要的表面力是压力和应力，不考虑其他源项力的情况下，重力是最重要的体积力。
压力
压力在这里指静压，用\(p\)来表示，其影响的是无穷小微团的压缩与膨胀，同时，压力作为正应力，其永远作用于无穷小微团面的法向。
以\(x\)方向为例，其左侧受压力为：
\[ p\mathrm{~d} y \mathrm{~d} z\]
其右侧受压力为：
\[ – p\mathrm{~d} y \mathrm{~d} z -\frac{\partial p}{\partial x} \mathrm{~d} x \mathrm{~d} y \mathrm{~d} z\]
则在\(x\)方向上，受压力为：
\[ -\frac{\partial p}{\partial x} \mathrm{~d} x \mathrm{~d} y \mathrm{~d} z\]
同理\(y\)方向上受压力为：
\[ -\frac{\partial p}{\partial y} \mathrm{~d} x \mathrm{~d} y \mathrm{~d} z\]
同理\(z\)方向上受压力为：
\[ -\frac{\partial p}{\partial z} \mathrm{~d} x \mathrm{~d} y \mathrm{~d} z\]
剪切应力
参考文章，暂时不打出具体公式及图解，目前根据剪切应力来看，可给出各个方向上的受力。
综上所述，无穷小微团受力为：
\[\mathbf{F}=\left[\begin{array}{l}
\left(-\frac{\partial p}{\partial x}+\frac{\partial \tau_{x x}}{\partial x}+\frac{\partial \tau_{y x}}{\partial y}+\frac{\partial \tau_{z x}}{\partial z}\right) \mathrm{d} x \mathrm{~d} y \mathrm{~d} z \
\left(-\frac{\partial p}{\partial y}+\frac{\partial \tau_{x y}}{\partial x}+\frac{\partial \tau_{y y}}{\partial y}+\frac{\partial \tau_{z y}}{\partial z}\right) \mathrm{d} x \mathrm{~d} y \mathrm{~d} z \
\left(-\frac{\partial p}{\partial z}+\frac{\partial \tau_{x z}}{\partial x}+\frac{\partial \tau_{y z}}{\partial y}+\frac{\partial \tau_{z z}}{\partial z}\right) \mathrm{d} x \mathrm{~d} y \mathrm{~d} z
\end{array}\right]\]
基于牛顿第二定律，已知无穷小流体微团质量为\(\mathrm{d} m=\rho \mathrm{d} x \mathrm{~d} y \mathrm{~d} z\)，流体加速度为\(\left[\frac{\mathrm{D} u}{\mathrm{D} t}, \frac{\mathrm{D} v}{\mathrm{D} t}, \frac{\mathrm{D} w}{\mathrm{D} t}\right]^{\mathrm{T}}\)，结合上述受力，可得到如下方程。
以\(x\)方向为例，在\(x\)方向上的动量方程为：
\[\rho \mathrm{d} x \mathrm{~d} y \mathrm{~d} z \frac{\mathrm{D} u}{\mathrm{D} t}=\left(-\frac{\partial p}{\partial x}+\frac{\partial \tau_{x x}}{\partial x}+\frac{\partial \tau_{y x}}{\partial y}+\frac{\partial \tau_{z x}}{\partial z}\right) \mathrm{d} x \mathrm{~d} y \mathrm{~d} z\]
此时有如下方程：
\[\rho \frac{\mathrm{D} u}{\mathrm{D} t}=-\frac{\partial p}{\partial x}+\frac{\partial \tau_{x x}}{\partial x}+\frac{\partial \tau_{y x}}{\partial y}+\frac{\partial \tau_{z x}}{\partial z}\]
\[ \rho \frac{\mathrm{D} v}{\mathrm{D} t}=-\frac{\partial p}{\partial y}+\frac{\partial \tau_{x y}}{\partial x}+\frac{\partial \tau_{y y}}{\partial y}+\frac{\partial \tau_{z y}}{\partial z}\]
\[ \rho \frac{\mathrm{D} w}{\mathrm{D} t}=-\frac{\partial p}{\partial z}+\frac{\partial \tau_{x z}}{\partial x}+\frac{\partial \tau_{y z}}{\partial y}+\frac{\partial \tau_{z z}}{\partial z}\]
结合散度、梯度等定义，方程可写为如下形式：
\[\rho \frac{\mathrm{DU}}{\mathrm{D} t}=-\nabla p+\nabla \cdot \tau\]
结合物质导数，方程可写为如下形式，这也是较为常见的表达：
\[\rho\left(\frac{\partial \mathbf{U}}{\partial t}+\mathbf{U} \cdot \nabla \mathbf{U}\right)=-\nabla p+\nabla \cdot \tau\]
非牛顿流体的引入
明日再写，嘿嘿