7.2 一阶封闭模型

7.2 一阶封闭模型

一阶封闭模型是近似计算 Reynolds／Favre 平均 Navier–Stokes 方程中 Reynolds 应力的最简便方法。它们基于 7.1.5、7.1.6 节分别讨论的 Boussinesq 模型或非线性涡黏性模型。因此，相应湍流模型的任务是计算涡黏度 $ \mu_{T} $。

在众多一阶封闭模型中，我们选取三种应用广泛、代表当前发展水平的方法。这三种模型均易于在结构网格和非结构网格上实现。首先讨论 Spalart 和 Allmaras 的单方程模型，然后介绍著名的 $ K-\varepsilon $ 双方程模型，最后讨论 Menter 提出的 $ K-\omega $ 剪切应力输运（SST）双方程模型。这些湍流模型在不同算例中的详细比较见文献 [46]。

下面的讨论中，密度和速度分量均应理解为 Reynolds／Favre 平均量；为简化书写，省略相应的平均记号。

7.2.1 Spalart–Allmaras 单方程模型

Spalart–Allmaras 单方程湍流模型 $ [47] $ 采用涡黏性变量 $ \tilde{\nu} $ 的输运方程。该模型依据经验、量纲分析和 Galilean 不变性建立，并利用二维混合层、尾流和平板边界层的结果进行标定。Spalart–Allmaras 模型对具有逆压梯度的湍流也能作出较为准确的预测，而且可以在用户指定的位置实现由层流向湍流的平滑过渡。该模型具有若干有利的数值特性。它是“局部的”，即某一点的方程不依赖于其他点的解，因此易于在多块结构网格或非结构网格上实现。此外，它具有较好的稳健性，向定常解收敛较快，对近壁区网格分辨率的要求也较适中。

微分形式

Spalart–Allmaras 湍流模型可用张量记号写为 $ [47] $

$$ \begin{aligned}\frac{\partial\tilde{\nu}}{\partial t}+\frac{\partial}{\partial x_{j}}(\tilde{\nu} v_{j})=&C_{b1}(1-f_{t2})\tilde{S}\tilde{\nu}\\&+\frac{1}{\sigma}\left\{\frac{\partial}{\partial x_{j}}\bigg[(\nu_{\mathrm{L}}+\tilde{\nu})\frac{\partial\tilde{\nu}}{\partial x_{j}}\bigg]+C_{b2}\frac{\partial\tilde{\nu}}{\partial x_{j}}\frac{\partial\tilde{\nu}}{\partial x_{j}}\right\}\\&-\left[C_{w1}f_{w}-\frac{C_{b1}}{\kappa^{2}}f_{t2}\right]\bigg(\frac{\tilde{\nu}}{d}\bigg)^{2}+f_{t1}\|\Delta\vec{v}\|_{2}^{2}.\end{aligned}\tag{7.36} $$

右端各项依次表示涡黏度生成、守恒扩散、非守恒扩散、近壁湍流破坏、转捩对生成项的抑制，以及转捩湍流源。此外，$ \nu_{L} = \mu_{L} / \rho $ 表示层流运动黏度，d 为到最近壁面的距离（壁面距离的计算见文献 [48]）。式（7.28）、（7.29）中的湍流涡黏度由下式得到：

$$ \mu_{\mathrm{T}}=f_{v1}\rho\tilde{\nu}.\tag{7.37} $$

生成项采用以下公式计算：

$$ \begin{aligned}\tilde{S}&=f_{v3}S+\frac{\tilde{\nu}}{\kappa^{2}d^{2}}f_{v2},\\f_{v1}&=\frac{\chi^{3}}{\chi^{3}+C_{v1}^{3}},\quad f_{v2}=\left(1+\frac{\chi}{C_{v2}}\right)^{-3},\\f_{v3}&=\frac{(1+\chi f_{v1})(1-f_{v2})}{\max(\chi,0.001)},\quad\chi=\frac{\tilde{\nu}}{\nu_{\mathrm{L}}}.\end{aligned}\tag{7.38} $$

其中，S 表示平均旋转率的大小，即

$$ S=\sqrt{2\Omega_{i j}\Omega_{i j}}, $$

