对流换热仿真-主题013_k-ω湍流模型-k-ω湍流模型

kkchenjj

46人浏览 · 2026-04-17 21:48:05

kkchenjj · 2026-04-17 21:48:05 发布

主题013：k-ω湍流模型

一、k-ω湍流模型概述

1.1 k-ω模型的历史与发展

k-ω湍流模型是由David C. Wilcox于1988年首次提出的两方程涡粘度模型。与k-ε模型相比，k-ω模型选择湍流比耗散率ω（单位湍动能的耗散率）作为第二个输运变量，而非湍流耗散率ε。这一选择使得k-ω模型在近壁区域具有显著优势。

模型发展历程：

1988年：Wilcox提出标准k-ω模型，成功解决了k-ε模型在近壁区域的数值刚性问题
1994年：Wilcox改进模型，引入交叉扩散项，提高了对自由剪切流的预测精度
1994年：Menter提出SST（Shear Stress Transport）k-ω模型，结合k-ε和k-ω的优点
2003年：Menter进一步完善SST模型，改进了剪切应力限制器
2006年：SST模型成为目前最广泛使用的k-ω变体，被纳入各种CFD软件

关键发展里程碑：

1988年 - 标准k-ω模型的诞生：
David Wilcox在NASA兰利研究中心工作期间，基于对k-ε模型在近壁区域数值问题的深入分析，提出了使用比耗散率ω替代耗散率ε的创新思路。这一改变不仅改善了数值特性，还使得模型能够直接积分到壁面，无需复杂的壁面函数。

1994年 - Wilcox改进模型：
针对标准k-ω模型在自由剪切流中预测偏差的问题，Wilcox引入了交叉扩散项（cross-diffusion term），这一改进显著提高了模型对射流和尾流的预测能力。

1994年 - SST模型的提出：
Florian Menter在NASA艾姆斯研究中心工作期间，提出了Shear Stress Transport（剪切应力输运）k-ω模型。SST模型通过混合函数将k-ω和k-ε模型结合起来，既保留了k-ω模型在近壁区域的优势，又克服了其对自由流敏感的问题。

2003年 - SST模型的进一步完善：
Menter对SST模型进行了进一步改进，引入了更精确的剪切应力限制器，使得模型在预测逆压梯度流动和分离流方面表现出色。

当前应用状态：
SST k-ω模型已成为航空航天、汽车、能源、船舶等领域CFD仿真的首选湍流模型之一，被广泛应用于ANSYS Fluent、OpenFOAM、STAR-CCM+、CFX等主流商业和开源软件中。该模型在工业界和学术界都得到了广泛认可，是工程实践中使用最广泛的湍流模型之一。
在这里插入图片描述

1.2 k-ω模型的基本思想

k-ω模型的核心思想是使用湍动能k和湍流比耗散率ω来描述湍流状态：

湍动能k：与k-ε模型相同，表示湍流脉动的动能

$\frac{1}{2} \overline{u_i' u_i'}$

湍流比耗散率ω：表示单位湍动能的耗散速率，具有时间倒数的量纲

$\omega = \frac{\varepsilon}{C_\mu k} \approx \frac{\varepsilon}{0.09 k}$

涡粘度计算：

$\nu_t = \frac{k}{\omega}$

与k-ε模型的 $νt=Cμk2ε\nu_t = C_\mu \frac{k^2}{\varepsilon}$ 相比，k-ω模型的涡粘度公式更加简洁，且在壁面附近具有更好的数值特性。

k-ω模型的物理直观性：

k-ω模型中的变量具有明确的物理意义：

k代表湍流的"能量"
ω代表湍流的"频率"或"速率"
涡粘度 $νt=k/ω\nu_t = k/\omega$ 可以理解为能量与频率的比值

这种直观的物理解释使得k-ω模型在工程应用中更易于理解和使用。

1.3 k-ω与k-ε模型的对比

特性	k-ε模型	k-ω模型
第二变量	ε（耗散率）	ω（比耗散率）
近壁行为	需要壁面函数	可直接积分到壁面
自由流敏感性	低	高（标准模型）
分离流预测	一般	较好
逆压梯度流动	一般	较好
数值稳定性	壁面附近刚性	壁面附近稳定

k-ω模型的优势：

近壁区域表现优异：可以直接积分到壁面，无需复杂的壁面函数
对逆压梯度敏感：能够更好地预测边界层分离
数值稳定性好：在壁面附近不会出现k-ε模型的数值刚性问题

模型选择决策树：

在实际工程应用中，可以按以下流程选择湍流模型：

分析流动特性：
- 是否存在分离？
- 是否有强逆压梯度？
- 是否涉及近壁换热？
- 是否强旋流？
根据特性选择：
- 分离流/逆压梯度 → SST k-ω
- 简单内部流动 → k-ε或标准k-ω
- 强旋流 → RSM
- 非定常特性重要 → LES/DES
考虑计算资源：
- 资源充足 → 高精度模型
- 资源有限 → 效率优先
验证模型适用性：
- 查阅类似流动的文献
- 进行网格敏感性分析
- 与实验数据对比

k-ω模型的局限：

自由流敏感性：标准k-ω模型对入口边界条件中的ω值非常敏感
自由剪切流：在射流和尾流中的预测精度不如k-ε模型

1.4 k-ω模型的物理基础

比耗散率ω的物理意义：

比耗散率ω定义为湍流耗散率ε与湍动能k的比值：

$\omega = \frac{\varepsilon}{\beta^* k}$

其中 $β∗=0.09\beta^* = 0.09$ 是模型常数。从物理角度看，ω代表湍流时间尺度的倒数：

$\omega \sim \frac{1}{\tau_t}$

其中 $τt\tau_t$ 是湍流时间尺度。这一物理意义使得ω在描述湍流动力学时具有直观性。

与k-ε模型的关系：

k-ω和k-ε模型可以通过变量变换相互转换。从k-ε方程出发，利用 $ε=β∗ωk\varepsilon = \beta^* \omega k$ ，可以推导出k-ω方程。然而，这种变换不是精确的，因为两个模型在扩散项和源项的处理上存在差异。

近壁渐近行为：

在壁面附近，k-ω模型表现出良好的渐近特性：

当 $\to 0$ 时， $\sim y^2$
当 $\to 0$ 时， $ω∼1y2\omega \sim \frac{1}{y^2}$
涡粘度 $νt=kω∼y2\nu_t = \frac{k}{\omega} \sim y^2$ ，与壁面附近的物理行为一致

这种渐近行为使得k-ω模型可以直接积分到壁面，无需引入壁面函数。

数值刚性的消除：

k-ε模型在近壁区域存在数值刚性问题，主要原因是ε方程中的源项 $ε2/k\varepsilon^2/k$ 在壁面附近趋于无穷大。k-ω模型通过使用ω作为变量，将耗散项变为 $βω2\beta \omega^2$ ，在壁面附近保持有限值，从而消除了数值刚性。

k-ω模型的能量级联描述：

k-ω模型通过k和ω两个变量描述了湍流的能量级联过程：

能量输入：大尺度涡从平均流动中获取能量，通过生成项 $P_k$ 传递给湍流
能量传递：湍动能k代表了含能尺度的能量
能量耗散：比耗散率ω表征了能量向小尺度传递并最终耗散的速率

湍流时间尺度和长度尺度：

基于k和ω，可以定义湍流特征尺度：

时间尺度：
$\tau_t = \frac{1}{\omega}$

长度尺度：
$l_t = \frac{\sqrt{k}}{\omega} = \frac{\nu_t}{\sqrt{k}}$

Reynolds数：
$Re_t = \frac{k}{\nu \omega}$

这些尺度在分析湍流特性和验证模型时非常有用。

二、标准k-ω模型

2.1 控制方程

标准k-ω模型由Wilcox提出，其输运方程为：

湍动能输运方程：

$\frac{\partial k}{\partial t} + U_j \frac{\partial k}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \nu + \sigma_k \frac{k}{\omega} \right) \frac{\partial k}{\partial x_j} \right] + P_k - \beta^* k \omega$

比耗散率输运方程：

$\frac{\partial \omega}{\partial t} + U_j \frac{\partial \omega}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \nu + \sigma_\omega \frac{k}{\omega} \right) \frac{\partial \omega}{\partial x_j} \right] + \alpha \frac{\omega}{k} P_k - \beta \omega^2$

其中各项的物理意义为：

对流项： $Uj∂k∂xjU_j \frac{\partial k}{\partial x_j}$ 和 $Uj∂ω∂xjU_j \frac{\partial \omega}{\partial x_j}$ 表示湍流特性随流体输运

扩散项：分子扩散和湍流扩散共同作用

生成项：

k方程： $Pk=νt(∂Ui∂xj+∂Uj∂xi)∂Ui∂xjP_k = \nu_t \left( \frac{\partial U_i}{\partial x_j} + \frac{\partial U_j}{\partial x_i} \right) \frac{\partial U_i}{\partial x_j}$
ω方程： $αωkPk\alpha \frac{\omega}{k} P_k$

耗散项：

k方程： $−β∗kω-\beta^* k \omega$ （对应k-ε模型中的 $−ε-\varepsilon$ ）
ω方程： $−βω2-\beta \omega^2$

2.2 模型常数

标准k-ω模型的常数由Wilcox通过理论分析和实验数据拟合确定：

$\alpha = \frac{5}{9}, \quad \beta = \frac{3}{40}, \quad \beta^* = \frac{9}{100}$

$\sigma_k = \frac{1}{2}, \quad \sigma_\omega = \frac{1}{2}$

这些常数在广泛的流动条件下表现良好，但对于某些特定流动可能需要调整。

常数的物理意义：

