第 17 章：湍流建模（Turbulence Modeling）
17.1 湍流建模（Turbulence Modeling）
本章开篇先澄清湍流作为流动状态的物理特征：层流是稳定的，而湍流是混沌的、扩散性的、能造成快速掺混的、随时间变化的，并涉及三维涡量脉动，具有宽广的时间和长度尺度（Tennekes & Lumley [1]）。湍流通常作为层流的不稳定性出现于某个临界雷诺数之上，这种不稳定性源于非线性惯性项对扰动的强烈放大。
目前被普遍接受的湍流理论是 Kolmogorov [2, 3] 提出的"能量级串"（energy cascade）概念。该理论把湍流看成由不同尺度涡旋组成的多尺度结构，每个涡旋都依其尺度携带一定能量。较大涡旋破裂并把能量传递给较小涡旋，后者又继续破裂并把能量再传递给更小涡旋。这一过程一直持续到最小可能尺度。在该尺度上分子粘性能够非常有效地把湍流动能耗散为热。
最小湍流涡旋的特征尺度由 Kolmogorov 微尺度 \(\eta\) 和微时间尺度 \(t_\eta\) 给出，其定义为 \(\eta = (\nu^3/\varepsilon)^{1/4}\)，\(t_\eta = (\nu/\varepsilon)^{1/2}\)（方程 17.1），其中 \(\nu\) 为分子运动粘度，\(\varepsilon\) 为湍流动能平均耗散率（后文给出定义）。最大涡旋的尺度（即积分尺度）则定义为与所涉及几何尺度成比例。基于能量级串概念，对湍流 Navier-Stokes 方程的直接数值求解（DNS）需要采用受 Courant 数 < 1 限制的极小时间步长与 \(\Delta x < \eta\) 的细网格，网格点数与 \(\mathrm{Re}^3\) 成正比。这一计算量极大的方法仅被少数研究者在简单算例中使用 [4–7]，目前尚无法用于工业问题。
为降低直接求解 Navier-Stokes 方程的计算成本，可采用统计方法简化湍流求解。湍流的随时间特性及其宽广的时间尺度范围提示可以对随机脉动使用统计平均技术。但时间平均会在运动方程的非线性项中产生事前未知的关联项，建模这些关联项即构成湍流建模中经典的封闭（closure）问题。统计方法中比 DNS 计算量较小的有：大涡模拟（LES）[8–11]，其大尺度湍流结构被直接求解而小尺度湍流结构被亚格子模型代替。LES 的核心思路是用空间统计滤波器（方程 17.2）过滤 Navier-Stokes 方程，决定哪些尺度保留、哪些尺度舍弃；滤波器 \(F\) 对大于滤波宽度 \(\Delta\) 的尺度上的 \(v\) 给出有效保留、对小尺度则有效为零。于是比滤波宽度大的湍流结构（包含大部分能量、且各向异性）被直接解析，而最小涡旋（位于各向同性且与流型及边界条件无关的尺度之下）则用模型代替。
目前工业湍流问题最常用的方法是求解 Reynolds 时均 Navier-Stokes 方程（RANS）[12]，其统计平均基于真实时间而非空间。该方法把流动变量分解为时均值与脉动值（图 17.1），代入原方程并作时间平均。即便名字上称为"Navier-Stokes"，分解和平均其实应用于所有控制方程。平均方程的两条主线分别是用于不可压缩湍流的标准 Reynolds 平均（推导 RANS）和用于可压缩湍流的 Favre 平均（推导 FANS）[13]。两者都对所有湍流尺度建模，因此网格限制不像 DNS/LES 那样严苛。本章后续即针对不可压缩流推导 RANS 方程。
17.2 Reynolds 平均（Reynolds Averaging）
令 \(\phi(x, t)\) 表示任一流场变量（速度、压力、能量、焓、温度、密度等）在时刻 \(t\)、位置 \(x\) 处的瞬时值。如图 17.1 所示，它被分解为均值分量 \(\overline{\phi}(x, t)\) 与脉动分量 \(\phi'(x, t)\) 之和（方程 17.3），其中均值由下述三种 Reynolds 平均方法之一 [12] 计算，其中时间平均最为常用。
17.2.1 时间平均（Time Averaging）
时间平均表示量在时间区间上的平均，适合稳态湍流（即平均意义上不随时间变化的流动）。若 \(T\) 是平均区间，则仅依赖于位置的 \(\overline{\phi}\) 由方程 17.4 计算。当 \(\phi\) 随时间的变化相比湍流脉动时间尺度缓慢时，上述方程用方程 17.5 代替。
17.2.2 空间平均（Spatial Averaging）
空间平均表示量在空间区间或体积 \(V\) 上的平均，适合均匀湍流。此时仅依赖于时间的 \(\overline{\phi}(t)\) 由方程 17.6 计算。
17.2.3 集合平均（Ensemble Averaging）
集合平均适合任意类型湍流（包括非稳态湍流），表示大量相同量在某一时刻的平均。设相同量的数量为 \(N\)，则该情况下作为空间与时间函数的 \(\overline{\phi}(x, t)\) 由方程 17.7 给出。
17.2.4 平均规则（Averaging Rules）
设 \(\phi\) 与 \(\varphi\) 为两个变量，\(\phi'\)、\(\varphi'\) 为其脉动分量，则推导 RANS 方程所需的平均规则如方程 17.8 所示，包括：脉动量的平均为零、均值的平均等于自身、均值的导数等于导数的均值、两个变量之和的平均等于均值之和、变量与常数乘积的平均等于均值与常数之积、均值与脉动量乘积的平均为零、变量乘积的平均等于均值乘积与脉动量乘积之和。
17.2.5 不可压缩 RANS 方程（Incompressible RANS Equations）
不可压缩 Reynolds 时均 Navier-Stokes 方程基于时均变量。将速度、压力与温度场分解（方程 17.9：\(v = \overline{v} + v'\)、\(p = \overline{p} + p'\)、\(T = \overline{T} + T'\)），代入方程 (3.13)、(3.39)、(3.78) 所表示的不可压缩连续性、动量、能量方程，假设牛顿流体并作时间平均，即得（方程 17.10–17.12）。
对其作 Reynolds 平均后得到方程 17.13–17.15：连续性、动量、能量方程。除方程右端因非线性项新增的脉动量乘积平均项外，上述方程与原守恒方程形式相同。这些项在动量方程中引入 6 个新未知量（张量 \(-\rho \overline{v'v'}\) 的各分量，即 Reynolds 应力张量 \(\sigma_R\)），并在能量方程中引入 3 个新未知脉动热通量 \(\dot{q}^R = -\rho c_p \overline{v'T'}\)。Reynolds 应力张量 \(\sigma_R\) 与脉动热通量矢量的展开形式由方程 17.16 给出。
由此 RANS 方程组不是封闭组：要求解它必须为 Reynolds 应力分量附加额外方程。由原守恒方程推导这些方程反而会引入更多未知量（如脉动量的三重乘积），使问题更加复杂。Reynolds 应力张量源自动量与温度方程的对流非线性项，再次说明湍流本身是非线性现象，对扰动高度敏感。因此，任何对方程的线性平均（如 Reynolds 平均方法）都不能降低问题阶数——复杂性递归：试图为三重乘积建立方程就出现四重乘积，依此类推。为克服此问题，任意湍流模型都必须在最后一步把非线性脉动应力分量仅用均值分量表达，如下文所述。
17.3 Boussinesq 假设（Boussinesq Hypothesis）
对 Reynolds 应力张量进行直接建模的基础是 Boussinesq 假设 [14–16]。与牛顿流动类比，它假设 Reynolds 应力是平均速度梯度的线性函数，即方程 17.17。对不可压缩流动，由于 \(\nabla \cdot v = 0\)，该式化简为方程 17.18。
式中（为简化记号，从此处开始省略平均量上的横线）\(k\) 为湍流动能，定义为 \(k = \frac{1}{2} v' \cdot v'\)（方程 17.19），\(\mu_t\) 为湍流涡粘度（与分子粘度类比）。它依赖流动而非流体。于是计算 Reynolds 应力分量的问题被转化为计算湍流动能与湍流粘度。对不可压缩流动，Reynolds 应力中的 \(-(2/3)\rho k I\) 项通常被合并到压力梯度项中，定义湍流压力 \(\widetilde{p} = p + (2/3)\rho k\)（方程 17.20），从而把未知量减至 \(\mu_t\) 一个，由各种湍流模型计算。
类似地，脉动热通量按 Fourier 定律类比建模（方程 17.21）：
，其中 \(k_t\) 为湍流热扩散率，后文给出计算方法。
17.4 湍流模型（Turbulence Models）
基于 Boussinesq 假设，已经发展出多种湍流模型，把湍流粘度 \(\mu_t\) 用速度尺度 \(\sqrt{k}\) 与长度尺度 \(\ell\) 表示，即 \(\mu_t = \rho \ell \sqrt{k}\)（方程 17.22）。这些模型分为四大类：
- 代数（零方程）模型
- 一方程模型
- 二方程模型
- 二阶封闭模型
已发展的模型中没有一个能普适地应用于所有流动条件；不过每一类都有各自的优势与长处。零方程模型 [17–19] 用代数方程计算 \(\mu_t\)，无需解任何微分方程。一方程模型 [20–22] 需求解一个输运微分方程以计算湍流涡粘度。二方程湍流模型 [23–35] 需要求解两个输运方程来计算 \(\mu_t\)。二阶封闭模型 [36–41] 是计算量最大的，因为需要为单个湍流通量（6 个方程）分别求解输运方程。
二方程湍流模型在工业应用模拟中使用最广：要求解两个输运方程，又能给出足够精确的预测。Jones 和 Launder [23] 的 \(k-\varepsilon\) 模型是最早的二方程模型之一，也是最常用的；Wilcox [28, 29] 的 \(k-\omega\) 模型紧随其后。两者都经历了大量改进 [26, 30]，极大扩展了其适用范围。
17.5.1 标准 \(k-\varepsilon\) 模型（Standard \(k-\varepsilon\) Model）
由 Jones 与 Launder [23] 提出的著名标准 \(k-\varepsilon\) 模型基于 Boussinesq 近似，把湍流粘度 \(\mu_t\) 与热扩散率 \(k_t\) 表示为 \(\mu_t = \rho C_\mu k^2/\varepsilon\)（方程 17.23）以及 \(k_t = c_p \mu_t/\mathrm{Pr}_t\)，其中 \(\varepsilon\) 是单位质量湍流动能的粘性耗散率（方程 17.24）。
模型中，湍流动能 \(k\) 与湍流动能耗散率 \(\varepsilon\) 通过方程 17.25 与 17.26 计算：方程 17.25 是 \(k\) 方程，\(\nabla \cdot (\mu_{\mathrm{eff},k} \nabla k)\) 为扩散项，\(S_k = P_k - \rho \varepsilon\) 为源项；方程 17.26 是 \(\varepsilon\) 方程，
。有效粘度 \(\mu_{\mathrm{eff},k} = \mu + \mu_t/\sigma_k\)， \(\mu_{\mathrm{eff},\varepsilon} = \mu + \mu_t/\sigma_\varepsilon\) （方程 17.27）。模型常数为 \(C_{\varepsilon 1} = 1.44\)、\(C_{\varepsilon 2} = 1.92\)、\(C_\mu = 0.09\)、\(\sigma_k = 1.0\)、\(\sigma_\varepsilon = 1.3\)、\(\mathrm{Pr}_t = 0.9\)。
湍流动能生成项的紧凑形式为 \(P_k = \sigma_R : \nabla v\)（方程 17.28），对不可压缩流动，将其乘以 \(\mu_t\) 即可由方程 (3.75) 得到其展开形式。在标准 \(k-\varepsilon\) 模型的推导中假设流动为完全湍流、分子粘性效应可忽略，因此它是仅对完全湍流的自由剪切流有效的高雷诺数模型，无法一直积分到壁面。
要模拟近壁面流动需要在细网格上积分这两个方程，以正确捕捉边界层内的湍流量及低雷诺数效应修正。能一直积分到壁面的湍流模型在文献中称为低雷诺数湍流模型或低雷诺数版本。多年来已提出多种所谓低雷诺数 \(k-\varepsilon\) 模型（综述见 Patel 等 [42] 与 Wilcox [29]）。其开发思路是通过阻尼函数使近壁面的湍流粘度随距壁面距离减小而趋于零；部分模型还会对湍流耗散方程中源项系数施加阻尼。所有模型共享相同的基本结构，差异仅在阻尼函数的调整与耗散方程中某些额外源项。
该规则唯一的例外是 Wilcox [28, 29] 的 \(k-\omega\) 湍流模型，它无需阻尼函数即可一直积分到壁面，但仍可作为高雷诺数模型使用。两方程模型的主要缺点之一是所谓驻点反常：在高应变率区域两方程模型倾向于高估湍流动能生成 \(P_k\)。该现象最初在驻点流动中被发现，但其实是更普遍的异常。问题根源是当中等大小的 \(k\) 受到大应变率作用时湍流动能生产过多。这可能源于汇项的低估与/或湍流粘度过高。这些思路可合并为对局部湍流时间尺度 \(t_s = k/\varepsilon\) 的约束（Medic 与 Durbin [43]）。第一步是把湍流粘度表达式改写为 \(\mu_t = \rho C_\mu k t_s\)（方程 17.29）。然后利用 \(t_s\)，\(\varepsilon\) 方程修正为方程 17.30 形式，其中源项中 \(1/t_s\) 与 \(\varepsilon/t_s\) 替换原本的 \(\varepsilon/k\) 与 \(\varepsilon^2/k\)。最后为约束 Reynolds 应力张量为正定，对 \(t_s\) 施加限制（方程 17.31）：\(t_s = \min\{k/\varepsilon, a/\sqrt{6 C_\mu S_t}\}\)，其中 \(a = 0.6\)，\(S_t = \sqrt{S_t \cdot S_t}\)，\(S_t = (1/2)(\nabla v + \nabla v^T)\)。其后果是在大应变率下 \(P_k\) 以 \(S_t\) 速率增长，而非 \(S_t^2\)。
17.5.2 \(k-\omega\) 模型（The \(k-\omega\) Model）
前文已经提到 \(k-\varepsilon\) 模型族在自由剪切流中表现良好，但倾向于在预测带逆压梯度的流动时失效。另一类模型将 \(\varepsilon\) 方程替换为 \(\omega\) 方程，其中 \(\omega\) 为单位体积、单位时间内湍流动能转化为内热的速率，能更好地预测分离流动。
该类中的第一个完整湍流模型由 Kolmogorov [27] 提出。除与 \(k-\varepsilon\) 相同的 \(k\) 方程外，Kolmogorov 还推导了第二个 \(\omega\) 方程。\(\omega\) 的倒数被用作局部湍流时间尺度 \(\sqrt{k}/\omega\)，湍流长度尺度为 \(\sqrt{k}/\omega\)。下述 \(k-\omega\) 模型由 Wilcox [29] 提出，作为 Wilcox [28] 中著名 \(k-\omega\) 模型的演进。
Wilcox [29] 的 \(k-\omega\) 模型在结构上与 \(k-\varepsilon\) 模型类似，也基于 Boussinesq 近似。两个输运方程被求解以确定两个（大）湍流尺度。特定湍流耗散率 \(\omega\) 由 \(\omega = \varepsilon/(C_\mu k)\)（方程 17.32）定义。用 \(\omega\) 方程代替 \(\varepsilon\) 方程的优点是：(i) 后者更易积分（更稳健）；(ii) 能在子层中积分而无需额外阻尼函数；(iii) 在弱逆压梯度流动中表现更好。守恒方程写作方程 17.33（\(k\) 方程）与 17.34（\(\omega\) 方程）。模型常数取 \(C_{a1} = 5/9\)、\(C_{b1} = 0.075\)、\(\beta^* = 0.09\)、\(\sigma_{k1} = 2\)、\(\sigma_{\omega 1} = 2\)、\(\mathrm{Pr}_t = 0.9\)，且有效粘度由方程 17.35 给出。
Wilcox 模型的主要缺点是对自由流 [30] 指定值的敏感性，使解对自由流 \(\omega\) 的任意指定有强依赖。这一依赖在 \(k-\varepsilon\) 模型中不存在。
17.5.3 基线（BSL）\(k-\omega\) 模型（The Baseline (BSL) \(k-\omega\) Model）
由 Menter [32] 提出的 Baseline（BSL）模型结合了 \(k-\varepsilon\) 与 \(k-\omega\) 模型，目的是利用它们各自的优势：\(k-\omega\) 模型因简单的低雷诺数公式在近壁面稳健，并能精确计算弱逆压梯度流动；\(k-\varepsilon\) 模型因对自由流值不敏感，在边界层外缘及远离壁面处表现更好。该方法的基础是把 \(k-\varepsilon\) 模型转换为 \(k-\omega\) 形式。这一转换除 \(k\) 与 \(\varepsilon\) 方程扩散系数差异造成的少量扩散项贡献外是精确的。\(k-\varepsilon\) 模型的 \(k-\omega\) 形式由方程 17.36（\(k\) 方程）与方程 17.37（\(\omega\) 方程）给出。
该形式与原始 \(k-\omega\) 模型之间的差异在于 \(\omega\) 方程中出现一项额外的交叉扩散项，以及模型常数不同：\(C_{a2} = 0.4404\)、\(C_{b2} = 0.0828\)、\(\sigma_{k2} = 1.0\)、\(\sigma_{\omega 2} = 0.856\)、\(\mathrm{Pr}_t = 0.9\)。
BSL \(k-\omega\) 模型由 \(k-\omega\)（方程 17.33、17.34）乘以混合函数 \(F_1\) 加上 \(k-\varepsilon\) 模型方程的 \(k-\omega\) 形式（方程 17.36、17.37）乘以 \((1 - F_1)\) 推导得到，从而给出 \(k\) 与 \(\omega\) 的方程（17.38、17.39，[32]）。这些方程在形式上与标准 \(k-\omega\) 模型方程非常相似，但其所有系数都依赖于混合函数 \(F_1\)，形如 \(\widetilde{\psi} = F_1 \psi_1 + (1 - F_1) \psi_2\)（方程 17.40），其中方程 17.39 中所用原始 \(k-\omega\) 模型常数为 \(C_{a1} = 0.5976\)、\(C_{b1} = 0.075\)、\(\beta^* = 0.09\)、\(\sigma_{k1} = 2\)、\(\sigma_{\omega 1} = 2\)、\(\mathrm{Pr}_t = 0.9\)。
混合函数 \(F_1\) 依赖于解变量和距最近壁面的法向距离 \(d_\perp\)，由方程 17.41 给出：\(F_1 = \tanh(\Phi_1^4)\)，其中
，以及
（方程 17.42）；湍流粘度与有效粘度由方程 17.43 给出。
BSL 模型对边界层流动的性能与 \(k-\omega\) 模型相当，对自由剪切流则几乎与 \(k-\varepsilon\) 模型相同。其稳健性接近 \(k-\omega\) 模型。
17.5.4 剪切应力输运（SST）\(k-\omega\) 模型（The Shear Stress Transport (SST) \(k-\omega\) Model）
对 BSL 模型作进一步修改得到剪切应力输运（SST）[32–35] 模型。与其他涡粘度模型相比，其在带逆压梯度流动中的性能得到改进。第一处修改与满足 Bradshaw 假设相关：Bradshaw 假设指出边界层中主剪切应力与湍流动能通过 \(\tau_{xy} = \rho a_1 k\)（方程 17.44）线性相关。
而对常规二方程湍流模型，主剪切应力可由方程 17.45 计算为
，其中 \(\Omega\) 为涡量。在带逆压梯度的流动中，湍流动能生成与耗散之比可远大于 1，从而严重违反 Bradshaw 假设。为在涡粘度模型框架内满足方程 17.44，Menter [35] 在 SST \(k-\omega\) 模型中对湍流粘度 \(\mu_t\) 按方程 17.46 进行约束： \(\mu_t = \rho a_1 k / \mathrm{Max}(a_1 \omega;\; 2 S_t F_2)\) ，其中 \(a_1 = 0.31\)，\(S_t\) 为方程 17.31 定义的应变率大小，\(F_2\) 由方程 17.47 给出：\(F_2 = \tanh(\Phi_2^2)\)，
。
为对自由剪切层保留原始涡粘度的形式，SST \(k-\omega\) 模型也采用与基线模型相同的混合函数方法。\(k\) 与 \(\omega\) 方程分别由方程 17.38 与 17.39 给出。第二处修改与 \(k\) 方程（方程 17.38）中湍流动能生成项 \(P_k\) 相关，被替换为 \(\widetilde{P}_k = \mathrm{Min}(P_k, c_1 \varepsilon)\) （方程 17.48），其中 \(\varepsilon\) 由方程 17.32 得到，混合函数 \(F_1\) 按 BSL 模型经方程 17.41 与 17.42 计算，系数由方程 17.40 计算，模型常数为 \(C_{a1} = 0.5532\)、\(C_{b1} = 0.075\)、\(\beta^* = 0.09\)、\(\sigma_{k1} = 2\)、\(\sigma_{\omega 1} = 2\)、\(c_1 = 10\)、\(C_{a2} = 0.4403\)、\(C_{b2} = 0.0828\)、\(\sigma_{k2} = 1.0\)、\(\sigma_{\omega 2} = 1.186\)、\(\mathrm{Pr}_t = 0.9\)。湍流热导率与有效湍流粘度由方程 17.49 计算。
17.6 不可压缩湍流方程汇总（Summary of Incompressible Turbulent Flow Equations）
不可压缩时均连续性方程、动量方程、能量方程、湍流动能方程、湍流耗散率方程与特定耗散率方程可分别写作方程 17.50–17.55。这些方程结构相似，可统一写成方程 (3.93) 的标准形式。因此，其离散化遵循前几章的一般流程。其离散形式的概述见下一节。
17.7.1 \(k\) 方程的离散形式（The Discretized Form of the k Equation）
湍流动能方程的最终代数形式可写为
（方程 17.56）。其中系数（方程 17.57）：
；
，其中源项 \(b_{k_C}^S\) 的取值视模型而定：\(k-\varepsilon\) 模型为 \(q_e^C V_C k - e_{model}\)，\(k-\omega\) 模型为 \(\beta^* q_C \omega_C V_C k - \omega_{models}\)；\(a_{k_C}^0 = q_C V_C/\Delta t\)；\(a_C = q_C V_C/\Delta t\)。
（SST \(k-\omega\) 模型）；\((P_k)_C V_C\)（其他情形）
。
湍流动能生成项 \(P_k\) 的离散化形式由方程 (16.27) 给出，但分子粘度被替换为湍流粘度。此外，与其他变量类似，湍流动能方程通常需要亚松弛。
17.7.2 \(\varepsilon\) 方程的离散形式（The Discretized Form of the \(\varepsilon\) Equation）
湍流耗散率方程的最终代数形式为
（方程 17.58）。其中系数（方程 17.59）：
；\(a_{\varepsilon_C}^0 = q_C V_C/\Delta t\)；\(a_C = q_C V_C/\Delta t\)；
；
。
此外，与其他变量类似，湍流耗散率方程通常需要亚松弛。
17.7.3 \(\omega\) 方程的离散形式（The Discretized Form of the \(\omega\) Equation）
特定湍流耗散率方程的最终代数形式为
（方程 17.60）。其中系数（方程 17.61）：
；
；\(a_{\omega_C}^0 = q_C V_C/\Delta t\)；\(a_C = q_C V_C/\Delta t\)；
。
\(a_{add,C}\) 与 \(b_{add,C}\) 的表达式取决于所用 \(k-\omega\) 模型的版本，对应表达式如下：
原始 \(k-\omega\) 模型（方程 17.62）：\(a_{add,C} = C_{b1} q_C \omega_C V_C \omega_C\)；\(b_{add,C} = C_{a1} (P_k)_C V_C / k_C\)。
\(k-\varepsilon\) 模型的 \(k-\omega\) 形式（方程 17.63）：
；
。
BSL 与 SST \(k-\omega\) 模型（方程 17.64）：
；
。
\((\nabla k \cdot \nabla \omega)_C\) 的离散形式按方程 17.65 计算。同样，特定湍流耗散率方程通常需要亚松弛。
17.8.1 壁面附近流动的建模（Modeling Flow Near the Wall）
当湍流流向壁面时，速度的均值与脉动分量以及随之而来的 \(k\) 都趋于零，形成大梯度。此外，远离壁面的湍流应力非常高，到近壁层减小至与粘性应力同量级的值。因此，若要解析近壁层，则需要相当数量的网格点。
低雷诺数湍流模型能够模拟壁面的阻尼效应，但代价是需要大量网格点。这是要求得到近壁区精确解所不可避免的成本。另一方面，高雷诺数湍流方法（如标准 \(k-\varepsilon\) 模型）通过使用壁面函数避免了对近壁层的解析。该方法假定从边界表面到第一个近壁节点之间的理论速度剖面并叠加于求解结果之上。与前一种方法相比，壁面函数大幅降低计算成本。该方法的主要缺点涉及这些剖面的有效性——剖面仅在近平衡边界层中已知且有效。该特殊处理的细节将在下面解释。
17.8.2 标准壁面函数（Standard Wall Functions）
壁面函数方法基于沿壁面边界层的通用流动剖面 [1]，可分为三个区域：粘性子层（\(0 < d^+ < 5\)）、缓冲层（\(5 < d^+ < 30\)）、惯性子层（\(30 < d^+ < 200\)），其规一化距离 \(d^+\) 由方程 17.66 定义：\(d^+ = d_\perp u_\tau / \nu = y^+\)，其中 \(d_\perp\) 为壁面法向距离，\(\nu\) 为运动粘度（\(=\mu/\rho\)），\(u_\tau\) 为摩擦速度，按壁面剪切应力 \(\tau_w\) 由方程 17.67 定义：\(u_\tau = \sqrt{|\tau_w|/\rho}\)。
测量与直接数值模拟都显示：粘性子层中湍流可忽略，惯性子层中粘性效应很小，而缓冲层中两种效应都重要 [44]；湍流生成最大值出现在约 \(d^+ = 12\) 处，位置略依赖雷诺数，使缓冲区的流动建模非常困难。因此，湍流模型通过将第一个内部网格点放置在粘性子层或惯性子层以避开缓冲区壁面附近的区域。把第一个网格点放在粘性子层中与低雷诺数湍流模型搭配使用，把另一个区域与高雷诺数湍流模型搭配使用。
粘性子层中适用的经验关系 [45] 为方程 17.68：\(u^+ = d^+\)、\(k^+ = 0.1 d^{+2}\)、\(\varepsilon^+ = 2 k^+/d^{+2} = 0.2\)、\(\omega^+ = 6/(C_{b1} d^{+2})\)。在有速度为 \(v_w\) 的运动壁面的一般情况下，\(u^+\)、\(k^+\)、\(\varepsilon^+\)、\(\omega^+\) 按方程 17.69 重新定义，\(|v - v_w|_{\parallel}\) 是平行于壁面的速度分量大小。把上述剖面与直接数值模拟数据比较显示：速度与耗散率在 \(d^+ = 10\) 以内保持很好的一致性，而湍流动能在 \(d^+ > 5\) 处被高估。最后在低雷诺数湍流公式中，基于 \(k-\omega\) 的模型只需满足模型渐近值的边界处理；而在基于 \(k-\varepsilon\) 的模型中，要为涡粘度方程加入一个模拟分子粘度对剪切应力直接效应的阻尼函数 [25]。在惯性子层中，速度剖面在零压梯度一维 Couette 流的假设下导出，湍流量的剖面则可针对特定湍流模型导出。对 \(k-\varepsilon\) 与 \(k-\omega\) 模型，这些剖面由方程 17.70 给出。
Von Karman 常数 \(\kappa\) 取 0.41，\(C_\mu = \beta^* = 0.09\)，\(B = 5.25\)。直接数值模拟（DNS）数据表明速度剖面非常一致，但湍流量准确性稍差。如图 17.2 所示，当求解湍流问题时，对靠近壁面的第一个内部点 \(C\) 处的控制体中守恒方程作修改。该位置的 \(d^+\) 值记为 \(d_C^+\)，首先计算以判断该点处于粘性子层还是惯性子层。\(d_C^+\) 按 \(d^+\) 定义与从壁面定律得到的 \(u_\tau\) 联合计算：先用 \(k_C^+ = k_C/u_\tau^2 = 1/\sqrt{C_\mu}\) 反解出 \(u_\tau = C_\mu^{1/4} \sqrt{k_C}\)，再代入 \(d_C^+ = d_\perp u_\tau / \nu\) 得 \(d_C^+ = C_\mu^{1/4} \sqrt{k_C} (d_\perp)_C / \nu\)（方程 17.71）。
从粘性层到惯性层的过渡假设发生在 \(d^+\) 的极限值 \(d_{\lim}^+\) 处。不同来源报告的 \(d_{\lim}^+\) 值不同，但都在 11 到 12 之间，本书采用 \(d_{\lim}^+ = 11.06\)。该极限值标记对数剖面与线性剖面的交点。若 \(d_C^+ \le d_{\lim}^+\)，则网格点位于粘性子层，否则位于惯性子层。
若第一个网格点位于粘性子层，则假设流动为层流，壁面粘度等于层流粘度 \(\mu\)，剪切应力按层流方式计算。在 \(k\) 方程中固定为零，并通过对数控制体上剪切应力为常数的假设修改 \(P_k\)，由方程 17.72 计算：
。
标准 \(k-\varepsilon\) 模型中，近壁第一个控制体形心处的湍流动能耗散率通过令层流粘度等于湍流粘度得到（方程 17.73）：\(\varepsilon_C = C_\mu \rho k_C^2/\mu\)。\(k-\omega\) 模型中，湍流频率 \(\omega_C\) 由粘性子层的解析解得到（方程 17.74）：\(\omega_C = 6 \nu/(C_{b1} (d_\perp)_C^2)\)。
若 \(d_C^+ \ge d_{\lim}^+\)，则网格点位于惯性子层，对第一个内部点 \(C\) 应用对数壁面函数。实施过程包括按方程 17.75 用对数壁面函数计算剪切应力：
，故 \(\tau_w = -\rho u_\tau / (1/\kappa \ln(d_C^+) + B) (v_C - v_w)_{\parallel}\) ，其中已用到方程 17.76： \(\tau_w = -|\tau_w| (v_C - v_w)_{\parallel}/|v_C - v_w|_{\parallel}\) 。该剪切应力可通过直接以源项形式代入（\(\tau_w S_b\)）或通过在壁面处的修改粘度 \(\mu_w\) 来求解动量方程（方程 17.77）：
，故 \(\mu_w = \rho u_\tau (d_\perp)_C / (1/\kappa \ln(d_C^+) + B)\) 。
壁面剪切应力的矢量形式如方程 17.78： \(\tau_w = -(\mu_w/(d_\perp)_C) (v_C - v_w)_{\parallel}\) 。两种情况下实现遵循第 15 章所述流程。
值得注意的是，剪切应力也可以用归一化量表示（方程 17.79）：
。
在求解湍流动能方程时，假设 \(k\) 值在整个控制体上不变（即对 \(k\) 用零梯度），在第一个内部点处的湍流动能生成项通过假设剪切应力在整个控制体上为常数并等于壁面值修改，其中速度梯度按方程 17.80 由壁面函数计算：
。因此，若第一个内部点位于惯性子层，\(k\) 方程中的生成项按方程 17.81 计算：\(P_k = |\tau_w| u_\tau / \kappa (d_\perp)_C\)。
把生成项 \(P_k\) 用归一化参数表示也是有用的。从方程 17.72 出发，可得（方程 17.82）：
。在 \(k-\varepsilon\) 模型中，\(\varepsilon\) 方程不在近壁第一个内部点处求解；其值通过要求湍流动能生成等于耗散率（方程 17.83）确定：
，\(u_\tau = C_\mu^{1/4} \sqrt{k_C}\)，故 \(\varepsilon_C = C_\mu^{3/4} k_C^{3/2}/\kappa (d_\perp)_C\) 。若 \(k-\omega\) 模型与壁面函数配合使用，则沿用相同流程；对 \(k\) 方程，由于 \(C_\mu = \beta^*\)，得到相同的修改生成项；\(\omega\) 方程不在第一个内部点处求解，其值同样按耗散等于生产给定（方程 17.84）： \(\rho \varepsilon_C = P_k = \rho C_\mu^{3/4} k_C^{3/2}/\kappa (d_\perp)_C\) ，\(\varepsilon_C/(\rho C_\mu k_C)\) 给出 \(\omega_C = \sqrt{k_C}/C_\mu^{1/4} \kappa (d_\perp)_C\) 。
17.8.3 改进壁面函数（Improved Wall Functions）
上述公式仅在局部平衡条件下有效，并在 \(\tau_w\) 为零（如分离点、再附点）时按方程 17.77 给出零粘度（因 \(u_\tau = \sqrt{\tau_w/\rho}\)）。Launder 与 Spalding [46] 对 \(k-\varepsilon\) 模型提出了对局部非平衡条件也适用的推广公式。在其工作中，\(\sqrt{k}\) 通过恒等式 \(u^* = C_\mu^{1/4} \sqrt{k}\)（方程 17.85）取代摩擦速度，被用作特征湍流速度尺度，从而壁面粘度与剪切应力变为（方程 17.86）： \(\mu_w = \rho u^* (d_\perp)_C / (1/\kappa \ln(d_C^*) + B)\) ，\(\tau_w = \rho u_\tau u^*\)， \(u_\tau = u^* |v_C - v_w|_{\parallel}/(1/\kappa \ln(d_C^*) + B)\) ，其中 \(d_C^*\) 由方程 17.87 定义：\(d_C^* = (d_\perp)_C u^*/\nu\)。方程 17.86 明确显示当 \(\tau_w = 0\) 时湍流粘度不为零。
与标准壁面函数类似，先计算 \(d_C^*\) 的值。若 \(d_C^* \le d_{\lim}^* = 11.06\)，则第一个点位于粘性子层，使用标准壁面函数所述流程。若 \(d_C^* > d_{\lim}^*\)，则第一个内部网格点位于惯性子层，壁面粘度按方程 17.86 计算。利用此粘度得到壁面剪切应力并按前述方法实施。为求 \(k_C\)，求解 \(k\) 守恒方程。标准壁面函数计算近壁控制体生成与耗散项的流程假设 \(k_C\) 与 \(\varepsilon_C\) 的值在整个控制体上不变。由于 \(P_k\) 与 \(\varepsilon\) 的值在近壁网格上变化剧烈，将其放在网格中心来离散 \(k\) 方程会带来不精确的近似。Launder 与 Spalding [46] 建议对其单元平均值作适当的近似。生产项从方程 17.72 出发，其平均值按方程 17.88 计算为
。引入方程 17.86 后，平均生产项按方程 17.89 计算：
。然后生产项的体积分由方程 17.90 给出：
。
为计算平均湍流耗散率，令耗散积分值等于生产率（方程 17.91）：
。该积分通过把剪切应力与速度梯度写为方程 17.92 形式来计算：\(|\tau_w| = \rho u_\tau^2\)，
。
由此得到方程 17.93：
。
平均湍流耗散率由方程 17.94 计算：
，其体积分由方程 17.95 给出：
。
与标准壁面函数一样，\(\varepsilon\) 方程不在近壁第一个内部点处求解，其值按方程 17.83 设置。此外，通过方程 17.32 把 \(\varepsilon\) 变换为 \(\omega\)，可对 \(k-\omega\) 模型推导类似方程。
17.8.4 可缩放壁面函数（Scalable Wall Functions）
壁面函数方法在近壁区第一个网格点位于惯性子层且归一化距离 \(d^* \ge d_{\mathrm{low}}^*\) 时最为准确，其中按数值格式不同 \(d_{\mathrm{low}}^* \approx 20\)。在边界层非常薄、不能用粗糙近壁网格分辨的情况下，这构成严重限制。[47, 48] 中提出的可缩放壁面函数方法通过对 \(u^+\) 的计算做轻微修改（按方程 17.96：\(u^+ = (1/\kappa) \ln \widetilde{d^*} + B\)）克服这一障碍，其中（方程 17.97）\(\widetilde{d^*} = \mathrm{Max}(d^*, d_{\lim}^*)\)，\(d_{\lim}^* \approx 11.06\)。按此公式，\(\widetilde{d^*}\) 的定义与网格间距无关，因它阻止第一个网格点位于粘性子层，从而对任意细化程度的网格都给出一致结果。所引入的模拟误差源于未对粘性子层建模——所有壁面函数公式都存在该误差。但应当澄清，对相对低雷诺数流动该误差可能很大。可缩放壁面函数的实现直接：上述流程不变，只需把 \(d^*\) 或 \(d^+\) 分别替换为 \(\widetilde{d^*}\) 或 \(\widetilde{d^+}\)。当边界层细节不关心时通常采用此方法。若用近壁细网格是为了检查边界层细节，则应采用下文介绍的低雷诺数模型。
17.8.5 低雷诺数模型的壁面边界条件（Wall Boundary Conditions for Low Reynolds Number Models）
对低雷诺数湍流模型，需要在近壁区使用细网格分辨粘性子层，以便施加层流流动边界条件。因此，动量方程中壁面边界条件为第 15 章所述无滑移条件，本章不再重复。\(k-\omega\) 与 \(k-\varepsilon\) 湍流模型及其变体均属此类。对所有这些模型，\(k\)、\(\varepsilon\) 与 \(\omega\) 的渐近壁面边界条件为（方程 17.98）：\(k_w \to 0\)，\(\varepsilon_w \to 2 \nu k/(d_\perp)_C^2\)，\(\omega_w \to 6 \nu/(C_{b1} (d_\perp)_C^2)\)。\(\omega_w\) 由方程 17.98 给出的值旨在施加于近壁第一个内部网格点，对位于 \(d^+\) 或 \(d^* < 2.5\) 的网格点有效，且在该区域需要 15–20 个节点以达到网格无关解。方程 17.98 中引入系数 10 是基于 Menter 的建议，它消除了在壁面附近点以外再为 \(\omega\) 指定值的需要。
17.8.6 自动近壁处理（Automatic Near-Wall Treatment）
使用壁面函数方法时，先计算到壁面的归一化距离，以判断第一个内部点位于粘性子层还是惯性子层。粘性子层中流动视为层流。由于标准 \(k-\varepsilon\) 模型在完全湍流区有效，并不包含模拟粘性效应的阻尼函数，它会把误差引入解。为消除此误差，应注意使第一个内部点位于对数区。但这又消除了粘性子层的影响，而该影响对解可能很重要。因此希望在不需要过细近壁网格的条件下得到粘性子层分辨的能力。这由 \(k-\omega\) 模型及其变体提供，因 \(\omega\) 方程可在无任何额外阻尼函数的条件下一直积分到壁面。其思想是开发一种方法，按到壁面的归一化距离值自动在低、高雷诺数公式之间切换。由于 \(\omega\) 方程在粘性子层与惯性子层有已知的解析解，该任务是可行的。为此，按 [34] 将两个解平滑混合（方程 17.99）： \(\omega = (\omega_{\mathrm{vis}}^2 + \omega_{\mathrm{log}}^2)^{0.5}\) ，其中（方程 17.100）\(\omega_{\mathrm{vis}} = 6 \nu/(C_{b1} d_\perp^2)\)，
。
动量方程中剪切应力按方程 17.101 计算： \(|\tau_w| = \rho u_\tau u^* = \rho (u^*/u^+) |v_C - v_w|_{\parallel}\) ，其中近壁区的 \(u_\tau\) 与 \(u^*\) 由方程 17.102 与 17.103 给出：
， \(u^* = \{u_{\mathrm{vis}}^{*4} + (u_{\mathrm{log}}^{*4})^{0.25}\}\) ， \(u_{\mathrm{vis}}^* = \sqrt{\mu |v_C - v_w|_{\parallel}/(\rho (d_\perp)_C)}\) ，\(u_{\mathrm{log}}^* = \beta^{*1/4} k^{1/2}\)。在求解 \(k\) 方程时，壁面梯度被置零，近壁单元中的生产项被修改为方程 17.104：
。
基于网格间距，该混合方法允许壁面处理从粘性子层到壁面函数的光滑切换。
17.8.7 近壁热传递（Near-Wall Heat Transfer）
与速度剖面类似，使用高雷诺数湍流模型时近壁温度剖面也应作修正。这些剖面由 Reynolds 类比从速度剖面得到，通过修改对数律以提供把近壁节点温度 \(T_C\)、壁面温度 \(T_w\) 与壁面热流 \(q_w\) 关联起来的方程。流程先按方程 17.105 定义归一化温度 \(T^+ = (T_w - T)/T^*\)，其中 \(T^*\) 由方程 17.106 给出：\(T^* = q_w/(\rho c_p u^*)\)。把上述两式合并，归一化温度方程变为方程 17.107：\(T^+ = \rho c_p u^* (T_w - T)/q_w\)。
对标准壁面函数公式，若第一个内部点位于粘性子层，则归一化温度按方程 17.108 计算：\(T^+ = \mathrm{Pr} \cdot d^*\)；若位于惯性子层，则按壁面定律（方程 17.109）计算：\(T^+ = 2.12 \ln(d^*) + b(\mathrm{Pr})\)，其中 \(b(\mathrm{Pr}) = (3.85 \mathrm{Pr}^{1/3} - 1.3)^2 + 2.12 \ln(\mathrm{Pr})\) （方程 17.110），\(\mathrm{Pr}\) 为层流 Prandtl 数。对可缩放壁面函数公式，同样的方程被使用，只是 \(d^+\) 被替换为 \(\widetilde{d^+}\)，如前所述。对自动近壁处理方法，采用 Kader [49] 给出的方程，它把粘性子层与壁面定律混合。事实上该方法也可与标准壁面函数与可缩放壁面函数方法联用。按此公式，归一化温度按方程 17.111 计算：
，其中混合函数 \(\Gamma\) 由方程 17.112 给出： \(\Gamma = 0.01 (\mathrm{Pr} \cdot d^*)^4/(1 + 5 \mathrm{Pr}^3 d^*)\) 。
此外，方程 17.107 可按温度的物理边界条件按两种方式解释。对给定的壁面热流边界条件，\(T_w\) 的数值条件是一个固定值并按方程 17.113 迭代更新：
。若温度被强加在边界上，则温度梯度必须保持数值固定并按方程 17.114 计算：
。
17.8.8 其他边界条件（Other Boundary Conditions）
除壁面外还需要其他边界条件，包括入口、出口与对称边界。在域的入口处，湍流动能与耗散率通常未知，是需要的。若值可由测量给出则应使用，否则必须估算。最简单的方法是直接给 \(k\) 与 \(\varepsilon\) 赋值。也可通过湍流强度 \(I = \sqrt{v' \cdot v'/\sqrt{v \cdot v}}\)（方程 17.115）来指定，湍流动能由 \(k = (1/2) I^2 (v \cdot v)\)（方程 17.116）得到。湍流强度通常使用 1%–10% 的值，其中 < 1% 视为低、> 10% 视为高。湍流耗散率 \(\varepsilon\) 与湍流频率 \(\omega\) 由方程 17.23 与 17.32 通过指定湍流长度尺度计算，湍流长度尺度的值取决于最大涡旋尺寸但通常设为剪切层宽度或域尺度的十分之一，表达为 \(\varepsilon = C_\mu k^{3/2}/\ell\)、\(\omega = \sqrt{k}/\ell\)（方程 17.117）。\(\varepsilon\) 与 \(\omega\) 的值也可由已知的 \(k\) 与湍流与层流粘度比 \(\mu_t/\mu\) 计算（方程 17.118）：\(\varepsilon = C_\mu q k^2 (\mu/\mu_t)\)，\(\omega = q k (\mu/\mu_t)\)。
在出口与对称边界，\(k\)、\(\varepsilon\) 与 \(\omega\) 的处理与前几章中一般标量变量 \(\phi\) 的处理类似，无需重复。
17.9 计算到壁面的法向距离（Calculating Normal Distance to the Wall）
在 BSL 与 SST 湍流模型中，整个域上都需要到最近壁面的法向距离 \(d_\perp\)，以确定 \(k-\varepsilon\) 与 \(k-\omega\) 计算区的界面，这一界面反映在 BSL 模型的方程 17.42 与 SST 模型的方程 17.47 中。在三维情况下，求解 \(d_\perp\) 的搜索过程即便对固定网格也是计算昂贵的；对运动网格情况更糟，每次时间步都必须重新搜索。这迫使研究者在 \(d_\perp\) 计算中引入带有较大误差的近似。
为避免昂贵的搜索过程，已经发展出基于求解 \(d_\perp\) 微分方程的技术。这些方法基于求解 Poisson、Eikonal 或 Hamilton-Jacobi 方程 [50–53]。这些方程的吸引力在于其组成涉及 CFD 解算器中现成的梯度与/或 Laplace 算子。这使得该方法易于实现、稳定、经济，尤其对运动网格。
采用 Poisson 类方法，求解如下变量 \(\phi\) 的微分方程（方程 17.119）：\(\nabla^2 \phi = -1\)，其约束为（方程 17.120）：\(\phi = 0\) 在壁面；\(\nabla \phi \cdot \mathbf{n} = 0\) 在其他位置。到最近壁面的法向距离按方程 17.121 用 \(\phi\) 的预测值与其梯度计算：
。
在求解湍流问题时，方程 17.119 用第 8 章所述方法在求解问题所用的网格上离散。在求解过程中，先对方程 17.119 获得收敛解，由此计算到壁面的法向距离，然后求解湍流问题。如图 17.3 所示，回忆 \(g_{D_{\mathrm{eff}}}\) 按方程 17.122 定义：\(g_{D_{\mathrm{eff}}} = E_f/d_{CF}\)。方程 17.119 的离散形式可写为
（方程 17.123），其中（方程 17.124）： \(a_F = \mathrm{Flux}_{F,f} = -g_{D_{\mathrm{eff}},f}\) ；
；
。
此外，方程 17.120 所述 Dirichlet 与 Von Neumann 边界条件按第 8 章所述方式处理。在对方程 17.119 求解后，域内所有单元形心处到壁面的法向距离用方程 17.121 计算，梯度按第 9 章详细方法计算。
17.10 计算指针（Computational Pointers）
OpenFOAM® [54] 为可压缩与不可压缩流动实现了若干 LES 与 RANS 湍流模型。所有模型的根目录记为"FOAM_SRC/turbulenceModels"。不可压缩湍流模型位于子目录"FOAM_SRC/turbulenceModels/incompressible"内，其下三个子目录分别为"LES"（指大涡模拟方法）、"RAS"（指 Reynolds 平均 Navier-Stokes 方法）与"turbulenceModel"。前两个子目录"LES"与"RAS"分别定义 LES 与 RAS 模型的特殊特性，"turbulenceModel"中定义不可压缩 RAS 与 LES 模型的抽象基类。基类定义了一系列必须在派生类中具体化的抽象虚函数，如 Listing 17.1 所示（返回湍流粘度 nut()、有效粘度 nuEff()、湍流动能 k()、耗散率 epsilon()、Reynolds 应力张量 R()）。
基类还定义了到壁面的法向距离，这是湍流模型中可被所有派生类使用的有用量。该语句如 Listing 17.2 所示：const nearWallDist& y() const { return y_; }。
为更好理解定义 OpenFOAM® 湍流模型的代码结构，下面以带 Spalding 壁面函数的 \(k-\varepsilon\) 模型与 SST \(k-\omega\) 模型为例。模型定义位于目录"FOAM_SRC/turbulenceModels/incompressible/RAS"，所有 RANS 湍流模型都放在此处。对 RANS 模型，OpenFOAM® 定义了一个额外的非虚基类 RASModel（由 turbulenceModel 派生），所有模型都由其派生（Listing 17.3）。该类主要作为包装类，用于限制有效粘度与湍流量值的范围，并提供字典定义。
17.10.1 \(k-\varepsilon\) 模型（The \(k-\varepsilon\) Model）
kEpsilon 类实现由方程 17.25 与 17.26 给出的标准 \(k-\varepsilon\) 模型。它定义模型所需的常数与变量集合，如 Listing 17.4 所示（受保护数据成员：Cmu_、C1_、C2_、sigmaEps_ 系数与 k_、epsilon_、nut_ 场）。它可从字典"RASProperties"中以"kEpsilon"名字（由类的 TypeName 定义）选择，如 Listing 17.5 所示：TypeName("kEpsilon")。
由于派生自基虚类 turbulenceModel，其所有虚基函数必须按下述定义（Listing 17.6）：返回包含层流应力的有效应力张量 devReff() 与返回动量方程的源项 divDevReff(volVectorField& U)。
divDevReff 函数返回动量方程中包括 Reynolds 应力在内的扩散贡献，即方程 17.125：
。该函数返回的数据类型（Listing 17.7）是 fvMatrix 形式，按 tmp<fvVectorMatrix> kEpsilon::divDevReff(volVectorField& U) const { return (- fvm::laplacian(nuEff(), U) - fvc::div(nuEff()*dev(T(fvc::grad(U))))); } 实现，其中速度场的 Laplace 算子被拆分为隐式与显式矩阵贡献。值得指出，divDevReff 项只依赖速度梯度场与总有效粘度 nuEff()。这意味着 OpenFOAM® 通过只修改壁面处的湍流粘度实现壁面剪切应力贡献，即实现方程 17.77。
湍流模型方程的组装与求解通过成员函数 correct() 定义，如 Listing 17.8。OpenFOAM® 先组装并求解 \(\varepsilon\) 方程，然后求解 \(k\) 方程，使用 Listing 17.9 中的脚本，其中 \(G\) 场定义方程 17.28 中的全局生产项 \(P_k\)。基于壁面函数方法，需要在求解 \(\varepsilon\) 方程前按方程 17.81 与 17.83 修改所有近壁网格形心处的 \(\varepsilon\) 值。在 OpenFOAM® 中，\(\varepsilon\) 与 \(P_k\) 在壁面处的值通过为 \(\varepsilon\) 场定义特殊壁面边界条件并通过 epsilon_.boundaryField().updateCoeffs(); 强制修改，\(\varepsilon\) 的专门壁面定义位于目录"FOAM_SRC/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction"。
基于壁面函数模型，该类必须更新两个变量 \(\varepsilon\) 与 \(P_k\) 的值。该操作由 calculate 函数执行（其中 \(G\) 在 Listing 17.10 中对应 \(P_k\) 生产项）。w 变量是考虑与多于一个壁面（即角点）相接触的边界网格的权重因子，对标准面其值为 1：epsilon[cellI] += w*Cmu75*pow(k[cellI], 1.5)/(kappa_*y[faceI])，G[cellI] += w * (nutw[faceI] + nuw[faceI]) * magGradUw[faceI] * Cmu25*sqrt(k[cellI])/(kappa_*y[faceI])。
值得注意的是，虽然 \(G\) 或 \(P_k\) 是源项，\(\varepsilon\) 是通过输运方程求解的场。为在网格中施加先前计算的值，必须在正确位置操控矩阵以返回正确的值。因此类 epsilonWallFunctionFvPatchScalarField 按 Listing 17.11 从 fixedInternalValueFvPatchField 派生。fixedInternalValueFvPatchField（Listing 17.12）是一个包装类，其中包含一个特殊函数 manipulateMatrix，通过 matrix.setValues(this->patch().faceCells(), this->patchInternalField()) 把边界内部场作为约束施加到矩阵中。该函数由 kEpsilon::correct() 在用 epsEqn().boundaryManipulate(epsilon_.boundaryField()) 组装矩阵后调用。
一旦 \(\varepsilon\) 方程被求解，OpenFOAM® 即继续组装并求解 \(k\) 方程，如 Listing 17.13 所示（kEqn() 矩阵由 fvm::ddt(k_) + fvm::div(phi_, k_) - fvm::laplacian(DkEff(), k_) == G - fvm::Sp(epsilon_/k_, k_) 构成）。此时由于 \(P_k\) 已在壁面被修改且 \(k\) 的边界条件只是零梯度（zeroGradient）类型，无需额外操控。
计算完 \(k\) 与 \(\varepsilon\) 的值后，湍流涡粘度 \(v_t\)（\(=\mu_t/\rho\)）被更新并在壁面边界按方程 17.77 修正，由 Listing 17.14 实现：nut_ = Cmu_*sqr(k_)/epsilon_，随后 nut_.correctBoundaryConditions()。
volScalarField 字段 nut（\(v_t\)）在 Listing 17.15 中 kEpsilon 类的构造函数中定义，其中使用函数 autoCreateNut("nut", mesh_)。该函数在文件"backwardsCompatibilityWallFunctions.C"中定义，置于目录"FOAM_SRC/src/turbulenceModels/incompressible/RAS/backwardsCompatibility/wallFunctions"。autoCreateNut 仅当工作目录中已有 nut 文件时，通过读取该文件及相关边界类型来创建 \(v_t\) 对象（Listing 17.16 中以 if (nutHeader.headerOk()) 检查）。
当 nut 文件未找到时，使用 Listing 17.17 中 if (isA<wallFvPatch>(bm[patchI])) 检查并通过边界补丁类定义 nutkWallFunctionFvPatchScalarField 施加标准 Spalding 壁面函数。Listing 17.18 中描述的类 nutkWallFunctionFvPatchScalarField 位于目录"FOAM_SRC/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/nutWallFunctions"，它继承自名为 nutWallFunctionFvPatchScalarField 的基类。Listing 17.19 中 nutWallFunctionFvPatchScalarField 是一个包装类，通过 updateCoeffs() 函数包装 \(v_t\) 边界值的修改：operator==(calcNut())，随后调用 fixedValueFvPatchScalarField::updateCoeffs()。其中 calcNut 被定义为纯虚，必须在派生类中具体化。派生类 nutkWallFunctionFvPatchScalarField 按 Spalding 假设和方程 17.77 实现 calcNut 函数（Listing 17.20）：yPlus = Cmu25*y[faceI]*sqrt(k[faceCellI])/nuw[faceI]，if (yPlus > yPlusLam_) nutw[faceI] = nuw[faceI]*(yPlus*kappa_/log(E_*yPlus) - 1.0);。
一旦基于壁面函数值更新了湍流涡粘度，动量剪切应力被正确计算。
17.10.2 SST \(k-\omega\) 模型（The SST \(k-\omega\) Model）
kOmegaSST 类实现由方程 17.38–17.49 描述的 SST \(k-\omega\) 模型版本。模型所用变量与常数由 Listing 17.21 中的脚本定义（受保护数据：alphaK1_、alphaK2_、alphaOmega1_、alphaOmega2_、gamma1_、gamma2_、beta1_、beta2_、betaStar_、a1_、b1_、c1_ 系数、F3_ 开关、到壁面距离 y_、k_、omega_、nut_ 字段）。还定义了额外的私有成员函数 blend、alphaK、alphaOmega、beta、gamma（按 blend(F1, alphaK1_, alphaK2_) = F1*(psi1 - psi2) + psi2 等）以表示方程 17.39 中出现的混合系数。函数 F1 与 F2 分别描述方程 17.41 与 17.47 给出的变量。函数 F3 不出现在 Menter 的原始模型中，是 Hellsten [55] 提出的针对粗糙壁面的修改。
如 Listing 17.22 所示，模型可从字典"RASProperties"中以"kOmegaSST"名字（由类的 TypeName 定义）选择：TypeName("kOmegaSST")。
kOmegaSST 类派生自基虚类 turbulenceModel，其所有虚基函数都相应地被具体化（Listing 17.23）。模型实现细节见"FOAM_SRC/turbulenceModels/incompressible/RAS/kOmegaSST/kOmegaSST.C"。如 kEpsilon 类，correct 函数求解 SST 湍流模型的全部方程。首先设置并求解 \(\omega\) 方程，Listing 17.24 给出了脚本：RASModel::correct() 调用后，构造 volScalarField S2(2*magSqr(symm(fvc::grad(U_)))) 与生产项 G = nut_*S2，随后用 omega_.boundaryField().updateCoeffs() 更新壁面值并构造 CDkOmega = (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_，再调用 F1(CDkOmega)。\(\omega\) 方程的矩阵按 fvm::ddt(omega_) + fvm::div(phi_, omega_) - fvm::laplacian(DomegaEff(F1), omega_) == gamma(F1)*S2 - fvm::Sp(beta(F1)*omega_, omega_) - fvm::SuSp((F1 - scalar(1))*CDkOmega/omega_, omega_) 构造。
该脚本中，\(G\) 场表示方程 17.28 定义的 \(P_k\) 生产项，而 CDkOmega 与 F1 场基于方程 17.41、17.42 与 17.47 计算。与 kEpsilon 类一样，\(\omega\) 与 \(P_k\) 的值按壁面函数方法在壁面被修改，修改通过 omega_.boundaryField().updateCoeffs() 强制施加。\(\omega\) 的边界定义位于目录"FOAM_SRC/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/omegaWallFunctions/omegaWallFunction"。基于壁面函数模型，该类必须更新两个变量 \(\omega\) 与 \(P_k\) 的值。Listing 17.25 给出了 calculate 函数（其中 G 对应 \(P_k\) 生产项）：omegaVis = 6.0*muw[faceI]/(rhow[faceI]*beta1_*sqr(y[faceI]))，omegaLog = sqrt(k[cellI])/(Cmu25*kappa_*y[faceI])，omega[cellI] += w*sqrt(sqr(omegaVis) + sqr(omegaLog))，并修改 G[cellI] += w*(mutw[faceI] + muw[faceI])*magGradUw[faceI]*Cmu25*sqrt(k[cellI])/(kappa_*y[faceI])，按方程 17.81、17.99、17.100 计算。
为施加先前计算的网格值，必须按与 kEpsilon 模型相同的方法操控矩阵。一旦 \(\omega\) 方程被求解，OpenFOAM® 继续组装并求解 \(k\) 方程，Listing 17.26（与 \(k-\varepsilon\) 模型中相同的脚本，但用 \(\omega\) 替换 \(\varepsilon\)）。
计算完 \(k\) 与 \(\omega\) 的值后，按方程 17.46 用 Listing 17.27 中的脚本更新湍流涡粘度 \(v_t\)（\(=\mu_t/\rho\)）：nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2))，随后 nut_.correctBoundaryConditions()。
Menter 可靠性约束被施加于涡粘度，常数 b1_ 取 1.0，而 Listing 17.28 中描述的 F23 函数默认返回方程 17.47 定义的 F2 函数；F3_ 为 true 时返回 F2() * F3()。最后按方程 17.77 在壁面边界处修正涡粘度，与 kEpsilon 类相同的流程，使用 Spalding 壁面函数由 nutkWallFunctionFvPatchScalarField 类施加。
17.10.3 simpleFoamTurbulent
simpleFoamTurbulent 解算器是 simpleFoamImproved 解算器（第 15 章）的湍流扩展。正如前文所述，湍流影响输运方程的扩散项。其效应通过有效粘度（即层流粘度与涡粘度之和）体现在方程中。
为在 OpenFOAM® 解算器中纳入湍流，调用基虚类 RASModel 并做少量修改。第一处修改是用变量湍流属性替换常量层流输运属性。因此在"createFields.H"文件中实例化 RASModel 类型的 turbulence 对象，如 Listing 17.29：singlePhaseTransportModel laminarTransport(U, mdotf) 后通过 incompressible::RASModel::New(U, mdotf, laminarTransport) 创建 autoPtr。
RASModel 的定义需要激活 singlePhaseTransportModel 类，它定义层流粘度的一般输运模型（可由字典设为常量或 Sutherland 模型那样的温度函数），以及速度与质量通量。如前所述，这些量是定义湍流量输运方程所必需的。在 Listing 17.29 中值得注意的是将 turbulence 对象定义为 autoPtr，基本上可当作 C++ 中的标准指针使用。
第二处主要修改（Listing 17.30）是对动量方程，其中扩散现在使用前述 divDevReff(U) 计算：fvVectorMatrix UEqn(fvm::ddt(U) + fvm::div(mdotf, U) + fvm::SuSp(-fvc::div(mdotf), U) + turbulence->divDevReff(U))。
最后一处主要修改（Listing 17.31）在主解算器文件中：#include "UEqn.H" 与 #include "ppEqn.H" 之后调用 turbulence->correct()。该语句在每次被调用时激活湍流模型方程的求解，从而计算出用于动量方程的涡粘度。
17.11 小结（Closure）
本章介绍了对不可压缩湍流建模所需的额外步骤，讨论了 \(k-\varepsilon\)、\(k-\omega\) 模型及其若干变体，详细说明了使用壁面函数对近壁区的建模。这结束了本书打算讨论的内容。下一章将讨论 OpenFOAM® 与 uFVM 中边界条件的实现。
17.12 习题（Exercises）
未覆盖的源文习题（17.12 共 14 道）按 skill 规则跳过，不在内概述中展开。Exercise 1–2 关于由湍流强度与长度尺度（及粘度比）反算 \(k\)、\(\varepsilon\)、\(\omega\)；Exercise 3 把方程 17.29–17.31 的可实现性约束引入标准 \(k-\varepsilon\)；Exercise 4 对标准 \(k-\omega\) 推导相应约束；Exercise 5 由 BSL 模型当 \(F_1 \equiv 1\) 时回退到标准 \(k-\varepsilon\)；Exercise 6 用 Newton-Raphson 线性化方程 17.63；Exercise 7 用 \(u^*\)、\(d^+\)、\(k\)、\(\nu\) 重写方程 17.83 与 17.84；Exercise 8 在 OpenFOAM® 中实现方程 17.119–17.121 的到壁面法向距离；Exercise 9 用等效粒径粗糙度 \(h_s^+\) 把粗糙壁对数律重写为经典形式；Exercise 10–14 通过 Doxygen 文档列出 RASModel 与 nutWallFunctionFvPatchScalarField 的派生类并对比 epsilonLowReWallFunction 与方程 17.98 的低雷诺数公式。
本章个人批注
本章把 Moukalled 等教材对湍流建模的处理集中讲清楚了一条主线：Boussinesq 假设把非线性 Reynolds 应力线性化后，问题被化简为求湍流粘度 \(\mu_t\)，而 \(\mu_t\) 由不同的二方程模型（\(k-\varepsilon\)、\(k-\omega\)）以及它们的混合版本（BSL、SST）给出。整个章节实际上是在建立一条逐步细化的"封闭关系"链：先解决 Reynolds 应力未知性（Boussinesq 假设），再解决 \(\mu_t\) 的微分方程（\(k\) 方程 + 第二个输运方程），最后解决近壁区积分（壁面函数 / 自动近壁处理）。
让我印象最深的是 17.5 节的两处缺陷修正：(1) 驻点反常（\(P_k\) 在高应变率下被高估），通过引入局部湍流时间尺度约束 \(t_s\) 把生成项增长率从 \(S_t^2\) 变为 \(S_t\)；(2) 边界层 Bradshaw 假设违反，通过 SST 中 \(\mu_t = \rho a_1 k/\mathrm{Max}(a_1 \omega, 2 S_t F_2)\) 把涡粘度限制住。这两处都反映出工业湍流模型的发展是问题驱动的——每个新模型都是对前一个模型在某种特定流型下失效的补丁。BSL 与 SST 的混合函数 \(F_1\)、\(F_2\) 是这种"组合式建模"思路的典型：通过对 \(k-\varepsilon\) 与 \(k-\omega\) 的线性混合，在不同区域（近壁 vs 远场）切换优势模型。
关于壁面函数：标准、改进、可缩放三套壁面函数（17.8.2–17.8.4）构成了一个"渐进修复"系列。标准壁面函数在局部平衡假设下成立但在分离点/再附点给出零粘度；改进壁面函数用 \(\sqrt{k}\) 替换 \(u_\tau\) 后允许非零粘度；可缩放壁面函数则通过 \(\widetilde{d^*} = \mathrm{Max}(d^*, d_{\lim}^*)\) 让网格无关性更稳健。三者并不互相替代，而是用于不同精度需求。
关于 17.10 节 OpenFOAM® 实现：作者用了 30 个 Listing 来展示 \(k-\varepsilon\) 与 SST \(k-\omega\) 模型的完整实现。divDevReff 通过只修改壁面湍流粘度来实现壁面剪切应力（即方程 17.77 的代码实现），这是 OpenFOAM® 框架的一个关键设计——把壁面处理从动量方程中分离出来作为边界条件而不是源项。epsilonWallFunctionFvPatchScalarField 通过继承 fixedInternalValueFvPatchField 在矩阵层面强加预计算的 \(\varepsilon\) 值，这是另一个值得记住的范式。SST 模型中的 CDkOmega、F1、F2、F23 完整复现了方程 17.41–17.47 的混合函数逻辑。
未覆盖的源文习题（17.12 共 14 道）按 skill 规则跳过。其中 Exercise 8（实现法向距离 Poisson 方程）与 Exercise 10–14（Doxygen 文档查阅）是直接对应 17.9 与 17.10 节的延伸；Exercise 9（粗糙壁等效 \(d_{\mathrm{eff}}^+\) 公式）则呼应 17.8.2 节的对数律公式。
与上下章的衔接（一段话）
第 17 章接续第 15 章（不可压缩流体流动计算）和第 16 章（可压缩流体流动计算）：前两章建立的压力基 SIMPLE 框架是 RANS 求解的载体，本章则把湍流这一未被前述封闭方程组显式建模的物理现象通过 Boussinesq 假设与二方程模型纳入框架。第 17 章结束本书的核心方程体系：从第 3 章的 Navier-Stokes 出发，到第 5 章的一般输运方程离散，第 7–12 章的线性解法与压力-速度耦合，第 13 章的非稳态项处理，第 14 章的源项线性化与松弛，第 15–16 章的可压缩/不可压缩压力基算法，再到本章的湍流模型。本章末段明确说"the next chapter will discuss the implementation of boundary conditions in OpenFOAM® and uFVM"，对应第 18 章 Boundary Conditions —— 是对第 5 章（一般输运方程边界条件）与第 15–17 章各具体边界条件处理（无滑移壁面、入流、出流、对称、入口湍流量指定）的总结式实现指南。从 Moukalled 团队的写作线索看，本书在第 17 章收尾、留第 18 章做边界条件实现是合理的——湍流模型使方程数加倍，是闭合方程组的最后一块拼图；而 OpenFOAM®/uFVM 的代码层面的边界条件实现则可以视为对前面所有边界条件讨论的一次"集中编译"。