其中，$ \Omega_{ij} $ 由式（7.4）给出。注意，$ \tilde{S} $ 与文献 [47] 中的原始定义不同。这一修正由 Spalart 提出，以防止 $ \tilde{S} $ 降为零（参见文献 [49]，第 155 页）。

控制涡黏度破坏的各项为

$$ \begin{aligned}f_{w}&=g\left(\frac{1+C_{w3}^{6}}{g^{6}+C_{w3}^{6}}\right)^{1/6},\\g&=r+C_{w2}(r^{6}-r),\quad r=\frac{\tilde{\nu}}{\tilde{S}\kappa^{2}d^{2}}.\end{aligned}\tag{7.39} $$

用于模拟层流—湍流转捩的函数为

$$ \begin{array}{l}f_{t1}=g_{t}C_{t1}\exp\left(-C_{t2}\frac{\omega_{t}^{2}}{\Delta U^{2}}\left(d^{2}+g_{t}^{2}d_{t}^{2}\right)\right),\\f_{t2}=C_{t3}\exp\left(-C_{t4}\chi^{2}\right),\quad g_{t}=\min\left[0.1,\|\Delta\vec{v}\|_{2}/(\omega_{t}\Delta x_{t})\right],\end{array}\tag{7.40} $$

其中，$ \omega_{t} $ 表示壁面触发点处的涡量（位置必须由用户指定），$ \|\Delta\vec{v}\|_{2} $ 表示触发点速度与当前场点速度之差的二范数，$ d_{t} $ 为到最近触发点的距离，$ \Delta x_{t} $ 表示触发点处沿壁面的网格间距。

最后，式（7.36）—（7.40）中的各常数定义为

$$ \begin{aligned}&C_{b1}=0.1355,\quad C_{b2}=0.622,\\&C_{v1}=7.1,\quad C_{v2}=5,\quad\sigma=2/3,\quad\kappa=0.41,\\&C_{w1}=C_{b1}/\kappa^{2}+(1+C_{b2})/\sigma,\quad C_{w2}=0.3,\quad c_{w3}=2,\\&C_{t1}=1,\quad C_{t2}=2,\quad C_{t3}=1.3,\quad C_{t4}=0.5.\\ \end{aligned}\tag{7.41} $$

为了计及非平衡效应对生成项的影响，有人提出将 $ C_{b1} $ 表示为应变率的函数 [50]。

正如文献 [47] 所指出，将式（7.36）中的非守恒扩散项，即

$$ \frac{1}{\sigma}\left\{\frac{\partial}{\partial x_{j}}\Biggl[(\nu_{\mathrm{L}}+\tilde{\nu})\frac{\partial\tilde{\nu}}{\partial x_{j}}\Biggr]+C_{b2}\frac{\partial\tilde{\nu}}{\partial x_{j}}\frac{\partial\tilde{\nu}}{\partial x_{j}}\right\} $$

替换为下列表达式较为方便：

$$ \frac{1+C_{b2}}{\sigma}\frac{\partial}{\partial x_{j}}\bigg[(\nu_{\mathrm{L}}+\tilde{\nu})\frac{\partial\tilde{\nu}}{\partial x_{j}}\bigg]-\frac{C_{b2}}{\sigma}(\nu_{\mathrm{L}}+\tilde{\nu})\nabla^{2}\tilde{\nu}.\tag{7.42} $$

这样就避免了项 $ (\partial\tilde{\nu}/\partial x_{j})^{2} $ 离散时的困难。

积分形式

将 Spalart–Allmaras 湍流模型（式（7.36））转换到有限体积框架后，得到

$$ \frac{\partial}{\partial t}\int_{\Omega}\tilde{\nu}\mathrm{d}\Omega+\oint_{\partial\Omega}\left(F_{\mathrm{c},T}-F_{\mathrm{v},T}\right)\mathrm{d}S=\int_{\Omega}Q_{T}\mathrm{d}\Omega,\tag{7.43} $$

其中，$ \Omega $ 表示控制体积，$ \partial\Omega $ 为其表面，dS 为 $ \Omega $ 上的面元。对流通量定义为

