主要参考了文献[4] ¹、[5] ²。但是要注意，很多模型有多种形式，在书中难免有出入，最终还是要以模拟软件采用的为准（对于FLUENT而言就是Fluent Theory Guide³）。

控制方程离散化

有限差分法（FDM）

主要思想为：差商代替微商，递推代替积分，从而完全使用有限的网格节点代替连续的求解域。例如对于一维情况，可设想\((x,t)\)构成网格，而物理量均离散化到各点上（差商代替微商），关键在于找出格点之间的关系（递推关系）。差分方程可以是双层/多层。可以是显式/隐式，但必须满足三项条件：

• 相容性：步长\(h\to0\)时，查分方程逼近微分方程
• 收敛性：查分方程的真解逼近微分方程的精确解
• 稳定性：差分方程的近似解与真解间的误差有界

例如，一维狭缝摩擦流的速度场满足抛物线形方程\(\frac{\partial u}{\partial t}=\nu\frac{\partial^2 u}{\partial y^2}\)，其中\(\nu\)为运动粘度。将\((x_i,t_n)\)处的格点速度记为\(u_i^n\)，若取空间步长为\(\Delta x\)，时间步长为\(\Delta t\)，则差分方程为 \[ \frac{u_i^{n+1}-u_i^n}{\Delta t}=\nu\frac{u_{i+1}^n-2u_i^n+u_{i-1}^n}{\Delta x^2} \] 已知初始条件即\(u_0\)，边界条件即\(u^0\)与\(u^N\)，则可得出递推公式： \[ \begin{pmatrix} u_1^{n+1}\\ u_2^{n+1}\\ \vdots\\ u_{I-2}^{n+1}\\ u_{I-1}^{n+1} \end{pmatrix}= \begin{pmatrix} \frac{\nu\Delta t}{\Delta x^2} & 1-\frac{2\nu\Delta t}{\Delta x^2} & \frac{\nu\Delta t}{\Delta x^2} & \cdots & 0 & 0 & 0\\ 0 & \frac{\nu\Delta t}{\Delta x^2} & 1-\frac{2\nu\Delta t}{\Delta x^2} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & \cdots & 1-\frac{2\nu\Delta t}{\Delta x^2} & \frac{\nu\Delta t}{\Delta x^2} & 0\\ 0 & 0 & 0 & \cdots & \frac{\nu\Delta t}{\Delta x^2} & 1-\frac{2\nu\Delta t}{\Delta x^2} & \frac{\nu\Delta t}{\Delta x^2} \end{pmatrix} \begin{pmatrix} u_0^{n}\\ u_1^{n}\\ u_2^{n}\\ \vdots\\ u_{I-2}^{n}\\ u_{I-1}^{n}\\ u_{I}^{n} \end{pmatrix} \] 该差分式的稳定性依赖于时间步长，在\(\Delta t\le\frac{\Delta x^2}{2\nu}\)时结果稳定。

有限元法（FEM）

主要思想为：同样划分网格，仅在网格节点储存值，但是使用基于网格节点值的插值函数（主要与形函数有关）代表连续的求解域，更为合理，在固体力学中广泛使用。详见有限元学习笔记。

有限体积法（FVM）

主要思想为：同样划分网格，仅在网格节点储存值，之后在每个格点周围划分控制体积，形成一种交错状的网格。在每个控制体积中，以控制方程（特定物理量的守恒方程）为基础推导出代数方程；控制体积之间的壁面上的情况则通过插值得出。因而有限体积法先天性满足守恒关系。

以一维稳态情况为例，对于某物理量\(\phi\)，其控制方程为（来源于N-S方程）： \[ \frac{\mathrm{d}(\rho u \phi)}{\mathrm{d}x}=\frac{\mathrm{d}}{\mathrm{d}x}\left(\Gamma\frac{\mathrm{d}\phi}{\mathrm{d}x}\right)+S \] 其中\(\phi\)可以是速度、温度、湍动能，在连续性方程中还可以是\(1\)；\(\Gamma,S\)为与之对应的广义扩散系数和广义源项。

