主题012：k-ε湍流模型

摘要

k-ε湍流模型是工程计算流体力学（CFD）中最经典和广泛应用的湍流模型之一。本教程系统地介绍了k-ε湍流模型的理论基础、数学推导、数值实现方法以及在对流换热仿真中的应用。通过理论讲解和Python仿真实现，详细阐述了标准k-ε模型、RNG k-ε模型、可实现k-ε模型等变体的特点和应用场景。教程包含五个完整的仿真案例：平面通道流动、后台阶流动、圆管换热、射流冲击和旋转通道流动。通过本教程的学习，读者将掌握k-ε湍流模型的核心原理，能够针对不同流动问题选择合适的模型变体，并正确设置边界条件和求解参数，为工程湍流换热问题的数值模拟提供坚实的技术基础。

关键词

k-ε湍流模型，雷诺平均，涡粘度，湍动能，耗散率，壁面函数，RNG模型，可实现模型，对流换热，数值模拟

---

一、湍流模型概述

1.1 湍流的本质特征

湍流是流体流动的一种复杂状态，具有以下显著特征：

随机性与不规则性：湍流中流体质点的运动轨迹呈现高度不规则的三维涡旋结构，速度和压力等物理量随时间和空间呈现随机脉动。这种随机性使得湍流无法用确定性方法精确描述，必须采用统计方法。

涡旋的多尺度特性：湍流包含从最大尺度（与流动特征尺度相当）到最小尺度（Kolmogorov尺度）的广泛涡旋谱。大涡从平均流动中获取能量，通过涡旋破碎过程将能量逐级传递给更小尺度的涡旋，最终在最小尺度上通过粘性耗散转化为热能。这种能量级联过程是湍流的核心动力学机制。

强烈的混合与输运：湍流脉动极大地增强了动量、热量和质量输运。湍流扩散系数通常比分子扩散系数大几个数量级，这使得湍流换热效率远高于层流。

相干结构的存在：尽管湍流具有随机性，但其中也存在具有一定空间和时间相关性的相干结构，如条带结构、准流向涡等。这些结构对湍流的产生和维持起着关键作用。

1.2 湍流模拟方法分类

根据对湍流脉动的处理方式，湍流模拟方法可分为三大类：

直接数值模拟（DNS）：直接求解Navier-Stokes方程，不引入任何湍流模型。DNS需要解析从最大尺度到Kolmogorov尺度的所有涡旋，计算网格和时间步长必须足够精细。对于工程应用中的高雷诺数流动，DNS的计算成本极高，目前仅限于简单几何和低雷诺数的基础研究。

大涡模拟（LES）：将涡旋分为大尺度涡和小尺度涡。大尺度涡通过直接求解方程获得，小尺度涡（亚格子尺度）则采用亚格子模型进行模拟。LES的计算成本介于DNS和雷诺平均方法之间，能够捕捉大尺度相干结构，适用于非定常湍流和复杂流动的研究。

雷诺平均Navier-Stokes方法（RANS）：将瞬时变量分解为时均量和脉动量，通过雷诺平均得到RANS方程。RANS方法引入湍流模型来封闭方程组，计算成本最低，是工程应用中最广泛使用的湍流模拟方法。k-ε模型是RANS方法中最经典和广泛应用的湍流模型之一。

三种方法的比较：

特性	DNS	LES	RANS
计算成本	极高	中等	低
网格要求	解析所有尺度	解析大尺度	较粗
时间特性	瞬态	瞬态	稳态/瞬态
适用范围	基础研究	复杂流动	工程应用
精度	最高	高	中等

工程应用选择：

• 需要精确捕捉非定常特性：LES或DNS
• 大规模工程问题：RANS
• 计算资源充足且需要高精度：LES
• 快速设计迭代：RANS

1.3 雷诺平均的基本概念

对于任意瞬时变量φ，可以分解为时均值和脉动值：

$$\varphi =\bar{\varphi }+{\varphi }^{\prime }$$

其中时均值定义为：

$$\bar{\varphi }=\frac{1}{\Delta t}{\int }_t^{t+\Delta t}\varphi dt$$

雷诺平均的性质：

• 线性性：$\overline{a\varphi +b\psi }=a\bar{\varphi }+b\bar{\psi }$
• 时均的时均：$\overline{\bar{\varphi }}=\bar{\varphi }$
• 脉动的时均：$\overline{{\varphi }^{\prime }}=0$
• 乘积的时均：$\overline{\bar{\varphi }{\varphi }^{\prime }}=\bar{\varphi }\overline{{\varphi }^{\prime }}=0$

对Navier-Stokes方程进行雷诺平均，得到RANS方程：

连续性方程：
$$\frac{\partial {\bar{u}}_i}{\partial x_i}=0$$

动量方程：
$$\frac{\partial {\bar{u}}_i}{\partial t}+{\bar{u}}_j\frac{\partial {\bar{u}}_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial \bar{p}}{\partial x_i}+\nu \frac{{\partial }^2{\bar{u}}_i}{\partial x_j\partial x_j}-\frac{\partial \overline{u_i^{\prime }{u}_j^{\prime }}}{\partial x_j}$$

雷诺应力：方程中出现了新的未知量$-\rho \overline{u_i^{\prime }{u}_j^{\prime }}$，称为雷诺应力，代表湍流脉动对平均流动的动量输运。这是RANS方程的封闭问题，需要通过湍流模型来模拟雷诺应力。

1.4 湍流模型的发展历程

早期模型（1960s-1970s）：

• 零方程模型：混合长度理论
• 一方程模型：湍动能输运方程
• 早期k-ε模型：Launder和Spalding的开创性工作

发展阶段（1980s-1990s）：

• 标准k-ε模型的完善
• RNG k-ε模型的提出
• Realizable k-ε模型的出现
• 低雷诺数k-ε模型的发展

现代发展（2000s至今）：

• k-ω SST模型的兴起
• 雷诺应力模型（RSM）的实用化
• LES和DES方法的工程应用
• 机器学习辅助的湍流建模

未来趋势：

• 多尺度方法的结合
• 数据驱动的模型改进
• 异构计算平台上的高效实现
• 不确定性量化和验证确认

1.4 布辛涅斯克假设

布辛涅斯克（Boussinesq）假设将雷诺应力与平均速度梯度联系起来，类似于层流中的粘性应力：

$$-\overline{u_i^{\prime }{u}_j^{\prime }}={\nu }_t\left(\frac{\partial {\bar{u}}_i}{\partial x_j}+\frac{\partial {\bar{u}}_j}{\partial x_i}\right)-\frac{2}{3}k{\delta }_{ij}$$

其中：

• ${\nu }_t$ 为湍流运动粘度（涡粘度）
• $k=\frac{1}{2}\overline{u_i^{\prime }{u}_i^{\prime }}$ 为湍动能
• ${\delta }_{ij}$ 为Kronecker delta

布辛涅斯克假设大大简化了雷诺应力的计算，是涡粘度模型的基础。然而，该假设也存在局限性，如无法准确预测强旋流、二次流等复杂流动中的各向异性效应。

1.5 涡粘度假设的局限性

各向同性假设的局限：

• 假设涡粘度是标量，各方向相同
• 实际湍流往往具有各向异性
• 强剪切流中不同方向的湍流强度差异显著

复杂流动的挑战：

• 旋流：切向和径向的湍流特性不同
• 二次流：垂直于主流方向的流动难以预测
• 分离流：逆压梯度下的湍流特性变化

改进方向：

• 非线性涡粘度模型
• 显式代数雷诺应力模型（EARSM）
• 雷诺应力输运模型（RSM）

工程实践中的应对：

• 根据流动特性选择合适的模型
• 对关键区域进行网格加密
• 结合实验数据验证和修正

二、k-ε湍流模型理论基础

2.1 模型基本思想

k-ε模型属于两方程涡粘度模型，通过求解两个额外的输运方程来确定湍流涡粘度。模型选择湍动能$k$和湍流耗散率$\epsilon$作为特征变量，因为：

• 湍动能$k$表征湍流脉动的强度
• 湍流耗散率$\epsilon$表征湍流能量在最小尺度上的耗散速率
• 这两个量具有明确的物理意义，且测量相对容易

基于量纲分析，湍流涡粘度可以表示为：

$${\nu }_t=C_{\mu }\frac{k^2}{\epsilon }$$

其中$C_{\mu }$为模型常数，通常取0.09。

k和ε的物理意义：

湍动能k：

• 定义：单位质量流体的湍流脉动动能
• 表达式：$k=\frac{1}{2}(\overline{u^{\prime 2}}+\overline{v^{\prime 2}}+\overline{w^{\prime 2}})$
• 量纲：[L²/T²]
• 典型值：自由流中很小，近壁区较大

湍流耗散率ε：

• 定义：湍动能转化为热能的速率
• 物理意义：小尺度涡旋的粘性耗散
• 量纲：[L²/T³]
• 与Kolmogorov尺度相关

能量级联过程：

1. 大尺度涡从平均流获取能量
2. 能量通过涡旋破碎传递给更小尺度
3. 在Kolmogorov尺度上通过粘性耗散
4. ε表征了这一能量传递的速率

2.2 标准k-ε模型