$$ F_{\mathrm{c},T}=\tilde{\nu}V\tag{7.44} $$

其中，V 为逆变速度（见式（2.22））。通常采用一阶迎风格式离散对流通量。黏性通量为

$$ F_{\mathrm{v},T}=n_{x}\tau_{x x}^{\mathrm{T}}+n_{ y}\tau_{ y y}^{\mathrm{T}}+n_{z}\tau_{z z}^{\mathrm{T}},\tag{7.45} $$

其中，$ n_{x} $、$ n_{ y} $、$ n_{z} $ 为单位法向量的分量。法向黏性应力为

$$ \begin{array}{c}{\tau_{x x}^{\mathrm{T}}=\frac{1}{\sigma}\left(\nu_{\mathrm{L}}+\tilde{\nu}\right)\frac{\partial\tilde{\nu}}{\partial x},\quad\tau_{ y y}^{\mathrm{T}}=\frac{1}{\sigma}\left(\nu_{\mathrm{L}}+\tilde{\nu}\right)\frac{\partial\tilde{\nu}}{\partial y},}\\ {\tau_{z z}^{\mathrm{T}}=\frac{1}{\sigma}\left(\nu_{\mathrm{L}}+\tilde{\nu}\right)\frac{\partial\tilde{\nu}}{\partial z}.}\end{array}\tag{7.46} $$

最后，式（7.43）中的源项为

$$ \begin{aligned}Q_{T}&=C_{b1}(1-f_{t2})\tilde{S}\tilde{\nu}+\frac{C_{b2}}{\sigma}\Bigg[\Bigg(\frac{\partial\tilde{\nu}}{\partial x}\Bigg)^{2}+\Bigg(\frac{\partial\tilde{\nu}}{\partial y}\Bigg)^{2}+\Bigg(\frac{\partial\tilde{\nu}}{\partial z}\Bigg)^{2}\Bigg]\\&\quad-\Bigg[C_{w1}f_{w}-\frac{C_{b1}}{\kappa^{2}}f_{t2}\Bigg]\Bigg(\frac{\tilde{\nu}}{d}\Bigg)^{2}+f_{t1}\|\Delta\vec{v}\|_{2}^{2}.\end{aligned}\tag{7.47} $$

也可以采用式（7.42）所建议的形式表示非守恒扩散。模型常数见式（7.41）。

初始条件与边界条件

$ \tilde{\nu} $ 的初始值通常取为 $ \tilde{\nu}=0.1\;\nu_{L} $，入口边界也指定相同的值。在出口边界，直接从计算域内部外推 $ \tilde{\nu} $。在固体壁面，适宜设置 $ \tilde{\nu}=0 $，从而有 $ \mu_{T}=0 $。

7.2.2 $ K - \varepsilon $ 双方程模型

$ K-\varepsilon $ 湍流模型可能是应用最广泛的双方程涡黏性模型。它通过求解湍动能 K 和湍流耗散率 $ \varepsilon $ 的方程构建。$ K-\varepsilon $ 模型的历史可以追溯到 Chou 的研究 [51]。20 世纪 70 年代，人们提出了该模型的多种形式，其中最重要的贡献来自 Jones 和 Launder [52,53]、Launder 和 Sharma [54]，以及 Launder 和 Spalding [55]。

为了使 $ K - \varepsilon $ 湍流模型在穿过黏性底层直到壁面的区域内仍然有效，必须加入所谓的阻尼函数。阻尼函数旨在保证 K 和 $ \varepsilon $ 在壁面附近具有正确的极限行为，即

$$ K\sim y^{2}\qquad\text{且}\quad\frac{\varepsilon}{K}\sim\frac{2\nu}{ y^{2}}\qquad\text{当}\quad y\rightarrow0,\tag{7.48} $$

其中，$ y $ 表示垂直于壁面的坐标。此外，可以证明 Reynolds 剪切应力具有如下行为（例如见 [23，第 138—139 页]）：

$$ \tau_{ij}^{R}\sim y^{3}\qquad\text{当}\quad y\rightarrow0,\ i\neq j.\tag{7.49} $$