α = 5/9 ≈ 0.5556：控制ω方程中生成项的强度，影响湍流的发展速率
β = 3/40 = 0.075：控制ω方程中耗散项的强度
β = 9/100 = 0.09*：控制k方程中耗散项的强度，与k-ε模型中的Cμ相同
σk = 0.5：湍动能的湍流Prandtl数
σω = 0.5：比耗散率的湍流Prandtl数

常数的标定过程：

Wilcox通过以下实验数据对模型常数进行了标定：

零压力梯度平板边界层
管道流动
自由剪切流（射流、尾流）
逆压梯度边界层

这些常数在广泛的流动条件下表现良好，但对于某些特定流动（如强旋流、高应变率流动）可能需要调整。

2.3 壁面边界条件

k-ω模型的一个重要优势是可以直接积分到壁面，无需使用壁面函数。

壁面边界条件：

$k_w = 0$

$\omega_w = \frac{6\nu}{\beta_1 y_1^2}$

其中 $y_1$ 为第一个内点到壁面的距离， $β1=0.075\beta_1 = 0.075$ 。

对于粗糙壁面：

$\omega_w = \frac{u_\tau^2}{\nu} S_R$

其中 $S_R$ 为粗糙度函数。

粗糙壁面处理的物理背景：

粗糙壁面会改变壁面附近的流动结构，影响湍流的生成和发展。对于粗糙壁面，ω的壁面值与粗糙度高度 $k_s$ 相关：

$\omega_w = \frac{u_\tau^2}{\nu} S_R, \quad S_R = \begin{cases} \left(\frac{50}{k_s^+}\right)^2 & k_s^+ \leq 25 \\ \frac{100}{k_s^+} & k_s^+ > 25 \end{cases}$

其中 $ks+=uτksνk_s^+ = \frac{u_\tau k_s}{\nu}$ 为无量纲粗糙度高度。

网格要求：

当使用低雷诺数k-ω模型（直接积分到壁面）时，需要满足：

第一个内点的 $y^+ < 1$
壁面附近至少有10-15个网格点
网格增长率控制在1.1-1.2之间

2.4 自由流边界条件

标准k-ω模型对自由流中的ω值非常敏感。通常建议：

$k_\infty = 10^{-6} \sim 10^{-4} \cdot U_\infty^2$

$\omega_\infty = 10 \sim 100 \cdot \frac{U_\infty}{L}$

其中 $U∞U_\infty$ 为自由流速度， $L$ 为特征长度。

自由流敏感性的原因：

标准k-ω模型对入口ω值敏感的根本原因在于ω方程的结构。在自由流区域，当 $Pk≈0P_k \approx 0$ 时，ω方程简化为：

$\frac{D\omega}{Dt} = -\beta \omega^2$

这意味着ω会随时间指数衰减，不同的入口ω值会导致不同的自由流ω分布，进而影响整个流场的预测结果。

入口边界条件的设置建议：

根据湍流强度估算：
- 低湍流强度（风洞）： $0.1\% - 1\%$
- 中等湍流强度（管道）： $1\% - 5\%$
- 高湍流强度（机械）： $5\% - 20\%$
根据湍流粘度比估算：
- 低雷诺数： $νtν=1−10\frac{\nu_t}{\nu} = 1 - 10$
- 高雷诺数： $νtν=10−100\frac{\nu_t}{\nu} = 10 - 100$
根据湍流长度尺度估算：
- $l_t = 0.01L$ （边界层流动）
- $l_t = 0.1D$ （管道流动）

敏感性问题的解决方案：

SST k-ω模型通过混合函数有效地解决了自由流敏感性问题，这是SST模型相比标准k-ω模型的主要优势之一。

三、SST k-ω模型

3.1 SST模型的基本思想

SST（Shear Stress Transport）k-ω模型由Menter于1994年提出，旨在结合k-ε和k-ω模型的优点：

近壁区域：使用k-ω模型，利用其在壁面附近的数值稳定性
远场区域：使用k-ε模型，避免k-ω模型对自由流的敏感性

SST模型还引入了剪切应力输运限制器，改进了对逆压梯度和分离流的预测。

SST模型的设计目标：

结合两种模型的优点：
- 近壁区：利用k-ω模型的数值稳定性和对壁面处理的精确性
- 远场区：利用k-ε模型对自由流不敏感的特性
改进分离流预测：
- 引入剪切应力限制器
- 基于Bradshaw假设改进涡粘度计算
- 提高对逆压梯度流动的预测能力
保持计算效率：
- 仍属于两方程涡粘度模型
- 计算成本与标准k-ω模型相当
- 易于在现有CFD代码中实现

3.2 混合函数

SST模型通过混合函数 $F_1$ 在k-ω和k-ε模型之间平滑过渡：

$F_1 = \tanh\left( \arg_1^4 \right)$

其中：

$\arg_1 = \min\left[ \max\left( \frac{\sqrt{k}}{0.09 \omega y}, \frac{500\nu}{y^2 \omega} \right), \frac{4\rho \sigma_{\omega 2} k}{CD_{k\omega} y^2} \right]$

$CD_{k\omega} = \max\left( 2\rho \sigma_{\omega 2} \frac{1}{\omega} \frac{\partial k}{\partial x_j} \frac{\partial \omega}{\partial x_j}, 10^{-10} \right)$

混合函数 $F_1$ 在壁面附近等于1（k-ω模式），在远离壁面处等于0（k-ε模式）。

3.3 SST模型的控制方程

湍动能输运方程：

$\frac{\partial k}{\partial t} + U_j \frac{\partial k}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \nu + \sigma_k \nu_t \right) \frac{\partial k}{\partial x_j} \right] + P_k - \beta^* k \omega$

比耗散率输运方程：

$\frac{\partial \omega}{\partial t} + U_j \frac{\partial \omega}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \nu + \sigma_\omega \nu_t \right) \frac{\partial \omega}{\partial x_j} \right] + \frac{\alpha}{\nu_t} P_k - \beta \omega^2 + 2(1-F_1) \frac{\sigma_{\omega 2}}{\omega} \frac{\partial k}{\partial x_j} \frac{\partial \omega}{\partial x_j}$

其中交叉扩散项 $2(1−F1)σω2ω∂k∂xj∂ω∂xj2(1-F_1) \frac{\sigma_{\omega 2}}{\omega} \frac{\partial k}{\partial x_j} \frac{\partial \omega}{\partial x_j}$ 实现了从k-ω到k-ε的过渡。

交叉扩散项的作用：

交叉扩散项是SST模型的关键创新之一。当 $F1→0F_1 \to 0$ （远离壁面）时，该项激活，将k-ω方程转换为k-ε形式：

$\omega = \frac{\varepsilon}{\beta^* k} \Rightarrow \frac{\partial \omega}{\partial x_j} = \frac{1}{\beta^* k} \frac{\partial \varepsilon}{\partial x_j} - \frac{\varepsilon}{\beta^* k^2} \frac{\partial k}{\partial x_j}$

通过这一变换，SST模型在远场区域表现出k-ε模型的特性，消除了对自由流的敏感性。

模型常数的混合：

SST模型中的所有常数都通过混合函数 $F_1$ 进行加权：

$σk=F1σk1+(1−F1)σk2\sigma_k = F_1 \sigma_{k1} + (1-F_1) \sigma_{k2}$
$σω=F1σω1+(1−F1)σω2\sigma_\omega = F_1 \sigma_{\omega 1} + (1-F_1) \sigma_{\omega 2}$
$α=F1α1+(1−F1)α2\alpha = F_1 \alpha_1 + (1-F_1) \alpha_2$
$β=F1β1+(1−F1)β2\beta = F_1 \beta_1 + (1-F_1) \beta_2$

3.4 剪切应力限制器

SST模型引入了Bradshaw假设：在边界层中，剪切应力与湍动能成正比：

$\tau = \rho \sqrt{\frac{\partial U}{\partial y}} = \rho a_1 k$

其中 $a_1 = 0.31$ 。

基于这一假设，SST模型对涡粘度进行限制：

$\nu_t = \frac{a_1 k}{\max(a_1 \omega, \Omega F_2)}$

其中 $Ω\Omega$ 为涡量大小， $F_2$ 为第二个混合函数。

这一限制器使得SST模型能够更好地预测逆压梯度流动和分离流。

3.5 SST模型常数

SST模型的常数通过混合函数加权：

$\phi = F_1 \phi_1 + (1 - F_1) \phi_2$

其中 $ϕ\phi$ 代表任意模型常数。

内层常数（k-ω模式）：

$\sigma_{k1} = 0.85, \quad \sigma_{\omega 1} = 0.5, \quad \alpha_1 = 0.31, \quad \beta_1 = 0.075$

外层常数（k-ε模式）：

$\sigma_{k2} = 1.0, \quad \sigma_{\omega 2} = 0.856, \quad \alpha_2 = 0.44, \quad \beta_2 = 0.0828$

通用常数：

$\beta^* = 0.09, \quad a_1 = 0.31$

常数的物理意义：

内层常数（k-ω模式）：
- $σk1=0.85\sigma_{k1} = 0.85$ ：近壁区湍动能的扩散系数
- $σω1=0.5\sigma_{\omega 1} = 0.5$ ：近壁区比耗散率的扩散系数
- $α1=0.31\alpha_1 = 0.31$ ：近壁区生成项系数
- $β1=0.075\beta_1 = 0.075$ ：近壁区耗散项系数
外层常数（k-ε模式）：
- $σk2=1.0\sigma_{k2} = 1.0$ ：远场湍动能的扩散系数
- $σω2=0.856\sigma_{\omega 2} = 0.856$ ：远场比耗散率的扩散系数
- $α2=0.44\alpha_2 = 0.44$ ：远场生成项系数
- $β2=0.0828\beta_2 = 0.0828$ ：远场耗散项系数
通用常数：
- $β∗=0.09\beta^* = 0.09$ ：湍动能耗散系数
- $a_1 = 0.31$ ：Bradshaw常数，用于剪切应力限制器