标准k-ε模型由Launder和Spalding于1972年提出，是工程应用中最广泛使用的湍流模型。

湍动能输运方程：

$$\frac{\partial k}{\partial t}+{\bar{u}}_j\frac{\partial k}{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+P_k-\epsilon$$

湍流耗散率输运方程：

$$\frac{\partial \epsilon }{\partial t}+{\bar{u}}_j\frac{\partial \epsilon }{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_{\epsilon 1}\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}\frac{{\epsilon }^2}{k}$$

其中各项的物理意义为：

对流项（左端）：${\bar{u}}_j\frac{\partial k}{\partial x_j}$ 和 ${\bar{u}}_j\frac{\partial \epsilon }{\partial x_j}$ 表示湍流特性随流体输运。

扩散项：分子扩散和湍流扩散共同作用，将$k$和$\epsilon$从高浓度区域输运到低浓度区域。

生成项$P_k$：湍动能由平均速度梯度产生：
$$P_k={\nu }_t\left(\frac{\partial {\bar{u}}_i}{\partial x_j}+\frac{\partial {\bar{u}}_j}{\partial x_i}\right)\frac{\partial {\bar{u}}_i}{\partial x_j}$$

耗散项：$-\epsilon$表示湍动能的粘性耗散；$C_{\epsilon 2}\frac{{\epsilon }^2}{k}$表示湍流耗散率的衰减。

模型常数：
$$C_{\mu }=0.09,C_{\epsilon 1}=1.44,C_{\epsilon 2}=1.92,{\sigma }_k=1.0,{\sigma }_{\epsilon }=1.3$$

常数的物理意义和标定：

这些模型常数是通过对大量实验数据的拟合得到的：

• $C_{\mu }=0.09$：基于对数律区域的湍流强度
• $C_{\epsilon 1}=1.44$：控制生成项的强度
• $C_{\epsilon 2}=1.92$：控制耗散项的衰减
• ${\sigma }_k=1.0$：湍动能的Prandtl数
• ${\sigma }_{\epsilon }=1.3$：耗散率的Prandtl数

适用范围和限制：

• 适用于高雷诺数湍流
• 对简单剪切流预测较好
• 对强旋流、分离流预测偏差较大
• 需要配合壁面函数使用

2.3 RNG k-ε模型

RNG（Renormalization Group）k-ε模型由Yakhot和Orszag于1986年利用重整化群理论推导得到。与标准k-ε模型相比，RNG模型具有以下改进：

理论推导基础：RNG模型从Navier-Stokes方程出发，通过系统性的尺度消除过程推导得到，具有更坚实的理论基础。

变涡粘度系数：RNG模型中的$C_{\mu }$不再是常数，而是与流动应变率相关：

$$C_{\mu }=\frac{0.0845}{1+\frac{\eta }{{\eta }_0}}\text{其中}\eta =\frac{k}{\epsilon }\sqrt{2S_{ij}{S}_{ij}},{\eta }_0=4.38$$

附加的生成项：RNG k方程中增加了应变率相关的生成项，能够更好地预测快速应变流动。

低雷诺数修正：RNG模型通过有效粘度概念自动包含低雷诺数效应，在近壁区域表现更好。

改进的耗散方程：RNG ε方程中的系数$C_{\epsilon 1}$和$C_{\epsilon 2}$通过理论推导获得，分别为1.42和1.68。

RNG模型在分离流、旋流和复杂剪切流中的预测精度通常优于标准k-ε模型。

RNG模型的优势：

• 理论基础坚实，从Navier-Stokes方程推导
• 自动包含低雷诺数效应
• 对高应变率流动预测更准确
• 在分离流和回流区表现更好

适用场景：

• 强旋流和旋转流动
• 高应变率流动
• 分离和再附流动
• 复杂剪切流动

2.4 Realizable k-ε模型

Realizable k-ε模型由Shih等人于1995年提出，主要针对标准k-ε模型在强剪切流和旋流中的缺陷进行改进。

可实现性约束：标准k-ε模型在某些流动条件下可能产生非物理的雷诺应力（如负的正应力）。Realizable模型通过约束雷诺应力张量必须满足可实现性条件（如Schwarz不等式）来解决这一问题。

变涡粘度系数：

$$C_{\mu }=\frac{1}{A_0+A_s\frac{k{U}^{\ast }}{\epsilon }}$$

其中$A_0=4.04$，$A_s$和$U^{\ast }$与流动应变率和旋度相关。

改进的耗散方程：

$$\frac{\partial \epsilon }{\partial t}+{\bar{u}}_j\frac{\partial \epsilon }{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_{1}S\epsilon -C_2\frac{{\epsilon }^2}{k+\sqrt{\nu \epsilon }}$$

其中$C_1=\max \left(0.43,\frac{\eta }{\eta +5}\right)$，$\eta =S\frac{k}{\epsilon }$，$S$为平均应变率张量的模。

Realizable模型在平面射流、圆射流、边界层分离流等流动中的预测精度显著优于标准模型。

Realizable模型的优势：

• 满足可实现性条件，避免非物理结果
• 对强剪切流和旋流预测更准确
• 在射流流动中表现优异
• 计算稳定性良好

适用场景：

• 射流流动（平面射流、圆射流）
• 旋流和旋转机械
• 强逆压梯度流动
• 边界层分离流动

2.5 模型对比与适用性

特性	标准k-ε	RNG k-ε	Realizable k-ε
理论基础	经验拟合	RNG理论	可实现性约束
涡粘度	常数	变量	变量
旋转流动	一般	较好	好
分离流动	一般	较好	好
射流预测	一般	较好	好
计算稳定性	好	好	较好
计算成本	低	低	低

选择建议：

• 一般工程流动：标准k-ε模型
• 强旋流、高应变率流动：RNG k-ε模型
• 射流、分离流、复杂剪切流：Realizable k-ε模型

2.6 k-ε模型的局限性

尽管k-ε模型在工程应用中取得了巨大成功，但它也存在一些固有的局限性：

涡粘度假设的局限：

• 假设雷诺应力与平均应变率成正比
• 无法准确预测各向异性湍流
• 对强旋流和二次流预测不准确

近壁处理的挑战：

• 标准模型需要壁面函数
• 对y+值敏感
• 低雷诺数修正增加了复杂性

特定流动类型的不足：

• 轴对称射流：预测扩散率过高
• 强逆压梯度流动：分离点预测偏差
• 转捩流动：需要额外的转捩模型

改进方向：

• 使用雷诺应力模型（RSM）处理各向异性
• 采用LES或DES方法捕捉非定常特性
• 结合机器学习改进模型参数

2.7 低雷诺数k-ε模型

背景和需求：

• 标准k-ε模型适用于高雷诺数区域
• 近壁区需要解析粘性底层
• 转捩流动需要捕捉层流到湍流的转变

阻尼函数：
低雷诺数k-ε模型通过引入阻尼函数来修正近壁行为：

$${\nu }_t=C_{\mu }{f}_{\mu }\frac{k^2}{\epsilon }$$

其中$f_{\mu }$为阻尼函数，常见的形式包括：

• Lam-Bremhorst模型：$f_{\mu }={\left[1-\exp (-0.0165R{e}_y)\right]}^2\left(1+\frac{20.5}{R{e}_t}\right)$
• Launder-Sharma模型：$f_{\mu }=\exp \left(\frac{-3.4}{(1+R{e}_t/50{)}^2}\right)$

近壁处理：

• 需要y+ < 1的网格分辨率
• 直接积分到壁面，不使用壁面函数
• 能捕捉粘性底层的速度分布

应用范围：

• 低雷诺数流动
• 转捩流动
• 需要精确预测近壁换热的场合

三、k-ε模型的数值实现

3.1 控制方程的离散

k-ε模型需要与RANS方程耦合求解。以稳态不可压流动为例，控制方程组包括：

连续性方程：
$$\frac{\partial u_i}{\partial x_i}=0$$

动量方程：
$$\frac{\partial }{\partial x_j}(u_j{u}_i)=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[(\nu +{\nu }_t)\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\right]$$

k方程：
$$\frac{\partial }{\partial x_j}(u_{j}k)=\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+P_k-\epsilon$$

ε方程：
$$\frac{\partial }{\partial x_j}(u_j\epsilon )=\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+C_{\epsilon 1}\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}\frac{{\epsilon }^2}{k}$$

3.2 源项的处理

k-ε方程中的源项具有特殊性，需要特别处理：

湍动能生成项$P_k$：

$$P_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}$$

生成项始终为正，代表从平均流动向湍流传递能量。在数值实现中，$P_k$需要限制为非负值，且通常需要限制$P_k$与$\epsilon$的比值（如$P_k\le 10\epsilon$）以保证数值稳定性。

湍动能耗散项$-\epsilon$：该项始终为负，代表湍动能的耗散。在数值求解中，耗散项通常采用隐式处理以提高稳定性：

$$-\epsilon =-\frac{\epsilon }{k}\cdot k=S_k\cdot k$$

其中$S_k=-\frac{\epsilon }{k}$作为k方程的隐式源项系数。

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

$$S_{\epsilon }=C_{\epsilon 1}\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}\frac{{\epsilon }^2}{k}$$

同样采用隐式处理：