包含阻尼函数的 $ K-\varepsilon $ 模型也称为低雷诺数模型。应用最广泛的阻尼函数形式由 Jones 和 Launder [52]、Launder 和 Sharma [54]、Lam 和 Bremhorst [56]，以及 Chien [57] 提出。七种不同低雷诺数 $ K-\varepsilon $ 模型的比较见文献 [58]。

在数值求解上，$ K-\varepsilon $ 湍流模型比前述 Spalart–Allmaras 模型（7.2.1 节）更困难。尤其是阻尼函数会使湍流方程具有刚性源项。这一点，加上近壁区为解析黏性底层而必须采用的高网格分辨率，要求至少使用点隐式时间推进格式，最好采用全隐式格式。文献 [59] 对 $ K-\varepsilon $ 方程的显式时间离散给出了有用建议。$ K-\varepsilon $ 模型在结构网格和非结构网格上的实现示例可参见文献 [60–74]。最后需要注意，对于逆压梯度流动，$ K-\varepsilon $ 模型的精度会下降 [23,58]。文献 [75] 研究了模型参数如何影响不同压力梯度条件下的精度。

微分形式

低雷诺数 $ K - \varepsilon $ 模型可写为

$$ \begin{gathered}\frac{\partial\rho K}{\partial t}+\frac{\partial}{\partial x_{j}}(\rho v_{j}K)=\frac{\partial}{\partial x_{j}}\left[\left(\mu_{\mathrm{L}}+\frac{\mu_{\mathrm{T}}}{\sigma_{K}}\right)\frac{\partial K}{\partial x_{j}}\right]+\tau_{i j}^{F}S_{i j}-\rho\varepsilon\\\begin{aligned}\frac{\partial\rho\varepsilon^{*}}{\partial t}+\frac{\partial}{\partial x_{j}}(\rho v_{j}\varepsilon^{*})=&\frac{\partial}{\partial x_{j}}\left[\left(\mu_{\mathrm{L}}+\frac{\mu_{\mathrm{T}}}{\sigma_{\varepsilon}}\right)\frac{\partial\varepsilon^{*}}{\partial x_{j}}\right]+C_{\varepsilon1}f_{\varepsilon1}\frac{\varepsilon^{*}}{K}\tau_{ij}^{F}S_{ij}\\&-C_{\varepsilon2}f_{\varepsilon2}\rho\frac{(\varepsilon^{*})^{2}}{K}+\phi_{\varepsilon}.\end{aligned}\end{gathered}\tag{7.50} $$

右端各项分别表示守恒扩散、涡黏度生成和耗散。此外，$ \phi_{\varepsilon} $ 表示所谓的显式壁面项。Favre 平均湍流应力 $ \tau_{ij}^{F} $ 由式（7.25）给出，应变率张量 $ S_{ij} $ 由式（7.3）得到。式（7.28）、（7.29）中的湍流涡黏度由下式确定：

$$ \mu_{\mathrm{T}}=C_{\mu}f_{\mu}\;\rho\;\frac{K^{2}}{\varepsilon^{*}}.\tag{7.51} $$

计算式（7.24）或（7.25）中的涡黏度时，也使用湍动能。量 $ \varepsilon^{*} $ 与湍流耗散率 $ \varepsilon $ 的关系为

$$ \varepsilon=\varepsilon_{w}+\varepsilon^{*}.\tag{7.52} $$

项 $ \varepsilon_{w} $ 是壁面处的耗散率。式（7.52）的定义大大简化了壁面边界条件的施加（见下文）。

不同 $ K - \varepsilon $ 模型采用的常数、近壁阻尼函数及壁面项各不相同。这里选用 Launder–Sharma 模型，因为它在广泛的应用中表现良好 [58]。对于 Launder–Sharma 模型，常数和湍流 Prandtl 数为 [54]

$$ \begin{aligned}&C_{\mu}=0.09,\quad C_{\varepsilon1}=1.44,\quad C_{\varepsilon2}=1.92,\\&\sigma_{K}=1.0,\quad\sigma_{\varepsilon}=1.3,\quad Pr_{T}=0.9.\\ \end{aligned}\tag{7.53} $$

此外，近壁阻尼函数为