记格点\(P​\)的两侧界面分别为\(w,e​\)，两侧格点分别为\(W,E​\)。则界面上的物理量可用两侧节点上的物理量插值得出，即\(\phi_w=f(\phi_W,\phi_P)​\)。为简化，设\(S​\)为线性，\(S=S_C+S_P \phi_P​\)，控制方程最终可化为： \[ a_P\phi_P=a_W\phi_W+a_E\phi_E+b \] 其中\(b=S_C\Delta x\)，而其他三个参数的形式取决于插值方式，有以下几种：

• 中心差分格式：\(\phi_w=\frac{\phi_W+\phi_P}{2}\)，条件稳定，扩散项中一般采用，不适合对流项

• 一阶迎风格式：界面物理量等于对流上游节点的物理量，即 \[ \left\{ \begin{aligned} \phi_w=\phi_W,\phi_e=\phi_P\qquad u_w>0\ \ \mathrm{and}\ \ u_e>0\\ \phi_w=\phi_P,\phi_e=\phi_E\qquad u_w<0\ \ \mathrm{and}\ \ u_e<0 \end{aligned} \right. \] 绝对稳定，但精度为一阶

• 混合格式、指数格式、乘方格式等等

• 高阶插值，需要使用更多节点，例如二阶迎风，在\(u_w>0\ \ \mathrm{and}\ \ u_e>0\)时， \[ \phi_w=\frac{3}{2}\phi_W-\frac{1}{2}\phi_{WW},\phi_e=\frac{3}{2}\phi_E-\frac{1}{2}\phi_{EE} \] 绝对稳定，精度较一阶高，但仍有假扩散。改进的QUICK格式更加优越，绝对稳定，假扩散小。

对于瞬态问题，要考虑时间积分，同样有插值问题（统一节点的心值与旧值之间），最常用的是全隐式，此时离散方程表达为： \[ a_P\phi_P=a_W\phi_W+a_E\phi_E+a_P^0\phi_P^0+b \] 需要对每个时间步列出线性方程组求解。

对于多维情况，控制方程为 \[ \frac{\partial(\rho \phi)}{\partial t}+\mathrm{div}(\rho\pmb u \phi)=\mathrm{div}(\Gamma\:\mathrm{grad}\:\phi)+S \] 离散方程通式为 \[ a_P\phi_P=\sum a_{nb}\phi_{nb}+b \] \(nb\)表示周围节点，对于瞬态问题\(b\)包括了\(a_P^0\phi_P^0\)，系数\(a\)取决于控制体积界面上的\(F=\rho u,D=\frac{\Gamma}{\delta x}\)以及\(\Delta V,\Delta t\)。

有限体积法有四条原则：

1. 控制体积界面上物理量表达式连续
2. 各系数\(a\)均为正数
3. \(S_P<0\)
4. \(b=0\)时，\(a_P=\sum a_{nb}\)

SIMPLE算法

Semi-Implicit Method for Pressure Linked Equations，求解压力耦合方程的半隐式方法，并不是什么图森破的算法。

在有限体积法中，若已知流场函数，其他物理量均可用离散控制方程得出；然而若将\(\pmb u\)本身带入离散控制方程，会得出三个动量方程和一个连续性方程，其中\(u,v,w,p\)相互耦合，采用耦合式解法联立求解耗时过长，因而产生了分离式解法，顺序、逐个的求解各变量代数方程组。其中包括消去\(p\)的非原始变量法，但是壁面条件难以给定；最常用的是基于“猜测-修正”迭代过程的SIMPLE算法⁴。

交错网格

使用普通网格计算流场时，压力梯度为 \[ \frac{\partial p}{\partial x}=\frac{p_E-p_W}{2\delta x} \] 这样无法检测出高频震荡的压力波。因此采用交错网格，即标量（包括\(p\)）使用主控制体积，\(u\)的控制体积在\(x\)方向上与其错位半格，\(v\)的控制体积在\(y\)方向上与其错位半格，这样\(\frac{\partial p}{\partial x}=\frac{p_E-p_P}{(\delta x)_e}\).

采用交错网格对动量方程离散时，标量储存在主控制体积节点上，速度分量储存在各自的错位节点上，因而\(F,D\)都要采用新的编号系统计算。另外边界上的错位控制体积会变大。

SIMPLE算法的步骤

以二维情况为例：

1. 假定一个速度分布，记为\(u^0,v^0\)，依此计算动量离散方程中的系数及常数项

2. 假定一个压力场\(p^*\)

3. 依次求解两个动量方程，得\(u^*,v^*\)，动量方程为 \[ \left\{ \begin{aligned} a_eu_e=\sum a_{nb}u_{nb}+b+(p_P-p_E)A_e\\ a_nv_n=\sum a_{nb}v_{nb}+b+(p_P-p_N)A_n\\ \end{aligned} \right. \]

然而我并没有看明白\(a_w,a_e,a_n,a_s\)该怎么求

4. 求解压力修正值方程，得\(p'\)。压力修正值方程为 \[ a_pp_P'=a_Ep_E'+a_Wp_W'+a_Np_N'+a_Sp_S'+b\\ \mathrm{where}\quad b=\frac{(\rho_P^0-\rho_P)\Delta x\Delta y}{\Delta t}+(\rho_wu_w^*-\rho_eu_e^*)\Delta y+(\rho_su_s^*-\rho_nu_n^*)\Delta x \]

5. 根据\(p'\)修正速度，得\(u,v\)，速度修正公式为 \[ \left\{ \begin{aligned} u_e=u_e^*+\frac{A_e}{a_e}(p_P'-p_E')\\ v_n=v_n^*+\frac{A_n}{a_n}(p_P'-p_N') \end{aligned} \right. \]

6. 若有通过源项等与速度场耦合的物理量\(\phi\)，则根据\(u,v\)求解\(\phi\)

7. 利用改进后的速度场重新计算动量离散方程的系数，并用改进后的压力场作为下一层迭代计算的初值，重复上述步骤，直到获得收敛的解。

改进

1. SIMPLEC（SIMPLE Consistent）算法⁵：将速度修正公式变为 \[ u_e=u_e^*+\frac{A_e}{a_e-\sum a_{nb}}(p_P'-p_E') \] 收敛特性远优于SIMPLE

2. 基于同位网格的SIMPLE算法：只需要一套网格，采用新的界面速度插值公式以克服2-\(\delta\)压差的问题

3. 基于非结构网格的SIMPLE算法：此时主要需要考虑节点连线与控制体积界面斜交的问题

三维湍流模型

雷诺数较高时，由层流转变为湍流，湍流瞬时速度\(u=\bar u+u'\)可认为是时均速度与随机脉动量的叠加。在湍流计算中，我们真正关心的不是随机脉动量，而是时均项，因而推导出时均形式的N-S方程，即Reynolds方程（采用张量的爱因斯坦求和约定，忘了？去看笔记 ）：

\[ \begin{aligned} \frac{\partial(\rho \overline{u}_i)}{\partial t}+ \frac{\partial(\rho \overline{u}_i\overline{u}_j)}{\partial x_j}&= -\frac{\partial\overline{p}}{\partial x_i}+ \frac{\partial}{\partial x_j}\left( \mu\frac{\partial\overline{u}_i}{\partial x_j} -\rho\overline{u_i'u_j'}\right)\\ \frac{\partial(\rho \overline{\phi})}{\partial t}+ \frac{\partial(\rho \overline{u}_j\overline{\phi})}{\partial x_j}&= \frac{\partial}{\partial x_j}\left( \Gamma\frac{\partial\overline{\phi}}{\partial x_j} -\rho\overline{u_j'\phi'}\right)+S \end{aligned} \]

由于引入了代表湍流脉动影响的Reynolds应力项\(-\rho\overline{u_i'u_j'}\)，仅靠上述方程无法封闭，需要附加关系式。可以通过对以上式子进行运算，再取时均，获得关于Reynolds应力的偏微分方程，这种方法称为Reynolds应力方程法。另一种方法不直接处理Reynolds应力，而是引入了湍流粘性系数\(\mu_t\)，称为湍流粘性系数法。仿照层流的本构方程，建立Reynolds应力模型（其中的物理量默认为时均量）：

\[ -\rho\overline{u_i'u_j'}= \mu_t\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right) -\frac23(\rho k+\mu_t\mathrm{div}\:\pmb u)\delta_{ij} \]

值得指出，虽然此式并无相应的物理基础，从纯学术的观点来看并不正确，但是以该式为基础的一些湍流模型仍能获得不少有实用意义的结果，因而工程计算中仍广为采用⁶。

其中\(k\)为单位质量流体湍流脉动动能： \[ k=\frac12(\overline{u'^2}+\overline{v'^2}+\overline{w'^2}) \] 对于其他物理量，类似于\(\mu_t​\)可以引入\(\Gamma_t​\)，于是计算湍流的关键变为将\(\mu_t​\)与湍流时均量联系起来。方法包括零方程模型（即混合长度模型，仅使用代数方程）、一方程模型（引入\(k​\)的输运方程），以及广泛使用的两方程模型（\(k\text{-}\varepsilon​\)模型）。

标准\(k\text{-}\varepsilon\)模型⁷

引入脉动动能耗散率\(\varepsilon\)，即小尺度涡的机械能转化为内能的速率： \[ \varepsilon=\frac\mu\rho\overline{\left(\frac{\partial u_i'}{\partial x_k}\right)\left(\frac{\partial u_i'}{\partial x_k}\right)} \] 则\(\mu_t=\rho C_\mu\frac{k^2}{\varepsilon}\)，其中\(C_\mu\)为经验常数。这样可以得出\(k\)与\(\varepsilon\)各自的控制方程，其中包含多个参数，FLUENT中默认为：

\[ C_{1\varepsilon}=1.44,\ C_{2\varepsilon}=1.92\ ,C_\mu=0.09,\ \sigma_k=1.0,\ \sigma_\varepsilon=1.3 \] 与其他物理量的控制方程相比，可发现所有控制方程都满足此形式： \[ \frac{\partial(\rho \phi)}{\partial t}+\mathrm{div}(\rho\pmb u \phi)=\mathrm{div}(\Gamma\:\mathrm{grad}\:\phi)+S \] 根据\(\phi\)的不同，扩散系数\(\Gamma\)和源项\(S\)不同，参见文献文献[5]P122⁸。

对于该模型的适用性，要注意以下几点：

• 对默认参数不要抱太高希望，这些常数也不是总是适合，有时需要修改
• 标准\(k\text{-}\varepsilon\)模型是针对充分发展的湍流（即高雷诺数）的，不适用于层流和壁面附近的流动
• 标准\(k\text{-}\varepsilon\)模型假定\(\mu_t\)为各向同性，因而对于强旋流、弯曲流线流等情况会出现失真

近壁面的处理

在壁面附近粘性占主导，分为三层：粘性底层、过渡层和湍流充分发展的对数律层，不适用标准\(k\text{-}\varepsilon\)模型，有两种方法处理：一是保持网格，采用半经验公式将壁面物理量与湍流核心区的物理量关联起来；另一种方法采用低雷诺数\(k\text{-}\varepsilon\)模型求解，这要求越靠近壁面网格越细。

1. 壁面函数法⁹：不考虑粘性底层和过渡层，确保壁面处第一个内节点在对数律层，此时无量纲的速度\(u^+\)（考虑了\(k\)和壁面切应力\(\tau_W\)）、无量纲的温度\(T^+\)（考虑了\(k\)和热流密度\(q_W\)）均关于广义坐标\(y^+\)满足对数律。

之后按照对数律即可确定\(\tau_W,q_W\)，结合壁面上的速度\(u_W\)以及温度\(T_W\)，可得到当量粘性系数\(\mu_t\)以及当量导热系数\(\lambda_t\)： \[ \tau_W=\mu_t\frac{u_P-u_W}{y_P},\ q_W=\lambda_t\frac{T_P-T_W}{y_P} \] 将\(\mu_t,\lambda_t\)带入\(k\)方程，边界条件一般取为\(\left(\frac{\partial k}{\partial y}\right)_W=0\)，即可求解第一个内节点出的\(k\)方程。而\(\varepsilon\)方程的边界条件难以确定，一般采用其他模型解出。

2. 低雷诺数\(k\text{-}\varepsilon​\)模型¹⁰¹¹：对原本的\(k\text{-}\varepsilon​\)模型进行了三方面修正，使控制方程中扩散系数同时包括湍流扩散系数和分子扩散系数两部分；使系数考虑了湍流雷诺数\(\mathit{Re}_t=\frac{\rho k^2}{\mu\varepsilon}​\)的影响；在\(k​\)方程中考虑到了壁面附近湍流脉动动能耗散的各向异性。固壁上的边界条件为\(k=\varepsilon=0​\)。根据不同的附加项构造，有多种模型。

改进的\(k\text{-}\varepsilon\)模型

1. RNG \(k\text{-}\varepsilon\)模型¹²：对脉动频谱进行滤波，在大尺度运动和修正后的黏度项体现小尺度的影响，从而使小尺度运动有系统地从控制方程中去除。该模型通过修正湍流粘性系数，考虑了平均流动中的旋转及旋流，并且引入了主流的时均应变率： \[ S_{ij}=\frac12\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right) \]

2. Realizable \(k\text{-}\varepsilon​\)模型¹³：时均应变率特别大时，标准\(k\text{-}\varepsilon​\)模型中会出现负的正应力，因此修正\(C_\mu​\)，使其不再是常数，而是与应变率联系起来： \[ C_\mu=\frac1{A_0+A_sU^*\frac k\varepsilon}\\ \mathrm{where} \left\{ \begin{aligned} A_0&=4.0\\ A_s&=\sqrt6\cos\left[\frac13\arccos(\sqrt6W)\right]\\ W&=\frac{S_{ij}S_{jk}S_{ki}}{\overline S^3}\\ \overline S&=\sqrt{S_{ij}S_{ij}}\\ U^*&=\sqrt{S_{ij}S_{ij}+\widetilde\Omega_{ij}\widetilde\Omega_{ij}}\\ \widetilde\Omega_{ij}&=\Omega_{ij}-2\varepsilon_{ijk}\omega_k\\ \Omega_{ij}&=\overline{\Omega}_{ij}-\varepsilon_{ijk}\omega_k \end{aligned} \right. \]

\(\sqrt{S_{ij}S_{ij}}\)之类的迷惑行为是原文写法，猜测只是\(|S_{ij}|\)的委婉表示方式……

其中\(\widetilde\Omega_{ij}\)是从角速度为\(\omega_k\)的参考系中观察到的时均转动速率，体现了旋转的影响，比标准\(k\text{-}\varepsilon​\)模型更为有效。

注意这两种改进模型也都是只适合于高雷诺数的。

Reynolds应力方程模型（RSM）

为了克服\(k\text{-}\varepsilon​\)模型中的问题，需要对\(-\rho\overline{u_i'u_j'}​\)直接建立控制方程求解。有多种形式的RSM，FLUENT中采用的是12方程模型。由于根据时均化法则\(\overline{u_i'u_j'}=\overline{u_iu_j}-\overline{u}_i\overline{u}_j​\)，因此对\(\overline{u_iu_j}​\)和\(\overline{u}_i\overline{u}_j​\)的控制方程进行操作，即可得到Reynolds应力的控制方程，但是同时又引入了三阶时均量\(\overline{u_i'u_j'u_k'}​\)以及压力脉动值\(p'​\)，必须补充模型将其与时均量以及低阶量联系起来。最终得到的方程中：

• 12个未知量：4个时均量（\(u,v,w,p\)）、6个Reynolds应力（\(\overline{u'^2},\overline{v'^2},\overline{w'^2},\overline{u'v'},\overline{v'w'},\overline{w'u'}\)）、湍流脉动动能\(k\)及其耗散率\(\varepsilon\)
• 12个方程：每个未知量都有自己的控制方程，他们相互耦合，联立之后可以求解

RSM也属于高雷诺数模型，在壁面需要用壁面函数或低雷诺数RSM。相较于\(k\text{-}\varepsilon\)模型，RSM能够考虑到一些各向异性效应，但对于一般回流流动并不优于\(k\text{-}\varepsilon\)模型且计算量大，不过还是有许多文献认为RSM颇具潜力。