NS 方程索引

流体流动类型众多，可以按照流体粘性种类划分，也可以按照流体压缩性划分。本文参考了 COMSOL 帮助文档，列出主要常用的各种流体流动方程，以作参考。

标准 NS 方程

NS 方程用于描述速度场 $\mathit{u}$ 和压力场 $p$ 的关系，可以写作 $$\begin{array}{rl}\underset{\text{惯性力}}{\underbrace{\rho (\frac{\partial \mathit{u}}{\partial t}+(\mathit{u}\cdot \mathrm{\nabla })\mathit{u})}}& =\underset{\text{压力}}{\underbrace{-\mathrm{\nabla }p}}+\underset{\text{粘性力}}{\underbrace{\mathrm{\nabla }\cdot \mathit{K}}}+\underset{\text{体积力}}{\underbrace{\mathit{F}}},\\ \frac{\partial \rho }{\partial t}+\mathrm{\nabla }\cdot (\rho \mathit{u})& =0.\end{array}$$ 其中 $\rho$ 是流体密度，粘性力项中张量 $\mathit{K}$ 与动力粘度 $\mu$ 和速度 $\mathit{u}$ 相关 $$\mathit{K}=\mu (\mathrm{\nabla }\mathit{u}+(\mathrm{\nabla }\mathit{u}{)}^{\mathsf{T}})-\frac{2}{3}\mu ((\mathrm{\nabla }\cdot \mathit{u})\mathit{I}).$$ 其中 $\mathit{I}$ 是单位矩阵，$\mathrm{\nabla }\mathit{u}+(\mathrm{\nabla }\mathit{u}{)}^{\mathsf{T}}=\mathit{D}$ 也被称为变形率张量。如果流体不可压，则有 $\text{d}\rho /\text{d}t=0$，代入方程 $\text{(???)}$ 得 $\mathrm{\nabla }\cdot \mathit{u}=0$。假设流体粘度恒定，此时粘性力项可简化为 $$\begin{array}{rl}\mathrm{\nabla }\cdot \mathit{K}& =\mathrm{\nabla }\cdot [\mu (\mathrm{\nabla }\mathit{u}+(\mathrm{\nabla }\mathit{u}{)}^{\mathsf{T}})-\frac{2}{3}\mu ((\mathrm{\nabla }\cdot \mathit{u})\mathit{I})]\\ & =\mathrm{\nabla }\cdot [\mu (\mathrm{\nabla }\mathit{u}+(\mathrm{\nabla }\mathit{u}{)}^{\mathsf{T}})]\\ & =\mu \mathrm{\nabla }\cdot \mathrm{\nabla }\mathit{u}=\mu {\mathrm{\nabla }}^2\mathit{u}.\end{array}$$ 则不可压、粘性、层流 NS 方程可写作 $$\boxed{{\displaystyle \begin{array}{rl}\rho (\frac{\partial \mathit{u}}{\partial t}+(\mathit{u}\cdot \mathrm{\nabla })\mathit{u})& =-\mathrm{\nabla }p+\mu {\mathrm{\nabla }}^2\mathit{u}+\mathit{F},\\ \mathrm{\nabla }\cdot \mathit{u}& =0.\end{array}}}$$ 流体的可压缩性可以通过无量纲数——马赫数 $$Ma=U/c$$ 来判定，其中 $U$ 是流体特征速度， $c$ 是流体声速。按马赫数大小可分为

• 不可压流动（$Ma=0$）
• 弱可压缩流动（$Ma<0.3$）
• 亚声速流动（$Ma=0.3～0.8$ 左右）
• 跨声速流动（$Ma=0.8～1.2$ 左右）
• 超声速流动（$Ma=1.2～5.0$ 左右）
• 高超声速流动（$Ma>5.0$）

不可压流动可通过方程 $\text{(???)}$ 求解，弱可压缩流动可通过方程 $\text{(???)}$ 和 $\text{(???)}$ 求解。当 $Ma>0.3$ 需要考虑湍流效应和热效应，需要额外方程来描述流体流动特征。

对于不可压流动，描述湍流效应还有一重要的无量纲数——雷诺数 $$Re=\frac{\rho UL}{\mu }=\frac{\text{惯性力}}{\text{粘性力}},$$ 其中 $L$ 是流动结构特征长度。尽管不同的流动结构产生湍流效应的临界雷诺数不同，一般可认为当 $Re>2000$ 时产生湍流效应，当 $Re<2000$ 流体流动可通过方程 $\text{(???)}$ 求解。如果惯性力相对于粘性力可以忽略（一般 $Re<<1$），流体流动及其缓慢，该种流动也被称为为蠕动流。在求解时，方程 $\text{(???)}$ 中惯性项中对流项 $(\mathit{u}\cdot \mathrm{\nabla })\mathit{u}$ 可以忽略 $$\boxed{{\displaystyle \begin{array}{rl}\rho \frac{\partial \mathit{u}}{\partial t}& =-\mathrm{\nabla }p+\mu {\mathrm{\nabla }}^2\mathit{u}+\mathit{F},\\ \mathrm{\nabla }\cdot \mathit{u}& =0.\end{array}}}$$

湍流

热效应

两相流

在工程上，有时需要追踪流体的界面，此时，表面张力在界面变形中起到重要作用。如最常见的水龙头接水时，液柱界面形状变化。考虑最普通的气液两相不可压流动，如下图所示