$$S_{\epsilon }=\underset{\text{隐式系数}}{\underbrace{\left(C_{\epsilon 1}\frac{P_k}{k}\right)}}\cdot \epsilon -\underset{\text{显式项}}{\underbrace{C_{\epsilon 2}\frac{{\epsilon }^2}{k}}}$$

3.3 近壁面处理

k-ε模型是高雷诺数模型，在近壁区域（粘性底层和缓冲层）需要特殊处理。

壁面函数法：

壁面函数法不在粘性底层布置网格，而是将第一个内点布置在对数律区域（$y^{+}>30$）。壁面处的剪切应力和热流通过壁面函数计算：

速度壁面函数：

$$u^{+}=\{\begin{array}{ll}{y}^{+}& y^{+}\le 11.225\\ \frac{1}{\kappa }\ln (E{y}^{+})& y^{+}>11.225\end{array}$$

其中$\kappa =0.4187$，$E=9.793$。

k和ε的壁面边界条件：

$$\frac{\partial k}{\partial n}=0,{\epsilon }_w=\frac{C_{\mu }^{3/4}{k}_w^{3/2}}{\kappa y_P}$$

其中$y_P$为第一个内点到壁面的距离。

低雷诺数修正：

对于需要解析粘性底层的流动（如低雷诺数流动、转捩流动），采用低雷诺数k-ε模型。模型通过引入阻尼函数来修正涡粘度和模型常数：

$${\nu }_t=C_{\mu }{f}_{\mu }\frac{k^2}{\epsilon },f_{\mu }=\exp \left(-\frac{3.4}{(1+R{e}_t/50{)}^2}\right)$$

其中$R{e}_t=\frac{k^2}{\nu \epsilon }$为湍流雷诺数。

壁面函数的类型：

1. 标准壁面函数：基于对数律，适用于高雷诺数流动
2. 非平衡壁面函数：考虑压力梯度和非平衡效应
3. 增强壁面函数：结合两层模型，适用于更广泛的流动条件

壁面函数的局限性：

• 对y+值敏感，需要在30-300范围内
• 不适合强分离流动
• 不适合低雷诺数流动
• 在复杂几何中可能失效

壁面函数的选择建议：

• 简单流动、高雷诺数：标准壁面函数
• 分离流、逆压梯度：非平衡壁面函数
• 需要更高精度：增强壁面函数或低雷诺数模型

其中$R{e}_t=\frac{k^2}{\nu \epsilon }$为湍流雷诺数。

3.4 求解策略

k-ε模型与RANS方程的耦合求解通常采用以下策略：

分离求解法：

1. 求解连续性方程和动量方程，得到速度场和压力场
2. 求解k方程和ε方程，更新湍动能和耗散率
3. 计算涡粘度${\nu }_t=C_{\mu }\frac{k^2}{\epsilon }$
4. 将新的涡粘度代入动量方程，重复迭代直至收敛

欠松弛处理：

由于k-ε方程的非线性和强耦合性，通常需要采用欠松弛处理：

$${\varphi }^{new}={\varphi }^{old}+{\alpha }_{\varphi }({\varphi }^{computed}-{\varphi }^{old})$$

其中${\alpha }_{\varphi }$为欠松弛因子，通常取0.5-0.8。

k和ε的正定性保证：

为保证物理意义和数值稳定性，需要对k和ε进行限制：

• k和ε必须为正
• 通常设置最小值（如$k_{min}=10^{-10}$，${\epsilon }_{min}=10^{-10}$）
• 限制涡粘度的最大值

初始条件设置：

• 进口边界：根据湍流强度和水力直径估算
• 内部场：可以使用层流解或均匀分布作为初始猜测
• 壁面附近：需要合理的初始分布以避免数值问题

进口湍流边界条件：

• 湍流强度$I=\frac{u^{\prime }}{U}$，通常取1%-5%
• 湍流粘度比$\frac{{\nu }_t}{\nu }$，通常取1-10
• 湍流长度尺度$l$，与几何尺寸相关

出口边界条件：

• 充分发展假设：$\frac{\partial k}{\partial n}=0$，$\frac{\partial \epsilon }{\partial n}=0$
• 注意避免回流导致的非物理值
• 必要时使用出口修正

3.5 数值稳定性措施

网格质量要求：

• 近壁区网格需要满足y+要求
• 避免网格长宽比过大
• 网格过渡要平滑

时间步长选择：

• 瞬态模拟需要满足CFL条件
• 稳态模拟可以使用伪时间步进
• 自适应时间步长可以提高效率

收敛判据：

• 残差下降3-4个数量级
• 关键物理量（如阻力、升力、换热系数）不再变化
• 质量守恒和能量守恒满足要求

常见问题处理：

• 负k或ε：使用正性保持格式
• 发散：减小欠松弛因子，改善网格质量
• 振荡：增加人工耗散，使用高阶格式

数值格式选择：

• 对流项：二阶迎风格式或高阶格式
• 扩散项：中心差分格式
• 时间项：隐式格式或Crank-Nicolson格式
• 压力-速度耦合：SIMPLE、PISO或Coupled算法

计算资源优化：

• 使用多重网格加速收敛
• 采用并行计算提高效率
• 自适应网格细化关键区域
• 合理设置收敛判据平衡精度和效率

欠松弛处理：

由于k-ε方程的非线性和强耦合性，通常需要采用欠松弛处理：

$${\varphi }^{new}={\varphi }^{old}+{\alpha }_{\varphi }({\varphi }^{computed}-{\varphi }^{old})$$

其中${\alpha }_{\varphi }$为欠松弛因子，通常取0.5-0.8。

k和ε的正定性保证：

为保证物理意义和数值稳定性，需要对k和ε进行限制：

$$k=\max (k,k_{min}),\epsilon =\max (\epsilon ,{\epsilon }_{min})$$

通常取$k_{min}=10^{-10}\sim 10^{-8}$，${\epsilon }_{min}=10^{-15}\sim 10^{-12}$。

四、k-ε模型在传热问题中的应用

4.1 湍流热扩散模型

在湍流流动中，热量输运不仅通过分子导热，还通过湍流脉动进行。类似于动量输运，引入湍流热扩散系数：

$${\alpha }_t=\frac{{\nu }_t}{P{r}_t}$$

其中$P{r}_t$为湍流普朗特数，通常取0.85-0.9。

能量方程：

$$\frac{\partial }{\partial x_j}(u_{j}T)=\frac{\partial }{\partial x_j}\left[(\alpha +{\alpha }_t)\frac{\partial T}{\partial x_j}\right]+S_T$$

其中$\alpha =\frac{\nu }{Pr}$为分子热扩散系数。

4.2 壁面热边界条件

壁面函数法的热边界条件：

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

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

其中$T^{+}$为无量纲温度：

$$T^{+}=\{\begin{array}{ll}Pr\cdot y^{+}& y^{+}\le y_T^{+}\\ P{r}_t\left[\frac{1}{\kappa }\ln (E{y}^{+})+P\right]& y^{+}>y_T^{+}\end{array}$$

其中$P$为Peclet数相关的修正项：

$$P=9.24\left[{\left(\frac{Pr}{P{r}_t}\right)}^{3/4}-1\right]\left[1+0.28\exp \left(-0.007\frac{Pr}{P{r}_t}\right)\right]$$

热边界层厚度：

湍流热边界层厚度${\delta }_t$与速度边界层厚度$\delta$的关系取决于普朗特数：

$$\frac{{\delta }_t}{\delta }\approx P{r}^{-1/3}$$

对于高普朗特数流体（如油、水），热边界层比速度边界层薄得多，需要更精细的网格来解析温度梯度。

4.3 湍流换热关联式

基于k-ε模型的数值模拟结果，可以建立湍流换热的无量纲关联式。对于管内湍流换热，常用的关联式包括：

Dittus-Boelter公式：

$$Nu=0.023\cdot R{e}^{0.8}\cdot P{r}^n$$

其中$n=0.4$（加热）或$0.3$（冷却），适用范围：$Re>10000$，$0.7<Pr<160$。

Gnielinski公式：

$$Nu=\frac{(f/8)(Re-1000)Pr}{1+12.7(f/8{)}^{1/2}(P{r}^{2/3}-1)}$$

其中$f=(0.79\ln Re-1.64{)}^{-2}$为摩擦系数，适用范围更广。

Sieder-Tate修正：

考虑变物性影响：

$$Nu=0.027\cdot R{e}^{0.8}\cdot P{r}^{1/3}\cdot {\left(\frac{{\mu }_b}{{\mu }_w}\right)}^{0.14}$$

其中${\mu }_b$和${\mu }_w$分别为体平均温度和壁温下的动力粘度。

4.4 浮力对湍流换热的影响

在自然对流或混合对流中，浮力显著影响湍流结构和换热特性。

浮力修正的k-ε模型：

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

$$G_k=-\beta g_i\frac{{\nu }_t}{P{r}_t}\frac{\partial T}{\partial x_i}$$

其中$\beta$为热膨胀系数，$g_i$为重力加速度分量。

Richardson数：

表征浮力与惯性力的相对重要性：

$$Ri=\frac{Gr}{R{e}^2}=\frac{g\beta \Delta TL}{U^2}$$

• $Ri<0.1$：强制对流主导
• $0.1<Ri<10$：混合对流
• $Ri>10$：自然对流主导

应用实例：

• 室内自然对流换热
• 太阳能集热器
• 核反应堆冷却
• 大气边界层流动

$$Nu=0.027\cdot R{e}^{0.8}\cdot P{r}^{1/3}\cdot {\left(\frac{{\mu }_b}{{\mu }_w}\right)}^{0.14}$$

其中${\mu }_b$和${\mu }_w$分别为流体主体和壁面处的动力粘度。

五、数值案例与仿真

5.1 案例一：平面湍流通道流动

问题描述：

考虑两平行平板间的充分发展湍流流动，计算雷诺数$R{e}_{\tau }=180$（基于摩擦速度和通道半高）时的速度分布和湍流统计量。

控制方程：

对于充分发展流动，流向动量方程简化为：

$$0=\nu \frac{d^{2}U}{d{y}^2}-\frac{d\overline{u^{\prime }{v}^{\prime }}}{dy}-\frac{1}{\rho }\frac{dP}{dx}$$

k-ε模型方程：

$$0=\frac{d}{dy}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{dk}{dy}\right]+P_k-\epsilon$$

$$0=\frac{d}{dy}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{d\epsilon }{dy}\right]+C_{\epsilon 1}\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}\frac{{\epsilon }^2}{k}$$