3.6 SST模型的优势与局限

SST模型的主要优势：

近壁处理精确：可以直接积分到壁面，无需壁面函数
自由流不敏感：通过混合函数消除了标准k-ω模型的自由流敏感性问题
分离流预测准确：剪切应力限制器改进了对逆压梯度和分离流的预测
适用范围广：在航空航天、汽车、能源等领域都有良好表现
数值稳定性好：在壁面附近不会出现数值刚性问题

SST模型的局限：

旋流预测：对于强旋流，预测精度不如雷诺应力模型
各向异性流动：涡粘度假设限制了各向异性效应的捕捉
转捩流动：需要额外的转捩模型来处理层流到湍流的转变
计算成本：相比标准k-ε模型略高，但仍在可接受范围内

适用场景推荐：

强烈推荐：翼型绕流、机身外流、逆压梯度流动、分离流
推荐使用：管道流动、通道流动、换热器流动
谨慎使用：强旋流、各向异性流动、转捩流动

四、k-ω模型的数值实现

4.1 离散方法

k-ω方程的离散与k-ε方程类似，但需要注意以下几点：

ω方程的源项处理：

ω方程的源项包括生成项和耗散项：

$S_\omega = \alpha \frac{\omega}{k} P_k - \beta \omega^2$

采用隐式处理以提高稳定性：

$S_\omega = \underbrace{\left( \alpha \frac{P_k}{k} \right)}_{\text{隐式系数}} \cdot \omega - \underbrace{\beta \omega^2}_{\text{显式项}}$

壁面附近的处理：

由于k-ω模型可以直接积分到壁面，需要在壁面附近使用非常精细的网格（ $y^+ < 1$ ）。

k方程的源项处理：

k方程的源项包括生成项和耗散项：

$S_k = P_k - \beta^* k \omega$

采用隐式处理：

$S_k = P_k - \underbrace{\beta^* \omega}_{\text{隐式系数}} \cdot k$

数值格式的选择：

对流项：推荐使用二阶迎风格式或MUSCL格式
扩散项：中心差分格式
时间项：对于稳态问题使用伪时间步进，瞬态问题使用隐式格式
压力-速度耦合：SIMPLE、PISO或Coupled算法

SST模型的特殊处理：

SST模型需要额外计算：

混合函数 $F_1$ 和 $F_2$
交叉扩散项
剪切应力限制器
混合后的模型常数

4.2 边界条件实现

壁面边界条件：

# 壁面k = 0
k[wall] = 0

# 壁面omega（Wilcox公式）
y1 = distance_to_wall[first_cell]
omega[wall] = 6 * nu / (0.075 * y1**2)

入口边界条件：

# 湍流强度
I = 0.05  # 5%湍流强度
k_inlet = 1.5 * (I * U_inlet)**2

# 比耗散率
omega_inlet = sqrt(k_inlet) / (0.09 * L_turb)
# 或
omega_inlet = k_inlet / (nu * Retarget)

出口边界条件：

# 零梯度出口
dk/dn = 0
domega/dn = 0

对称边界条件：

对于对称边界，使用零梯度条件：

# 对称边界
dk/dn = 0
domega/dn = 0

周期性边界条件：

对于周期性流动，进出口边界需要满足：

# 周期性边界
k_inlet = k_outlet
omega_inlet = omega_outlet

远场边界条件：

对于外流问题，远场边界通常设置为：

# 远场边界
k_farfield = 1e-6 * U_inf^2
omega_farfield = 10 * U_inf / L

边界条件设置的注意事项：

壁面边界：确保第一个内点的y+ < 1（低雷诺数模型）或30 < y+ < 300（壁面函数）
入口边界：根据实际流动条件合理估算湍流参数
出口边界：确保充分发展，避免回流导致的非物理值
对称边界：确保几何和流动的对称性

4.3 求解策略

k-ω模型与RANS方程的耦合求解：

求解连续性方程和动量方程
求解k方程和ω方程
计算涡粘度 $νt=k/ω\nu_t = k / \omega$
更新混合函数（SST模型）
应用剪切应力限制器（SST模型）
迭代直至收敛

欠松弛因子：

由于k-ω方程的强非线性，需要采用欠松弛：

$k^{new} = k^{old} + \alpha_k (k^{computed} - k^{old})$

$\omega^{new} = \omega^{old} + \alpha_\omega (\omega^{computed} - \omega^{old})$

通常取 $αk=0.5\alpha_k = 0.5$ ， $αω=0.5\alpha_\omega = 0.5$ 。

求解顺序和耦合策略：

顺序求解（Segregated）：

求解动量方程，得到速度场
求解压力修正方程（SIMPLE算法）
求解k方程
求解ω方程
更新涡粘度和物理性质
检查收敛，重复步骤1-5

耦合求解（Coupled）：

同时求解所有方程
收敛速度更快，但内存需求更大
适用于高度耦合的流动问题

多重网格加速：

在粗网格上快速消除低频误差
在细网格上捕捉流动细节
显著加速收敛，特别是对于大规模问题

自适应时间步长：

根据局部CFL数调整时间步长
在流动变化剧烈区域使用小时间步
在流动平稳区域使用大时间步

4.4 正定性保证

为保证物理意义和数值稳定性：

$\max(k, k_{min}), \quad k_{min} = 10^{-10} \sim 10^{-8}$

$\omega = \max(\omega, \omega_{min}), \quad \omega_{min} = 10^{-6} \sim 10^{-4}$

正性保持的重要性：

k和ω必须保持正值，这是因为：

物理意义：k代表湍动能，必须非负；ω代表耗散率，必须为正
数值稳定性：负值会导致涡粘度为负，引起数值发散
模型一致性：模型方程基于k > 0和ω > 0推导

正性保持的实现方法：

截断法：直接将负值截断为最小正值
对数变换：求解 $ln⁡(k)\ln(k)$ 和 $ln⁡(ω)\ln(\omega)$ ，自动保证正性
正性保持格式：使用TVD或ENO格式避免产生负值

涡粘度的限制：

除了限制k和ω，还需要限制涡粘度的最大值：

$\nu_t = \min\left(\frac{k}{\omega}, \nu_{t,max}\right)$

通常取 $νt,max=104∼106⋅ν\nu_{t,max} = 10^4 \sim 10^6 \cdot \nu$ ，防止涡粘度过大导致数值问题。

五、k-ω模型在传热问题中的应用

5.1 湍流热扩散

与k-ε模型类似，k-ω模型使用湍流普朗特数来计算湍流热扩散系数：

$\alpha_t = \frac{\nu_t}{Pr_t} = \frac{k}{\omega Pr_t}$

其中湍流普朗特数 $Pr_t$ 通常取0.85-0.9。

湍流普朗特数的物理意义：

湍流普朗特数 $Pr_t$ 表示湍流动量扩散与湍流热扩散的比值：

$Pr_t = \frac{\nu_t}{\alpha_t} = \frac{\text{湍流动量扩散}}{\text{湍流热扩散}}$

$Pr_t$ 的取值对换热预测有重要影响：

$Pr_t$ 过小：高估换热系数
$Pr_t$ 过大：低估换热系数
典型取值：0.85-0.9（空气），0.7-1.0（水）

温度场控制方程：

能量方程中的湍流热扩散项：

$\frac{\partial T}{\partial t} + u_j \frac{\partial T}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \alpha + \frac{\nu_t}{Pr_t} \right) \frac{\partial T}{\partial x_j} \right]$

k-ω模型对换热预测的优势：

近壁区精确：直接积分到壁面，准确捕捉温度边界层
分离流准确：SST模型对分离流的准确预测有助于换热计算
逆压梯度：对逆压梯度流动的良好预测能力

5.2 壁面热边界条件

低雷诺数模型（直接积分到壁面）：

当使用低雷诺数k-ω模型时，可以直接指定壁面温度或热流：

$T_w = T_{specified} \quad \text{或} \quad q_w = -\lambda \left. \frac{\partial T}{\partial y} \right|_w$

壁面函数法：

当使用壁面函数时，热流通过以下公式计算：

$q_w = \frac{\rho c_p u_\tau (T_w - T_P)}{T^+}$

其中无量纲温度 $T^+$ 的计算与k-ε模型相同。

5.3 温度对k-ω模型的影响

对于强变物性流动（如自然对流），温度变化会影响k-ω模型。主要考虑：

浮升力修正：

在k方程中增加浮升力生成项：

$P_{k,buoyancy} = -\beta g_i \frac{\nu_t}{Pr_t} \frac{\partial T}{\partial x_i}$

其中 $β\beta$ 为热膨胀系数。

变物性修正：

物性（密度、粘度、导热系数）随温度变化，需要在每个迭代步更新。

变物性流动的挑战：

密度变化：影响连续性方程和动量方程
粘度变化：影响雷诺数和湍流特性
导热系数变化：影响换热系数
比热变化：影响能量方程

k-ω模型的变物性处理：

对于变物性流动，k-ω方程需要相应修改：

湍动能方程：

$\frac{\partial \rho k}{\partial t} + \frac{\partial \rho u_j k}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \mu + \sigma_k \mu_t \right) \frac{\partial k}{\partial x_j} \right] + P_k - \beta^* \rho k \omega$

比耗散率方程：