$$ \begin{array}{l}f_{\mu}=\exp\left(\frac{-3.4}{(1+0.02\mathrm{Re}_{\mathrm{T}})^{2}}\right)\\f_{\varepsilon1}=1\\f_{\varepsilon2}=1-0.3\exp\left(\mathrm{Re}_{\mathrm{T}}^{2}\right)\end{array}\tag{7.54} $$

校注：原书最后一个阻尼函数的指数写为正的 $\mathrm{Re}_T^2$。该写法在雷诺数增大时会使阻尼函数成为很大的负值；用于模型实现前应核对指数的负号。译文保留原式。

其中，$ Re_{\mathrm{T}} = \rho K^{2} / (\varepsilon^{*}\mu_{\mathrm{L}}) $ 为湍流雷诺数。最后，显式壁面项 $ \phi_{\varepsilon} $ 及量 $ \varepsilon_{w} $ 定义为

$$ \phi_{\varepsilon}=2\mu_{\mathrm{T}}\frac{\mu_{\mathrm{L}}}{\rho}\left(\frac{\partial^{2} v_{s}}{\partial y_{n}^{2}}\right)^{2}\quad\text{且}\varepsilon_{w}=\frac{2\mu_{\mathrm{L}}}{\rho}\left(\frac{\partial\sqrt{K}}{\partial y_{n}}\right)^{2},\tag{7.55} $$

其中，$ v_{s} $ 表示平行于壁面的速度，$ y_{n} $ 为垂直于壁面的坐标。为避免显式确定壁面距离和方向，通常采用以下笛卡尔张量形式计算壁面项和 $ \varepsilon_{w} $ [66,76]：

$$ \phi_{\varepsilon}=2\mu_{\mathrm{T}}\frac{\mu_{\mathrm{L}}}{\rho}\left(\frac{\partial^{2} v_{i}}{\partial x_{j}\partial x_{k}}\right)^{2}\quad\text{且}\varepsilon_{w}=\frac{2\mu_{\mathrm{L}}}{\rho}\left(\frac{\partial\sqrt{K}}{\partial x_{j}}\right)^{2}\tag{7.56} $$

以替代式（7.55）。

积分形式

对于具有面元 dS 的控制体 $ \Omega $，低雷诺数 $ K - \varepsilon $ 湍流模型的非定常积分形式为

$$ \frac{\partial}{\partial t}\int_{\Omega}\vec{W}_{T}\mathrm{d}\Omega+\oint_{\partial\Omega}(\vec{F}_{\mathrm{c},T}-\vec{F}_{\mathrm{v},T})\mathrm{d}S=\int_{\Omega}\vec{Q}_{T}\mathrm{d}\Omega.\tag{7.57} $$

守恒变量矢量为

$$ \vec{W}_{T}=\begin{bmatrix}\rho K\\ \rho\varepsilon^{*}\end{bmatrix}.\tag{7.58} $$

对流通量矢量定义为

$$ \vec{F}_{c,T}=\begin{bmatrix}\rho K V\\ \rho\varepsilon^{*}V\end{bmatrix},\tag{7.59} $$

其中，V 表示逆变速度（见式（2.22））。黏性通量矢量为

$$ \vec{F}_{\mathrm{v},T}=\begin{bmatrix}n_{x}\tau_{xx}^{K}+n_{ y}\tau_{ y y}^{K}+n_{z}\tau_{zz}^{K}\\ n_{x}\tau_{xx}^{\varepsilon}+n_{ y}\tau_{ y y}^{\varepsilon}+n_{z}\tau_{zz}^{\varepsilon}\end{bmatrix}\tag{7.60} $$

其中，法向湍流黏性应力为

$$ \begin{gathered}\tau_{xx}^{K}=\left(\mu_{\mathrm{L}}+\frac{\mu_{\mathrm{T}}}{\sigma_{K}}\right)\frac{\partial K}{\partial x},\cdots\\\tau_{x x}^{\varepsilon}=\left(\mu_{\mathrm{L}}+\frac{\mu_{\mathrm{T}}}{\sigma_{\varepsilon}}\right)\frac{\partial\varepsilon^{*}}{\partial x},\cdots\end{gathered}\tag{7.61} $$