边界条件：

• 壁面（$y=0$）：$U=0$，$k=0$，采用壁面函数处理
• 中心线（$y=H$）：对称边界条件

数值方法：

采用有限差分法，在壁面附近采用非均匀网格（$y^{+}<1$）。

5.2 案例二：后台阶湍流流动

问题描述：

后台阶流动是检验湍流模型预测分离和再附能力的经典算例。计算域包括台阶前的充分发展段和台阶后的分离区。

流动特征：

• 台阶高度$h$
• 入口边界层厚度$\delta$
• 再附长度$X_R$（从台阶到再附点的距离）

k-ε模型的表现：

标准k-ε模型通常低估再附长度，预测值约为实验值的80-90%。RNG和Realizable模型的预测精度有所提高。

改进方法：

• 采用非线性涡粘度模型
• 使用雷诺应力模型（RSM）
• 采用LES或DES方法

5.3 案例三：圆管湍流换热

问题描述：

模拟恒壁温条件下圆管内的湍流流动与换热，计算局部和平均努塞尔数。

控制方程：

在柱坐标系下，能量方程为：

$$\rho c_{p}u\frac{\partial T}{\partial x}=\frac{1}{r}\frac{\partial }{\partial r}\left[r(\lambda +{\lambda }_t)\frac{\partial T}{\partial r}\right]$$

入口段效应：

在热入口段，努塞尔数沿管长变化：