若空间内气液界面 $\mathit{s}=\mathit{s}(\mathit{\xi };t)$ 用两个参数 $\mathit{\xi }=({\xi }_1,{\xi }_2)$ 的参数方程表示，则界面的运动方程满足 $$(\frac{\partial \mathit{s}(\mathit{\xi };t)}{\partial t}-\mathit{u}(\mathit{\xi };t))\cdot \mathit{n}=0,$$ 其中 $\mathit{u}(\mathit{\xi };t)$ 表示位置 $\mathit{\xi }=({\xi }_1,{\xi }_2)$ 的速度，该方程通常也简写为动力学边界条件 $\mathit{n}\cdot (\mathit{u}-{\mathit{u}}_S)=0$，${\mathit{u}}_S$ 是气液界面运动速度。如果用等值面 $S(\mathit{s}(\mathit{\xi };t);t)=S_0$ 来表示气液界面，上式可写作 $$\begin{array}{r}\frac{\partial S}{\partial t}+\mathit{u}\cdot \mathrm{\nabla }S=0.\end{array}$$ 在气相和液相 ${\mathrm{\Omega }}_1\cup {\mathrm{\Omega }}_2$ 中，流体流动均满足方程 $\text{(???)}$ 。在气液界面 $S$ 上，两侧的应力差驱动界面流动，满足界面应力平衡方程 $$\begin{array}{r}\boxed{{\displaystyle \mathit{n}\cdot [{\mathit{T}}_i{]}_2^1=\gamma \kappa \mathit{n},}}\end{array}$$ 其中 $[x_i{]}_2^1=x_1-x_2$ 表示物理量跨越界面时的差值，${\mathit{T}}_i=-p_i\mathit{I}\mathbf{+}{\mu }_i{\mathit{D}}_i$ 称为水动力应力张量（hydrodynamic stress tensor），$\kappa =-\mathrm{\nabla }\cdot \mathit{n}$ 是界面平均曲率。如果气液两相可当作无粘流体或者均处于静止状态，则两侧压力差为 $$\mathrm{\Delta }p=p_1-p_2=\gamma \kappa ,$$ 该压力差也被称为拉普拉斯压力。注意NS 方程是体方程，而方程 $\text{(???)}$ 是表面方程。追踪界面方法有流体体积（Volume of fluid， VOF）法、水平集（Level set，LS）法和任意拉格朗日（Arbitrary Lagrangian Eulerian，ALE）法等。其中 VOF 和 LS 法将界面处理成一定厚度的过渡层。

VOF 方法

VOF 法引入相体积分数 $c$ 来实现对计算域内界面进行追踪，同时将双流体方程（气液相各一个）转化为单流体方程。流体体积分数一般设置为 $$\{\begin{array}{rlr}& c=0& \text{气相}\\ & 0<c<1& \text{界面}\\ & c=1& \text{液相}\end{array}.$$

通过求出整个计算域内各网格单元的相分数，从而可以构建出界面

在 VOF 法中，所有物理量经过界面都是渐变的，如转化为单流体方程中密度和粘度可表示为 $$\rho (c)=(1-c){\rho }_1+c{\rho }_2,\mu (c)=(1-c){\mu }_1+c{\mu }_2.$$ 此时求解的单流体动量方程为 $$\begin{array}{r}\rho (c)(\frac{\partial \mathit{u}}{\partial t}+(\mathit{u}\cdot \mathrm{\nabla })\mathit{u})=\mathrm{\nabla }\cdot \mathit{T}+{\mathit{F}}_{\mathit{S}},\end{array}$$ 此处我们用 $\mathit{T}=-p\mathit{I}\mathbf{+}\mu (c)\mathit{D}$ 代替了 $\text{(???)}$ 中的压力项和粘性项，${\mathit{F}}_{\mathit{S}}$ 为表面张力导出的体积力。单流体仍然当作是不可压流体，有 $$\frac{\text{d}\rho }{\text{d}t}=0\to \frac{\text{d}((1-c){\rho }_1+c{\rho }_2)}{\text{d}t}=({\rho }_2-{\rho }_1)\frac{\text{d}c}{\text{d}t}=0,$$ 即 $$\boxed{{\displaystyle \frac{\partial c}{\partial t}+\mathit{u}\cdot \mathrm{\nabla }c=0.}}$$ 该方程是不可压缩 VOF 模型中的相方程，也是最常见的对流方程，也与方程 $\text{(???)}$ 一致。VOF 模型中，界面法向向量和曲率可以写为 $$\mathit{n}=\frac{\mathrm{\nabla }c}{\left|\mathrm{\nabla }c\right|},\kappa =-\mathrm{\nabla }\cdot \left(\frac{\mathrm{\nabla }c}{\left|\mathrm{\nabla }c\right|}\right).$$ 为了将表面方程 $\text{(???)}$ 考虑到 VOF 方法中，需要用到连续力（Continuum surface force，CSF）模型，该模型将间断的压力处理成连续过渡的体积力，即（不考虑马兰戈尼效应） $$\mathit{n}\cdot ({\mathit{T}}_1-{\mathit{T}}_2)=\gamma \kappa \mathit{n}\to \mathrm{\nabla }\cdot \mathit{T}=\mathrm{\nabla }(\gamma \kappa )=\gamma \kappa \mathrm{\nabla }c,$$ 结合方程 $\text{(???)}$ 不可压缩 VOF 模型中动量方程 $$\boxed{{\displaystyle \rho (c)(\frac{\partial \mathit{u}}{\partial t}+(\mathit{u}\cdot \mathrm{\nabla })\mathit{u})=-\mathrm{\nabla }p+\mathrm{\nabla }\cdot (\mu (c)\mathit{D})-\gamma \mathrm{\nabla }\cdot \left(\frac{\mathrm{\nabla }c}{\left|\mathrm{\nabla }c\right|}\right)\mathrm{\nabla }c.}}$$

Level-set 方法

ALE 方法