式（7.60）中，$ n_{x} $、$ n_{ y} $、$ n_{z} $ 表示表面 $ \partial\Omega $ 向外的单位法向量的各分量。源项由下式计算：

$$ \vec{Q}_{T}=\left[\begin{array}{c}{P-\rho\varepsilon}\\ {\left(C_{\varepsilon1}f_{\varepsilon1}P-C_{\varepsilon2}f_{\varepsilon2}\rho\varepsilon^{*}\right)\frac{\varepsilon^{*}}{K}+\phi_{\varepsilon}}\end{array}\right],\tag{7.62} $$

其中，P 表示湍动能生成项，定义为

$$ \begin{aligned}P=&\tau_{xx}^{F}\frac{\partial u}{\partial x}+\tau_{ y y}^{F}\frac{\partial v}{\partial y}+\tau_{zz}^{F}\frac{\partial w}{\partial z}\\&+\tau_{x y}^{F}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)+\tau_{xz}^{F}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)+\tau_{ y z}^{F}\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)\end{aligned}\tag{7.63} $$

其中，Favre 平均湍流应力 $ \tau_{ij}^{F} $ 由式（7.25）给出。Launder–Sharma 模型的常数、近壁阻尼函数和壁面项由式（7.52）—（7.56）确定。湍流涡黏度 $ \mu_{T} $ 由式（7.51）得到。

初始条件与边界条件

最简单的方法是用自由来流值初始化 K 和 $ \varepsilon^{*} $。更好的方法是在固体壁面附近指定 K 和 $ \varepsilon^{*} $ 的剖面，这些剖面可以通过与湍流平板边界层的类比得到 [60]。不过，这需要知道壁面距离；对于非结构网格等情况，该信息可能不易获得。

如果采用式（7.52）的变换，固体壁面处适宜的边界条件为 K = 0 和 $ \varepsilon^{*} = 0 $，同时意味着壁面处有 $ \mu_{T} = 0 $。在入口边界，K 和 $ \varepsilon^{*} $ 可由湍流强度及长度尺度的关系计算，即

$$ (T_{u})_{\infty}=\frac{\sqrt{\frac{2}{3}K_{\infty}}}{\|\vec{v}_{\infty}\|_{2}},\quad(l_{\mathrm{T}})_{\infty}=\frac{C_{\mu}K_{\infty}^{3/2}}{\varepsilon_{\infty}^{*}},\tag{7.64} $$

这里假设 $ \varepsilon_{\infty}^{*}=\varepsilon_{\infty} $。在叶轮机械中，$ (l_{\mathrm{T}})_{\infty} $ 取平均径向叶片间距的 $ 10^{-3} $ 至 $ 10^{-2} $ 倍 [77]。在出口边界，从内部外推 K 和 $ \varepsilon^{*} $ 的值。

壁面函数

如前所述，低雷诺数模型要求壁面附近的网格非常细。标准条件是第一个节点（或单元形心）距壁面的无量纲距离为 $ y^{+} \leq 1 $。为了降低湍流方程的刚性，并减少网格点或单元数，常采用满足 $ 10 \leq y^{+} \leq 100 $ 的较粗网格。此时，$ K - \varepsilon $ 模型（式（7.50）或（7.57））不再使用阻尼函数（$ f_{\mu} = f_{\varepsilon1} = f_{\varepsilon2} = 1 $；$ \varepsilon_{w} = 0 $）和壁面项（$ \phi_{\varepsilon} = 0 $），即采用高雷诺数湍流模型。显然，第一节点（单元形心）与壁面之间的区域必须通过壁面函数衔接。壁面函数给出紧邻壁面节点（单元形心）处的 K 和 $ \varepsilon^{*} $ 值；在壁面本身及第一层节点（单元）上不求解湍流方程。壁面函数有多种形式，通常基于对数壁面律。Spalding 函数 [78] 就是一个例子，它对黏性底层、过渡区以及对数层均进行建模。高雷诺数模型的实现可参见文献 [79–85] 或 [67]。