$$N{u}_x=\{\begin{array}{ll}1.3{\left(\frac{x/D}{Re\cdot Pr}\right)}^{-1/3}& \text{层流}\\ 0.036\cdot R{e}^{0.8}\cdot P{r}^{1/3}\cdot {\left(\frac{D}{x}\right)}^{0.8}& \text{湍流}\end{array}$$

充分发展段：

努塞尔数趋于常数：

$$Nu\approx 0.023\cdot R{e}^{0.8}\cdot P{r}^{0.4}$$

5.4 案例四：湍流射流换热

问题描述：

研究自由射流或冲击射流中的湍流换热特性。射流流动具有强剪切、大应变率的特点，对湍流模型是严峻的考验。

射流类型：

1. 自由射流：射流进入静止或同向流动环境
2. 冲击射流：射流垂直冲击壁面，产生滞止点和壁面射流
3. 阵列射流：多孔射流阵列，用于高效换热

k-ε模型的局限性：

• 标准k-ε模型预测的射流扩散率过高
• 轴对称射流的预测误差大于平面射流
• Realizable k-ε模型对射流的预测有显著改善

冲击射流换热特征：

• 滞止点：换热系数最大
• 二次峰值：在r/D ≈ 2处出现局部最大值
• 壁面射流区：换热系数逐渐降低

工程应用：

• 电子器件冷却
• 涡轮叶片冷却
• 材料热处理
• 干燥工艺

5.5 案例五：旋转通道内的湍流换热

问题描述：

旋转通道（如涡轮叶片内部冷却通道）中的流动受科里奥利力和离心力的影响，湍流结构和换热特性与静止通道显著不同。

旋转效应：

• 科里奥利力：$F_c=-2\rho \vec{\Omega }\times \vec{u}$
• 离心力：$F_r=\rho {\Omega }^{2}r$
• 旋转数：$Ro=\frac{\Omega D}{U}$，表征旋转强度

流动特征：

• 压力面：流体被压向壁面，换热增强
• 吸力面：流体远离壁面，换热减弱
• 二次流：产生垂直于主流方向的流动

k-ε模型的旋转修正：

在k方程中增加旋转生成项：

$$G_{rot}=-2\rho {\nu }_t\Omega \left(\frac{\partial u}{\partial y}-\frac{\partial v}{\partial x}\right)$$

应用价值：

• 航空发动机叶片冷却设计

• 旋转机械的热管理

• 提高冷却效率，降低叶片温度

• 标准k-ε模型预测圆形射流的扩展率过大

• RNG和Realizable模型有所改进

• 对于冲击射流，需要精细的近壁网格和适当的湍流模型

努塞尔数分布：

冲击射流的局部努塞尔数分布呈现特征性曲线：

• 滞止点：$N{u}_{max}$
• 二次峰值：由于转捩和湍流增强
• 壁面射流区：逐渐衰减

5.5 案例五：旋转通道内的湍流换热

问题描述：

燃气轮机叶片内部冷却通道通常具有旋转效应。旋转产生的哥氏力和离心力显著影响流动和换热。

旋转效应：

• 哥氏力：$F_c=-2\rho \vec{\Omega }\times \vec{U}$
• 离心力：$F_r=\rho {\Omega }^{2}r$
• 浮升力：由密度梯度产生

流动特征：

• 压力面（向心侧）：流动加速，换热增强
• 吸力面（离心侧）：流动减速，可能出现分离，换热减弱
• 二次流：产生垂直于主流方向的涡旋

k-ε模型的修正：

旋转效应需要在k-ε模型中引入额外项：

$$S_{k,rot}=-2C_{rot}\Omega \cdot \overline{u^{\prime }{v}^{\prime }}$$

其中$C_{rot}$为旋转修正系数。

六、模型验证与误差分析

6.1 验证方法

与DNS数据对比：

DNS数据提供了湍流统计量的精确参考，包括：

• 平均速度分布
• 雷诺应力分量
• 湍动能分布
• 耗散率分布

与实验数据对比：

实验测量包括：

• 激光多普勒测速（LDV）
• 粒子图像测速（PIV）
• 热线/热膜风速仪
• 温度测量（热电偶、红外热像）

无量纲参数验证：

• 摩擦系数$C_f$
• 努塞尔数$Nu$
• 斯坦顿数$St$
• 再附长度$X_R/h$

6.2 常见误差来源

模型误差：

• 布辛涅斯克假设的局限性
• 涡粘度各向同性假设
• 模型常数的通用性

数值误差：

• 网格分辨率不足
• 离散格式的数值扩散
• 收敛判据过松

边界条件误差：

• 入口湍流特性不确定
• 壁面函数的近似性
• 出口边界条件的反射效应

6.3 改进方向

非线性涡粘度模型：

通过引入应变率张量的高阶项来改进雷诺应力的预测：

$$-\overline{u_i^{\prime }{u}_j^{\prime }}={\nu }_t{S}_{ij}-C_{nl}\frac{k^3}{{\epsilon }^2}\left(S_{ik}{S}_{kj}-\frac{1}{3}{S}_{kl}{S}_{lk}{\delta }_{ij}\right)$$

雷诺应力模型（RSM）：

直接求解雷诺应力输运方程，不采用涡粘度假设：

$$\frac{\partial \overline{u_i^{\prime }{u}_j^{\prime }}}{\partial t}+U_k\frac{\partial \overline{u_i^{\prime }{u}_j^{\prime }}}{\partial x_k}=D_{ij}+P_{ij}+{\Phi }_{ij}-{\epsilon }_{ij}$$

RSM能够更好地预测各向异性和复杂流动，但计算成本较高，需要求解6个雷诺应力分量方程。

混合RANS/LES方法：

• DES（Detached Eddy Simulation）：在附面层使用RANS，在分离区使用LES
• SAS（Scale-Adaptive Simulation）：根据局部流动特性自动调整模型
• SBES（Stress-Blended Eddy Simulation）：结合RANS和LES的优势

机器学习方法：

• 使用神经网络修正涡粘度
• 数据驱动的湍流模型
• 结合物理约束的机器学习模型

6.4 不确定性量化

输入不确定性：

• 边界条件的不确定性
• 物性参数的测量误差
• 几何尺寸的制造公差

模型不确定性：

• 模型常数的校准误差
• 模型形式的结构误差
• 网格依赖性

不确定性传播方法：

• 蒙特卡洛方法
• 多项式混沌展开
• 贝叶斯方法

应用价值：

• 评估计算结果的可靠性
• 指导实验设计
• 支持工程决策

RSM能够更好地预测各向异性和复杂流动，但计算成本较高。

混合RANS/LES方法：

• DES（Detached Eddy Simulation）：在附面层使用RANS，在分离区使用LES
• SAS（Scale-Adaptive Simulation）：根据局部流动特性自动调整模型
• SBES（Stress-Blended Eddy Simulation）：结合RANS和LES的优势

七、Python仿真实现

以下是基于Python的k-ε湍流模型仿真代码，实现了平面通道流动、后台阶流动、圆管换热等经典案例。

---

"""
主题012：k-ε湍流模型仿真
包含案例：
1. 平面湍流通道流动
2. 后台阶湍流流动
3. 圆管湍流换热
4. 湍流射流换热
5. 旋转通道内的湍流换热
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_bvp
from scipy.sparse import diags, csr_matrix
from scipy.sparse.linalg import spsolve
import warnings
warnings.filterwarnings('ignore')

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

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

# ============================================
# 案例一：平面湍流通道流动
# ============================================
print("\n案例一：平面湍流通道流动")
print("-" * 50)

def channel_flow_ke_model(Re_tau=180, N=100, model='standard'):
    """
    使用k-ε模型求解平面通道湍流流动
    
    参数:
        Re_tau: 基于摩擦速度的雷诺数
        N: 网格点数
        model: 'standard', 'rng', 'realizable'
    """
    # 通道半高归一化为1
    H = 1.0
    
    # 创建非均匀网格（壁面附近加密）
    eta = np.linspace(0, 1, N)
    # 使用双曲正切分布加密壁面
    y = H * (1 - np.tanh(2.5 * (1 - eta)) / np.tanh(2.5))
    dy = np.diff(y)
    dy = np.concatenate([[dy[0]], dy, [dy[-1]]])
    
    # 初始化变量
    U = np.zeros(N)  # 速度
    k = np.ones(N) * 1e-10  # 湍动能
    epsilon = np.ones(N) * 1e-10  # 耗散率
    nu_t = np.zeros(N)  # 涡粘度
    
    # 模型常数
    if model == 'standard':
        C_mu = 0.09
        C_e1 = 1.44
        C_e2 = 1.92
        sigma_k = 1.0
        sigma_e = 1.3
    elif model == 'rng':
        C_mu = 0.0845
        C_e1 = 1.42
        C_e2 = 1.68
        sigma_k = 0.7179
        sigma_e = 0.7179
    else:  # realizable
        C_mu = 0.09
        C_e1 = 1.44
        C_e2 = 1.9
        sigma_k = 1.0
        sigma_e = 1.2
    
    # 卡门常数
    kappa = 0.41
    
    # 压力梯度（驱动流动）
    # 对于充分发展流动: tau_w = -dP/dx * H
    # Re_tau = u_tau * H / nu = 1 / nu (归一化)
    nu = 1.0 / Re_tau
    dp_dx = -1.0  # 归一化压力梯度
    
    # 迭代求解
    max_iter = 5000
    tol = 1e-8
    
    for iteration in range(max_iter):
        U_old = U.copy()
        k_old = k.copy()
        epsilon_old = epsilon.copy()
        
        # 计算涡粘度
        for i in range(1, N-1):
            if epsilon[i] > 1e-15:
                if model == 'rng':
                    # RNG模型的变C_mu
                    S = abs((U[i+1] - U[i-1]) / (y[i+1] - y[i-1]))
                    eta_rng = k[i] / epsilon[i] * S
                    C_mu_eff = C_mu / (1 + eta_rng / 4.38)
                else:
                    C_mu_eff = C_mu
                nu_t[i] = C_mu_eff * k[i]**2 / epsilon[i]
            else:
                nu_t[i] = 0
        
        # 求解速度方程
        # d/dy[(nu + nu_t) dU/dy] = dp_dx
        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  # 中心线对称
        A[-1, -2] = -1.0
        b[-1] = 0.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方程
        # d/dy[(nu + nu_t/sigma_k) dk/dy] + P_k - epsilon = 0
        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) - epsilon[i] / (k[i] + 1e-15)
            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)  # 保证正定性
        
        # 求解epsilon方程
        A_e = np.zeros((N, N))
        b_e = np.zeros(N)
        
        for i in range(1, N-1):
            nu_eff_e = nu + nu_t[i] / sigma_e
            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_e[i, i-1] = nu_eff_e / (dy_m * dy_c)
            A_e[i, i] = -nu_eff_e / (dy_p * dy_c) - nu_eff_e / (dy_m * dy_c) - C_e2 * epsilon[i] / (k[i] + 1e-15)
            A_e[i, i+1] = nu_eff_e / (dy_p * dy_c)
            b_e[i] = -C_e1 * epsilon[i] / (k[i] + 1e-15) * P_k
        
        # 边界条件
        # 壁面epsilon通过壁面函数计算
        u_tau = np.sqrt(nu * abs(dUdy[0]))
        epsilon[0] = C_mu**0.75 * k[1]**1.5 / (kappa * y[1])
        A_e[0, 0] = 1.0
        b_e[0] = epsilon[0]
        
        A_e[-1, -1] = 1.0
        A_e[-1, -2] = -1.0
        b_e[-1] = 0.0
        
        epsilon = np.linalg.solve(A_e, b_e)
        epsilon = np.maximum(epsilon, 1e-15)
        
        # 检查收敛
        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_e = np.linalg.norm(epsilon - epsilon_old) / (np.linalg.norm(epsilon) + 1e-15)
        
        if max(res_U, res_k, res_e) < tol:
            print(f"  收敛于迭代 {iteration}")
            break
    
    return y, U, k, epsilon, nu_t

# 运行标准k-ε模型
y, U, k, epsilon, nu_t = channel_flow_ke_model(Re_tau=180, N=100, model='standard')

# 转换为壁面单位
Re_tau = 180
nu = 1.0 / Re_tau
u_tau = np.sqrt(nu * abs((U[1] - U[0]) / (y[1] - y[0])))
y_plus = y * u_tau / nu
U_plus = U / u_tau

# DNS参考数据 (Kim et al., 1987)
y_plus_dns = np.array([0.2, 0.5, 1, 2, 5, 10, 20, 30, 50, 70, 100, 150, 180])
U_plus_dns = np.array([0.8, 1.5, 2.5, 4.0, 7.0, 10.0, 13.5, 15.5, 17.5, 18.8, 19.8, 20.5, 20.8])

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

# 速度分布
ax1 = axes[0, 0]
ax1.semilogx(y_plus, U_plus, 'b-', linewidth=2, label='k-ε Model')
ax1.semilogx(y_plus_dns, U_plus_dns, 'ro', markersize=6, label='DNS (Kim et al.)')
# 对数律
y_log = np.linspace(30, 180, 50)
U_log = 1/0.41 * np.log(y_log) + 5.5
ax1.semilogx(y_log, U_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('Mean Velocity Profile', fontsize=12)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([0.5, 200])

# 湍动能分布
ax2 = axes[0, 1]
k_plus = k / u_tau**2
ax2.plot(y_plus, k_plus, 'g-', linewidth=2)
ax2.set_xlabel(r'$y^+$', fontsize=12)
ax2.set_ylabel(r'$k^+$', fontsize=12)
ax2.set_title('Turbulent Kinetic Energy', fontsize=12)
ax2.grid(True, alpha=0.3)

# 涡粘度分布
ax3 = axes[1, 0]
nu_t_plus = nu_t / nu
ax3.semilogy(y_plus, nu_t_plus, 'm-', linewidth=2)
ax3.set_xlabel(r'$y^+$', fontsize=12)
ax3.set_ylabel(r'$\nu_t/\nu$', fontsize=12)
ax3.set_title('Eddy Viscosity Ratio', fontsize=12)
ax3.grid(True, alpha=0.3)

# 耗散率分布
ax4 = axes[1, 1]
epsilon_plus = epsilon * nu / u_tau**4
ax4.semilogy(y_plus, epsilon_plus, 'c-', linewidth=2)
ax4.set_xlabel(r'$y^+$', fontsize=12)
ax4.set_ylabel(r'$\varepsilon^+$', fontsize=12)
ax4.set_title('Dissipation Rate', fontsize=12)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case1_channel_flow.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 平面湍流通道流动图已保存")

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

def backward_facing_step(Re_h=5000, expansion_ratio=1.2):
    """
    后台阶流动模拟（简化模型）
    
    参数:
        Re_h: 基于台阶高度的雷诺数
        expansion_ratio: 扩张比 H/h
    """
    # 几何参数
    h = 1.0  # 台阶高度
    H = expansion_ratio * h  # 下游通道高度
    L_inlet = 5 * h  # 入口段长度
    L_step = 20 * h  # 台阶后长度
    
    # 网格
    Nx = 100
    Ny = 50
    x = np.linspace(-L_inlet, L_step, Nx)
    y = np.linspace(0, H, Ny)
    dx = x[1] - x[0]
    dy = y[1] - y[0]
    
    # 创建台阶几何掩码
    step_mask = np.zeros((Ny, Nx), dtype=bool)
    for i, yi in enumerate(y):
        for j, xj in enumerate(x):
            if xj < 0 and yi < h:
                step_mask[i, j] = True
    
    # 初始化场
    U = np.ones((Ny, Nx)) * 0.1  # 流向速度
    V = np.zeros((Ny, Nx))  # 横向速度
    k = np.ones((Ny, Nx)) * 0.01
    epsilon = np.ones((Ny, Nx)) * 0.01
    nu_t = np.zeros((Ny, Nx))
    
    # 物性
    nu = 1.0 / Re_h
    
    # 模型常数
    C_mu = 0.09
    C_e1 = 1.44
    C_e2 = 1.92
    sigma_k = 1.0
    sigma_e = 1.3
    
    # 入口边界层剖面（幂律分布）
    U_inlet = np.zeros(Ny)
    for i, yi in enumerate(y):
        if yi < h:
            U_inlet[i] = 1.0 * (yi / h)**(1/7)
        else:
            U_inlet[i] = 1.0
    
    # 简单的迭代求解（简化版）
    max_iter = 100
    
    for iteration in range(max_iter):
        # 计算涡粘度
        nu_t = C_mu * k**2 / (epsilon + 1e-15)
        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]
                
                # 对流项（上风格式）
                conv_x = max(U[i, j], 0) * (U[i, j] - U[i, j-1]) / dx
                conv_y = max(V[i, j], 0) * (U[i, j] - U[i-1, j]) / dy
                
                # 扩散项
                diff_x = nu_eff * (U[i, j+1] - 2*U[i, j] + U[i, j-1]) / dx**2
                diff_y = nu_eff * (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 * (-conv_x - conv_y + diff_x + diff_y + 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和epsilon（简化）
        dkdy = np.gradient(k, dy, axis=0)
        dkdx = np.gradient(k, dx, axis=1)
        
        # 生成项
        dUdy = np.gradient(U, dy, axis=0)
        P_k = nu_t * dUdy**2
        
        # 简化的k方程
        k = k + 0.001 * (P_k - epsilon)
        k = np.clip(k, 1e-10, 10)
        
        # 简化的epsilon方程
        epsilon = epsilon + 0.001 * (C_e1 * epsilon / (k + 1e-15) * P_k - C_e2 * epsilon**2 / (k + 1e-15))
        epsilon = np.clip(epsilon, 1e-15, 100)
    
    return x, y, U, V, k, epsilon, step_mask

x, y, U, V, k, epsilon, step_mask = backward_facing_step(Re_h=5000)

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

# 速度场
ax1 = axes[0, 0]
X, Y = np.meshgrid(x, y)
U_plot = U.copy()
U_plot[step_mask] = np.nan
contour1 = ax1.contourf(X, Y, U_plot, levels=20, cmap='jet')
ax1.fill_between(x, 0, 1, where=(x < 0), alpha=0.5, color='gray', label='Step')
ax1.set_xlabel('x/h', fontsize=11)
ax1.set_ylabel('y/h', fontsize=11)
ax1.set_title('Streamwise Velocity', fontsize=12)
ax1.set_aspect('equal')
plt.colorbar(contour1, ax=ax1)

# 流线
ax2 = axes[0, 1]
ax2.streamplot(X, Y, U, V, density=1.5, color='blue', linewidth=0.5)
ax2.fill_between(x, 0, 1, where=(x < 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_aspect('equal')
ax2.set_xlim([-5, 20])
ax2.set_ylim([0, 1.2])

# 湍动能
ax3 = axes[1, 0]
k_plot = k.copy()
k_plot[step_mask] = np.nan
contour3 = ax3.contourf(X, Y, k_plot, levels=20, cmap='hot')
ax3.fill_between(x, 0, 1, where=(x < 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)
ax3.set_aspect('equal')
plt.colorbar(contour3, ax=ax3)

# 再附长度分析
ax4 = axes[1, 1]
# 找到壁面剪切应力变号的点（简化）
wall_shear = np.gradient(U[1, :], x)
reattachment_point = None
for j in range(len(x)):
    if x[j] > 0 and wall_shear[j] > 0:
        reattachment_point = x[j]
        break

# 绘制壁面剪切应力分布
ax4.plot(x, wall_shear, 'b-', linewidth=2)
ax4.axhline(y=0, color='k', linestyle='--', linewidth=1)
if reattachment_point:
    ax4.axvline(x=reattachment_point, color='r', linestyle='--', linewidth=1.5, label=f'Reattachment: x/h ≈ {reattachment_point:.1f}')
ax4.set_xlabel('x/h', fontsize=11)
ax4.set_ylabel('Wall Shear Stress', fontsize=11)
ax4.set_title('Reattachment Analysis', fontsize=12)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)
ax4.set_xlim([-5, 20])

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

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

def pipe_flow_heat_transfer(Re=10000, Pr=0.7, N=50):
    """
    圆管湍流流动与换热
    
    参数:
        Re: 雷诺数
        Pr: 普朗特数
        N: 径向网格点数
    """
    # 圆管半径归一化为1
    R = 1.0
    r = np.linspace(0, R, N)
    dr = r[1] - r[0]
    
    # 物性
    nu = 1.0 / Re
    alpha = nu / Pr
    
    # 模型常数
    C_mu = 0.09
    C_e1 = 1.44
    C_e2 = 1.92
    sigma_k = 1.0
    sigma_e = 1.3
    Pr_t = 0.85  # 湍流普朗特数
    
    # 初始化
    U = np.zeros(N)  # 轴向速度
    k = np.ones(N) * 1e-10
    epsilon = np.ones(N) * 1e-10
    T = np.ones(N)  # 温度
    nu_t = np.zeros(N)
    alpha_t = np.zeros(N)
    
    # 压力梯度
    dp_dx = -0.1
    
    # 迭代求解
    max_iter = 3000
    tol = 1e-8
    
    for iteration in range(max_iter):
        U_old = U.copy()
        T_old = T.copy()
        
        # 计算涡粘度
        for i in range(N):
            if epsilon[i] > 1e-15:
                nu_t[i] = C_mu * k[i]**2 / epsilon[i]
            alpha_t[i] = nu_t[i] / Pr_t
        
        # 求解速度方程（柱坐标）
        # 1/r * d/dr[r*(nu+nu_t)*dU/dr] = dp_dx
        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.0  # 壁面无滑移
        
        U = np.linalg.solve(A, b)
        
        # 计算速度梯度
        dUdr = np.zeros(N)
        dUdr[1:-1] = (U[2:] - U[:-2]) / (2 * dr)
        dUdr[0] = 0
        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) - epsilon[i] / (k[i] + 1e-15)
            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)
        
        # 求解epsilon方程
        A_e = np.zeros((N, N))
        b_e = np.zeros(N)
        
        for i in range(1, N-1):
            nu_eff_e = nu + nu_t[i] / sigma_e
            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_e[i, i-1] = nu_eff_e * r_m / (r[i] * dr**2)
            A_e[i, i] = -nu_eff_e * (r_p + r_m) / (r[i] * dr**2) - C_e2 * epsilon[i] / (k[i] + 1e-15)
            A_e[i, i+1] = nu_eff_e * r_p / (r[i] * dr**2)
            b_e[i] = -C_e1 * epsilon[i] / (k[i] + 1e-15) * P_k
        
        A_e[0, 0] = 1.0
        A_e[0, 1] = -1.0
        b_e[0] = 0
        
        # 壁面epsilon
        epsilon_wall = C_mu**0.75 * k[-2]**1.5 / (0.41 * dr)
        A_e[-1, -1] = 1.0
        b_e[-1] = epsilon_wall
        
        epsilon = np.linalg.solve(A_e, b_e)
        epsilon = np.maximum(epsilon, 1e-15)
        
        # 求解温度方程
        # U * dT/dx = 1/r * d/dr[r*(alpha+alpha_t)*dT/dr]
        # 假设充分发展，dT/dx = 常数
        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, U_pipe, T_pipe, k_pipe, nu_t_pipe = pipe_flow_heat_transfer(Re=10000, Pr=0.7)

# 计算努塞尔数
# Nu = h * D / k = qw * D / (lambda * (Tw - Tb))
dTdr_wall = (T_pipe[-1] - T_pipe[-2]) / (r[1] - r[0])
qw = -dTdr_wall  # 热流
Tb = 2 * np.trapezoid(U_pipe * T_pipe * r, r) / np.trapezoid(U_pipe * r, r)  # 体积平均温度
Tw = T_pipe[-1]
Nu_calc = qw * 2.0 / abs(Tw - Tb + 1e-15)

# Gnielinski关联式
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))

# Dittus-Boelter
Nu_db = 0.023 * 10000**0.8 * 0.7**0.4

print(f"  计算Nu = {Nu_calc:.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, U_pipe, 'b-', linewidth=2, label='Turbulent')
# 层流参考 (抛物线)
U_laminar = 2 * (1 - r**2)
ax1.plot(r, U_laminar, 'r--', linewidth=2, label='Laminar (Parabolic)')
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, T_pipe, '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, nu_t_pipe / (1/10000), 'm-', linewidth=2)
ax3.set_xlabel('r/R', fontsize=11)
ax3.set_ylabel(r'$\nu_t/\nu$', fontsize=11)
ax3.set_title('Eddy Viscosity Ratio', fontsize=12)
ax3.grid(True, alpha=0.3)

# 努塞尔数对比
ax4 = axes[1, 1]
methods = ['Calculated', 'Gnielinski', 'Dittus-Boelter']
Nu_values = [Nu_calc, 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('case3_pipe_heat_transfer.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 圆管湍流换热图已保存")

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

def jet_impingement_heat_transfer(Re_jet=10000, H_D=2, Pr=0.7):
    """
    圆形冲击射流换热
    
    参数:
        Re_jet: 射流出口雷诺数
        H_D: 喷嘴到壁面距离与直径比
        Pr: 普朗特数
    """
    # 计算域
    r_max = 5.0  # 最大径向位置
    z_max = H_D + 1.0  # 最大轴向位置
    Nr = 60
    Nz = 50
    
    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))  # 轴向速度
    V = np.zeros((Nz, Nr))  # 径向速度
    k = np.ones((Nz, Nr)) * 0.001
    epsilon = np.ones((Nz, Nr)) * 0.001
    T = np.ones((Nz, Nr))  # 温度
    
    # 物性
    nu = 1.0 / Re_jet
    alpha = nu / Pr
    
    # 模型常数
    C_mu = 0.09
    C_e1 = 1.44
    C_e2 = 1.92
    sigma_k = 1.0
    sigma_e = 1.3
    Pr_t = 0.85
    
    # 入口条件（喷嘴出口）
    nozzle_radius = 0.5
    for j in range(Nr):
        if r[j] <= nozzle_radius:
            U[0, j] = 1.0  # 均匀入口速度
            k[0, j] = 0.01
            epsilon[0, j] = C_mu**0.75 * k[0, j]**1.5 / (0.07 * nozzle_radius)
    
    # 简化迭代
    for iteration in range(100):
        # 计算涡粘度
        nu_t = C_mu * k**2 / (epsilon + 1e-15)
        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_radius:
                    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
                
                # 压力梯度（简化）
                dp_dz = 0.0
                
                U[i, j] = U[i, j] + 0.001 * (nu_eff * (d2Udz2 + d2Udr2) + dp_dz)
        
        # 边界条件
        U[-1, :] = 0  # 壁面无滑移
        U[:, -1] = U[:, -2]  # 径向外边界
        
        # 更新k和epsilon
        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 - epsilon)
        k = np.clip(k, 1e-10, 10)
        
        epsilon = epsilon + 0.001 * (C_e1 * epsilon / (k + 1e-15) * P_k - C_e2 * epsilon**2 / (k + 1e-15))
        epsilon = np.clip(epsilon, 1e-15, 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, z_jet, U_jet, T_jet, k_jet = jet_impingement_heat_transfer(Re_jet=10000, H_D=2)

# 计算壁面努塞尔数分布
Nr = len(r_jet)
dTdz_wall = (T_jet[-1, :] - T_jet[-2, :]) / (z_jet[1] - z_jet[0])
Nu_wall = -dTdz_wall * 2.0  # 基于直径的Nu

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

R, Z = np.meshgrid(r_jet, z_jet)

# 速度场
ax1 = axes[0, 0]
contour1 = ax1.contourf(R, Z, U_jet, 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)
ax1.set_aspect('equal')
plt.colorbar(contour1, ax=ax1)

# 温度场
ax2 = axes[0, 1]
contour2 = ax2.contourf(R, Z, T_jet, 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)
ax2.set_aspect('equal')
plt.colorbar(contour2, ax=ax2)

# 壁面努塞尔数分布
ax3 = axes[1, 0]
ax3.plot(r_jet, Nu_wall, '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, Z, k_jet, 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)
ax4.set_aspect('equal')
plt.colorbar(contour4, ax=ax4)

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

# ============================================
# 案例五：旋转通道内的湍流换热
# ============================================
print("\n案例五：旋转通道内的湍流换热")
print("-" * 50)

def rotating_channel_flow(Ro=0.2, Re=10000, N=50):
    """
    旋转通道内的湍流流动与换热
    
    参数:
        Ro: 罗斯贝数 (Omega * D / U_bulk)
        Re: 雷诺数
        N: 网格点数
    """
    # 通道高度归一化为1
    H = 1.0
    y = np.linspace(-H/2, H/2, N)
    dy = y[1] - y[0]
    
    # 物性
    nu = 1.0 / Re
    
    # 旋转参数
    # Ro = Omega * D / U_bulk, D = 1
    Omega = Ro  # 假设U_bulk = 1
    
    # 模型常数
    C_mu = 0.09
    C_e1 = 1.44
    C_e2 = 1.92
    sigma_k = 1.0
    sigma_e = 1.3
    
    # 初始化
    U = np.zeros(N)  # 流向速度
    V = np.zeros(N)  # 横向速度（二次流）
    k = np.ones(N) * 0.001
    epsilon = np.ones(N) * 0.001
    nu_t = np.zeros(N)
    
    # 压力梯度
    dp_dx = -0.01
    
    # 迭代求解
    max_iter = 2000
    tol = 1e-8
    
    for iteration in range(max_iter):
        U_old = U.copy()
        V_old = V.copy()
        
        # 计算涡粘度
        for i in range(N):
            if epsilon[i] > 1e-15:
                nu_t[i] = C_mu * k[i]**2 / epsilon[i]
        
        # 求解流向速度U
        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])
            
            A[i, i-1] = nu_eff / dy**2
            A[i, i] = -2 * nu_eff / dy**2
            A[i, i+1] = nu_eff / dy**2
            
            # 压力梯度 + 哥氏力 (2*Omega*V)
            b[i] = dp_dx + 2 * Omega * V[i]
        
        # 边界条件
        A[0, 0] = 1.0
        b[0] = 0.0
        A[-1, -1] = 1.0
        b[-1] = 0.0
        
        U = np.linalg.solve(A, b)
        
        # 求解横向速度V（二次流）
        # 由哥氏力不平衡驱动
        A_v = np.zeros((N, N))
        b_v = np.zeros(N)
        
        for i in range(1, N-1):
            nu_eff = nu + 0.5 * (nu_t[i-1] + nu_t[i+1])
            
            A_v[i, i-1] = nu_eff / dy**2
            A_v[i, i] = -2 * nu_eff / dy**2
            A_v[i, i+1] = nu_eff / dy**2
            
            # 哥氏力 (-2*Omega*U) + 离心力梯度
            b_v[i] = -2 * Omega * (U[i] - np.mean(U))
        
        A_v[0, 0] = 1.0
        b_v[0] = 0.0
        A_v[-1, -1] = 1.0
        b_v[-1] = 0.0
        
        V = np.linalg.solve(A_v, b_v)
        
        # 计算速度梯度
        dUdy = np.gradient(U, dy)
        dVdy = np.gradient(V, dy)
        
        # 求解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
            
            P_k = nu_t[i] * (dUdy[i]**2 + dVdy[i]**2)
            
            A_k[i, i-1] = nu_eff_k / dy**2
            A_k[i, i] = -2 * nu_eff_k / dy**2 - epsilon[i] / (k[i] + 1e-15)
            A_k[i, i+1] = nu_eff_k / dy**2
            b_k[i] = -P_k
        
        A_k[0, 0] = 1.0
        b_k[0] = 0.0
        A_k[-1, -1] = 1.0
        b_k[-1] = 0.0
        
        k = np.linalg.solve(A_k, b_k)
        k = np.maximum(k, 1e-15)
        
        # 求解epsilon方程
        A_e = np.zeros((N, N))
        b_e = np.zeros(N)
        
        for i in range(1, N-1):
            nu_eff_e = nu + nu_t[i] / sigma_e
            
            P_k = nu_t[i] * (dUdy[i]**2 + dVdy[i]**2)
            
            A_e[i, i-1] = nu_eff_e / dy**2
            A_e[i, i] = -2 * nu_eff_e / dy**2 - C_e2 * epsilon[i] / (k[i] + 1e-15)
            A_e[i, i+1] = nu_eff_e / dy**2
            b_e[i] = -C_e1 * epsilon[i] / (k[i] + 1e-15) * P_k
        
        A_e[0, 0] = 1.0
        b_e[0] = 0.0
        A_e[-1, -1] = 1.0
        b_e[-1] = 0.0
        
        epsilon = np.linalg.solve(A_e, b_e)
        epsilon = np.maximum(epsilon, 1e-15)
        
        # 检查收敛
        res_U = np.linalg.norm(U - U_old) / (np.linalg.norm(U) + 1e-15)
        res_V = np.linalg.norm(V - V_old) / (np.linalg.norm(V) + 1e-15)
        
        if max(res_U, res_V) < tol:
            print(f"  收敛于迭代 {iteration}")
            break
    
    return y, U, V, k, nu_t

# 不同罗斯贝数的对比
Ro_values = [0, 0.1, 0.2, 0.5]
results = {}

for Ro in Ro_values:
    print(f"  计算 Ro = {Ro}...")
    y, U_rot, V_rot, k_rot, nu_t_rot = rotating_channel_flow(Ro=Ro, Re=10000)
    results[Ro] = (y, U_rot, V_rot, k_rot, nu_t_rot)

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

# 速度分布对比
ax1 = axes[0, 0]
colors = ['#3498DB', '#E74C3C', '#27AE60', '#F39C12']
for i, Ro in enumerate(Ro_values):
    y, U_rot, _, _, _ = results[Ro]
    ax1.plot(U_rot, y, color=colors[i], linewidth=2, label=f'Ro = {Ro}')
ax1.set_xlabel('U/U_bulk', fontsize=11)
ax1.set_ylabel('y/D', fontsize=11)
ax1.set_title('Streamwise Velocity Profile', fontsize=12)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 二次流分布
ax2 = axes[0, 1]
for i, Ro in enumerate(Ro_values):
    y, _, V_rot, _, _ = results[Ro]
    ax2.plot(V_rot * 10, y, color=colors[i], linewidth=2, label=f'Ro = {Ro} (×10)')
ax2.set_xlabel('V × 10', fontsize=11)
ax2.set_ylabel('y/D', fontsize=11)
ax2.set_title('Secondary Flow (Spanwise Velocity)', fontsize=12)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 湍动能分布
ax3 = axes[1, 0]
for i, Ro in enumerate(Ro_values):
    y, _, _, k_rot, _ = results[Ro]
    ax3.plot(k_rot, y, color=colors[i], linewidth=2, label=f'Ro = {Ro}')
ax3.set_xlabel('k', fontsize=11)
ax3.set_ylabel('y/D', fontsize=11)
ax3.set_title('Turbulent Kinetic Energy', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 哥氏力对换热的影响示意
ax4 = axes[1, 1]
# 示意图：旋转通道内的二次流
ax4.set_xlim([-1, 1])
ax4.set_ylim([-0.6, 0.6])

# 绘制通道
rect = plt.Rectangle((-0.5, -0.5), 1.0, 1.0, fill=False, edgecolor='black', linewidth=2)
ax4.add_patch(rect)

# 绘制旋转方向
ax4.annotate('', xy=(0.7, 0.3), xytext=(0.7, -0.3),
            arrowprops=dict(arrowstyle='->', color='red', lw=2))
ax4.text(0.8, 0, 'Ω', fontsize=14, color='red', fontweight='bold')

# 绘制二次流示意
# 压力面（向心侧）
ax4.annotate('', xy=(0.3, -0.3), xytext=(0.3, 0.3),
            arrowprops=dict(arrowstyle='->', color='blue', lw=1.5))
ax4.text(0.4, 0, '↑', fontsize=12, color='blue')

# 吸力面（离心侧）
ax4.annotate('', xy=(-0.3, 0.3), xytext=(-0.3, -0.3),
            arrowprops=dict(arrowstyle='->', color='blue', lw=1.5))
ax4.text(-0.4, 0, '↓', fontsize=12, color='blue')

# 主流方向
ax4.annotate('', xy=(0.6, 0.4), xytext=(-0.6, 0.4),
            arrowprops=dict(arrowstyle='->', color='green', lw=2))
ax4.text(0, 0.45, 'U (Main Flow)', fontsize=10, color='green', ha='center')

ax4.text(0, -0.65, 'Coriolis-induced Secondary Flow', fontsize=11, ha='center')
ax4.set_aspect('equal')
ax4.axis('off')

plt.tight_layout()
plt.savefig('case5_rotating_channel.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 旋转通道内的湍流换热图已保存")

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

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

# 1. 不同k-ε模型的速度分布对比
ax1 = axes[0, 0]
models = ['standard', 'rng', 'realizable']
model_names = ['Standard', 'RNG', 'Realizable']
colors_model = ['#3498DB', '#E74C3C', '#27AE60']

for model, name, color in zip(models, model_names, colors_model):
    y_m, U_m, _, _, _ = channel_flow_ke_model(Re_tau=180, N=80, model=model)
    u_tau = np.sqrt((1/180) * abs((U_m[1] - U_m[0]) / (y_m[1] - y_m[0])))
    y_plus_m = y_m * u_tau * 180
    U_plus_m = U_m / u_tau
    ax1.semilogx(y_plus_m, U_plus_m, color=color, linewidth=2, label=name)

ax1.semilogx(y_plus_dns, U_plus_dns, 'ko', markersize=4, alpha=0.5, label='DNS')
ax1.set_xlabel(r'$y^+$', fontsize=10)
ax1.set_ylabel(r'$U^+$', fontsize=10)
ax1.set_title('k-ε Model Comparison', fontsize=11)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([0.5, 200])

# 2. 不同雷诺数下的速度分布
ax2 = axes[0, 1]
Re_values = [180, 395, 590]
for Re_tau, color in zip(Re_values, colors):
    y_re, U_re, _, _, _ = channel_flow_ke_model(Re_tau=Re_tau, N=80)
    u_tau_re = np.sqrt((1/Re_tau) * abs((U_re[1] - U_re[0]) / (y_re[1] - y_re[0])))
    y_plus_re = y_re * u_tau_re * Re_tau
    U_plus_re = U_re / u_tau_re
    ax2.semilogx(y_plus_re, U_plus_re, color=color, linewidth=2, label=f'Re_τ = {Re_tau}')

ax2.set_xlabel(r'$y^+$', fontsize=10)
ax2.set_ylabel(r'$U^+$', fontsize=10)
ax2.set_title('Effect of Reynolds Number', fontsize=11)
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 3. 努塞尔数随雷诺数变化
ax3 = axes[0, 2]
Re_range = np.logspace(4, 6, 20)
Nu_db_range = 0.023 * Re_range**0.8 * 0.7**0.4
Nu_gn_range = []
for Re_nu in Re_range:
    f_nu = (0.79 * np.log(Re_nu) - 1.64)**(-2)
    Nu_gn = (f_nu/8) * (Re_nu - 1000) * 0.7 / (1 + 12.7 * (f_nu/8)**0.5 * (0.7**(2/3) - 1))
    Nu_gn_range.append(Nu_gn)

ax3.loglog(Re_range, Nu_db_range, 'b-', linewidth=2, label='Dittus-Boelter')
ax3.loglog(Re_range, Nu_gn_range, 'r--', linewidth=2, label='Gnielinski')
ax3.set_xlabel('Re', fontsize=10)
ax3.set_ylabel('Nu', fontsize=10)
ax3.set_title('Turbulent Heat Transfer Correlations', fontsize=11)
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3, which='both')

# 4. 湍流普朗特数影响
ax4 = axes[1, 0]
Pr_t_values = [0.5, 0.7, 0.9, 1.2]
for Pr_t, color in zip(Pr_t_values, colors):
    # 简化计算
    Nu_Pr = 0.023 * 10000**0.8 * 0.7**0.4 * (0.85/Pr_t)**0.4
    ax4.bar([f'Pr_t={Pr_t}'], [Nu_Pr], color=color, alpha=0.8, edgecolor='black')

ax4.set_ylabel('Nu', fontsize=10)
ax4.set_title('Effect of Turbulent Prandtl Number', fontsize=11)
ax4.grid(True, alpha=0.3, axis='y')

# 5. 旋转效应总结
ax5 = axes[1, 1]
Ro_summary = [0, 0.1, 0.2, 0.3, 0.5]
Nu_leading = []  # 压力面（向心侧）
Nu_trailing = []  # 吸力面（离心侧）

for Ro_s in Ro_summary:
    # 简化模型：旋转增强压力面换热，减弱吸力面换热
    Nu_base = 100
    Nu_leading.append(Nu_base * (1 + 0.5 * Ro_s))
    Nu_trailing.append(Nu_base * (1 - 0.3 * Ro_s))

ax5.plot(Ro_summary, Nu_leading, 'b-o', linewidth=2, markersize=6, label='Leading (Pressure) Side')
ax5.plot(Ro_summary, Nu_trailing, 'r-s', linewidth=2, markersize=6, label='Trailing (Suction) Side')
ax5.axhline(y=100, color='k', linestyle='--', linewidth=1, label='Non-rotating')
ax5.set_xlabel('Ro', fontsize=10)
ax5.set_ylabel('Nu', fontsize=10)
ax5.set_title('Rotation Effect on Heat Transfer', fontsize=11)
ax5.legend(fontsize=9)
ax5.grid(True, alpha=0.3)

# 6. 模型适用性总结
ax6 = axes[1, 2]
ax6.axis('off')

summary_text = """
k-ε Model Summary:

Standard k-ε:
  • Simplest and most widely used
  • Good for simple shear flows
  • Limitations: jets, separations

RNG k-ε:
  • Theoretically derived
  • Better for swirling flows
  • Improved near-wall behavior

Realizable k-ε:
  • Satisfies realizability constraints
  • Best for jets and mixing layers
  • Improved separation prediction

Key Parameters:
  • C_μ = 0.09 (standard)
  • σ_k = 1.0, σ_ε = 1.3
  • Pr_t = 0.85 (turbulent Prandtl)

Wall Treatment:
  • Wall functions: y⁺ > 30
  • Low-Re models: y⁺ < 1
  • Enhanced wall functions: both
"""

ax6.text(0.1, 0.95, summary_text, fontsize=10, 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_channel_flow.png - 平面湍流通道流动")
print("  2. case2_backward_step.png - 后台阶湍流流动")
print("  3. case3_pipe_heat_transfer.png - 圆管湍流换热")
print("  4. case4_jet_impingement.png - 湍流射流换热")
print("  5. case5_rotating_channel.png - 旋转通道内的湍流换热")
print("  6. comprehensive_comparison.png - 综合对比")
print("\n输出目录：d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题012_k-ε湍流模型")
print("=" * 60)