常用的湍流模型

2021-02-09

沃纳·海森堡曾说：“当我见到上帝后，我一定要问他两个问题——什么是相对论，什么是湍流。但我相信他应该只对第一个问题有了答案。”

流体力学中为便于研究，对流体特性引入假设，而将流动分成两大类：无粘流动（有势流动,简称“势流”）和粘性流动（简称“粘流”，包括层流和湍流）。目前势流计算相对成熟，它具有简易、计算量小等特点而成为切实的工程方法。就CFD技术本身功能发展而言，目前主要是处理粘性流动问题，特别是湍流，理论上来说，对于任何空间点，湍流中流体运动都满足非线性N-S方程，所以至少可以用数值计算机做数值运算，然而这一工作量是十分巨大的。数值求解这一方程的方法有三大类：直接数值模拟方法(DNS)，大尺度涡动模拟方法(LES)，雷诺平均方法(RANS)，由于计算机硬件资源的限制，前二种方法暂还没有得到广泛的应用，而雷诺平均方法在工程上得到了大量的应用。而雷诺平均方法(RANS)是一个湍流模型大家族，包含了很多的子模型，因而本文将主要对RANS模型下的一些子模型进行简要介绍。

控制方程

在介绍常用的湍流模型之前，先把描述流体流动的控制方程进行罗列。

连续方程 \[ \begin{equation} \frac{\partial \rho }{\partial t}+\boldsymbol{\nabla}\cdot \left ( \rho \vec{\boldsymbol{U}} \right )= 0 \end{equation} \] 动量方程 \[ \begin{equation} \frac{\partial\left ( \rho \vec{\boldsymbol{U}}\right ) }{\partial t}+\boldsymbol{\nabla}\cdot \left ( \rho \vec{\boldsymbol{U}} \vec{\boldsymbol{U}}\right )= -\boldsymbol{\nabla}p+\boldsymbol{\nabla}\cdot\tau +\rho \vec{f} \end{equation} \] 能量方程 \[ \begin{equation} \begin{split} \frac{\partial }{\partial t}\left [ \rho \left ( e+\frac{U^{2}}{2}\right ) \right ]+\boldsymbol{\nabla}\cdot \left [ \rho \left ( e+\frac{U^{2}}{2}\right )\vec{\boldsymbol{U}} \right ]&=\rho \dot{q}+\frac{\partial }{\partial x}\left ( K\frac{\partial T}{\partial x} \right )+\frac{\partial }{\partial y}\left ( K\frac{\partial T}{\partial y} \right )+\frac{\partial }{\partial z}\left ( K\frac{\partial T}{\partial z} \right )-\frac{\partial \left ( up \right )}{\partial x}\\ &-\frac{\partial \left ( \upsilon p \right )}{\partial y}-\frac{\partial \left ( \omega p \right )}{\partial z}+\frac{\partial \left ( u\tau _{xx} \right )}{\partial x}+\frac{\partial \left ( u\tau _{yx} \right )}{\partial y}+\frac{\partial \left ( u\tau _{zx} \right )}{\partial z}+\frac{\partial \left ( \upsilon \tau _{xy} \right )}{\partial x}\\ &+\frac{\partial \left ( \upsilon \tau _{yy} \right )}{\partial y}+\frac{\partial \left ( \upsilon \tau _{zy} \right )}{\partial z}+\frac{\partial \left ( \omega \tau _{xz} \right )}{\partial x}+\frac{\partial \left ( \omega \tau _{yz} \right )}{\partial y}+\frac{\partial \left ( \omega \tau _{zz} \right )}{\partial z}+\rho \boldsymbol{f}\cdot \vec{\boldsymbol{U}} \end{split} \end{equation} \]

雷诺平均方法(RANS)