只要网格不太粗，采用壁面函数就能对附着边界层给出较为准确的结果，也允许使用纯显式时间推进格式。然而，对于分离流动，壁面函数的适用性存在很大疑问。

7.2.3 Menter 的 SST 双方程模型

Menter 的 $ K-\omega $ SST 湍流模型 [86,87] 将 Wilcox 的 $ K-\omega $ 模型 [23,88] 与高雷诺数 $ K-\varepsilon $ 模型相结合，后者被转换为 $ K-\omega $ 形式。SST 模型旨在结合两者的优点。因此，边界层底层采用 $ K-\omega $ 方法，因为 $ K-\omega $ 模型不需要阻尼函数，在精度相近的情况下，其数值稳定性显著高于 $ K-\varepsilon $ 模型。此外，边界层的对数区也采用 $ K-\omega $ 模型，因为在逆压梯度流动和可压缩流动中，它优于 $ K-\varepsilon $ 方法。另一方面，边界层尾迹区采用 $ K-\varepsilon $ 模型，因为 $ K-\omega $ 模型对自由来流 $ \omega $ 的值高度敏感 [89]。自由剪切层也采用 $ K-\varepsilon $ 方法，因为它在尾流、射流和混合层的精度之间提供了较好的折中。

SST 湍流模型的一个显著特点是修正后的湍流涡黏度函数，其目的是提高对强逆压梯度流动及压力诱导边界层分离的预测精度。这一修正计及湍流剪切应力的输运，依据是 Bradshaw 关于主要剪切应力与湍动能成正比的观察。

SST 模型的一个不足是必须显式知道到最近壁面的距离，因此在多块结构网格或非结构网格上需要专门处理。壁面距离的计算可参见文献 $ [48] $；SST 湍流模型的应用示例见文献 $ [90-97] $。

微分形式

湍动能与湍流比耗散率的输运方程，其微分形式为 $ [86] $

$$ \begin{aligned}\frac{\partial\rho K}{\partial t}+\frac{\partial}{\partial x_{j}}(\rho v_{j} K)&=\frac{\partial}{\partial x_{j}}\left[(\mu_{\mathrm{L}}+\sigma_{K}\mu_{\mathrm{T}})\frac{\partial K}{\partial x_{j}}\right]+\tau_{ij}^{F} S_{ij}-\beta^{*}\rho\omega K,\\\frac{\partial\rho\omega}{\partial t}+\frac{\partial}{\partial x_{j}}(\rho v_{j}\omega)&=\frac{\partial}{\partial x_{j}}\left[(\mu_{\mathrm{L}}+\sigma_\omega\mu_{\mathrm{T}})\frac{\partial\omega}{\partial x_{j}}\right]+\frac{C_\omega\rho}{\mu_{\mathrm{T}}}\tau_{ij}^{F} S_{ij}\\&\quad-\beta\rho\omega^{2}+2(1-f_{1})\frac{\rho\sigma_{\omega2}}{\omega}\frac{\partial K}{\partial x_{j}}\frac{\partial\omega}{\partial x_{j}}.\end{aligned}\tag{7.65} $$

式（7.65）右端各项依次表示守恒扩散、涡黏度生成和耗散。此外，$ \omega $ 方程的最后一项表示交叉扩散。Favre 平均湍流应力 $ \tau_{ij}^{F} $ 由式（7.25）给出，应变率张量 $ S_{ij} $ 由式（7.3）得到。式（7.28）、（7.29）中的湍流涡黏度采用文献 [86] 的关系计算。

$$ \mu_{\mathrm{T}}=\frac{a_{1}\rho K}{\max(a_{1}\omega,f_{2}\|\mathrm{curl}\vec{v}\|_{2})}.\tag{7.66} $$

这一湍流黏度定义保证：在湍动能 K 的生成大于其耗散 $ \omega $（因而有 $ a_{1}\omega < \|curl\vec{v}\|_{2} $）的逆压梯度边界层中，满足 Bradshaw 假设，即 $ \tau = a_{1}\rho K $（剪切应力与湍动能成正比）。

式（7.65）中的函数 $ f_{1} $ 用于混合边界层内的 $ K-\omega $ 模型系数，以及自由剪切层和自由来流区内转换后的 $ K-\varepsilon $ 模型系数，其定义为

