一些前置知识可以参考笔者的先前的文章：
不可压雷诺方程与涡黏模型流体力学基本方程（不可压、含部分细节以及推导）湍流瞬时量的时均化法则湍流的基本方程组利用雷诺分解，将瞬时量分为时均量和脉动量之和，\bm{u}^{*}=\bm{u}+\bm{u}^{\prime}，p^{*}=p+p^{\prime}，这里假设外力是线性的，得到时均运动方程组，即\begin{align*}
\frac{\partial u_{i}}{\partial x_i}&=0
\\
{\frac{\partial u_{i}}{\partial t}}
+
{u_j \frac{\partial u_{i}}{\partial x_j}}
&=
- \frac{1}{\rho}
{\frac{\partial p}{\partial x_i}}
+ \frac{\partial}{\partial x_j}\left(\nu \frac{\partial u_{i}}{\partial x_j} \right)
-\frac{\partial}{\partial x_j}\left(\overline{{u}_{i}^{\prime} {u}_{j}^{\prime}} \right)
+ F_i
\end{align*}在时均运动的雷诺方程组中出现的雷诺应力项（脉动速度二阶相关项）是未知的，从而导致雷诺方程组不封闭。
如果继续建立二阶相关项的输运方程，则会引出三阶相关项的未知量，方程仍然不封闭，以此类推，三阶方程会出现四阶相关项未知量，四阶方程会出现五阶相关项未知量，……，方程组永不能封闭，这就是著名的湍流方程组不封闭问题。为了实现工程湍流的计算，需要使用一些经验假设建立各种脉动量相关的统计值补充非常，特别是雷诺应力的补充方程。
实用湍流理论的实质是以雷诺时均运动方程和有关脉动量输运方程为基础，依靠理论和经验的结合，引进一系列模化假设，建立一组湍流时均运动的封闭方程组来解决湍流工程计算问题。湍涡粘性基于唯象学原理，1877 年 法国 力学家 Boussinesq 首先将湍流脉动产生的附加切应力（后称为雷诺应力）与粘性应力比拟，提出著名的涡粘假设，建立了雷诺应力与时均速度梯度之间的比拟关系。
虽然涡粘性的概念早于雷诺方程组的出现，但是却为后来的工程湍流计算奠定了基础。对于简单的近壁区时均二维流动，湍动切应力（雷诺应力） 可表达为\tau_t = -\rho \overline{u^{\prime} v^{\prime}} = \rho \nu_t \frac{\partial u}{\partial y}其中， \nu_t 为涡粘性系数 (Turbulent or Eddy viscosity)。与之相对应，时均流动产生的粘性切应力为\tau_l =
\rho \nu \frac{\partial u}{\partial y}得到作用于流层之间的总切应力为\tau_0 = \nu +\nu_t = \rho (\nu+\nu_t )\frac{\partial u}{\partial y}这里，涡粘性系数 \nu_t 不是流体的物理属性，而是湍流运动状态的函数。后经 Hinze 等推广到三维流动中，即-\overline{u^{\prime}_i u^{\prime}_j}
=-\frac{2}{3}k \delta_{ij} + \nu_t \left(\frac{\partial u_j}{\partial x_i} + \frac{\partial u_i}{\partial x_j} \right)这样，湍流方程的封闭问题，就变成了如何确定 \nu_t 的大小和分布。如果吧湍流流场看作由一系列大小不同流体微团碰撞和动量交换的结果，则从唯象学的观点出发，类似分子粘性系数（分子粘性系数正比于分子运动的平均速度 C 和平均自由程 l ，\nu \propto C \,l）认为：涡粘性系数
正比于表征大尺度涡（载能涡）运动的特征速度尺度 V 和特征长度尺度 l_t 的乘积，即\nu_t \propto V l_t补充：这里，三个正应力分量为：\begin{align*}
-\overline{u^{\prime}_1 u^{\prime}_1}
=\nu_t \left[2 \frac{\partial u_1}{\partial x_1} - \frac{1}{3}\left(\frac{\partial u_1}{\partial x_1}+\frac{\partial u_2}{\partial x_2}+\frac{\partial u_3}{\partial x_3} \right)\right]- \frac{2}{3}k
\\
-\overline{u^{\prime}_2 u^{\prime}_2}
=\nu_t \left[2 \frac{\partial u_2}{\partial x_2} - \frac{1}{3}\left(\frac{\partial u_1}{\partial x_1}+\frac{\partial u_2}{\partial x_2}+\frac{\partial u_3}{\partial x_3} \right)\right]- \frac{2}{3}k
\\
-\overline{u^{\prime}_3 u^{\prime}_3}
=\nu_t \left[2 \frac{\partial u_3}{\partial x_3} - \frac{1}{3}\left(\frac{\partial u_1}{\partial x_1}+\frac{\partial u_2}{\partial x_2}+\frac{\partial u_3}{\partial x_3} \right)\right]- \frac{2}{3}k \\
\end{align*}或\begin{align*}
-\overline{u^{\prime}_i u^{\prime}_j}
=\nu_t \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{1}{3}\frac{\partial u_k}{\partial x_k}\delta_{ij}\right)-\frac{2}{3}k\delta_{ij}
\end{align*}写成向量形式：\begin{align*}
-\rho \overline{\bm{u}^{\prime} \bm{u}^{\prime}}
=\rho \nu_t \left((\nabla \bm{u})+(\nabla \bm{u})^{T}-\frac{1}{3}(\nabla \cdot \bm{u})\bm{I}\right)-\frac{2}{3}\rho k \bm{I}
\end{align*}应力张量S_{ij}=\frac{1}{2} \left(\frac{\partial u_j}{\partial x_i} + \frac{\partial u_i}{\partial x_j} \right)写成更紧致的形式，其非对角量为S_{ij}^{*}=\frac{1}{2} \left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} -\frac{1}{3} \frac{\partial u_k}{\partial x_k}\delta_{ij}\right)许多学者使用非对角分量使得方程更紧致，但是并没有特殊物理意义。零方程模式（混合长度理论）根据湍流运动特性，\nu_t 的大小和分布在流场中发生明显的变化。按照量纲分析和湍流研究结果，涡粘性 \nu_t 由载能涡的特征长度尺度和特征速度尺度决定，即\nu_t \propto l_t V_t雷诺应力与粘性应力的比值为：\frac{\tau_t}{\tau_l} = \frac{\rho \overline{u^{\prime} v^{\prime}}}{\rho \nu \frac{\partial u}{\partial y}}
=\frac{\rho \nu_t \frac{\partial u}{\partial y}}{\rho \nu \frac{\partial u}{\partial y}}=\frac{\nu_t}{\nu}
=\frac{l_t V_t}{\nu}=Re_t式中， Re_t 表征大尺度湍涡运动特性雷诺数，称为湍动雷诺数。一般取 Re_t \sim 10^3 \sim 10^51925 年普朗特基于分子运动论的比拟，提出混合长度理论 (Mixing-Length Hypothesis)，认为：对于剪切湍流，湍流涡体的特征速度 V_t 正比于时均速度梯度和混合长度的乘积（混合长度——流体质点受湍涡的作用发生自由掺混的平均尺度，与湍涡的平均尺度同量级），即：V_t \propto l_m \left \vert \frac{\partial u}{\partial y} \right \vert所以得到\begin{align*}
\tau_t=-\rho \overline{u^{\prime} v^{\prime}}
=\rho l_m^2 \frac{\partial u}{\partial y} \left \vert \frac{\partial u}{\partial y} \right \vert
,\quad
\nu_t = l_m^2 \left \vert \frac{\partial u}{\partial y} \right \vert
\end{align*}几点说明：在近壁湍流中，靠近壁面附近受壁面影响，脉动速度很小，湍流切应力也很小，但流速梯度很大，粘性切应力起主导作用，速度分布是线性的，这一层称为粘性底层；在粘性底层外区是湍流核心区，此时湍动切应力起主导作用，速度分布符合对数或幂次分布；在湍流核心区和粘性底层区之间为过渡区；粘性底层不是层流，也不是完全湍流，在该层内存在湍斑；粘性底层厚度与壁面粗糙度直接影响沿程损失。在近壁湍流区，假设湍动切应力近似等于壁面切应力 \tau_w，以及混合长度与流体质点到壁面的距离 y 成正比，即 l_m=\kappa y（k 为卡门 (Karman) 常数，k\approx 0.4），得到\begin{align*}
\frac{\tau_w}{\rho}=\kappa^2 y^2 (\frac{\mathrm{d} u}{\mathrm{d} y})^2
\end{align*}积分上式得到著名的近壁区时均速度对数分布曲线：\begin{align*}
\frac{u}{u^{*}}=\frac{1}{\kappa} \ln \frac{u^{*}y}{\nu} + C
\end{align*}其中，C 为常数，u^{*}=\sqrt{\frac{\tau_w}{\rho}} 为摩阻速度。这一模式 建立了涡粘性系数和当地时均速度梯度的关系，现在的问题归结到如何确定未知参数 l_m。对于一些常见的剪切层流动问题，混合长度由比较简单的经验关系式确定。若用 b 表示剪切层的厚度（对于壁面边界层为边界层厚度 \delta，对于自由射流为射流的半厚度），则 l_m 与
b 的关系由下表给出：
1942 年，普朗特通过对自由剪切层（混合层、射流和尾流）流动问题的研究，进一步提出了适用于这类流动的更为简单的涡粘性模式。普朗特认为在 自由剪切层 流动问题中，因为没有壁面的限制和影响，确定涡粘性系数 \nu_t 不能用壁面率而应采用尾迹率，即 \nu_t 在剪切层的任一横断面上是常数，湍动涡体的特征长度尺度正比于剪切层的厚度 b，特征速度尺度正比于剪切层横断面上的最大速度差，\begin{align*}
\nu_t= \alpha b \left \vert u_{max}-u_{min} \right \vert
\end{align*}其中， \alpha 为经验常数。对于静止环境中的自由射流，上式可以写成\nu_t= \alpha b u_{m}关于零方程模型的一些说明：与混合长度理论平行的还有 Taylor 的涡量输运理论、卡门的相似性理论等。由于这些模式只考虑了湍动应力和时均速度梯度的关系，未引入表征湍流高阶统计量的微分方程，故被称为半经验理论或零方程模式(Zero-Equation Model)。半经验理论成功地解决了诸如湍流边界层、湍射流、湍动尾流和管道湍流等一些剪切层流动问题时均物理量的湍流计算，可以说在计算机出现之前是湍流工程计算的主要依据，甚至目前仍然是工程上广泛采用的一类模式。一方程模式虽然混合长度模式较为成功地解决了一些剪切层流动的问题，但由于其在模化过程中仅着眼于涡粘性与时均量的关系，而未考虑湍流扩散和对流输运，即假定湍流处于局部平衡状态，认为流场中任一点处湍动能量的产生和耗散是相等的，这意味着任一点的湍动量不可能通过湍动输运而影响流场中的其他点，存在不合理之处，导致混合长度模式存在下列若干问题：混合长度模式汇总的经验常数缺乏通用性，不同流动问题的经验常数取值不同；混合长度模式不适用于处理那些湍动输运过程起主要作用的流动问题，比如快速发展的流动、弱剪切流动和分离流动等；对于较为复杂的流动问题，混合长度 l_m 不易确定。湍动能 k 方程模式为了克服混合长度模式的缺点，Kolmogorov 和 Prandtl 首先提出了一方程模式 (One-Equation Model)。他们的基本思想是用湍动能
k=\overline{{u}_{i}^{\prime} {u}_{j}^{\prime}}/2 来代替湍流速度尺度 V，即V = \sqrt{k}湍涡粘性系数 \nu_t （\nu_t \propto l_t V_t） 写成\nu_t = C_{\mu}^{\prime} \sqrt{k} L这就是著名的
表达式，式中，C_{\mu}^{\prime} 为通用常数，L 为湍流特征长度尺度，由实验确定。湍动能 k 方程模式由湍动能的精确方程通过模化得到。湍动能 k 的方程为\begin{align*}
{\frac{\partial k}{\partial t}}
+
{u_j \frac{\partial k}{\partial x_j}}
=
\frac{\partial}{\partial x_j} \left( -\overline{\frac{{u}_{i}^{\prime} {u}_{i}^{\prime}}{2} {u}_{j}^{\prime}}
- \frac{ \overline{p^{\prime} {u}_{j}^{\prime}} }{\rho} + \nu \frac{\partial k}{\partial x_j} \right)
-\overline{u^{\prime}_{i} u^{\prime}_{j}} \frac{\partial u_{i}}{\partial x_j}
-\nu \overline{\frac{\partial u^{\prime}_{i}
}{\partial x_j} \frac{\partial u^{\prime}_{i}
}{\partial x_j}}
\end{align*}上式中，需要模化的有湍动能扩散项 \frac{\partial}{\partial x_j} \left( -\overline{\frac{{u}_{i}^{\prime} {u}_{i}^{\prime}}{2} {u}_{j}^{\prime}}
- \frac{ \overline{p^{\prime} {u}_{j}^{\prime}} }{\rho} \right) 和湍动能耗散率 \epsilon。对于扩散项，通常假定与 k 的梯度成正比（梯度型假设），即-\overline{\frac{{u}_{i}^{\prime} {u}_{i}^{\prime}}{2} {u}_{j}^{\prime}}
- \frac{ \overline{p^{\prime} {u}_{j}^{\prime}} }{\rho} = \frac{\nu_t}{\sigma_k}\frac{\partial k}{\partial x_j}其中，\sigma_k 为经验系数。
湍动能耗散虽然主要发生在粘性起作用的小尺度涡范围内，但湍动能耗散 \epsilon 是源自大尺度涡提供的能量，这些大尺度的涡可用 V(=\sqrt{k}) 和 L 来表征。对于湍涡而言，在涡体的运动过程中单位质量所消耗的湍动能，根于量纲分析，得到\left[\epsilon \right]=
\left[\nu \overline{\frac{\partial u^{\prime}_{i}
}{\partial x_j} \frac{\partial u^{\prime}_{i}
}{\partial x_j}}\right]
=\left[LV \frac{V}{L} \frac{V}{L} \right] =
\left[\frac{V^3}{L}\right]这样，\epsilon 的常用的一个模化式为:\epsilon=C_D \frac{k^{3/2}}{L}式中，C_D 是另一个经验系数。使用上述模化假定，湍动能 k 方程模式为\begin{align*}
{\frac{\partial k}{\partial t}}
+
{u_j \frac{\partial k}{\partial x_j}}
=
\frac{\partial}{\partial x_j} \left[ \left(\frac{\nu_t}{\sigma_k}+\nu
\right) \frac{\partial k}{\partial x_j} \right]
-\overline{u^{\prime}_{i} u^{\prime}_{j}} \frac{\partial u_{i}}{\partial x_j}
-C_D \frac{k^{3/2}}{L}
\end{align*}一方程模式中的经验系数，推荐值为 C_{\mu}^{\prime}C_D \approx 0.08，\sigma_k=1.0。应当指出，对湍流场起作用的是经验系数 C_{\mu}^{\prime} 和 C_D 的乘积，而不是它们独立的取值。一个特例：如果流动为定常流动，且不计 k 方程中对流输运和扩散输运，即得到湍动能生成等于耗散，即湍流处于局部平衡状态，也就是0=-\overline{u^{\prime}_{i} u^{\prime}_{j}} \frac{\partial u_{i}}{\partial x_j}
-C_D \frac{k^{3/2}}{L}对于薄剪切层流动 (\tau_t = -\rho \overline{u^{\prime} v^{\prime}} = \rho \nu_t \frac{\partial u}{\partial y})，上式可以进一步简化为\nu_t \left(\frac{\partial u }{\partial y}\right)^2=C_D \frac{k^{3/2}}{L}利用上式，并联立 \nu_t = C_{\mu}^{\prime} \sqrt{k} L 中，消去 k，得到关于 \nu_t 的表达式\nu_t \left(\frac{\partial u }{\partial y}\right)^2=C_D \frac{\left({\frac{\nu_t}{C_{\mu}^{\prime} L}}\right)^{3}}{L}
\Rightarrow
\nu_t= \left({\frac{C_{\mu}^{\prime3}}{C_D}}\right)^{1/2} L^2 \left \vert \frac{\partial u }{\partial y} \right \vert这个公式即为混合长度模式。这个推导清楚地表明，混合长度仅适用于湍流处于局部平衡状态的流动。湍动剪切应力输运方程模式除了用湍动能 k 方程给出的一方程模式外，Bradshaw 针对剪切流动问题，放弃涡粘性的概念，提出用湍动切应力取代湍动能 k，发展了一个用湍动切应力 -\rho \overline{u^{\prime} v^{\prime}} 表征的一方程模式。
则对于二维薄剪切层流动，湍动能 k 的精确方程可简化为\begin{align*}
{\frac{\partial k}{\partial t}}
+{u \frac{\partial k}{\partial x}}
+{v \frac{\partial k}{\partial y}}
=
\frac{\partial}{\partial y} \left[ -\overline{\frac{u^{\prime}_{i} u^{\prime}_{i}}{2} v^{\prime}} -\frac{\overline{p^{\prime }v^{\prime}}}{\rho} \right]
-\overline{u^{\prime}
v^{\prime} } \frac{\partial u}{\partial y}
-\epsilon
\end{align*}为了封闭这个方程，Bradshaw 定义了下列三个量，即\begin{align*}
a_1 &=\frac{-\overline{u^{\prime}
v^{\prime}}}{k} \quad (\text{无量纲数})
\\
L &= \frac{(-\overline{u^{\prime}
v^{\prime}})^{3/2}}{\epsilon}
\quad (\text{长度尺度})
\\
G &=\frac{\left(-\overline{\frac{u^{\prime}_{i} u^{\prime}_{i}}{2} v^{\prime}} -\frac{\overline{p^{\prime }v^{\prime}}}{\rho}\right)}{(-\overline{u^{\prime}
v^{\prime}})(-\overline{u^{\prime}
v^{\prime}})_{\text{max}}^{1/2}}
\quad (\text{无量纲数})
\end{align*}将上面三式代入湍动能方程，得到剪切应力输运模式为:\begin{align*}
\frac{\mathrm{D} }{\mathrm{D}t} \left(\frac{\overline{u^{\prime}
v^{\prime}}}{a_1}\right)
=- \left(\overline{u^{\prime} v^{\prime}}\right) \frac{\partial u}{\partial y} - \frac{\partial}{\partial y} \left[G (-\overline{u^{\prime}
v^{\prime}})(-\overline{u^{\prime}
v^{\prime}})_{\text{max}}^{1/2}\right]
-\frac{(-\overline{u^{\prime}
v^{\prime}})^{3/2}}{L}
\end{align*}这个模式中出现的三个量 a_1，G，L 在求解方程之前必须事先给定。
在壁面边界层流动问题中，Bradshaw 利用大量的实验数据给出了这些量的经验关系式，其中系数 a_1 近似等于0.3，G 和 L 是离壁面无量纲距离 \frac{y}{\delta} 的函数。另外，如果不计这个方程中的对流输运和扩散输运，同样可以退化到混合长度模式，即\begin{align*}
0=-\left(\overline{u^{\prime} v^{\prime}}\right) \frac{\partial u}{\partial y}
-\frac{(-\overline{u^{\prime}
v^{\prime}})^{3/2}}{L}
\end{align*}Bradshaw 的剪应力模式，在许多壁面边界层流动中得到成功的应用，且精度相当高。目前，这一模式已被引申到有换热的不可压和可压缩三维流动的计算中。二方程模式k-\epsilon 模式的建立一方程模式虽然比零方程模式前进了一大步，引入了表征湍流脉动场的 k 方程，但是这类模式仍然含有需要通过实验确定的特征长度尺度 L，对问题的依赖性较强，并没有得到广泛应用。为避免由实验方法确定 L，人们期望用一个微分输运方程来确定 L，由此出现了二方程模式 (Two-Equation Models)。
二方程模式的主要特点是把表征大尺度涡运动的 特征速度尺度 V 和 特征长度尺度 L 均用微分输运方程来描述。
一般，V=\sqrt{k} 仍然用湍动能 k 方程表述；
此外，为了方便推导方程，常常不直接用 L 作为未知量，而是以 Z = k^mL^n 作为未知量，来间接确定 L。
一些学者采用不同的 Z 变量来表征 L 的输运方程，虽然不同的 Z 变量所表达的物理过程不同，但是最终的模化结果是相似的。关于 Z 变量输运方程的一个通用模式是\begin{align*}
{\frac{\partial Z}{\partial t}}
+\underbrace{{u_j \frac{\partial Z}{\partial x_j}} }_{\text{对流输运项}}
=
\underbrace{\frac{\partial}{\partial x_j} \left(\frac{\sqrt{k}L}{\sigma_{Z}}
\right) \frac{\partial Z}{\partial x_j} }_{\text{扩散输运项}}
+\underbrace{C_{Z1} \frac{Z}{k}P}_{\text{生成项}}
-\underbrace{C_{Z2}\frac{\sqrt{k}}{L}Z}_{\text{破毁项}}
+\underbrace{S}_{\text{源项}}
\end{align*}其中，P=-\overline{u^{\prime}_{i} u^{\prime}_{j}} \frac{\partial u_{i}}{\partial x_j} 为湍动能生成项，\sigma_{Z}、 C_{Z1}、 C_{Z2} 为经验系数；S 是针对不同的 Z 所出现的二次源项（主要是在近壁区起作用），在湍动能耗散率 \epsilon 不需要此项。
与湍动能输运方程类似，此方程也具有 对流输运、扩散输运、湍动产生 和 耗散过程。在诸多有关 L 的输运方程中，目前以 \epsilon=\frac{k^{3/2}}{L} 作为未知量应用最广。其主要原因是：\epsilon 作为未知量具有明确的物理意义，表示湍动能耗散率 \left(\nu \overline{\frac{\partial u^{\prime}_{i}
}{\partial x_j} \frac{\partial u^{\prime}_{i}
}{\partial x_j}}\right)，不仅容易推导精确输运方程，而且作为未知量的 \epsilon 直接出现在 k 方程中，不需要再在 k 方程中建立模式。这样对于二方程模式，涡粘性 \nu_t 被看作 k, \epsilon 的函数。由于 L \sim \frac{k^{3/2}}{\epsilon}，由涡粘性系数的定义，有\nu_t \propto \sqrt{k}L = \sqrt{k} \frac{k^{3/2}}{\epsilon} = \frac{k^2}{\epsilon}写成等式关系为\nu_t = C_{\mu}\frac{k^{2}}{\epsilon}式中，C_{\mu} 为经验常数。精确的 k 和 \epsilon 方程为\begin{align*}
{\frac{\partial k}{\partial t}}
+
{u_j \frac{\partial k}{\partial x_j}}
=
\frac{\partial}{\partial x_j} \left( -\overline{\frac{{u}_{i}^{\prime} {u}_{i}^{\prime}}{2} {u}_{j}^{\prime}}
- \frac{ \overline{p^{\prime} {u}_{j}^{\prime}} }{\rho} + \nu \frac{\partial k}{\partial x_j} \right)
-\overline{u^{\prime}_{i} u^{\prime}_{j}} \frac{\partial u_{i}}{\partial x_j}
-\nu \overline{\frac{\partial u^{\prime}_{i}
}{\partial x_j} \frac{\partial u^{\prime}_{i}
}{\partial x_j}}
\end{align*}\begin{align*}
{\frac{\partial \epsilon}{\partial t}}
+u_{j}{\frac{\partial \epsilon}{\partial x_j}}
=&{\frac{\partial
}{\partial x_k}} \left(-\overline{\epsilon^{\prime} u_{k}^{\prime}} - \frac{2}{\rho} \nu \overline{\frac{\partial u_{k}^{\prime} }{\partial x_j} \frac{\partial p^{\prime} }{\partial x_j}} + \nu \frac{\partial \epsilon }{\partial x_k} \right)
-2 \nu \; \frac{\partial^2 u_{i}
}{\partial x_k \partial x_j} \overline{u_{k}^{\prime} \frac{\partial u_{i}^{\prime} }{\partial x_j}}
\\&
-2\nu \; \frac{\partial u_{i}
}{\partial x_j} \left ( \overline{\frac{\partial u^{\prime}_{k}
}{\partial x_i} \frac{\partial u^{\prime}_{k}
}{\partial x_j}} + \overline{\frac{\partial u^{\prime}_{i}
}{\partial x_k} \frac{\partial u^{\prime}_{j}
}{\partial x_k}}\right)
-2\nu \; \overline{\frac{\partial u_{i}^{\prime}
}{\partial x_j} \frac{\partial u_{i}^{\prime}
}{\partial x_k} \frac{\partial u_{j}^{\prime}
}{\partial x_k}}
-2 \; \overline{ \left(\nu \frac{\partial^2 u_{i}^{\prime}
}{\partial x_j \partial x_k} \right)^2}
\end{align*}为了封闭精确的 k 方程和 \epsilon 方程，现对这些方程中出现的未知量给出如下模化 。湍动扩散输运项的梯度型假设：
k 和 \epsilon 方程中的湍动扩散输运项均采用梯度型假设，即（注意表示湍动能的符号 k 和张量角标 x_k中的k 不同）\begin{align*}
\mathrm{diff} (k) &= -\overline{\frac{{u}_{i}^{\prime} {u}_{i}^{\prime}}{2} {u}_{j}^{\prime}}
- \frac{ \overline{p^{\prime} {u}_{j}^{\prime}} }{\rho}
&= \frac{\nu_t}{\sigma_{K}} \frac{\partial k}{\partial x_j}
\\
\mathrm{diff}
(\epsilon) &= -\overline{\epsilon^{\prime} u_{k}^{\prime}} - \frac{2}{\rho} \nu \overline{\frac{\partial u_{k}^{\prime} }{\partial x_j} \frac{\partial p^{\prime} }{\partial x_j}}
&= \frac{\nu_t}{\sigma_{\epsilon}} \frac{\partial \epsilon}{\partial x_k}
\end{align*}对于非各向同性的扩散，可考虑采用如下模式\begin{align*}
\mathrm{diff} ( k )
&= C_{k} \frac{k}{\epsilon} \overline{{u}_{j}^{\prime} {u}_{k}^{\prime} } \frac{\partial k}{\partial x_k}
\\
\mathrm{diff}
(\epsilon )
&= C_{k} \frac{k}{\epsilon} \overline{{u}_{j}^{\prime} {u}_{k}^{\prime} } \frac{\partial \epsilon}{\partial x_k}
\end{align*}湍动耗散的各向同性假设：
根据 Kolmogorov 的局部各向同性假设，在高雷诺数下，湍流大尺度涡决定的性质不受粘性的影响，而小尺度涡结构在统计上则与时均运动和大尺度涡运动无关，是各向同性的。考虑到湍动耗散过程主要决定于各向同性的小尺度运动，因此可以认为湍动耗散也是各向同性的。从各向同性湍流的统计理论出发，得到\nu \overline{\frac{\partial u^{\prime}_{i}
}{\partial x_j} \frac{\partial u^{\prime}_{i}
}{\partial x_j}}
=\frac{\epsilon}{30}(4 \delta_{ij}\delta_{kl}-\delta_{ik}\delta_{jl}-\delta_{il}\delta_{jk})当指标 k=l 时，由上式可得\nu \overline{\frac{\partial u^{\prime}_{i}
}{\partial x_j} \frac{\partial u^{\prime}_{i}
}{\partial x_j}}=\frac{1}{3} \delta_{ij} \epsilon在 \epsilon 方程中产生项-2\nu \; \frac{\partial u_{i}
}{\partial x_j} \left ( \overline{\frac{\partial u^{\prime}_{k}
}{\partial x_i} \frac{\partial u^{\prime}_{k}
}{\partial x_j}} + \overline{\frac{\partial u^{\prime}_{i}
}{\partial x_k} \frac{\partial u^{\prime}_{j}
}{\partial x_k}}\right)
=-2\frac{\partial u_{i}
}{\partial x_j} \left(\frac{2}{3}\delta_{ij}\right) \epsilon
=-\frac{4}{3}\epsilon \frac{\partial u_{j}
}{\partial x_j} = 0可略去不计。
此外，通过量级比较， \epsilon 方程中另一个产生项-2 \nu \; \frac{\partial^2 u_{i}
}{\partial x_k \partial x_j} \overline{u_{k}^{\prime} \frac{\partial u_{i}^{\prime} }{\partial x_j}} = o(\frac{1}{Re_t}) \to 0与湍动雷诺数的倒数同量级。因此，在高雷诺数下，该项也可以略去。小涡拉伸引起的生成项和粘性耗散项的模化I=-2\nu \; \overline{\frac{\partial u_{i}^{\prime}
}{\partial x_j} \frac{\partial u_{i}^{\prime}
}{\partial x_k} \frac{\partial u_{j}^{\prime}
}{\partial x_k}}
-2 \; \overline{ \left(\nu \frac{\partial^2 u_{i}^{\prime}
}{\partial x_j \partial x_k} \right)^2}式中，第一项表示小涡拉伸引起的生成项，相当于湍动能耗散率 \epsilon 的一个源项，从物理角度看，正比于湍动能生成项 P，这是由于 P 的增大会引起湍动能增加，从而引起 \epsilon 也增加。第二项是粘性耗散项，为 \epsilon 方程的一个汇项，应与湍动能耗散率 \epsilon 成正比。
由于这两项之差对
\epsilon
的发展起主要作用，故通常这两项不分开模化。如 Lumley 根据湍流在局部平衡状态下 P=\epsilon 的条件，认为这两项之差应与 \left(\frac{P}{\epsilon} - 1 \right) 成正比，通过量纲分析，可得I \propto \frac{\epsilon^2}{k} \left(\frac{P}{\epsilon} - 1 \right)写成等式关系为I = C_{\epsilon 1}\frac{\epsilon}{k} P - C_{\epsilon 2}\frac{\epsilon^2}{k} = \frac{\epsilon}{k} \left( C_{\epsilon 1} P -
C_{\epsilon 2} \epsilon \right)将以上各模化项代入精确的 k 和 \epsilon 的精确方程式，得到标准的 k-\epsilon 模式：\begin{align*}
{\frac{\partial k}{\partial t}}
+
{u_j \frac{\partial k}{\partial x_j}}
&=
\frac{\partial}{\partial x_j} \left[\left(\frac{\nu_t}{\sigma_{K}} + \nu \right) \frac{\partial k}{\partial x_j} \right] + P - \epsilon
\\
{\frac{\partial \epsilon}{\partial t}}
+u_{j}{\frac{\partial \epsilon}{\partial x_j}}
&=
\frac{\partial}{\partial x_j} \left[\left(\frac{\nu_t}{\sigma_{\epsilon}} + \nu \right) \frac{\partial \epsilon}{\partial x_j} \right] +
C_{\epsilon 1}\frac{\epsilon}{k} P - C_{\epsilon 2}\frac{\epsilon^2}{k}
\end{align*}其中，P=-\overline{u^{\prime}_{i} u^{\prime}_{j}} \frac{\partial u_{i}}{\partial x_j} 为湍动能生成项。k-\epsilon 模式中系数的确定C_\mu 的确定 (\nu_t = C_{\mu}\frac{k^{2}}{\epsilon})在近壁剪切湍流中，湍流处于局部平衡区。有湍动能方程式有P=\epsilon, \quad
-\overline{u^{\prime} v^{\prime}} \frac{\partial u}{\partial y} = \epsilon, \quad
-\overline{u^{\prime} v^{\prime}}=\nu_t \frac{\partial u}{\partial y}注C_{\mu}=\nu_t \frac{\epsilon}{k^2}=\frac{-\overline{u^{\prime} v^{\prime}}}{\frac{\partial u}{\partial y}} \frac{\epsilon}{k^2} = \frac{-\overline{u^{\prime} v^{\prime}}}{\frac{\epsilon}{-\overline{u^{\prime} v^{\prime}}}} \frac{\epsilon}{k^2}=\frac{(-\overline{u^{\prime} v^{\prime}})^2}{k^2}
由此可以得到C_{\mu} = \left(\frac{-\overline{u^{\prime} v^{\prime}}}{k} \right)^2在湍流边界层的常应力区，实验结果表明 -\overline{u^{\prime} v^{\prime}}/k=0.3，得到 C_{\mu}=0.09。C_{\epsilon 2} 的取值根据网格湍流的衰变指数，确定 C_{\epsilon 2} 的值。对于栅格后均匀衰变 （没有生成项）的湍流，由 k 和 \epsilon 得到\begin{align*}
{\frac{\partial k}{\partial t}}&=-\epsilon \\
{\frac{\partial \epsilon}{\partial t}}&=-C_{\epsilon 2}\frac{\epsilon^2}{k}
\end{align*}对于均匀衰变的湍流，湍动能 k \sim t^{-m}，代入 k 方程得到 \epsilon \sim t^{-(m+1)}，从而得到\begin{align*}
-(m+1)t^{-(m+2)}
& \sim C_{\epsilon 2}\frac{\epsilon^2}{k}
\sim -C_{\epsilon 2} t^{-(m+1)} \frac{m t^{-(m+1)}}{t^{-m}}
\\
C_{\epsilon 2} &= \frac{m+1}{m}
\end{align*}对于网格后衰变的湍流，测得的衰变指数为 m=1.1 \sim 1.3， 代入上式，得到 C_{\epsilon 2}=1.92。C_{\epsilon 1} 和 \sigma_{\epsilon} 的取值利用壁湍流和均匀剪切湍流的结果，可以近似确定 C_{\epsilon 1}和 \sigma_{\epsilon} 的值。
在壁湍流等应力区，平均速度是对数分布 u =u_{\tau} \frac{1}{\kappa} \ln {y} + C。雷诺应力表达式为{-\overline{u^{\prime} v^{\prime}}}=\frac{\tau_0}{\rho}=u_{\tau}^2=\nu_t \frac{\partial u}{\partial y} = C_{\mu} \frac{k^2}{\epsilon} \frac{u_\tau}{\kappa y},通过上式得到耗散率的表达式为\epsilon = C_{\mu} \frac{k^2}{\kappa y u_{\tau}}另一方面，在等应力区中，湍动能生成和耗散平衡，有以下等式：\begin{align*}
\epsilon=P
\Rightarrow
-\overline{u^{\prime} v^{\prime}} \frac{\partial u}{\partial y}
=\nu_t\left(\frac{\partial u}{\partial y}\right)^2
=C_{\mu} \frac{k^2}{\epsilon} \left(\frac{u_{\tau}}{\kappa y}\right)^2
\end{align*}由上式得到耗散率的另一个表达式是\epsilon = \frac{C_{\mu}^{1/2} \kappa u_{\tau}}{\kappa y}对以上两个耗散率的表达式做乘积，可以消去 u_{\tau}，得到\epsilon =\frac{ C_{\mu}^{3/4} k^{3/2}}{\kappa y}在对数层中，雷诺应力和湍动能近似为常数，湍动能耗散的流向梯度近似为零，生成项和耗散项平衡 (P=\epsilon)，于是有湍动能耗散的输运方程为\begin{align*}
0&=\frac{\partial}{\partial x_j} \left[\left(\frac{\nu_t}{\sigma_{\epsilon}} + \nu \right) \frac{\partial \epsilon}{\partial x_j} \right] +
C_{\epsilon 1}\frac{\epsilon}{k} P - C_{\epsilon 2}\frac{\epsilon^2}{k}
\\
\Rightarrow
0&=\frac{\mathrm{d}}{\mathrm{d} y} \left(\frac{\nu_{\tau}}{\sigma_\epsilon} \frac{\mathrm{d}\epsilon}{\mathrm{d} y} \right) + C_{\epsilon 1} \frac{\epsilon^2}{\kappa} -
C_{\epsilon 2} \frac{\epsilon^2}{\kappa}
\end{align*}将 \epsilon 代入上式（ \nu_t 也是 \epsilon 的函数，需要挨个代入，两次求导化简，最后得到 C_{\epsilon 2}-C_{\epsilon 1} 的表达式 ），经简化后可得C_{\epsilon 2}-C_{\epsilon 1}=\frac{\kappa^2}{\sigma_\epsilon C_{\mu}^{1/2}}此外，还需要一个方程来确定 C_{\epsilon 1} 和 \sigma_{\epsilon}。可以利用均匀剪切湍流的实验结果来确定另一个关系。
在均匀剪切湍流流场中，湍流统计量的空间导数等于零，因此它的 k 和 \epsilon 方程可以简化为\begin{align*}
{\frac{\partial k}{\partial t}}
&= P - \epsilon
\\
{\frac{\partial \epsilon}{\partial t}}
&= C_{\epsilon 1}\frac{\epsilon}{k} P - C_{\epsilon 2}\frac{\epsilon^2}{k}=\frac{\epsilon^2}{k} \left(C_{\epsilon 1} \frac{P}{\epsilon} - C_{\epsilon 2} \right)
\end{align*}实验和直接数值模拟都证实，均匀剪切湍流中 k/\epsilon 趋向于常数，即 \left[\frac{\mathrm{d}}{\mathrm{d} t} \left(\frac{k}{\epsilon} \right) \right]_{t\to \infty}=0，将上面两式代入可得 （将 \left[\frac{\mathrm{d}}{\mathrm{d} t} \left(\frac{k}{\epsilon} \right) \right] 分部求导展开，然后将 {\frac{\partial k}{\partial t}} 和 {\frac{\partial \epsilon}{\partial t}} 代入，得到最终关系式）\begin{align*}
C_{\epsilon 1}=1+\frac{C_{\epsilon 2} - 1}{P/\epsilon}
\end{align*}在均匀剪切湍流中，P/\epsilon=1.4，可以估算 \sigma_\epsilon 和 C_{\epsilon 1}，建议采用C_{\epsilon 1}=1.44,\quad \sigma_\epsilon=1.3\sigma_{k} 的取值湍动能扩散发生在非均匀湍流场中，很难在简单的湍流中确定湍动能扩散系数 \sigma_{k}，假定湍动能传输和动量传输以相同的机制进行，最简单的模型假定是\sigma_{k}=1.0经过大量计算和修正，目前多数学者推荐的各常数取值为 C_{\mu}=0.09，\sigma_{k}=1.0，\sigma_{\epsilon}=1.3，C_{\epsilon 1}=1.44，C_{\epsilon 2}=1.92。
最后，将以上的模式常数汇总（由 Launder 和 Spalding (1974) 最早建议），合并称之为标准的 k-\epsilon 模式常数：
k-\epsilon 模式被广泛应用于湍流的工程计算中，现已经得到许多成功的算例，如壁湍流、湍射流、突扩分离流和其他的一些剪切流动问题。标量的 k-\epsilon 模式定义标量能量 k_{\theta}=\overline{\theta^{\prime}\theta^{\prime}}/2，可建立标量 \theta^{*}=\theta+\theta^{\prime} 的能量 k_{\theta} 的微分输运方程：\begin{align*}
{\frac{\partial k_{\theta}}{\partial t}}
+
{u_j \frac{\partial k_{\theta}}{\partial x_j}}
&=
\kappa_t \frac{\partial \theta}{\partial x_i} \frac{\partial \theta}{\partial x_i}-
\frac{\partial}{\partial x_j} \left[\left(\frac{\kappa_t}{\sigma_{\theta}} + \kappa \right) \frac{\partial k_{\theta}}{\partial x_j} \right]
-\epsilon_{\theta}
\end{align*}其中，\kappa_t 是涡扩散系数，在 k_{\theta}-\epsilon_{\theta} 模式中，它等于\kappa_t = C_t k \sqrt{\left(\frac{k}{\epsilon} \frac{k_{\theta}}{\epsilon_{\theta}}\right)}\epsilon_{\theta}=\nu \frac{\partial \theta^{\prime}}{\partial x_i} \frac{\partial \theta^{\prime}}{\partial x_i} 是标量输运的耗散项，类似于湍动能耗散的模型方程，它的输运方程为\begin{align*}
{\frac{\partial \epsilon_{\theta}}{\partial t}}
+u_{j}{\frac{\partial \epsilon_{\theta}}{\partial x_j}}
=&
-C_{\theta 1} \frac{\epsilon_\theta}{k_\theta} \overline{u_i^{\prime} \theta^{\prime}} \frac{\partial \theta}{\partial x_i}
-\frac{\partial}{\partial x_j} \left[\left(\frac{\kappa_t}{\sigma_{\epsilon_\theta}} + \kappa \right) \frac{\partial \epsilon_\theta}{\partial x_j} \right] \\
&-C_{\theta 2}\frac{\epsilon_\theta^2}{k_\theta}
-C_{\theta 3}\overline{u_i^{\prime} u_j^{\prime}} \frac{\partial u_i}{\partial x_j} \frac{\epsilon_\theta}{k}
+C_{\theta 4}\frac{\epsilon \epsilon_\theta}{k}
+\kappa \frac{\partial^2 \epsilon_\theta}{\partial x_j \partial x_j}
\end{align*}上式中标量通量 -\overline{u_i^{\prime} \theta^{\prime}}=\kappa_t \frac{\partial \theta}{\partial x_i}。
与 k-\epsilon
模式类似，k_{\theta}-\epsilon_{\theta} 模式常数用简单湍流场的标量输运实验确定（ Nagano 和 Kim(1988) ），推荐
关于标准 k-\epsilon 模式的评价标准的 k-\epsilon 模式可以计算比较复杂的湍流，但是在定量结果方面并没有比代数涡粘模式有明显的优势。
标准 k-\epsilon
模式的主要缺点是：标准 k-\epsilon
模式假定雷诺应力和当时当地的平均切变率成正比，所以它不能准确反映雷诺应力沿着流向 的 历史效应；标准 k-\epsilon
模式的涡粘系数是标量，不能反映雷诺正应力的 各向异性，尤其是近壁湍流，雷诺正应力具有明显的各向异性。（例如，方管湍流中的二次流是由于雷诺正应力之差产生的，标准
模式不能正确地表达雷诺正应力，因此不能预测到方管湍流的二次流）标准 k-\epsilon
模式不能反映平均涡量的影响，而平均涡量对雷诺应力的分布有影响，特别是在湍流分离流中，这种影响是十分重要的。如果能够克服标准
k-\epsilon
模式的缺点，他将有更好的预测结果。此后又在 标准
k-\epsilon
模式之上提出了诸多改进，可参阅相关文献。参考文献[1] 刘沛清. 湍流模式理论[M]. 北京：科学出版社,2020.[2]张兆顺,崔桂香,许春晓.湍流理论与模拟[M].清华大学出版社,2005.[3]Fluid Mechanics 101—Reynolds-Averaged Navier Stokes (RANS)本文使用 Zhihu On VSCode 创作并发布}