在雷诺平均法中，把瞬时速度u分解为时均速度\(\overline{u_{i}}\) 和脉动速度\({u_{i}}'\) ，则对于速度分量：

\[ \begin{equation} u_{i}=\overline{u_{i}}+{u_{i}}' \end{equation} \]

同样，对于压力和能量、相分数或其他标量也采用这种方式进行分解。

\[ \begin{align} \phi(x_i, t) = \underbrace{\overline{\phi}(x_i)}_{时间平均项} + \underbrace{\phi'(x_i, t)}_{脉动项} \end{align} \] 对于稳态问题，方程中的时间平均项可以表示为

\[ \begin{equation} \overline{\phi}(x_i) = \lim_{T \to \infty} \frac{1}{T} \int_0^T \phi(x_i, t) \mathrm{d}t \end{equation} \]

设\(\phi\)及\(\xi\)是两个瞬时值，\({\phi}'\)及\({\xi}'\)是相应的脉动值，则有以下关系成立：

\[ \begin{equation} \overline{\overline{\phi}} = \overline{\phi}, \quad \overline{\phi'} = 0\\ \overline{\overline{\phi}+{\phi'}}=\overline{\phi}, \quad \overline{\overline{\phi}\overline{\xi}}=\overline{\phi \xi} \nonumber \end{equation} \]

\[ \begin{equation} \begin{aligned} \overline{\phi \xi} &= \overline{(\overline{\phi}+\phi')(\overline{\xi}+\xi')} \\ &= \overline{(\overline{\phi}\overline{\xi} + \overline{\xi} \phi' + \overline{\phi}\xi' + \xi' \xi')} \\ &= \overline{\overline{\phi}\overline{\xi}} + \overline{\overline{\xi} \phi'} + \overline{\overline{\phi}\xi'} + \overline{\phi' \xi'} \\ &= \overline{\phi}\overline{\xi} + \overline{\phi' \xi'} \end{aligned} \nonumber \end{equation} \]

\[ \begin{equation} \overline{\alpha \phi + \beta \xi} = \alpha \overline{\phi} + \beta \overline{\xi} \quad (\alpha 和 \beta 为常数) \nonumber \end{equation} \]

\[ \begin{equation} \frac{\overline{\partial \phi}}{\partial t} = \frac{\partial \overline{\phi}}{\partial t}, \quad \frac{\overline{\partial \phi}}{\partial x_i} = \frac{\partial \overline{\phi}}{\partial x_i}\\ \frac{\overline{\partial^2 \phi} }{\partial x_{i}^2}=\frac{\partial^2 \overline{\phi} }{\partial x_{i}^2}, \quad \frac{\overline{\partial \phi'}}{\partial x_i} =0, \quad \frac{\overline{\partial^2 \phi'} }{\partial x_{i}^2}=0 \nonumber \end{equation} \]

对于连续性方程，将三个坐标方向的瞬时速度表示为时均速度\(\overline{u_{i}}\)和脉动速度\({u_{i}}'\)之和并代入，再对该式作时均运算，得到：

\[ \begin{equation} \overline{\frac{\partial\left(\overline{u}+u^\prime\right)}{\partial x}}+\overline{\frac{\partial\left(\overline{v}+v^\prime\right)}{\partial y}}+\overline{\frac{\partial\left(\overline{\omega}+\omega^\prime\right)}{\partial z}}=\frac{\partial\overline{u}}{\partial x}+\frac{\partial\overline{v}}{\partial y}+\frac{\partial\overline{\omega}}{\partial z}+\frac{\partial\overline{u^\prime}}{\partial x}+\frac{\partial\overline{v^\prime}}{\partial y}+\frac{\partial\overline{\omega^\prime}}{\partial z}=0 \end{equation} \] 显然可有

\[ \begin{equation} \frac{\partial\overline{u}}{\partial x}+\frac{\partial\overline{v}}{\partial y}+\frac{\partial\overline{\omega}}{\partial z}=0 \end{equation} \]

\[ \begin{equation} \frac{\partial\overline{u^\prime}}{\partial x}+\frac{\partial\overline{v^\prime}}{\partial y}+\frac{\partial\overline{\omega^\prime}}{\partial z}=0 \end{equation} \]

这表明湍流速度的时均值和脉动值仍然满足连续性方程。

对于动量方程，以x方向为例，作类似于上面的处理可得

\[ \begin{align} \overline{\frac{\partial\left(\overline{u}+u^\prime\right)}{\partial t}}+\overline{\frac{\partial\left(\overline{u}+u^\prime\right)^2}{\partial x}}+\overline{\frac{\partial\left(\overline{u}+u^\prime\right)\partial\left(\overline{v}+v^\prime\right)}{\partial y}}+\overline{\frac{\partial\left(\overline{u}+u^\prime\right)\partial\left(\overline{\omega}+\omega^\prime\right)}{\partial z}}=\\ -\frac{1}{\rho}\overline{\frac{\partial\left(\overline{p}+p^\prime\right)}{\partial x}}+\nu\left[\overline{\frac{\partial^2\left(\overline{u}+u^\prime\right)}{\partial x^2}}+\overline{\frac{\partial^2\left(\overline{v}+v^\prime\right)}{\partial y^2}}+\overline{\frac{\partial^2\left(\overline{\omega}+\omega^\prime\right)}{\partial z^2}}\right] \nonumber \end{align} \]

利用上面给出的关系式可得

\[ \begin{align} \frac{\partial\overline{u}}{\partial t}+\frac{\partial\left({\overline{u}}^2\right)}{\partial x}+\frac{\partial\left(\overline{uv}\right)}{\partial y}+\frac{\partial\left(\overline{u\omega}\right)}{\partial z}+\frac{\overline{\partial\left(u^\prime\right)^2}}{\partial x}+\frac{\overline{\partial\left(u^\prime v^\prime\right)}}{\partial y}+\frac{\overline{\partial\left(u^\prime\omega^\prime\right)}}{\partial z}=-\frac{1}{\rho}\frac{\partial\overline{p}}{\partial x}+\nu\left[\frac{\partial^2\overline{u}}{\partial x^2}+\frac{\partial^2\overline{v}}{\partial y^2}+\frac{\partial^2\overline{\omega}}{\partial z^2}\right] \end{align} \]

将左端脉动分量乘积的时均值放到右端粘性项中，即

\[ \begin{align} \frac{\partial\overline{u}}{\partial t}+\frac{\partial\left({\overline{u}}^2\right)}{\partial x}+\frac{\partial\left(\overline{uv}\right)}{\partial y}+\frac{\partial\left(\overline{u\omega}\right)}{\partial z}=-\frac{1}{\rho}\frac{\partial\overline{p}}{\partial x}+\frac{\partial}{\partial x}\left[\nu\frac{\partial\overline{u}}{\partial x}-\overline{\left(u^\prime\right)^2}\right]+\frac{\partial}{\partial y}\left[\nu\frac{\partial\overline{v}}{\partial y}-\overline{\left(u^\prime v^\prime\right)}\right]+\frac{\partial}{\partial z}\left[\nu\frac{\partial\overline{\omega}}{\partial z}-\overline{\left(u^\prime\omega^\prime\right)}\right] \end{align} \]

对其他两个方向也可以做类似的推导，对不可压缩流体，常写作张量符号形式，可得到下列时均形式的N-S方程，即Reynolds方程

\[ \begin{equation} \frac {\partial \overline{u}_i}{\partial t} + \frac{\partial}{\partial x_j}\left( \overline{u}_i \overline{u}_j \right) - \frac{\partial}{\partial x_j}\left[\nu\left( \frac{\partial \overline{u}_i}{\partial x_j} + \frac{\partial \overline{u}_j}{\partial x_i} \right) \underbrace{- \overline{u_i' u_j'}}_{雷诺应力} \right] = - \frac{1}{\rho}\frac{\partial \overline{p}}{\partial x_i} \end{equation} \]

上式的雷诺应力通常记为\(R_{ij}\)。写作分量形式如下

\[ \begin{equation} R_{ij} \equiv - \overline{u_i' u_j'} = -\left ( \begin{matrix} \overline{u' u'} & \overline{u' v'} & \overline{u' w'} \\ \overline{v' u'} & \overline{v' v'} & \overline{v' w'} \\ \overline{w' u'} & \overline{w' v'} & \overline{w' w'} \\ \end{matrix} \right ) \end{equation} \]

\(R_{ij}\)是对角对称的，包含六个未知量：三个正应力分量（对角）和三个切应力分量（非对角）。

要求解 RANS 方程，直观的做法是求解各应力分量的对流输运方程。对应力分量的对流输运方程进行雷诺平均，会得到更高阶的未知量 \(\overline{u_i' u_j' u_k'}\)。实际上这样下去会得到无限多的未知量，整个方程组永远无法封闭。这就需要采用一些近似的模型来封闭方程组，比如用某些平均量来近似表示雷诺应力。这些不同的近似方法则被统称为湍流模型（turbulence models）。所谓湍流模型就是把湍流的脉动值附加项与时均值联系起来的一些特定关系式。

从数学形式上看，雷诺应力很像二阶应力张量。因此 Boussinesq 仿照分子粘性应力与速度变形率的关系，提出了关于涡粘性（eddy viscosity）的 Boussinesq 假设，认为雷诺应力与平均速度变形率成线性关系：

\[ \begin{equation} -\overline{u_i'u_j'} = \nu_t \left( \frac{\partial \overline{u}_i}{\partial x_j} + \frac{\partial \overline{u}_j}{\partial x_i} - \frac{2}{3} \frac{\partial \overline{u}_k}{\partial x_k} \delta_{ij} \right) - \frac{2}{3} k \delta_{ij} \end{equation} \] 式 (14) 中的 \(k\)为湍动能（turbulence kinetic energy），定义为\(k=\frac{1}{2}(\overline{u'u'}+\overline{v'v'}+\overline{w'w'})\)，\(\delta_{ij}\)为 Kronecker delta。\(\nu_t\)被称为湍流粘度或涡粘系数，湍流粘度和一般流体粘度不同，它不是流体的固有物理属性，而和当地的湍流强度有很大关系。对于流场中不同位置的点，\(\nu_t\)的值都不尽相同，也就是说它是空间坐标的函数。湍流对流换热的研究归结为确定\(\nu_t\)，确定湍流黏性系数所需微分方程的个数成为湍流工程计算模型的名称。

基于以上假设的模型被称为涡粘模型（Eddy Viscosity Model, EVM）。Boussinesq 假设被用于Spalart-Allmaras单方程模型和双方程模型，其优点是与求解湍流粘性系数有关的计算花费时间较少，但缺点是认为湍流粘性系数是各向同性标量，对于一些复杂的流动条件下并不严格成立，所以具有其应用限制性。工业界常用的涡粘模型有一方程 Spalart-Allmaras 模型，两方程 k−ε 及 k−ω 模型等。

如果湍流场各向异性很明显，如强旋流动以及应力驱动的二次流等流动中，这时可采用雷诺应力模型（Reynolds Stress Model, RSM）在更高阶的应力上封闭 RANS 方程。具体做法通常是求解六个未知应力分量的方程及其对应的耗散率方程。雷诺应力模型可以较好的模拟具有强烈非线性特征的复杂流动，然而由于它需要求解许多额外的方程，因此计算量相比涡粘模型大得多，对计算机内存也有更高的要求，因此在工程中应用不多。常见的雷诺应力模型有 LRR （Launder-Reece-Rodi）和 SSG（Speziale-Sarkar-Gatski）模型。

零方程模型

所谓零方程模型是指确定湍流黏性系数不需要微分方程的模型，最简单的零方程模型是常系数模型，对自由剪切层流动，Prandtl提出在同一截面上\(\nu_t\)为常数。 \[ \nu_t=C\delta\left|{u_{min}-u}_{max}\right| \] 其中，\(\delta\)为剪切层厚度（对称轴到1%速度点的距离），\(u_{max}\)与\(u_{min}\)为同一截面上的最大和最小流速。

在二维坐标系中，湍流切应力表示成为 \[ \tau_t=-\rho{u_i'}{u_j'}=\rho\ l^2\left(\frac{\partial U}{\partial y}\right)^2=\rho\ l^2\left|\frac{\partial U}{\partial y}\right|\frac{\partial U}{\partial y} \] 则湍流黏性系数可表示为 \[ \nu_t=\rho\ l^2\left|\frac{\partial U}{\partial y}\right| \] 其中\(U\)为主流的时均速度，\(y\)是与主流方向相垂直的坐标，\(l\)为混合长度，是模型中需要加以确定的参数，确定合适的混合长度是零方程模型的关键。

混合长度理论适用于一些比较简单的流动，如边界层类型流动；平直通道内的流动；回流较弱接近于边界层类型流动等。但是混合长度理论在物理概念上存在一些不足如，在管道中心线处速度梯度为零，但实际\(\nu_t\)不为零；不考虑湍流的历史（上游情况）的影响；不考虑湍流强度的影响。因此需要引入偏微分方程来确定湍流黏性系数来适应这方面的需要。

一方程模型

以Spalart-Allmaras单方程模型为例，它具有良好的鲁棒性和数值收敛性，可以很好的模拟绝大部分的附着流动和薄层自由剪切流动。Spalart-Allmaras单方程模型基于线性涡粘假设，湍流粘度按以下公式计算 \[ \begin{equation} \nu_t = \tilde \nu f_{v1} \end{equation} \] 其中，\(\tilde{\nu}\)为需要通过输运方程求解的量，\(f_{v1}\)是一个关于湍流粘度比\(\chi\)的函数，定义如下 \[ \begin{equation} f_{v1} = \frac{\chi^3}{\chi^3 + C_{v1}^3} \end{equation} \] 湍流粘度比\(\chi\)定义如下： \[ \begin{equation} \chi =\tilde{\nu}/\nu \end{equation} \] 在A one-equation turbulence model for aerodynamic flows这篇原始文献中，关于\(\tilde{\nu}\)的输运方程表示为以下形式： \[ \begin{equation} \frac{\partial \tilde{\nu}}{\partial t} = M(\tilde{\nu}) + P(\tilde{\nu}) + D(\tilde{\nu}) + T \end{equation} \] 其中，\(M(\tilde{\nu})\)表示 advection+diffusion： \[ \begin{equation} M(\tilde{\nu}) = - \nabla \cdot ({\tilde{\nu} \mathbf u}) + \frac{1}{\sigma} \left[ \nabla \cdot \left[ (\nu+\tilde{\nu}) \nabla \tilde{\nu} \right] + C_{b2} (\nabla \tilde{\nu})^2 \right] \end{equation} \] \(P(\tilde{\nu})\)表示 production： \[ \begin{equation} P(\tilde{\nu}) = C_{b1} (1 - f_{t2}) \tilde{S} \tilde{\nu} \end{equation} \] \(D(\tilde{\nu})\)表示 wall destruction： \[ \begin{equation} D(\tilde{\nu}) = \left(C_{w1} f_w - \frac{C_{b1}}{\kappa^2} f_{t2}\right) \left( \frac{\tilde{\nu}}{d} \right)^2 \end{equation} \] \(T\)表示 trip function： \[ \begin{equation} T = f_{t1} \Delta \mathbf{u}^2 \end{equation} \] 于是，关于\(\tilde{\nu}\)的输运方程可以简化为： \[ \begin{equation} \frac{\partial \tilde{\nu}}{\partial t} + \underbrace{\nabla \cdot ({\tilde{\nu} \mathbf u})}_{\mathrm{advection}} = \underbrace{\frac{1}{\sigma} \left[ \nabla \cdot \left[ (\nu+\tilde{\nu}) \nabla \tilde{\nu} \right] + C_{b2} (\nabla \tilde{\nu})^2 \right]}_{\mathrm{diffusion}} + \underbrace{C_{b1} \tilde{S} \tilde{\nu}}_{\mathrm{production}} - \underbrace{\left(C_{w1} f_w \right) \left( \frac{\tilde{\nu}}{d} \right)^2}_{\mathrm{destruction}} \end{equation} \] \(\tilde{S}\)的定义如下： \[ \tilde{S} \equiv S + \frac{\tilde \nu}{\kappa^2 d^2} f_{v2} \] \(S\)的定义如下： \[ \begin{equation} S=\sqrt{2} |\mathbf S|, \quad \mathbf S = \frac{1}{2}\left( \nabla \mathbf u + (\nabla \mathbf u)^T\right) \end{equation} \] 计算\(\tilde{S}\)时用到的\(f_{v2}\)定义如下： \[ \begin{equation} f_{v2} = 1 - \frac{\chi}{1 + \chi f_{v1}} \end{equation} \] wall destruction term 中的\(f_{w}\)定义如下 \[ \begin{align*} f_{w} = g \left( \frac{1 + C_{w3}^6}{g^6 + C_{w3}^6} \right)^{1/6}, \quad g = r + C_{w2}(r^6 - r), \quad r = \min \left( \frac{\tilde \nu}{\tilde S \kappa^2 d^2}, 10 \right) \end{align*} \] 各系数的默认值如下：

\(\sigma\)	\(\kappa\)	\(C_{b1}\)	\(C_{b2}\)	\(C_{w1}\)	\(C_{w2}\)	\(C_{w3}\)	\(C_{v1}\)	\(C_{s}\)
0.66666	0.41	0.1355	0.622	\(\displaystyle\frac{C_{b1}}{\kappa^2}+\frac{1+C_{b2}}{\sigma}\)	0.3	2.0	7.1	0.3

二方程模型

k−ε模型

在一方程模型中，湍流长度的标尺\(l\)是由经验公式给出的。文献中广泛采用形如\(z=k^ml^n\)的公式来选择与湍流脉动的长度标尺有关的量。己经采用的\(z\)变量的主要形式列于表中。Rodi指出，所有这些\(z\)变量的微分方程形式均类似，但对靠近壁面地区来说\(\varepsilon\)方程计算最为方便。因而，在湍流的工程计算中，\(k-\varepsilon\)两方程模型使用最宽广。

\(z\)变量	\(\frac{k^{1/2}}{l}\)	\(\frac{k^{3/2}}{l}\)	\(\frac{k}{l}\)	\(\frac{k}{l^{2}}\)
提出者	Kolmogorov	Chou(周培源)	Rodi, Spalding	Spalding
符号	\(f\)	\(\varepsilon\)	\(kl\)	\(W\)
物理意义	涡旋频率	能量的耗散	能量与标尺之积	涡量脉动的时均方值

耗散率\(\varepsilon\)的严格定义： \[ \varepsilon=\nu\overline{\left(\frac{\partial u_i'}{\partial x_k}\right)\left(\frac{\partial u_i'}{\partial x_k}\right)} \]

耗散率\(\varepsilon\)的模拟定义： \[ \varepsilon=c_D\rho\frac{k^{3/2}}{l} \] 它可看作是单位质量流体脉动动能耗散率，\(c_D\)为经验常数。这一模拟定义来源可以如下的解：从较大的涡向较小的涡传递能量的速率对单位体积流体正比于\(\rho k\)，而反比于传递时间。传递时间与湍流长度标尺\(l\)成正比，而与脉动速度\(k^{1/2}\)成反比。

标准\(k-\varepsilon\)模型需要求解湍动能及其耗散率方程。湍动能输运方程是通过精确的方程推导得到，但耗散率方程是通过物理推理，数学上模拟相似原形方程得到的。该模型假设流动为完全湍流，分子粘性的影响可以忽略。因此，标准\(k-\varepsilon\)模型只适合完全湍流的流动过程模拟。

标准\(k-\varepsilon\)模型的湍动能\(k\)及其耗散率\(\varepsilon\)方程为如下形式 \[ \rho\frac{Dk}{Dt}=\frac{\partial}{\partial x_i}\left[\left(\mu+\frac{\mu_t}{\sigma_k}\right)\frac{\partial k}{\partial x_i}\right]+G_k+G_b-\rho\varepsilon-Y_M \]

\[ \rho\frac{D\varepsilon}{Dt}=\frac{\partial}{\partial x_i}\left[\left(\mu+\frac{\mu_t}{\sigma_\varepsilon}\right)\frac{\partial\varepsilon}{\partial x_i}\right]+C_{1\varepsilon}\frac{\varepsilon}{k}\left(G_k+C_{3\varepsilon}G_b\right)-C_{2\varepsilon}\rho\frac{\varepsilon^2}{k} \]

在上述方程中，\(G_k\)表示由于平均速度梯度引起的湍动能产生，\(G_b\)是用于浮力影响引起的湍动能产生，\(Y_M\)可压速湍流脉动膨胀对总的耗散率的影响。湍流粘性系数\(\mu_t=\rho\ C_\mu\frac{k^2}{\varepsilon}\)

在FLUENT中，作为默认值常数，\(C_{1\varepsilon}=1.44\)，\(C_{2\varepsilon}=1.92\)，\(C_\mu=0.09\)，湍动能\(k\)与耗散率\(\varepsilon\)的湍流普朗特数分别为\(\sigma_k=1.0\)，\(\sigma_\varepsilon=1.3\)可以通过调节“粘性模型”面板来调节这些常数值。

k−ω模型

湍动能\(k\)和比耗散率\(\omega\)由以下输运方程得出： \[ \frac{\partial}{\partial t}\left(\rho k\right)+\frac{\partial}{\partial x_i}\left(\rho k\ u_i\right)=\frac{\partial}{\partial x_j}\left(\Gamma_k\frac{\partial k}{\partial x_j}\right)+G_k-Y_k+S_k+G_b \]

\[ \frac{\partial}{\partial t}\left(\rho\omega\right)+\frac{\partial}{\partial x_i}\left(\rho\omega u_i\right)=\frac{\partial}{\partial x_j}\left(\Gamma_\omega\frac{\partial\omega}{\partial x_j}\right)+G_\omega-Y_\omega+S_\omega+G_\omega \]

其中\(G_k\)表示平均速度梯度产生的湍流动能；\(G_w\)表示w的生成；\(\Gamma_k\)和 \(\Gamma_\omega\)分别代表\(k\)和\(\omega\)的有效扩散系数；\(Y_k\)和\(Y_w\)表示\(k\)和\(\omega\)在湍流作用下的耗散；\(S_k\)和\(S_w\)是用户定义的源项。

\(k-\omega\)模型的有效扩散系数由下式给出： \[ \Gamma_k=\mu+\frac{\mu_t}{\sigma_k} \]

\[ \Gamma_\omega=\mu+\frac{\mu_t}{\sigma_\omega} \]

\(\sigma_k\)和\(\sigma_\omega\)分别代表\(k\)和\(\omega\)的湍流普朗特数，结合\(k\)和\(\omega\)计算湍流粘度\(\mu_t\): \[ \mu_t=\alpha^\ast\frac{\rho k}{\omega} \] 系数\(\alpha^\ast\)抑制湍流粘度，所以要对低雷诺数修正，它是由下式给出： \[ \alpha^\ast=\alpha_\infty^\ast\left(\frac{\alpha_0^\ast+{Re}_t{/}R_k}{1+{Re}_t{/}R_k}\right) \] 其中， \[ {Re}_t{=}\frac{\rho k}{\mu\omega}, \quad R_k=6, \quad \alpha_0^\ast=\frac{\beta_i}{3}, \quad \beta_i=0.072 \] 注意,在\(k-\omega\)模型的高雷诺数形式下,\(\alpha^\ast=\alpha_\infty^\ast=1\)。

\(G_k\)表示平均速度梯度产生的湍流动能。由\(k\)的输运方程可知，这一项可以定义为： \[ G_k=-\rho\overline{u_i' u_j'}\frac{\partial u_j}{\partial x_i} \] 用符合Boussinesq假设的方法计算\(G_k\)， \[ G_k=\mu_tS^2 \] 其中\(S\)是平均应变率张量的模量，定义方式与\(k-\varepsilon\)模型相同。

\(\omega\)的产生由下式给出： \[ G_\omega=\alpha\frac{\omega}{k}G_k \] 系数\(\alpha\)是 \[ \alpha=\frac{\alpha_\infty^\ast}{\alpha^\ast}\left(\frac{\alpha_0^\ast+{Re}_t{/}R_k}{1+{Re}_t{/}R_k}\right) \] 其中\(R_ω = 2.95\)。注意,在\(k-\omega\)模型的高雷诺数形式下,\(\alpha^\ast=\alpha_\infty^\ast=1\)。

\(k\)的耗散为： \[ Y_k=\rho\beta^\ast f_{\beta^\ast}k\omega \] 其中， \[ f_{\beta^\ast}=\left\{\begin{matrix}\begin{matrix}1&\chi_k\le0\\\end{matrix}\\\begin{matrix}\frac{1+680\chi_k^2}{1+400\chi_k^2}&\chi_k>0\\\end{matrix}\\\end{matrix}\right. \]

\[ \chi_k=\frac{1}{\omega^3}\frac{\partial k}{\partial x_j}\frac{\partial\omega}{\partial x_j}, \quad \beta^\ast=\beta_i^\ast\left[1+\zeta^\ast F\left(M_t\right)\right] \]

\[ \beta_i^\ast=\beta_\infty^\ast\left(\frac{4/15+\left({Re}_t{/}R_\beta\right)^4}{1+\left({Re}_t{/}R_\beta\right)^4}\right), \quad \zeta^\ast=1.5, \quad R_\beta=8, \quad \beta_\infty^\ast=0.09 \]

\(\omega\)的耗散由下式给出： \[ Y_k=\rho\beta f_{\beta}\omega^2 \] 其中， \[ f_\beta=\frac{1+70\chi_\omega}{1+80\chi_\omega}, \quad \chi_\omega=\left|\frac{\Omega_{ij}\Omega_{jk}S_{ki}}{\left(\beta_\infty^\ast\omega\right)^3}\right|, \quad \Omega_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}-\frac{\partial u_j}{\partial x_i}\right) \]

\[ \beta=\beta_i\left[1-\frac{\beta_i^\ast}{\beta_i}\zeta^\ast F\left(M_t\right)\right] \]

可压缩函数\(F\left(M_t\right)\)由下式给出： \[ F\left(M_t\right)=\left\{\begin{matrix}0&M_t\le M_{t0}\\M_t^2-M_{t0}^2&M_t>M_{t0}\\\end{matrix}\right. \] 其中， \[ M_t^2=\frac{2k}{\alpha^2}, \quad M_{t0}=0.25, \quad \alpha=\sqrt{\gamma RT} \] 注意,在\(k-\omega\)模型的高雷诺数形式下，\(\beta_i^\ast=\beta_\infty^\ast\)，不可压缩形式\(\beta^\ast=\beta_i^\ast\)。

压缩效应在非常有限的自由剪切流实验中进行了校准，不推荐普遍使用。默认情况下是禁用的。

模型常数如下： \[ \alpha_\infty^\ast=1, \quad \alpha_\infty=0.52, \quad \alpha_0=\frac{1}{9}, \quad \beta_\infty^\ast=0.09, \quad \beta_i=0.072, \quad R_\beta=8 \]

\[ R_k=6, \quad R_\omega=2.95, \quad \zeta^\ast=1.5, \quad M_{t0}=0.25, \quad \sigma_k=2.0, \quad \sigma_\omega=2.0 \]