$$ \begin{array}{c}{f_{1}=\tanh(\mathrm{arg}_{1}^{4}),}\\ {\mathrm{arg}_{1}=\min\left[\max\left(\frac{\sqrt{K}}{0.09\omega d},\frac{500\mu_{\mathrm{L}}}{\rho\omega d^{2}}\right),\frac{4\rho\sigma_{\omega2}K}{\mathrm{CD}_{K\omega}d^{2}}\right],}\\ \end{array}\tag{7.67} $$

其中，d 表示到最近壁面的距离，$ CD_{K\omega} $ 为式（7.65）中交叉扩散项的正值部分，即

$$ \mathrm{CD}_{K\omega}=\max\left(2\frac{\rho\sigma_{\omega2}}{\omega}\frac{\partial K}{\partial x_{j}}\frac{\partial\omega}{\partial x_{j}},10^{-20}\right).\tag{7.68} $$

式（7.66）中的辅助函数 $ f_{2} $ 为

$$ \begin{aligned}f_{2}&=\tanh(\mathrm{arg}_{2}^{2}),\\\mathrm{arg}_{2}&=\max\left(\frac{2\sqrt{K}}{0.09\omega d},\frac{500\mu_{L}}{\rho\omega d^{2}}\right).\end{aligned}\tag{7.69} $$

模型常数如下：

$$ a_{1}=0.31,\quad\beta^{*}=0.09,\quad\kappa=0.41.\tag{7.70} $$

最后，SST 湍流模型的系数 $ \beta $、$ C_{\omega} $、$ \sigma_{K} $ 和 $ \sigma_{\omega} $，通过混合 $ K-\omega $ 模型的系数（记为 $ \phi_{1} $）与转换后的 $ K-\varepsilon $ 模型系数 $ (\phi_{2}) $ 得到，相应关系为

$$ \phi=f_{1}\phi_{1}+(1-f_{1})\phi_{2}.\tag{7.71} $$

内部模型 $ (K - \omega) $ 的系数为

$$ \begin{aligned}&\sigma_{K1}=0.85,\quad\sigma_{\omega1}=0.5,\quad\beta_{1}=0.075,\\&C_{\omega1}=\beta_{1}/\beta^{*}-\sigma_{\omega1}\kappa^{2}/\sqrt{\beta^{*}}=0.533.\\ \end{aligned}\tag{7.72} $$

外部模型 $ (K - \varepsilon) $ 的系数定义为

$$ \begin{aligned}&\sigma_{K2}=1.0,\quad\sigma_{\omega2}=0.856,\quad\beta_{2}=0.0828,\\&C_{\omega2}=\beta_{2}/\beta^{*}-\sigma_{\omega2}\kappa^{2}/\sqrt{\beta^{*}}=0.440.\\ \end{aligned}\tag{7.73} $$

SST 湍流模型的积分形式在原则上与 7.2.2 节的 $ K-\varepsilon $ 模型相同，因此不再重复。

边界条件

固体壁面处湍动能和比耗散率的边界条件为

$$ K=0\qquad\text{且}\quad\omega=10\frac{6\mu_{\mathrm{L}}}{\rho\beta_{1}(d_{1})^{2}}\tag{7.74} $$

其中，$ d_{1} $ 为第一个节点（单元形心）到壁面的距离。必须加密网格，以满足 $ y^{+} < 3 $。

对于入口边界，建议采用以下自由来流值：

$$ \omega_{\infty}=C_{1}\frac{\|\vec{v}_{\infty}\|_{2}}{L},\quad(\mu_{\mathrm{T}})_{\infty}=(\mu_{\mathrm{L}})_{\infty}10^{-C_{2}},\quad K_{\infty}=\frac{(\mu_{\mathrm{T}})_{\infty}}{\rho_{\infty}}\omega_{\infty},\tag{7.75} $$

其中，L 为计算域的长度，并且分别有 $ 1 \leq C_{1} \leq 10 $ 和 $ 2 \leq C_{2} \leq 5 $。在出口边界，从内部外推 K 和 $ \omega $ 的值。