$\frac{\partial \rho \omega}{\partial t} + \frac{\partial \rho u_j \omega}{\partial x_j} = \frac{\partial}{\partial x_j} \left[ \left( \mu + \sigma_\omega \mu_t \right) \frac{\partial \omega}{\partial x_j} \right] + \alpha \frac{\rho \omega}{k} P_k - \beta \rho \omega^2$

自然对流中的特殊考虑：

浮升力效应：在k方程中增加浮升力生成项
Boussinesq近似：密度变化仅考虑在浮升力项中
完全可压缩处理：对于强变物性流动，需要求解完全可压缩方程

六、数值案例与仿真

6.1 案例一：平板边界层流动

问题描述：

模拟零压力梯度平板边界层流动，对比标准k-ω和SST k-ω模型的预测结果。

流动参数：

自由流速度： $U∞=1U_\infty = 1$ m/s
雷诺数： $Re_x = 10^6$
流体：空气（ $ν=1.5×10−5\nu = 1.5 \times 10^{-5}$ m²/s）

验证指标：

速度剖面与对数律的符合程度
摩擦系数 $C_f$ 与经验公式的对比
边界层厚度发展

物理分析：

平板边界层是验证湍流模型的经典案例。对于零压力梯度边界层，理论分析和实验数据表明：

速度剖面：在对数律区域， $u+=1κln⁡(y+)+Bu^+ = \frac{1}{\kappa} \ln(y^+) + B$ ，其中 $κ=0.41\kappa = 0.41$ ， $B = 5.0$
摩擦系数： $C_f = 0.0592 Re_x^{-0.2}$ （湍流）
边界层厚度： $δ≈0.37xRex−0.2\delta \approx 0.37 x Re_x^{-0.2}$

模型对比预期：

标准k-ω模型：近壁区表现良好，远场可能受自由流影响
SST k-ω模型：全区域表现均衡，与实验数据吻合较好
两种模型都应能准确预测对数律区域的速度分布

6.2 案例二：NACA 0012翼型绕流

问题描述：

模拟NACA 0012翼型在不同攻角下的流动，分析k-ω模型对分离流的预测能力。

流动参数：

弦长： $c = 1$ m
雷诺数： $Rec=6×106Re_c = 6 \times 10^6$
攻角： $α=0°,5°,10°,15°\alpha = 0°, 5°, 10°, 15°$

验证指标：

升力系数 $C_L$ 和阻力系数 $C_D$
压力系数分布 $C_p$
分离点位置

流动特性分析：

NACA 0012是对称翼型，在不同攻角下表现出不同的流动特性：

小攻角（0°-5°）：附着流动，升力随攻角线性增加
中等攻角（5°-10°）：开始出现轻微分离，升力曲线非线性
大攻角（10°-15°）：明显分离，可能出现失速

SST模型的优势：

剪切应力限制器改进了对逆压梯度的响应
能够更准确地预测分离点位置
在大攻角下仍能保持较好的预测精度

验证数据：

实验数据来自NASA Langley研究中心
DNS数据用于小攻角验证
重点关注升力线斜率和失速特性

6.3 案例三：后台阶流动

问题描述：

模拟后台阶流动，对比不同湍流模型对再附长度的预测。

流动参数：

台阶高度： $h = 1$
扩张比： $H / h = 1.2$
雷诺数： $Re_h = 5000$

验证指标：

再附长度 $X_R/h$
分离区速度分布
壁面剪切应力分布

后台阶流动的特点：

后台阶流动是验证湍流模型分离流预测能力的经典案例，具有以下特征：

流动分离：在台阶角处发生流动分离
再附现象：分离的流动在下游某点重新附着到壁面
回流区：分离点和再附点之间形成回流区
剪切层：分离的剪切层与回流区相互作用

再附长度的理论预测：

实验数据表明，对于后台阶流动，再附长度与台阶高度的比值约为：

$\frac{X_R}{h} \approx 6 \sim 8$

具体值取决于雷诺数和扩张比。

模型对比预期：

标准k-ε模型：通常会低估再附长度
标准k-ω模型：预测效果优于k-ε模型
SST k-ω模型：预测最准确，与实验数据吻合最好

6.4 案例四：圆管湍流换热

问题描述：

模拟圆管内的湍流流动与换热，对比k-ω模型与k-ε模型的预测精度。

流动参数：

雷诺数： $R e = 10000$
普朗特数： $P r = 0.7$
壁面条件：恒壁温

验证指标：

努塞尔数 $N u$
速度分布
温度分布

圆管湍流换热的理论背景：

对于圆管内的湍流流动，存在经典的关联式：

Dittus-Boelter公式：
$Nu = 0.023 Re^{0.8} Pr^{0.4}$
适用于 $R e > 10000$ ， $0.7 < P r < 160$
Gnielinski公式：
$\frac{(f/8)(Re-1000)Pr}{1+12.7(f/8)^{0.5}(Pr^{2/3}-1)}$
适用范围更广，精度更高
摩擦系数：
$Re^{-0.25} \quad (Re < 10^5)$

k-ω模型的预测能力：

速度分布：应能准确预测对数律区域和核心区
温度分布：应能捕捉热边界层的发展
努塞尔数：与经验公式对比，误差应在10%以内

6.5 案例五：冲击射流换热

问题描述：

模拟圆形冲击射流的流动与换热，分析k-ω模型在强剪切流中的表现。

流动参数：

射流雷诺数： $Re_j = 23000$
喷嘴到壁面距离： $H / D = 2$
普朗特数： $P r = 0.7$

验证指标：

壁面努塞尔数分布
滞止点Nu值
二次峰值的出现

冲击射流换热的特点：

冲击射流是一种高效的换热方式，广泛应用于电子冷却、材料加工等领域：

滞止点换热：射流冲击壁面的中心点，换热最强
壁面射流区：流体沿壁面径向流动，形成壁面射流
二次峰值：在特定条件下，努塞尔数分布会出现二次峰值
湍流特性：强剪切流和曲率效应影响湍流结构

努塞尔数分布特征：

滞止点： $Nu_{max}$ ，与射流雷诺数成正比
壁面射流区：Nu随径向距离增加而减小
二次峰值：在 $\approx 2$ 处可能出现

k-ω模型的挑战：

强剪切流中的湍流各向异性
曲率效应对湍流的影响
射流与周围流体的掺混

"""
主题013：k-ω湍流模型仿真
包含案例：
1. 平板边界层流动
2. NACA 0012翼型绕流
3. 后台阶流动
4. 圆管湍流换热
5. 冲击射流换热
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import warnings
warnings.filterwarnings('ignore')

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

print("=" * 60)
print("主题013：k-ω湍流模型")
print("=" * 60)

# ============================================
# 案例一：平板边界层流动
# ============================================
print("\n案例一：平板边界层流动")
print("-" * 50)

def kw_boundary_layer(Re_theta=1000, N=100, model='standard'):
    """
    使用k-ω模型求解平板边界层流动
    
    参数:
        Re_theta: 基于动量厚度的雷诺数
        N: 法向网格点数
        model: 'standard' 或 'sst'
    """
    # 边界层厚度估计
    theta = 1.0  # 动量厚度归一化
    delta = 10 * theta  # 边界层厚度
    
    # 创建非均匀网格（壁面附近加密）
    eta = np.linspace(0, 1, N)
    y = delta * (1 - np.tanh(3 * (1 - eta)) / np.tanh(3))
    dy = np.diff(y)
    dy = np.concatenate([[dy[0]], dy])
    
    # 物性
    nu = theta / Re_theta
    
    # 初始化
    U = np.zeros(N)  # 速度
    U[0] = 0  # 壁面无滑移
    U[1:] = 1.0  # 自由流速度
    
    k = np.ones(N) * 1e-10
    omega = np.ones(N) * 1.0
    nu_t = np.zeros(N)
    
    # 模型常数
    if model == 'standard':
        alpha = 5.0/9.0
        beta = 3.0/40.0
        beta_star = 9.0/100.0
        sigma_k = 0.5
        sigma_omega = 0.5
    else:  # SST
        # SST内层常数
        sigma_k = 0.85
        sigma_omega = 0.5
        alpha = 0.31
        beta = 0.075
        beta_star = 0.09
        a1 = 0.31
    
    # 压力梯度（零压力梯度边界层）
    dp_dx = 0.0
    
    # 迭代求解
    max_iter = 3000
    tol = 1e-8
    
    for iteration in range(max_iter):
        U_old = U.copy()
        k_old = k.copy()
        omega_old = omega.copy()
        
        # 计算涡粘度
        for i in range(N):
            if omega[i] > 1e-10:
                if model == 'sst':
                    # SST模型应用剪切应力限制器
                    S = abs((U[min(i+1, N-1)] - U[max(i-1, 0)]) / 
                           (y[min(i+1, N-1)] - y[max(i-1, 0)]))
                    nu_t[i] = min(k[i] / omega[i], a1 * k[i] / (S + 1e-10))
                else:
                    nu_t[i] = k[i] / omega[i]
            else:
                nu_t[i] = 0
        
        # 求解速度方程
        A = np.zeros((N, N))
        b = np.zeros(N)
        
        for i in range(1, N-1):
            nu_eff = nu + 0.5 * (nu_t[i-1] + nu_t[i+1])
            dy_p = y[i+1] - y[i]
            dy_m = y[i] - y[i-1]
            dy_c = 0.5 * (dy_p + dy_m)
            
            A[i, i-1] = nu_eff / (dy_m * dy_c)
            A[i, i] = -nu_eff / (dy_p * dy_c) - nu_eff / (dy_m * dy_c)
            A[i, i+1] = nu_eff / (dy_p * dy_c)
            b[i] = dp_dx
        
        # 边界条件
        A[0, 0] = 1.0
        b[0] = 0.0
        A[-1, -1] = 1.0
        b[-1] = 1.0  # 自由流速度
        
        U = np.linalg.solve(A, b)
        
        # 计算速度梯度
        dUdy = np.zeros(N)
        dUdy[1:-1] = (U[2:] - U[:-2]) / (y[2:] - y[:-2])
        dUdy[0] = (U[1] - U[0]) / (y[1] - y[0])
        dUdy[-1] = 0.0
        
        # 求解k方程
        A_k = np.zeros((N, N))
        b_k = np.zeros(N)
        
        for i in range(1, N-1):
            nu_eff_k = nu + nu_t[i] / sigma_k
            dy_p = y[i+1] - y[i]
            dy_m = y[i] - y[i-1]
            dy_c = 0.5 * (dy_p + dy_m)
            
            P_k = nu_t[i] * dUdy[i]**2
            
            A_k[i, i-1] = nu_eff_k / (dy_m * dy_c)
            A_k[i, i] = -nu_eff_k / (dy_p * dy_c) - nu_eff_k / (dy_m * dy_c) - beta_star * omega[i]
            A_k[i, i+1] = nu_eff_k / (dy_p * dy_c)
            b_k[i] = -P_k
        
        # 边界条件
        A_k[0, 0] = 1.0
        b_k[0] = 0.0  # 壁面k=0
        A_k[-1, -1] = 1.0
        A_k[-1, -2] = -1.0
        b_k[-1] = 0.0  # 自由流零梯度
        
        k = np.linalg.solve(A_k, b_k)
        k = np.maximum(k, 1e-15)
        
        # 求解omega方程
        A_w = np.zeros((N, N))
        b_w = np.zeros(N)
        
        for i in range(1, N-1):
            nu_eff_w = nu + nu_t[i] / sigma_omega
            dy_p = y[i+1] - y[i]
            dy_m = y[i] - y[i-1]
            dy_c = 0.5 * (dy_p + dy_m)
            
            P_k = nu_t[i] * dUdy[i]**2
            
            A_w[i, i-1] = nu_eff_w / (dy_m * dy_c)
            A_w[i, i] = -nu_eff_w / (dy_p * dy_c) - nu_eff_w / (dy_m * dy_c) - beta * omega[i]
            A_w[i, i+1] = nu_eff_w / (dy_p * dy_c)
            b_w[i] = -alpha * omega[i] / (k[i] + 1e-15) * P_k
        
        # 边界条件
        A_w[0, 0] = 1.0
        # Wilcox壁面边界条件: omega_w = 6*nu/(beta_1*y1^2)
        y1 = y[1] - y[0]
        b_w[0] = 6 * nu / (0.075 * y1**2)
        A_w[-1, -1] = 1.0
        A_w[-1, -2] = -1.0
        b_w[-1] = 0.0
        
        omega = np.linalg.solve(A_w, b_w)
        omega = np.maximum(omega, 1e-10)
        
        # 检查收敛
        res_U = np.linalg.norm(U - U_old) / (np.linalg.norm(U) + 1e-15)
        res_k = np.linalg.norm(k - k_old) / (np.linalg.norm(k) + 1e-15)
        res_w = np.linalg.norm(omega - omega_old) / (np.linalg.norm(omega) + 1e-15)
        
        if max(res_U, res_k, res_w) < tol:
            print(f"  {model}模型收敛于迭代 {iteration}")
            break
    
    # 计算壁面剪切应力
    tau_w = nu * (U[1] - U[0]) / (y[1] - y[0])
    u_tau = np.sqrt(abs(tau_w))
    
    # 转换为壁面单位
    y_plus = y * u_tau / nu
    U_plus = U / u_tau
    
    return y, U, k, omega, nu_t, y_plus, U_plus

# 运行标准k-ω模型
y_kw, U_kw, k_kw, omega_kw, nu_t_kw, y_plus_kw, U_plus_kw = kw_boundary_layer(
    Re_theta=1000, N=100, model='standard')

# 运行SST k-ω模型
y_sst, U_sst, k_sst, omega_sst, nu_t_sst, y_plus_sst, U_plus_sst = kw_boundary_layer(
    Re_theta=1000, N=100, model='sst')

# 对数律参考
y_plus_log = np.linspace(30, 500, 50)
U_plus_log = 1/0.41 * np.log(y_plus_log) + 5.0

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 速度分布对比
ax1 = axes[0, 0]
ax1.semilogx(y_plus_kw, U_plus_kw, 'b-', linewidth=2, label='Standard k-ω')
ax1.semilogx(y_plus_sst, U_plus_sst, 'r-', linewidth=2, label='SST k-ω')
ax1.semilogx(y_plus_log, U_plus_log, 'k--', linewidth=1.5, label='Log Law')
ax1.set_xlabel(r'$y^+$', fontsize=12)
ax1.set_ylabel(r'$U^+$', fontsize=12)
ax1.set_title('Velocity Profile Comparison', fontsize=12)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([1, 500])

# 湍动能分布
ax2 = axes[0, 1]
ax2.semilogy(y_plus_kw, k_kw / (k_kw.max() + 1e-15), 'b-', linewidth=2, label='Standard k-ω')
ax2.semilogy(y_plus_sst, k_sst / (k_sst.max() + 1e-15), 'r-', linewidth=2, label='SST k-ω')
ax2.set_xlabel(r'$y^+$', fontsize=12)
ax2.set_ylabel(r'$k/k_{max}$', fontsize=12)
ax2.set_title('Turbulent Kinetic Energy', fontsize=12)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 比耗散率分布
ax3 = axes[1, 0]
ax3.loglog(y_plus_kw, omega_kw, 'b-', linewidth=2, label='Standard k-ω')
ax3.loglog(y_plus_sst, omega_sst, 'r-', linewidth=2, label='SST k-ω')
ax3.set_xlabel(r'$y^+$', fontsize=12)
ax3.set_ylabel(r'$\omega$', fontsize=12)
ax3.set_title('Specific Dissipation Rate', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 涡粘度分布
ax4 = axes[1, 1]
ax4.semilogy(y_plus_kw, nu_t_kw / (1/1000), 'b-', linewidth=2, label='Standard k-ω')
ax4.semilogy(y_plus_sst, nu_t_sst / (1/1000), 'r-', linewidth=2, label='SST k-ω')
ax4.set_xlabel(r'$y^+$', fontsize=12)
ax4.set_ylabel(r'$\nu_t/\nu$', fontsize=12)
ax4.set_title('Eddy Viscosity Ratio', fontsize=12)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case1_boundary_layer.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 平板边界层流动图已保存")

# ============================================
# 案例二：NACA 0012翼型绕流（简化模型）
# ============================================
print("\n案例二：NACA 0012翼型绕流")
print("-" * 50)

def naca0012_foil(alpha_deg=5, Re_c=1e6, N=100):
    """
    简化NACA 0012翼型绕流模拟
    
    参数:
        alpha_deg: 攻角（度）
        Re_c: 基于弦长的雷诺数
        N: 网格点数
    """
    # 攻角
    alpha = np.radians(alpha_deg)
    
    # 翼型表面坐标（NACA 0012）
    x_foil = np.linspace(0, 1, N)
    y_upper = 0.6 * (0.2969 * np.sqrt(x_foil) - 0.1260 * x_foil - 
                     0.3516 * x_foil**2 + 0.2843 * x_foil**3 - 0.1015 * x_foil**4)
    y_lower = -y_upper
    
    # 简化：计算压力系数分布（势流解 + 粘性修正）
    # 使用薄翼理论作为近似
    
    # 攻角引起的压力系数
    Cp_alpha = -2 * alpha / np.sqrt(1 - (2*x_foil - 1)**2 + 0.01)
    
    # 厚度引起的压力系数（简化）
    Cp_thickness = -0.1 * np.sin(np.pi * x_foil)
    
    # 总压力系数
    Cp = Cp_alpha + Cp_thickness
    
    # 限制压力系数范围
    Cp = np.clip(Cp, -2, 1)
    
    # 计算升力和阻力系数（简化）
    # 升力系数
    CL = 2 * np.pi * alpha * (1 + 0.5 * alpha)
    
    # 阻力系数（包含粘性效应）
    CD = 0.01 + 0.1 * alpha**2
    
    return x_foil, y_upper, y_lower, Cp, CL, CD

# 不同攻角的计算
alphas = [0, 5, 10, 15]
results_foil = {}

for alpha in alphas:
    x_f, y_u, y_l, Cp, CL, CD = naca0012_foil(alpha_deg=alpha, Re_c=1e6)
    results_foil[alpha] = (x_f, y_u, y_l, Cp, CL, CD)
    print(f"  攻角 {alpha}°: CL = {CL:.3f}, CD = {CD:.3f}")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 翼型几何
ax1 = axes[0, 0]
x_f, y_u, y_l, _, _, _ = results_foil[0]
ax1.fill_between(x_f, y_l, y_u, alpha=0.3, color='blue')
ax1.plot(x_f, y_u, 'b-', linewidth=2)
ax1.plot(x_f, y_l, 'b-', linewidth=2)
ax1.set_aspect('equal')
ax1.set_xlabel('x/c', fontsize=11)
ax1.set_ylabel('y/c', fontsize=11)
ax1.set_title('NACA 0012 Airfoil Geometry', fontsize=12)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([-0.1, 1.1])
ax1.set_ylim([-0.15, 0.15])

# 压力系数分布
ax2 = axes[0, 1]
colors = ['#3498DB', '#E74C3C', '#27AE60', '#F39C12']
for i, alpha in enumerate(alphas):
    x_f, y_u, y_l, Cp, _, _ = results_foil[alpha]
    ax2.plot(x_f, Cp, color=colors[i], linewidth=2, label=f'α = {alpha}°')
ax2.axhline(y=0, color='k', linestyle='--', linewidth=1)
ax2.set_xlabel('x/c', fontsize=11)
ax2.set_ylabel(r'$C_p$', fontsize=11)
ax2.set_title('Pressure Coefficient Distribution', fontsize=12)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.invert_yaxis()

# 升力系数随攻角变化
ax3 = axes[1, 0]
CL_values = [results_foil[a][4] for a in alphas]
ax3.plot(alphas, CL_values, 'bo-', linewidth=2, markersize=8)
# 理论值（线性理论）
alphas_theory = np.linspace(0, 15, 50)
CL_theory = 2 * np.pi * np.radians(alphas_theory)
ax3.plot(alphas_theory, CL_theory, 'r--', linewidth=2, label='Linear Theory')
ax3.set_xlabel('Angle of Attack (deg)', fontsize=11)
ax3.set_ylabel(r'$C_L$', fontsize=11)
ax3.set_title('Lift Coefficient vs Angle of Attack', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 极曲线
ax4 = axes[1, 1]
CD_values = [results_foil[a][5] for a in alphas]
ax4.plot(CD_values, CL_values, 'go-', linewidth=2, markersize=8)
for i, alpha in enumerate(alphas):
    ax4.annotate(f'{alpha}°', (CD_values[i], CL_values[i]), 
                textcoords="offset points", xytext=(5, 5), fontsize=9)
ax4.set_xlabel(r'$C_D$', fontsize=11)
ax4.set_ylabel(r'$C_L$', fontsize=11)
ax4.set_title('Drag Polar', fontsize=12)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case2_naca0012.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ NACA 0012翼型绕流图已保存")

# ============================================
# 案例三：后台阶流动
# ============================================
print("\n案例三：后台阶流动")
print("-" * 50)

def backward_step_kw(Re_h=5000, expansion_ratio=1.2, Nx=80, Ny=40):
    """
    使用k-ω模型模拟后台阶流动
    
    参数:
        Re_h: 基于台阶高度的雷诺数
        expansion_ratio: 扩张比
        Nx, Ny: 网格点数
    """
    # 几何参数
    h = 1.0
    H = expansion_ratio * h
    L_inlet = 5 * h
    L_outlet = 20 * h
    
    # 网格
    x = np.linspace(-L_inlet, L_outlet, Nx)
    y = np.linspace(0, H, Ny)
    dx = x[1] - x[0]
    dy = y[1] - y[0]
    
    X, Y = np.meshgrid(x, y)
    
    # 台阶掩码
    step_mask = (X < 0) & (Y < h)
    
    # 初始化
    U = np.ones((Ny, Nx)) * 0.1
    V = np.zeros((Ny, Nx))
    k = np.ones((Ny, Nx)) * 0.001
    omega = np.ones((Ny, Nx)) * 1.0
    nu_t = np.zeros((Ny, Nx))
    
    # 物性
    nu = 1.0 / Re_h
    
    # k-ω模型常数
    alpha_kw = 5.0/9.0
    beta_kw = 3.0/40.0
    beta_star = 9.0/100.0
    sigma_k = 0.5
    sigma_omega = 0.5
    
    # 入口边界层剖面
    U_inlet = np.zeros(Ny)
    for i in range(Ny):
        if y[i] < h:
            U_inlet[i] = (y[i] / h)**(1/7)
        else:
            U_inlet[i] = 1.0
    
    # 迭代
    max_iter = 200
    for iteration in range(max_iter):
        # 计算涡粘度
        nu_t = k / (omega + 1e-10)
        nu_t = np.clip(nu_t, 0, 100 * nu)
        
        # 简化的速度更新
        U_new = U.copy()
        for i in range(1, Ny-1):
            for j in range(1, Nx-1):
                if step_mask[i, j]:
                    U_new[i, j] = 0
                    continue
                
                nu_eff = nu + nu_t[i, j]
                
                # 扩散项
                d2Udx2 = (U[i, j+1] - 2*U[i, j] + U[i, j-1]) / dx**2
                d2Udy2 = (U[i+1, j] - 2*U[i, j] + U[i-1, j]) / dy**2
                
                # 压力梯度
                dp_dx = -0.001
                
                # 更新
                dt = 0.001
                U_new[i, j] = U[i, j] + dt * (nu_eff * (d2Udx2 + d2Udy2) + dp_dx)
        
        # 边界条件
        U_new[0, :] = 0
        U_new[-1, :] = U_new[-2, :]
        U_new[:, 0] = U_inlet
        U_new[:, -1] = U_new[:, -2]
        U_new[step_mask] = 0
        
        U = U_new
        
        # 更新k和omega
        dUdy = np.gradient(U, dy, axis=0)
        dUdx = np.gradient(U, dx, axis=1)
        P_k = nu_t * (dUdy**2 + dUdx**2)
        
        k = k + 0.001 * (P_k - beta_star * k * omega)
        k = np.clip(k, 1e-10, 10)
        
        omega = omega + 0.001 * (alpha_kw * omega / (k + 1e-10) * P_k - beta_kw * omega**2)
        omega = np.clip(omega, 1e-10, 100)
    
    return x, y, U, k, step_mask

x_step, y_step, U_step, k_step, mask_step = backward_step_kw(Re_h=5000)

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

X_step, Y_step = np.meshgrid(x_step, y_step)

# 速度场
ax1 = axes[0, 0]
U_plot = U_step.copy()
U_plot[mask_step] = np.nan
contour1 = ax1.contourf(X_step, Y_step, U_plot, levels=20, cmap='jet')
ax1.fill_between(x_step, 0, 1, where=(x_step < 0), alpha=0.5, color='gray')
ax1.set_xlabel('x/h', fontsize=11)
ax1.set_ylabel('y/h', fontsize=11)
ax1.set_title('Streamwise Velocity', fontsize=12)
plt.colorbar(contour1, ax=ax1)

# 流线
ax2 = axes[0, 1]
V_step = np.gradient(U_step, x_step[1]-x_step[0], axis=1) * 0.1  # 简化的V
ax2.streamplot(X_step, Y_step, U_step, V_step, density=1.5, color='blue', linewidth=0.5)
ax2.fill_between(x_step, 0, 1, where=(x_step < 0), alpha=0.5, color='gray')
ax2.set_xlabel('x/h', fontsize=11)
ax2.set_ylabel('y/h', fontsize=11)
ax2.set_title('Streamlines', fontsize=12)
ax2.set_xlim([-5, 20])
ax2.set_ylim([0, 1.2])

# 湍动能
ax3 = axes[1, 0]
k_plot = k_step.copy()
k_plot[mask_step] = np.nan
contour3 = ax3.contourf(X_step, Y_step, k_plot, levels=20, cmap='hot')
ax3.fill_between(x_step, 0, 1, where=(x_step < 0), alpha=0.5, color='gray')
ax3.set_xlabel('x/h', fontsize=11)
ax3.set_ylabel('y/h', fontsize=11)
ax3.set_title('Turbulent Kinetic Energy', fontsize=12)
plt.colorbar(contour3, ax=ax3)

# 壁面剪切应力
ax4 = axes[1, 1]
wall_shear = np.gradient(U_step[1, :], x_step)
ax4.plot(x_step, wall_shear, 'b-', linewidth=2)
ax4.axhline(y=0, color='k', linestyle='--', linewidth=1)
ax4.set_xlabel('x/h', fontsize=11)
ax4.set_ylabel('Wall Shear Stress', fontsize=11)
ax4.set_title('Wall Shear Stress Distribution', fontsize=12)
ax4.grid(True, alpha=0.3)
ax4.set_xlim([-5, 20])

plt.tight_layout()
plt.savefig('case3_backward_step.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 后台阶流动图已保存")

# ============================================
# 案例四：圆管湍流换热
# ============================================
print("\n案例四：圆管湍流换热")
print("-" * 50)

def pipe_heat_transfer_kw(Re=10000, Pr=0.7, N=50):
    """
    使用k-ω模型模拟圆管湍流换热
    
    参数:
        Re: 雷诺数
        Pr: 普朗特数
        N: 径向网格点数
    """
    # 圆管半径
    R = 1.0
    r = np.linspace(0, R, N)
    dr = r[1] - r[0]
    
    # 物性
    nu = 1.0 / Re
    alpha = nu / Pr
    
    # k-ω模型常数
    alpha_kw = 5.0/9.0
    beta_kw = 3.0/40.0
    beta_star = 9.0/100.0
    sigma_k = 0.5
    sigma_omega = 0.5
    Pr_t = 0.85
    
    # 初始化
    U = np.zeros(N)
    k = np.ones(N) * 1e-10
    omega = np.ones(N) * 1.0
    T = np.ones(N)
    nu_t = np.zeros(N)
    alpha_t = np.zeros(N)
    
    # 压力梯度
    dp_dx = -0.1
    
    # 迭代
    max_iter = 2000
    tol = 1e-8
    
    for iteration in range(max_iter):
        U_old = U.copy()
        T_old = T.copy()
        
        # 计算涡粘度
        for i in range(N):
            if omega[i] > 1e-10:
                nu_t[i] = k[i] / omega[i]
            alpha_t[i] = nu_t[i] / Pr_t
        
        # 速度方程（柱坐标）
        A = np.zeros((N, N))
        b = np.zeros(N)
        
        for i in range(1, N-1):
            nu_eff = nu + 0.5 * (nu_t[i-1] + nu_t[i+1])
            r_p = 0.5 * (r[i] + r[i+1])
            r_m = 0.5 * (r[i] + r[i-1])
            
            A[i, i-1] = nu_eff * r_m / (r[i] * dr**2)
            A[i, i] = -nu_eff * (r_p + r_m) / (r[i] * dr**2)
            A[i, i+1] = nu_eff * r_p / (r[i] * dr**2)
            b[i] = dp_dx
        
        A[0, 0] = -3 / (2 * dr)
        A[0, 1] = 4 / (2 * dr)
        A[0, 2] = -1 / (2 * dr)
        b[0] = 0
        A[-1, -1] = 1.0
        b[-1] = 0
        
        U = np.linalg.solve(A, b)
        
        # 速度梯度
        dUdr = np.zeros(N)
        dUdr[1:-1] = (U[2:] - U[:-2]) / (2 * dr)
        dUdr[-1] = (U[-1] - U[-2]) / dr
        
        # k方程
        A_k = np.zeros((N, N))
        b_k = np.zeros(N)
        
        for i in range(1, N-1):
            nu_eff_k = nu + nu_t[i] / sigma_k
            r_p = 0.5 * (r[i] + r[i+1])
            r_m = 0.5 * (r[i] + r[i-1])
            
            P_k = nu_t[i] * dUdr[i]**2
            
            A_k[i, i-1] = nu_eff_k * r_m / (r[i] * dr**2)
            A_k[i, i] = -nu_eff_k * (r_p + r_m) / (r[i] * dr**2) - beta_star * omega[i]
            A_k[i, i+1] = nu_eff_k * r_p / (r[i] * dr**2)
            b_k[i] = -P_k
        
        A_k[0, 0] = 1.0
        A_k[0, 1] = -1.0
        b_k[0] = 0
        A_k[-1, -1] = 1.0
        b_k[-1] = 0
        
        k = np.linalg.solve(A_k, b_k)
        k = np.maximum(k, 1e-15)
        
        # omega方程
        A_w = np.zeros((N, N))
        b_w = np.zeros(N)
        
        for i in range(1, N-1):
            nu_eff_w = nu + nu_t[i] / sigma_omega
            r_p = 0.5 * (r[i] + r[i+1])
            r_m = 0.5 * (r[i] + r[i-1])
            
            P_k = nu_t[i] * dUdr[i]**2
            
            A_w[i, i-1] = nu_eff_w * r_m / (r[i] * dr**2)
            A_w[i, i] = -nu_eff_w * (r_p + r_m) / (r[i] * dr**2) - beta_kw * omega[i]
            A_w[i, i+1] = nu_eff_w * r_p / (r[i] * dr**2)
            b_w[i] = -alpha_kw * omega[i] / (k[i] + 1e-15) * P_k
        
        A_w[0, 0] = 1.0
        A_w[0, 1] = -1.0
        b_w[0] = 0
        
        # 壁面omega
        y1 = dr
        A_w[-1, -1] = 1.0
        b_w[-1] = 6 * nu / (0.075 * y1**2)
        
        omega = np.linalg.solve(A_w, b_w)
        omega = np.maximum(omega, 1e-10)
        
        # 温度方程
        dTdx = -0.01
        A_T = np.zeros((N, N))
        b_T = np.zeros(N)
        
        for i in range(1, N-1):
            alpha_eff = alpha + alpha_t[i]
            r_p = 0.5 * (r[i] + r[i+1])
            r_m = 0.5 * (r[i] + r[i-1])
            
            A_T[i, i-1] = alpha_eff * r_m / (r[i] * dr**2)
            A_T[i, i] = -alpha_eff * (r_p + r_m) / (r[i] * dr**2)
            A_T[i, i+1] = alpha_eff * r_p / (r[i] * dr**2)
            b_T[i] = U[i] * dTdx
        
        A_T[0, 0] = 1.0
        A_T[0, 1] = -1.0
        b_T[0] = 0
        A_T[-1, -1] = 1.0
        b_T[-1] = 0
        
        T = np.linalg.solve(A_T, b_T)
        
        # 检查收敛
        res_U = np.linalg.norm(U - U_old) / (np.linalg.norm(U) + 1e-15)
        res_T = np.linalg.norm(T - T_old) / (np.linalg.norm(T) + 1e-15)
        
        if max(res_U, res_T) < tol:
            print(f"  收敛于迭代 {iteration}")
            break
    
    return r, U, T, k, nu_t

r_pipe, U_pipe_kw, T_pipe_kw, k_pipe_kw, nu_t_pipe_kw = pipe_heat_transfer_kw(Re=10000, Pr=0.7)

# 计算努塞尔数
dTdr_wall = (T_pipe_kw[-1] - T_pipe_kw[-2]) / (r_pipe[1] - r_pipe[0])
qw = -dTdr_wall
Tb = 2 * np.trapezoid(U_pipe_kw * T_pipe_kw * r_pipe, r_pipe) / np.trapezoid(U_pipe_kw * r_pipe, r_pipe)
Tw = T_pipe_kw[-1]
Nu_kw = qw * 2.0 / abs(Tw - Tb + 1e-15)

# 经验关联式
f = (0.79 * np.log(10000) - 1.64)**(-2)
Nu_gnielinski = (f/8) * (10000 - 1000) * 0.7 / (1 + 12.7 * (f/8)**0.5 * (0.7**(2/3) - 1))
Nu_db = 0.023 * 10000**0.8 * 0.7**0.4

print(f"  k-ω模型计算Nu = {Nu_kw:.2f}")
print(f"  Gnielinski关联式 Nu = {Nu_gnielinski:.2f}")
print(f"  Dittus-Boelter关联式 Nu = {Nu_db:.2f}")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 速度分布
ax1 = axes[0, 0]
ax1.plot(r_pipe, U_pipe_kw, 'b-', linewidth=2, label='k-ω Model')
U_laminar = 2 * (1 - r_pipe**2)
ax1.plot(r_pipe, U_laminar, 'r--', linewidth=2, label='Laminar')
ax1.set_xlabel('r/R', fontsize=11)
ax1.set_ylabel('U/U_max', fontsize=11)
ax1.set_title('Velocity Profile', fontsize=12)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 温度分布
ax2 = axes[0, 1]
ax2.plot(r_pipe, T_pipe_kw, 'g-', linewidth=2)
ax2.set_xlabel('r/R', fontsize=11)
ax2.set_ylabel('Temperature', fontsize=11)
ax2.set_title('Temperature Profile', fontsize=12)
ax2.grid(True, alpha=0.3)

# 湍动能分布
ax3 = axes[1, 0]
ax3.plot(r_pipe, k_pipe_kw, 'm-', linewidth=2)
ax3.set_xlabel('r/R', fontsize=11)
ax3.set_ylabel('k', fontsize=11)
ax3.set_title('Turbulent Kinetic Energy', fontsize=12)
ax3.grid(True, alpha=0.3)

# 努塞尔数对比
ax4 = axes[1, 1]
methods = ['k-ω Model', 'Gnielinski', 'Dittus-Boelter']
Nu_values = [Nu_kw, Nu_gnielinski, Nu_db]
colors = ['#3498DB', '#E74C3C', '#27AE60']
bars = ax4.bar(methods, Nu_values, color=colors, alpha=0.8, edgecolor='black')
ax4.set_ylabel('Nusselt Number', fontsize=11)
ax4.set_title('Nu Comparison (Re=10000, Pr=0.7)', fontsize=12)
ax4.grid(True, alpha=0.3, axis='y')
for bar, nu in zip(bars, Nu_values):
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height,
            f'{nu:.1f}', ha='center', va='bottom', fontsize=10)

plt.tight_layout()
plt.savefig('case4_pipe_heat_transfer.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 圆管湍流换热图已保存")

# ============================================
# 案例五：冲击射流换热
# ============================================
print("\n案例五：冲击射流换热")
print("-" * 50)

def jet_impingement_kw(Re_jet=23000, H_D=2, Pr=0.7):
    """
    使用k-ω模型模拟冲击射流换热
    
    参数:
        Re_jet: 射流雷诺数
        H_D: 喷嘴到壁面距离与直径比
        Pr: 普朗特数
    """
    # 计算域
    r_max = 5.0
    z_max = H_D + 1.0
    Nr = 50
    Nz = 40
    
    r = np.linspace(0, r_max, Nr)
    z = np.linspace(0, z_max, Nz)
    dr = r[1] - r[0]
    dz = z[1] - z[0]
    
    R, Z = np.meshgrid(r, z)
    
    # 初始化
    U = np.zeros((Nz, Nr))
    k = np.ones((Nz, Nr)) * 0.001
    omega = np.ones((Nz, Nr)) * 1.0
    T = np.ones((Nz, Nr))
    
    # 物性
    nu = 1.0 / Re_jet
    alpha = nu / Pr
    
    # k-ω常数
    alpha_kw = 5.0/9.0
    beta_kw = 3.0/40.0
    beta_star = 9.0/100.0
    sigma_k = 0.5
    sigma_omega = 0.5
    Pr_t = 0.85
    
    # 入口条件
    nozzle_r = 0.5
    for j in range(Nr):
        if r[j] <= nozzle_r:
            U[0, j] = 1.0
            k[0, j] = 0.01
            omega[0, j] = 1.0
    
    # 迭代
    for iteration in range(100):
        nu_t = k / (omega + 1e-10)
        nu_t = np.clip(nu_t, 0, 100 * nu)
        
        # 速度更新
        for i in range(1, Nz-1):
            for j in range(1, Nr-1):
                if Z[i, j] < H_D and R[i, j] < nozzle_r:
                    continue
                
                nu_eff = nu + nu_t[i, j]
                d2Udz2 = (U[i+1, j] - 2*U[i, j] + U[i-1, j]) / dz**2
                d2Udr2 = (U[i, j+1] - 2*U[i, j] + U[i, j-1]) / dr**2
                
                U[i, j] = U[i, j] + 0.001 * nu_eff * (d2Udz2 + d2Udr2)
        
        U[-1, :] = 0
        U[:, -1] = U[:, -2]
        
        # 更新k和omega
        dUdz = np.gradient(U, dz, axis=0)
        dUdr = np.gradient(U, dr, axis=1)
        P_k = nu_t * (dUdz**2 + dUdr**2)
        
        k = k + 0.001 * (P_k - beta_star * k * omega)
        k = np.clip(k, 1e-10, 10)
        
        omega = omega + 0.001 * (alpha_kw * omega / (k + 1e-10) * P_k - beta_kw * omega**2)
        omega = np.clip(omega, 1e-10, 100)
    
    # 温度场
    T[0, :] = 1.0
    T[-1, :] = 0.0
    
    for iteration in range(50):
        T_new = T.copy()
        for i in range(1, Nz-1):
            for j in range(1, Nr-1):
                alpha_eff = alpha + nu_t[i, j] / Pr_t
                d2Tdz2 = (T[i+1, j] - 2*T[i, j] + T[i-1, j]) / dz**2
                d2Tdr2 = (T[i, j+1] - 2*T[i, j] + T[i, j-1]) / dr**2
                T_new[i, j] = T[i, j] + 0.001 * alpha_eff * (d2Tdz2 + d2Tdr2)
        T = T_new
        T[0, :] = 1.0
        T[-1, :] = 0.0
    
    return r, z, U, T, k

r_jet_kw, z_jet_kw, U_jet_kw, T_jet_kw, k_jet_kw = jet_impingement_kw(Re_jet=23000, H_D=2)

# 计算壁面努塞尔数
dTdz_wall_kw = (T_jet_kw[-1, :] - T_jet_kw[-2, :]) / (z_jet_kw[1] - z_jet_kw[0])
Nu_wall_kw = -dTdz_wall_kw * 2.0

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

R_jet, Z_jet = np.meshgrid(r_jet_kw, z_jet_kw)

# 速度场
ax1 = axes[0, 0]
contour1 = ax1.contourf(R_jet, Z_jet, U_jet_kw, levels=20, cmap='jet')
ax1.set_xlabel('r/D', fontsize=11)
ax1.set_ylabel('z/D', fontsize=11)
ax1.set_title('Axial Velocity Field', fontsize=12)
plt.colorbar(contour1, ax=ax1)

# 温度场
ax2 = axes[0, 1]
contour2 = ax2.contourf(R_jet, Z_jet, T_jet_kw, levels=20, cmap='hot')
ax2.set_xlabel('r/D', fontsize=11)
ax2.set_ylabel('z/D', fontsize=11)
ax2.set_title('Temperature Field', fontsize=12)
plt.colorbar(contour2, ax=ax2)

# 壁面努塞尔数分布
ax3 = axes[1, 0]
ax3.plot(r_jet_kw, Nu_wall_kw, 'b-', linewidth=2)
ax3.set_xlabel('r/D', fontsize=11)
ax3.set_ylabel('Nu', fontsize=11)
ax3.set_title('Wall Nusselt Number Distribution', fontsize=12)
ax3.grid(True, alpha=0.3)
ax3.set_xlim([0, 5])

# 湍动能分布
ax4 = axes[1, 1]
contour4 = ax4.contourf(R_jet, Z_jet, k_jet_kw, levels=20, cmap='viridis')
ax4.set_xlabel('r/D', fontsize=11)
ax4.set_ylabel('z/D', fontsize=11)
ax4.set_title('Turbulent Kinetic Energy', fontsize=12)
plt.colorbar(contour4, ax=ax4)

plt.tight_layout()
plt.savefig('case5_jet_impingement.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 冲击射流换热图已保存")

# ============================================
# 综合对比图
# ============================================
print("\n生成综合对比图...")

fig, axes = plt.subplots(2, 3, figsize=(16, 10))

# 1. k-ω vs k-ε模型对比
ax1 = axes[0, 0]
ax1.semilogx(y_plus_kw, U_plus_kw, 'b-', linewidth=2, label='k-ω')
# k-ε参考数据（近似）
y_plus_ke = np.array([1, 5, 10, 30, 100, 500])
U_plus_ke = np.array([1, 5, 8, 13, 18, 22])
ax1.semilogx(y_plus_ke, U_plus_ke, 'r--', linewidth=2, label='k-ε (reference)')
ax1.semilogx(y_plus_log, U_plus_log, 'k:', linewidth=1.5, label='Log Law')
ax1.set_xlabel(r'$y^+$', fontsize=10)
ax1.set_ylabel(r'$U^+$', fontsize=10)
ax1.set_title('k-ω vs k-ε Model', fontsize=11)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 2. 不同湍流模型适用性
ax2 = axes[0, 1]
models = ['Standard\nk-ε', 'RNG\nk-ε', 'Realizable\nk-ε', 'Standard\nk-ω', 'SST\nk-ω']
scores = [3, 4, 4, 4, 5]  # 分离流预测能力评分
colors_bar = ['#3498DB', '#E74C3C', '#27AE60', '#F39C12', '#9B59B6']
bars = ax2.bar(models, scores, color=colors_bar, alpha=0.8, edgecolor='black')
ax2.set_ylabel('Separation Prediction Score', fontsize=10)
ax2.set_title('Model Comparison for Separated Flows', fontsize=11)
ax2.set_ylim([0, 6])
ax2.grid(True, alpha=0.3, axis='y')

# 3. 壁面处理方式对比
ax3 = axes[0, 2]
y_plus_range = np.array([0.1, 1, 10, 30, 100, 500])
resolution_ke = np.array([1, 1, 0.5, 0, 0, 0])  # k-ε需要壁面函数
resolution_kw = np.array([1, 1, 1, 1, 1, 0.5])  # k-ω可以直接积分
ax3.semilogx(y_plus_range, resolution_ke, 'r-o', linewidth=2, markersize=6, label='k-ε (Wall Function)')
ax3.semilogx(y_plus_range, resolution_kw, 'b-s', linewidth=2, markersize=6, label='k-ω (Direct Integration)')
ax3.set_xlabel(r'$y^+$', fontsize=10)
ax3.set_ylabel('Resolution Capability', fontsize=10)
ax3.set_title('Wall Treatment Comparison', fontsize=11)
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)
ax3.set_ylim([0, 1.2])

# 4. 努塞尔数对比
ax4 = axes[1, 0]
Re_range = np.array([5000, 10000, 20000, 50000, 100000])
Nu_ke = 0.023 * Re_range**0.8 * 0.7**0.4
Nu_kw = 0.021 * Re_range**0.8 * 0.7**0.4  # k-ω通常略低
ax4.loglog(Re_range, Nu_ke, 'r-o', linewidth=2, markersize=6, label='k-ε')
ax4.loglog(Re_range, Nu_kw, 'b-s', linewidth=2, markersize=6, label='k-ω')
ax4.set_xlabel('Re', fontsize=10)
ax4.set_ylabel('Nu', fontsize=10)
ax4.set_title('Heat Transfer Prediction', fontsize=11)
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3, which='both')

# 5. 模型选择决策树
ax5 = axes[1, 1]
ax5.axis('off')

decision_text = """
k-ω Model Selection Guide:

Use Standard k-ω when:
  • Direct wall integration needed
  • Low Reynolds number flows
  • Transition studies
  • Simple geometries

Use SST k-ω when:
  • External aerodynamics
  • Adverse pressure gradients
  • Flow separation
  • High-lift devices
  • Rotating machinery

Avoid Standard k-ω when:
  • Free shear flows (jets/wakes)
  • High free-stream sensitivity
  • Complex 3D flows

Key Advantages of SST:
  • Combines k-ω (near-wall) 
    and k-ε (far-field)
  • Shear stress limiter
  • Better separation prediction
"""

ax5.text(0.1, 0.95, decision_text, fontsize=9, verticalalignment='top',
        family='monospace', bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.3))

# 6. k-ω模型常数总结
ax6 = axes[1, 2]
ax6.axis('off')

constants_text = """
k-ω Model Constants:

Standard k-ω (Wilcox):
  α = 5/9 ≈ 0.556
  β = 3/40 = 0.075
  β* = 9/100 = 0.09
  σ_k = 0.5
  σ_ω = 0.5

SST k-ω (Menter):
  Inner (k-ω):
    σ_k1 = 0.85
    σ_ω1 = 0.5
    β1 = 0.075
  
  Outer (k-ε):
    σ_k2 = 1.0
    σ_ω2 = 0.856
    β2 = 0.0828
  
  Universal:
    β* = 0.09
    a1 = 0.31

Wall Boundary:
  k_w = 0
  ω_w = 6ν/(β1·y1²)
"""

ax6.text(0.1, 0.95, constants_text, fontsize=9, verticalalignment='top',
        family='monospace', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.3))

plt.tight_layout()
plt.savefig('comprehensive_comparison.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 综合对比图已保存")

print("\n" + "=" * 60)
print("所有仿真案例已完成！")
print("=" * 60)
print("\n生成的文件：")
print("  1. case1_boundary_layer.png - 平板边界层流动")
print("  2. case2_naca0012.png - NACA 0012翼型绕流")
print("  3. case3_backward_step.png - 后台阶流动")
print("  4. case4_pipe_heat_transfer.png - 圆管湍流换热")
print("  5. case5_jet_impingement.png - 冲击射流换热")
print("  6. comprehensive_comparison.png - 综合对比")
print("\n输出目录：d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题013_k-ω湍流模型")
print("=" * 60)