主题014：近壁面处理方法

近壁面处理是计算流体力学（CFD）中最重要的环节之一，直接影响壁面摩擦阻力、传热系数和流动分离的预测精度。本主题将详细介绍近壁面区域的流动特性、壁面函数法、低雷诺数模型、增强壁面处理方法，以及工程实践中的最佳实践。

一、近壁面流动概述

1.1 近壁面区域的重要性

在壁面附近的流动区域，即边界层内，存在着复杂的物理现象，包括：

• 粘性主导效应：壁面无滑移条件导致速度梯度极大
• 湍流各向异性：垂直于壁面方向的脉动受到抑制
• 能量耗散集中：湍动能的生成和耗散在近壁区达到峰值
• 传热关键区域：温度边界层与速度边界层在此发展，壁面换热系数直接取决于近壁面温度梯度

准确模拟近壁面区域对于预测壁面摩擦阻力、传热系数和流动分离至关重要，是CFD仿真成功的关键。

工程应用中的重要性：

1. 航空航天：机翼和机身的摩擦阻力占总阻力的50%以上
2. 能源动力：换热器的换热效率直接取决于壁面传热系数
3. 汽车工程：车身气动性能与边界层特性密切相关
4. 化工过程：反应器内的传质传热受壁面效应显著影响
5. 海洋工程：船舶阻力与船体表面边界层特性密切相关
6. 环境工程：大气边界层模拟需要准确的壁面处理

近壁面区域的流动特征：

• 速度梯度极大：从壁面的零速度迅速增加到自由流速度
• 湍流结构复杂：存在条带结构、准流向涡等相干结构
• 能量转换剧烈：机械能通过粘性耗散转化为热能
• 压力梯度敏感：逆压梯度容易导致流动分离
  
  
  
  
  
  

1.2 边界层结构

近壁面区域可分为三个子层：

粘性底层（Viscous Sublayer）：

• 范围：$y^{+}<5$
• 特征：粘性力主导，速度呈线性分布
• 速度分布：$U^{+}=y^{+}$
• 物理机制：分子粘性主导动量输运

缓冲层（Buffer Layer）：

• 范围：$5<y^{+}<30$
• 特征：粘性和湍流效应同等重要
• 最复杂的区域，难以用简单公式描述
• 物理机制：湍流生成和耗散达到平衡

对数律层（Logarithmic Layer）：

• 范围：$30<y^{+}<500$
• 特征：湍流效应主导，速度呈对数分布
• 速度分布：$U^{+}=\frac{1}{\kappa }\ln (y^{+})+B$
• 物理机制：湍流输运主导动量传递

外层（Outer Layer）：

• 范围：$y^{+}>500$或$y/\delta >0.2$
• 特征：受外流条件影响
• 存在尾迹律（Wake Law）
• 速度亏损律适用

其中$y^{+}$为无量纲壁面距离，定义为：

$$y^{+}=\frac{y{u}_{\tau }}{\nu }$$

$u_{\tau }=\sqrt{{\tau }_w/\rho }$为摩擦速度，${\tau }_w$为壁面剪切应力。

边界层厚度定义：

• 名义厚度$\delta$：速度达到99%自由流速度的位置
• **位移厚度${\delta }^{\ast }$：${\delta }^{\ast }={\int }_0^{\infty }(1-U/U_{\infty })dy$
• 动量厚度$\theta$：$\theta ={\int }_0^{\infty }(U/U_{\infty })(1-U/U_{\infty })dy$
• 形状因子$H$：$H={\delta }^{\ast }/\theta$

边界层发展的影响因素：

1. 压力梯度：顺压梯度加速，逆压梯度减速
2. 表面粗糙度：增加湍流，提前转捩
3. 壁面曲率：凸曲率稳定，凹曲率不稳定
4. 自由流湍流度：影响转捩位置
5. 可压缩性：高速流动中密度和温度变化显著

1.3 近壁面处理的挑战

1. 网格分辨率要求：直接模拟近壁面需要极细的网格
2. 数值刚性：壁面附近变量梯度大，导致数值不稳定
3. 湍流模型适用性：高雷诺数湍流模型在壁面附近失效
4. 计算成本：精细网格显著增加计算量

网格分辨率要求的定量分析：

对于直接数值模拟（DNS）近壁面区域：

• 流向分辨率：$\Delta x^{+}\approx 10-20$
• 法向分辨率：$\Delta y_{min}^{+}<1$，至少10-15个网格点
• 展向分辨率：$\Delta z^{+}\approx 5-10$

对于大涡模拟（LES）：

• 流向分辨率：$\Delta x^{+}\approx 50-100$
• 法向分辨率：$\Delta y_{min}^{+}<1$
• 展向分辨率：$\Delta z^{+}\approx 20-40$

对于雷诺平均（RANS）配合壁面函数：

• 第一个内点：$30<y^{+}<300$
• 壁面法向：5-10个网格点
• 计算成本降低1-2个数量级

数值刚性的来源：

1. 速度梯度：壁面附近$\partial U/\partial y$极大
2. 湍动能梯度：k在壁面附近变化剧烈
3. 耗散率奇异：ε在壁面附近趋于无穷大
4. 时间步长限制：CFL条件要求极小的时间步长

湍流模型的近壁行为：

高雷诺数湍流模型（如标准k-ε模型）在壁面附近存在以下问题：

• 涡粘度计算失效
• 源项出现奇异性
• 边界条件难以确定
• 需要低雷诺数修正或壁面函数

二、壁面函数法

2.1 壁面函数的基本思想

壁面函数法的核心思想是：

1. 粗化壁面网格：第一个内点放置在对数律层（$30<y^{+}<300$）
2. 解析解替代：用半经验公式代替粘性底层的详细计算
3. 边界条件修正：通过壁面函数计算壁面剪切应力和热流

2.2 标准壁面函数

速度壁面函数：

对于$30<y^{+}<500$：

$$U^{+}=\frac{1}{\kappa }\ln (E{y}^{+})$$

其中$\kappa =0.41$为卡门常数，$E=9.8$为粗糙度参数。

对于$y^{+}<5$（低雷诺数修正）：

$$U^{+}=y^{+}$$

对于$5<y^{+}<30$（缓冲层，使用混合公式）：

$$U^{+}=\frac{1}{\kappa }\ln (E{y}^{+})-\frac{1}{\kappa }\ln (1+0.3y^{+})$$

Spalding统一壁面函数：

Spalding提出了一个适用于全$y^{+}$范围的统一公式：

$$y^{+}=U^{+}+\frac{1}{E}\left[e^{\kappa U^{+}}-1-\kappa U^{+}-\frac{(\kappa U^{+}{)}^2}{2}-\frac{(\kappa U^{+}{)}^3}{6}\right]$$

这个公式在$y^{+}<5$时退化为$U^{+}=y^{+}$，在$y^{+}>30$时退化为对数律。

壁面函数的物理基础：

壁面函数基于以下关键假设：

1. 常应力层假设：在近壁区，总剪切应力近似为常数
2. 局部平衡假设：湍流生成率等于耗散率
3. 混合长度理论：涡粘度与距离成正比

壁面函数的验证：

壁面函数的有效性已通过大量实验数据验证：

• 零压力梯度边界层：误差<5%
• 管道流动：误差<3%
• 通道流动：误差<4%

但对于以下情况需要特殊处理：

• 强压力梯度流动
• 分离和再附流动
• 三维复杂流动
• 高马赫数流动

温度壁面函数：

对于恒壁温条件：

$$T^{+}=\frac{(T_w-T_P)\rho c_p{u}_{\tau }}{q_w}={\sigma }_t\left(U^{+}+P\right)$$

其中$P$为Jayatilleke函数：

$$P=9.24\left[{\left(\frac{{\sigma }_l}{{\sigma }_t}\right)}^{0.75}-1\right]\left[1+0.28\exp \left(-0.007\frac{{\sigma }_l}{{\sigma }_t}\right)\right]$$

${\sigma }_l$为层流普朗特数，${\sigma }_t$为湍流普朗特数。

壁面函数的理论基础：

壁面函数基于以下理论假设：

1. 局部平衡假设：在近壁区域，湍流生成与耗散近似平衡
2. 相似性假设：无量纲速度剖面具有普适性
3. 常应力层假设：在$30<y^{+}<500$区域，剪切应力近似为常数

壁面函数的适用范围：

• 适用条件：

  • 零压力梯度或弱压力梯度
  • 高雷诺数湍流
  • 无分离或轻微分离
  • 壁面粗糙度较小

• 不适用情况：

  • 强逆压梯度流动
  • 大分离流动
  • 低雷诺数流动
  • 强三维效应
  • 转捩流动

壁面函数的改进方向：

1. 非平衡壁面函数：考虑压力梯度和非平衡效应
2. 增强壁面函数：扩展适用范围，提高精度
3. 自适应壁面函数：根据流动条件自动调整

2.3 壁面剪切应力计算

从壁面函数可以推导壁面剪切应力：

$${\tau }_w=\frac{\rho \kappa U_P{u}_{\tau }}{\ln (E{y}_P^{+})}$$

其中$U_P$为第一个内点的速度，$y_P^{+}$为该点的无量纲距离。

迭代求解过程：

1. 假设初始$u_{\tau }$
2. 计算$y_P^{+}=y_P{u}_{\tau }/\nu$
3. 从壁面函数计算$U^{+}$
4. 更新$u_{\tau }=U_P/U^{+}$
5. 重复直到收敛

2.4 壁面热流计算

对于恒壁温边界：

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

对于恒热流边界：

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

热流计算的迭代过程：

1. 计算壁面剪切应力${\tau }_w$和摩擦速度$u_{\tau }$
2. 计算无量纲温度$T^{+}$
3. 根据边界条件类型计算$q_w$或$T_w$
4. 更新壁面边界条件
5. 重复直到收敛

努塞尔数计算：

$$Nu=\frac{hL}{k}=\frac{q_{w}L}{k(T_w-T_{ref})}$$

其中$h$为对流换热系数，$L$为特征长度，$k$为导热系数。

斯坦顿数计算：

$$St=\frac{h}{\rho c_p{U}_{\infty }}=\frac{Nu}{Re\cdot Pr}$$

2.5 可压缩流壁面函数

对于高速可压缩流动，需要考虑密度变化和动能效应：

Van Driest变换：

$$U_{VD}={\int }_0^{U^{+}}\sqrt{\frac{\rho }{{\rho }_w}}d{U}^{+}$$

可压缩流的速度剖面：

$$U_{VD}^{+}=\frac{1}{\kappa }\ln (y^{+})+C$$

其中$U_{VD}^{+}=U_{VD}/u_{\tau }$为变换后的无量纲速度。

温度-速度关系（Crocco-Busemann关系）：

$$T=T_w+(T_r-T_w)\frac{U}{U_{\infty }}-r\frac{U^2}{2c_p}$$

其中$T_r$为恢复温度，$r$为恢复因子（层流$r=P{r}^{1/2}$，湍流$r=P{r}^{1/3}$）。

绝热壁面温度：

$$T_{aw}=T_{\infty }\left(1+r\frac{\gamma -1}{2}{M}_{\infty }^2\right)$$

可压缩修正的壁面函数：

对于高马赫数流动，需要对标准壁面函数进行修正：

$$U^{+}=\frac{1}{\kappa }\ln (E{y}^{+})+\frac{2\Pi }{\kappa }{\sin }^2\left(\frac{\pi y}{2\delta }\right)$$

其中$\Pi$为尾迹参数，$\delta$为边界层厚度。

对于高速可压缩流动，需要考虑密度变化和动能效应：

Van Driest变换：

$$U_{VD}={\int }_0^{U^{+}}\sqrt{\frac{\rho }{{\rho }_w}}d{U}^{+}$$

温度-速度关系：

$$T=T_w+\frac{T_r-T_w}{U_{\infty }}U-\frac{r{U}^2}{2c_p}$$

其中$T_r$为恢复温度，$r$为恢复因子。

2.6 非平衡壁面函数

标准壁面函数假设局部平衡（生成=耗散），在非平衡流动中（强压力梯度、分离流）需要修正：

Shih等人的非平衡壁面函数：

$$U^{+}=\frac{1}{\kappa }\ln (E{y}^{+})+\frac{1}{\kappa }\left[2\Pi {\sin }^2\left(\frac{\pi y}{2\delta }\right)+\frac{y}{\delta }(1-\cos (2\pi y/\delta ))\right]$$

其中$\Pi$为压力梯度参数，$\delta$为边界层厚度。

2.7 粗糙壁面壁面函数

对于粗糙壁面，需要对标准壁面函数进行修正：

等效砂粒粗糙度$k_s$：

粗糙壁面的速度分布：

$$U^{+}=\frac{1}{\kappa }\ln \left(\frac{y}{k_s}\right)+8.5$$

或写成：

$$U^{+}=\frac{1}{\kappa }\ln (E{y}^{+})-\Delta U^{+}$$

其中$\Delta U^{+}$为粗糙度修正项：

$$\Delta U^{+}=\frac{1}{\kappa }\ln (1+0.3k_s^{+})$$

粗糙度区域分类：

• 水力光滑：$k_s^{+}<5$，粗糙度影响可忽略
• 过渡粗糙：$5<k_s^{+}<70$，部分粗糙元突出粘性底层
• 完全粗糙：$k_s^{+}>70$，粘性底层完全被破坏

粗糙度对传热的影响：

粗糙度不仅增加摩擦阻力，还增强传热：

$$\frac{N{u}_{rough}}{N{u}_{smooth}}={\left(\frac{f_{rough}}{f_{smooth}}\right)}^n$$

其中$n\approx 0.5$对于充分发展的湍流。

等效砂粒粗糙度的确定：

• 对于均匀砂粒粗糙面：直接使用砂粒直径
• 对于工业粗糙面：需要通过实验标定
• 对于涂层表面：考虑涂层厚度和分布

2.8 壁面函数的数值实现

实现步骤：

1. 计算$y^{+}$：

   y_plus = y * u_tau / nu
   

2. 确定壁面函数区域：

   • $y^{+}<5$：粘性底层
   • $5<y^{+}<30$：缓冲层
   • $30<y^{+}<500$：对数律层

3. 计算无量纲速度$U^{+}$：

   if y_plus < 5:
       U_plus = y_plus
   elif y_plus < 30:
       U_plus = 5 * log(y_plus) + 5.5  # 近似
   else:
       U_plus = (1/kappa) * log(E * y_plus)
   

4. 计算壁面剪切应力：

   tau_w = rho * u_tau**2
   u_tau = U_P / U_plus
   

5. 迭代求解：直到$u_{\tau }$收敛

收敛判据：

• 相对变化小于$10^{-6}$
• 通常3-5次迭代即可收敛

三、低雷诺数模型

3.1 低雷诺数模型的基本思想

低雷诺数模型通过以下方式直接模拟近壁面区域：

1. 阻尼函数：在壁面附近抑制湍流
2. 修正的输运方程：考虑粘性效应
3. 精细网格：第一个网格点放置在$y^{+}<1$

与壁面函数法的对比：

特性	壁面函数法	低雷诺数模型
网格要求	$y^{+}>30$	$y^{+}<1$
计算成本	低	高
精度	中等	高
适用范围	高雷诺数	全雷诺数
分离流	一般	好
近壁细节	无	有

低雷诺数模型的优势：

1. 直接模拟近壁区：可以捕捉粘性底层的详细流动结构
2. 适用于分离流：对逆压梯度和分离流动的预测更准确
3. 全雷诺数适用：从层流到湍流都可以使用
4. 物理细节丰富：可以提供近壁区的详细流动信息

低雷诺数模型的劣势：

1. 计算成本高：需要更多的网格点和计算时间
2. 网格要求严格：第一个点必须在$y^{+}<1$
3. 数值稳定性差：壁面附近梯度大，收敛困难
4. 内存需求大：精细网格需要更多的存储空间

3.2 低雷诺数k-ε模型

Lam-Bremhorst模型：

在标准k-ε模型中引入阻尼函数：

$$f_{\mu }={\left[1-\exp (-0.0165R_y)\right]}^2\left(1+\frac{20.5}{R_t}\right)$$

$$f_1=1+{\left(\frac{0.05}{f_{\mu }}\right)}^3$$

$$f_2=1-\exp (-R_t^2)$$

其中：

$$R_y=\frac{y\sqrt{k}}{\nu },R_t=\frac{k^2}{\nu \epsilon }$$

修正的涡粘度：

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

Lam-Bremhorst模型的特点：

1. 阻尼函数$f_{\mu }$：在壁面附近抑制涡粘度
2. 修正函数$f_1$和$f_2$：调整生成项和耗散项
3. 适用范围：可以模拟$y^{+}<1$的近壁区域

其他低雷诺数k-ε模型：

• Launder-Sharma模型：引入修正的耗散率$\tilde{\epsilon }$
• Chien模型：使用不同的阻尼函数形式
• Yang-Shih模型：针对旋转流动优化

3.3 低雷诺数k-ω模型

标准k-ω模型本身就是低雷诺数模型，可以直接积分到壁面。

SST k-ω模型的壁面处理：

• 在壁面附近自动切换到k-ω模式
• 不需要额外的阻尼函数
• 壁面边界条件：$k_w=0$，${\omega }_w=\frac{6\nu }{{\beta }_1{y}_1^2}$

k-ω模型的壁面边界条件推导：

在壁面附近，k的渐近行为为：
$$k\sim y^2$$

因此壁面边界条件为：
$$k_w=0$$

对于ω，Wilcox建议的壁面边界条件为：
$${\omega }_w=\frac{6\nu }{{\beta }_1{y}_1^2}$$

其中$y_1$为第一个内点到壁面的距离，${\beta }_1=0.075$。

数值实现注意事项：

1. 避免在壁面直接计算$\omega$，使用理论值
2. 第一个内点的$y^{+}$应小于1
3. 网格在壁面法向应足够精细
4. 使用高阶离散格式提高精度
5. 采用隐式时间推进保证稳定性

k-ω模型与k-ε模型的近壁处理对比：

特性	k-ε模型	k-ω模型
阻尼函数	需要	不需要
壁面边界条件	复杂	简单
数值稳定性	较差	较好
自由流敏感性	低	高（标准模型）
分离流预测	一般	好

3.4 雷诺应力模型（RSM）的近壁面处理

对于RSM模型，需要在壁面附近对雷诺应力进行修正：

压力-应变项的壁面反射效应：

$${\varphi }_{ij}^w=C_1^w\frac{\epsilon }{k}\left(\overline{u_k{u}_m}{n}_k{n}_m{\delta }_{ij}-\frac{3}{2}\overline{u_i{u}_k}{n}_k{n}_j-\frac{3}{2}\overline{u_j{u}_k}{n}_k{n}_i\right)\frac{k^{3/2}}{C_{\mu }\epsilon y}$$

其中$n_i$为壁面法向单位向量。

3.5 近壁面网格要求

壁面函数法：

• 第一个内点：$30<y^{+}<300$
• 壁面法向增长率：1.1-1.3
• 边界层内至少10-20个网格点

低雷诺数模型：

• 第一个内点：$y^{+}<1$
• 壁面法向增长率：1.1-1.2
• 边界层内至少15-30个网格点
• 粘性底层内至少5-10个网格点

网格质量检查清单：

1. y+检查：

   • 壁面函数法：$30<y^{+}<100$
   • 低雷诺数模型：$y^{+}<1$
   • 避免$5<y^{+}<30$区域

2. 增长率检查：

   • 壁面法向增长率：1.1-1.3
   • 避免增长率突变
   • 确保网格平滑过渡

3. 网格正交性：

   • 壁面法向与壁面正交
   • 避免高度倾斜的网格
   • 检查网格扭曲度

4. 边界层覆盖：

   • 边界层内足够网格点
   • 粘性底层内足够网格点
   • 外层网格适当粗化

四、增强壁面处理（EWT）

4.1 EWT的基本概念

增强壁面处理（Enhanced Wall Treatment）结合了壁面函数和低雷诺数模型的优点：

1. 混合方法：根据$y^{+}$自动选择处理方式
2. 双层模型：在近壁区使用一方程模型
3. 连续过渡：在缓冲层平滑过渡

4.2 双层模型

在近壁区域（$y^{+}<200$）使用一方程湍流模型：

$$\frac{\partial k}{\partial t}+U_j\frac{\partial k}{\partial x_j}=\frac{\partial }{\partial x_j}\left[(\nu +{\nu }_t)\frac{\partial k}{\partial x_j}\right]+P_k-\frac{k^{3/2}}{l_{\epsilon }}$$

其中混合长度：

$$l_{\epsilon }=C_{l}y\left[1-\exp \left(-\frac{y^{+}}{A_{\epsilon }}\right)\right]$$

涡粘度：

$${\nu }_t=C_{\mu }{l}_{\mu }\sqrt{k}$$

4.3 混合函数

EWT使用混合函数在双层模型和核心区域湍流模型之间过渡：

$$\varphi ={\varphi }_{inner}{F}_{blend}+{\varphi }_{outer}(1-F_{blend})$$

其中混合函数：

$$F_{blend}=\tanh {\left(\frac{y^{+}}{200}\right)}^4$$

EWT的优势：

1. 适用范围广：可以处理$y^{+}$从1到300的各种情况
2. 自动适应：不需要用户手动选择壁面处理方式
3. 精度较高：在近壁区使用更精确的模型
4. 鲁棒性好：对网格质量要求相对较低

EWT的局限性：

1. 计算成本：比标准壁面函数高
2. 复杂性：实现和调试更复杂
3. 验证需求：需要更多的验证工作
4. 适用范围：对某些复杂流动可能不适用

EWT的适用场景：

• 网格质量不均匀
• 流动条件复杂多变
• 无法确定合适的$y^{+}$范围
• 需要自动适应不同区域

EWT的实现要点：

1. 混合函数设计：确保平滑过渡
2. 内层模型选择：一方程或双方程模型
3. 外层模型耦合：与核心区域湍流模型匹配
4. 边界条件处理：统一处理各种边界条件

商业软件中的EWT：

• ANSYS Fluent：Enhanced Wall Treatment
• STAR-CCM+：All y+ Wall Treatment
• OpenFOAM：continuous wall function
• CFX：Automatic Wall Function

1. 计算成本：比标准壁面函数高
2. 复杂性：实现和调试更复杂
3. 验证需求：需要更多的验证工作

五、近壁面传热模型

5.1 热边界层结构

与速度边界层类似，热边界层也可以分为三个区域：

热粘性底层：

• 范围：$y^{+}<5$
• 特征：热传导主导
• 温度分布：$T^{+}=Pr\cdot y^{+}$

热缓冲层：

• 范围：$5<y^{+}<30$
• 特征：热传导和湍流热扩散同等重要

热对数律层：

• 范围：$30<y^{+}<500$
• 特征：湍流热扩散主导
• 温度分布：$T^{+}=\frac{P{r}_t}{\kappa }\ln (y^{+})+C(Pr)$

5.2 温度壁面函数

恒壁温边界条件：

对于湍流普朗特数$P{r}_t=0.85$：

$$T^{+}=\frac{(T_w-T_P)\rho c_p{u}_{\tau }}{q_w}=P{r}_t\left(U^{+}+P\right)$$

其中$P$为Jayatilleke函数：

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

不同普朗特数的修正：

• 对于气体（$Pr\approx 0.7$）：标准壁面函数适用
• 对于水（$Pr\approx 7$）：需要修正
• 对于油（$Pr\approx 100$）：需要特殊处理
• 对于液态金属（$Pr\ll 1$）：分子热传导主导

温度壁面函数的适用范围：

温度壁面函数适用于：

• 恒壁温或恒热流边界条件
• 中等普朗特数（$0.5<Pr<10$）
• 无强烈变物性效应
• 弱压力梯度流动

不适用于：

• 大温差流动（强变物性）
• 高普朗特数液体（$Pr>50$）
• 低普朗特数液态金属（$Pr<0.1$）
• 强压力梯度流动

温度壁面函数的改进：

1. Kader修正：考虑压力梯度效应
2. Kays修正：适用于宽普朗特数范围
3. Spalding统一公式：适用于全$y^{+}$范围

• 对于水（$Pr\approx 7$）：需要修正
• 对于油（$Pr\approx 100$）：需要特殊处理
• 对于液态金属（$Pr\ll 1$）：分子热传导主导

5.3 热流计算

壁面热流密度：

$$q_w=-k{\frac{\partial T}{\partial y}|}_{y=0}$$

使用壁面函数计算：

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

努塞尔数计算：

$$Nu=\frac{q_{w}L}{k(T_w-T_{\infty })}=\frac{St\cdot Re\cdot Pr}{1}$$

其中$St$为斯坦顿数，$Re$为雷诺数。

5.4 变物性效应

在高温或低温条件下，物性参数随温度变化，需要考虑变物性效应：

参考温度法：

使用参考温度计算物性：
$$T_{ref}=T_w-0.5(T_w-T_{\infty })$$

或使用膜温度：
$$T_f=\frac{T_w+T_{\infty }}{2}$$

物性比修正：

对于气体，粘度随温度变化：
$$\frac{\mu }{{\mu }_{\infty }}={\left(\frac{T}{T_{\infty }}\right)}^n$$

其中$n\approx 0.7$对于空气。

变物性对传热的影响：

1. 大温差气体流动：

   • 密度变化显著
   • 浮力效应可能重要
   • 需要使用可压缩流壁面函数

2. 液体加热：

   • 粘度随温度降低
   • 边界层变薄
   • 传热增强

3. 液体冷却：

   • 粘度随温度增加
   • 边界层增厚
   • 传热减弱

变物性修正公式：

对于管内流动：
$$Nu=N{u}_{const}{\left(\frac{{\mu }_b}{{\mu }_w}\right)}^{0.14}$$

其中${\mu }_b$为体平均温度下的粘度，${\mu }_w$为壁面温度下的粘度。

六、近壁面处理的数值方法

6.1 网格生成策略

边界层网格要求：

1. 壁面法向网格：

   • 壁面函数法：第一个点$y^{+}=30-100$
   • 低雷诺数模型：第一个点$y^{+}<1$
   • 增长率：1.1-1.3

2. 流向网格：

   • 边界层内：$\Delta x^{+}\approx 50-100$
   • 分离区附近：需要加密

3. 展向网格：

   • 三维流动：$\Delta z^{+}\approx 20-50$

网格质量检查：

• 检查$y^{+}$分布
• 检查网格正交性
• 检查网格长宽比
• 检查网格扭曲度

6.2 离散格式选择

对流项离散：

• 壁面函数法：二阶迎风格式或混合格式
• 低雷诺数模型：需要更高精度，推荐三阶格式

扩散项离散：

• 始终使用中心差分格式
• 确保二阶精度

时间离散：

• 稳态问题：隐式格式
• 非稳态问题：二阶隐式或Crank-Nicolson格式
• 时间步长选择：满足CFL条件

离散格式的稳定性分析：

1. CFL条件：
   $$\Delta t<\frac{\Delta x}{|U|+c}$$
   其中$c$为声速。

2. 扩散稳定性：
   $$\Delta t<\frac{\Delta y^2}{2\nu }$$

3. 湍流方程稳定性：

   • k方程：需要较小的时间步长
   • ε或ω方程：可能需要更小的时间步长

高阶格式的应用：

• MUSCL格式：用于捕捉激波和间断

• WENO格式：高分辨率，适用于复杂流动

• 谱方法：高精度，适用于简单几何

• 稳态问题：隐式格式

• 非稳态问题：二阶隐式或Crank-Nicolson格式

6.3 求解策略

耦合求解 vs 分离求解：

• 耦合求解：适用于高速流动和强耦合问题
• 分离求解：适用于低速流动，计算成本低

欠松弛因子：

• 速度：0.7-0.9
• 压力：0.1-0.3
• 湍流变量：0.5-0.8
• 温度：0.9-1.0

收敛判据：

• 连续性方程残差：$10^{-3}$或更小
• 动量方程残差：$10^{-4}$或更小
• 能量方程残差：$10^{-6}$或更小
• 壁面剪切应力和热流变化小于0.1%

加速收敛的方法：

1. 多重网格：

   • 在粗网格上快速消除低频误差
   • 在细网格上修正高频误差
   • 显著加速收敛

2. 预处理：

   • 块ILU预处理
   • 不完全Cholesky分解
   • 代数多重网格预处理

3. 初始猜测：

   • 使用粗网格解作为初始猜测
   • 使用势流解作为初始猜测
   • 使用相似流动解作为初始猜测

收敛监控：

• 监控残差下降曲线
• 监控关键物理量（如阻力、升力）
• 监控质量流量平衡
• 监控能量平衡

七、工程应用案例

7.1 平板边界层

问题描述：

• 来流速度：$U_{\infty }=50$ m/s
• 平板长度：$L=1$ m
• 流体：空气，$R{e}_L=3.3\times 10^6$

网格设置：

• 壁面函数法：$y^{+}\approx 50$，100个网格点
• 低雷诺数模型：$y^{+}<1$，200个网格点

结果对比：

• 摩擦系数与理论值误差：壁面函数法5%，低雷诺数模型2%
• 计算时间：壁面函数法1小时，低雷诺数模型5小时

7.2 翼型绕流

问题描述：

• 翼型：NACA 0012
• 攻角：$\alpha =5°$
• 雷诺数：$Re=2\times 10^6$

近壁面处理选择：

• 前缘：低雷诺数模型（曲率大）
• 中部：壁面函数法
• 后缘：低雷诺数模型（可能分离）

关键结果：

• 升力系数：与实验数据误差<2%
• 阻力系数：与实验数据误差<5%
• 分离点位置：误差<5%弦长

7.3 管内流动换热

问题描述：

• 圆管直径：$D=0.05$ m
• 雷诺数：$Re=10,000$
• 恒壁温边界条件

近壁面处理：

• 使用低雷诺数模型
• $y^{+}<1$，边界层内20个网格点

结果验证：

• 努塞尔数：与Gnielinski公式误差<8%
• 摩擦系数：与Blasius公式误差<5%
• 温度分布：与解析解吻合良好

案例总结：

案例	推荐方法	关键考虑因素	典型误差
平板边界层	壁面函数法	简单流动，高雷诺数	3-5%
翼型绕流	混合方法	曲率、分离	2-5%
管内流动	低雷诺数模型	传热预测	5-8%
后台阶流动	低雷诺数模型	分离再附	10-15%

工程实践建议：

1. 初步设计阶段：使用壁面函数法，快速获得结果
2. 详细设计阶段：使用低雷诺数模型，提高精度
3. 关键区域：对重要区域使用低雷诺数模型
4. 验证阶段：与实验数据对比，评估模型适用性

八、常见问题与解决方法

8.1 $y^{+}$计算错误

问题表现：

• 壁面剪切应力计算不准确
• 收敛困难
• 结果与实验偏差大

解决方法：

1. 检查网格质量
2. 确保第一个内点位置正确
3. 检查物性参数
4. 使用正确的参考速度

8.2 数值不稳定

问题表现：

• 残差震荡
• 解发散
• 物理量出现非物理值

解决方法：

1. 降低欠松弛因子
2. 改善网格质量
3. 使用更稳健的离散格式
4. 检查边界条件设置

8.3 分离流预测不准

问题表现：

• 分离点位置偏差大
• 再附长度不准确
• 压力分布错误

解决方法：

1. 使用低雷诺数模型代替壁面函数
2. 加密分离区网格
3. 使用SST k-ω等适合分离流的模型
4. 考虑使用LES或DES

8.4 传热预测偏差

问题表现：

• 努塞尔数计算误差大
• 壁面温度分布不正确
• 热流密度预测偏差

解决方法：

1. 检查温度壁面函数适用性
2. 考虑变物性效应
3. 使用低雷诺数模型提高精度
4. 检查普朗特数设置

8.5 收敛困难

问题表现：

• 残差停滞不降
• 物理量震荡
• 计算时间过长

解决方法：

1. 改善初始猜测
2. 调整欠松弛因子
3. 使用多重网格加速
4. 检查网格质量
5. 简化物理模型

8.6 网格相关问题

问题表现：

• y+分布不均匀
• 网格扭曲严重
• 边界层覆盖不足

解决方法：

1. 重新生成网格
2. 调整增长率
3. 增加边界层网格点
4. 使用结构化网格
5. 检查网格正交性

问题表现：

• 分离点位置偏差大
• 再附长度不准确
• 压力分布错误

解决方法：

1. 使用低雷诺数模型代替壁面函数
2. 加密分离区网格
3. 使用SST k-ω等适合分离流的模型
4. 考虑使用LES或DES

九、Python仿真实现

9.1 代码概述

本节提供了基于Python的近壁面处理仿真代码，包括：

1. 壁面函数实现：标准壁面函数及其应用
2. 边界层网格生成：自动生成满足y+要求的网格
3. 近壁面流动可视化：边界层结构和速度剖面

代码特点：

• 使用NumPy进行数值计算
• 使用Matplotlib进行可视化
• 模块化设计，便于扩展
• 详细的注释和文档字符串

9.2 壁面函数实现

import numpy as np

def wall_function(y_plus, kappa=0.41, E=9.8):
    """
    标准壁面函数
    
    参数:
        y_plus: 无量纲壁面距离
        kappa: 卡门常数
        E: 粗糙度参数
    
    返回:
        U_plus: 无量纲速度
    """
    if y_plus < 5:
        # 粘性底层
        U_plus = y_plus
    elif y_plus < 30:
        # 缓冲层 - 使用混合公式
        U_viscous = y_plus
        U_log = (1/kappa) * np.log(E * y_plus)
        # 混合
        blend = (y_plus - 5) / 25
        U_plus = (1 - blend) * U_viscous + blend * U_log
    else:
        # 对数律层
        U_plus = (1/kappa) * np.log(E * y_plus)
    
    return U_plus

def compute_tau_wall(U_P, y_P, rho, nu, max_iter=100, tol=1e-6):
    """
    计算壁面剪切应力
    
    参数:
        U_P: 第一个内点速度
        y_P: 第一个内点到壁面距离
        rho: 密度
        nu: 运动粘度
        max_iter: 最大迭代次数
        tol: 收敛容差
    
    返回:
        tau_w: 壁面剪切应力
        u_tau: 摩擦速度
        y_plus: 无量纲距离
    """
    # 初始猜测
    u_tau = np.sqrt(0.01 * U_P**2)  # 初始猜测
    
    for i in range(max_iter):
        y_plus = y_P * u_tau / nu
        U_plus = wall_function(y_plus)
        
        # 更新u_tau
        u_tau_new = U_P / U_plus
        
        # 检查收敛
        if abs(u_tau_new - u_tau) / u_tau < tol:
            break
        
        u_tau = u_tau_new
    
    tau_w = rho * u_tau**2
    
    return tau_w, u_tau, y_plus

# 示例使用
if __name__ == "__main__":
    # 流动参数
    U_P = 20.0  # m/s
    y_P = 1e-4  # m
    rho = 1.225  # kg/m3
    nu = 1.5e-5  # m2/s
    
    tau_w, u_tau, y_plus = compute_tau_wall(U_P, y_P, rho, nu)
    
    print(f"壁面剪切应力: {tau_w:.4f} Pa")
    print(f"摩擦速度: {u_tau:.4f} m/s")
    print(f"y+: {y_plus:.2f}")
    print(f"无量纲速度U+: {wall_function(y_plus):.2f}")

9.2 边界层网格生成

import numpy as np
import matplotlib.pyplot as plt

def generate_boundary_layer_mesh(y_max, n_points, y_plus_first, growth_rate=1.2):
    """
    生成边界层网格
    
    参数:
        y_max: 边界层外缘位置
        n_points: 网格点数
        y_plus_first: 第一个点的y+值
        growth_rate: 网格增长率
    
    返回:
        y: 网格点坐标
        dy: 网格间距
    """
    # 使用指数分布生成网格
    y = np.zeros(n_points)
    
    # 第一个点
    y[0] = y_plus_first
    
    # 生成网格
    for i in range(1, n_points):
        dy = y[i-1] * (growth_rate - 1)
        y[i] = y[i-1] + dy
        
        if y[i] > y_max:
            y[i] = y_max
            y = y[:i+1]
            break
    
    # 归一化到y_max
    y = y / y[-1] * y_max
    
    dy = np.diff(y)
    
    return y, dy

def check_y_plus(y, u_tau, nu):
    """
    检查y+分布
    
    参数:
        y: 网格点坐标
        u_tau: 摩擦速度
        nu: 运动粘度
    
    返回:
        y_plus: 无量纲距离
    """
    y_plus = y * u_tau / nu
    return y_plus

# 示例
if __name__ == "__main__":
    # 参数
    y_max = 0.01  # m
    n_points = 50
    u_tau = 1.0  # m/s
    nu = 1.5e-5  # m2/s
    
    # 目标y+ = 1
    y_plus_target = 1.0
    y_first = y_plus_target * nu / u_tau
    
    # 生成网格
    y, dy = generate_boundary_layer_mesh(y_max, n_points, y_first, growth_rate=1.15)
    
    # 检查y+
    y_plus = check_y_plus(y, u_tau, nu)
    
    print(f"网格点数: {len(y)}")
    print(f"第一个点y+: {y_plus[0]:.2f}")
    print(f"最大y+: {y_plus[-1]:.2f}")
    print(f"平均增长率: {np.mean(dy[1:]/dy[:-1]):.3f}")
    
    # 可视化
    plt.figure(figsize=(10, 5))
    
    plt.subplot(1, 2, 1)
    plt.semilogy(y_plus, 'bo-', markersize=4)
    plt.xlabel('Grid Index')
    plt.ylabel('y+')
    plt.title('y+ Distribution')
    plt.grid(True)
    
    plt.subplot(1, 2, 2)
    plt.plot(dy, 'ro-', markersize=4)
    plt.xlabel('Grid Index')
    plt.ylabel('dy')
    plt.title('Grid Spacing')
    plt.grid(True)
    
    plt.tight_layout()
    plt.savefig('boundary_layer_mesh.png', dpi=150)
    plt.show()

网格生成代码说明：

该代码实现了以下功能：

1. 指数分布网格：使用指数增长生成边界层网格
2. y+检查：计算并检查网格点的无量纲距离
3. 可视化输出：绘制y+分布和网格间距

使用建议：

• 根据目标y+调整第一个网格点高度
• 控制增长率在1.1-1.3之间
• 确保边界层内足够网格点

9.4 近壁面流动可视化

import numpy as np
import matplotlib.pyplot as plt

def velocity_profile_y_plus(y_plus_max=500, n_points=100):
    """
    计算边界层速度剖面
    """
    y_plus = np.linspace(0.1, y_plus_max, n_points)
    U_plus = np.zeros_like(y_plus)
    
    kappa = 0.41
    E = 9.8
    
    for i, yp in enumerate(y_plus):
        if yp < 5:
            U_plus[i] = yp
        elif yp < 30:
            # 缓冲层
            U_viscous = yp
            U_log = (1/kappa) * np.log(E * yp)
            blend = (yp - 5) / 25
            U_plus[i] = (1 - blend) * U_viscous + blend * U_log
        else:
            U_plus[i] = (1/kappa) * np.log(E * yp)
    
    return y_plus, U_plus

def plot_boundary_layer_regions():
    """
    绘制边界层各区域
    """
    y_plus, U_plus = velocity_profile_y_plus()
    
    plt.figure(figsize=(10, 6))
    
    # 绘制速度剖面
    plt.semilogx(y_plus, U_plus, 'b-', linewidth=2, label='Velocity Profile')
    
    # 标记各区域
    plt.axvline(x=5, color='r', linestyle='--', alpha=0.5, label='Viscous Sublayer End')
    plt.axvline(x=30, color='g', linestyle='--', alpha=0.5, label='Buffer Layer End')
    
    # 填充区域
    plt.fill_between([0.1, 5], 0, 25, alpha=0.2, color='red', label='Viscous Sublayer')
    plt.fill_between([5, 30], 0, 25, alpha=0.2, color='yellow', label='Buffer Layer')
    plt.fill_between([30, 500], 0, 25, alpha=0.2, color='green', label='Logarithmic Layer')
    
    plt.xlabel('y+', fontsize=12)
    plt.ylabel('U+', fontsize=12)
    plt.title('Turbulent Boundary Layer Structure', fontsize=14)
    plt.legend(loc='upper left')
    plt.grid(True, alpha=0.3)
    plt.xlim([0.1, 500])
    plt.ylim([0, 25])
    
    plt.tight_layout()
    plt.savefig('boundary_layer_regions.png', dpi=150)
    plt.show()

if __name__ == "__main__":
    plot_boundary_layer_regions()

可视化代码说明：

该代码绘制了湍流边界层的速度剖面，并标记了三个区域：

1. 粘性底层（红色）：$y^{+}<5$
2. 缓冲层（黄色）：$5<y^{+}<30$
3. 对数律层（绿色）：$30<y^{+}<500$

应用价值：

• 理解边界层结构
• 验证壁面函数
• 教学演示

代码扩展建议：

1. 添加温度剖面：绘制温度边界层分布
2. 添加粗糙度效应：考虑粗糙壁面的修正
3. 添加压力梯度效应：考虑顺压和逆压梯度
4. 添加可压缩性效应：考虑高速流动的修正

九、Python仿真实现

以下是基于Python的近壁面处理方法仿真代码，实现了y+计算、壁面函数应用和网格敏感性分析。

"""
主题014：近壁面处理方法仿真
包含案例：
1. 平板边界层y+分布分析
2. 壁面函数vs低雷诺数模型
3. 圆管流动网格敏感性
4. 后台阶流动的近壁面处理
5. y+计算器与网格优化
"""

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("主题014：近壁面处理方法")
print("=" * 60)

# ============================================
# 案例一：平板边界层y+分布分析
# ============================================
print("\n案例一：平板边界层y+分布分析")
print("-" * 50)

def calculate_yplus(y, U_inf, nu, L=1.0):
    """
    计算y+值
    
    参数:
        y: 壁面距离
        U_inf: 自由流速度
        nu: 运动粘度
        L: 特征长度
    """
    # 估算摩擦系数（基于平板湍流边界层）
    Re_L = U_inf * L / nu
    Cf = 0.0592 * Re_L**(-0.2)  # 湍流平板摩擦系数
    
    # 壁面剪切应力
    tau_w = 0.5 * Cf * U_inf**2
    u_tau = np.sqrt(tau_w)
    
    # y+
    y_plus = y * u_tau / nu
    
    return y_plus, u_tau, tau_w

def create_grid_yplus(target_yplus, U_inf, nu, N=50, growth_rate=1.2):
    """
    创建满足目标y+的网格
    
    参数:
        target_yplus: 目标y+值
        U_inf: 自由流速度
        nu: 运动粘度
        N: 网格点数
        growth_rate: 增长率
    """
    # 计算第一层网格高度
    Re_L = U_inf * 1.0 / nu
    Cf = 0.0592 * Re_L**(-0.2)
    u_tau = np.sqrt(0.5 * Cf * U_inf**2)
    y1 = target_yplus * nu / u_tau
    
    # 创建几何增长网格
    y = np.zeros(N)
    y[0] = y1
    for i in range(1, N):
        y[i] = y[i-1] * growth_rate
    
    # 计算实际的y+分布
    y_plus = y * u_tau / nu
    
    return y, y_plus

# 流动参数
U_inf = 1.0
nu = 1e-5
Re_L = U_inf * 1.0 / nu

print(f"  雷诺数 Re_L = {Re_L:.2e}")

# 三种网格配置
configs = [
    {'name': 'Coarse (Wall Function)', 'target_yplus': 50, 'N': 30, 'growth': 1.3},
    {'name': 'Medium (Transition)', 'target_yplus': 15, 'N': 40, 'growth': 1.2},
    {'name': 'Fine (Low-Re)', 'target_yplus': 0.5, 'N': 50, 'growth': 1.1}
]

results_yplus = {}
for config in configs:
    y, y_plus = create_grid_yplus(
        config['target_yplus'], U_inf, nu, 
        N=config['N'], growth_rate=config['growth']
    )
    results_yplus[config['name']] = {'y': y, 'y_plus': y_plus}
    print(f"  {config['name']}: y1+ = {y_plus[0]:.2f}, y_max+ = {y_plus[-1]:.2f}")

# 绘制y+分布
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# y+分布对比
ax1 = axes[0, 0]
colors = ['#E74C3C', '#F39C12', '#27AE60']
for i, (name, data) in enumerate(results_yplus.items()):
    ax1.semilogy(range(len(data['y_plus'])), data['y_plus'], 
                'o-', color=colors[i], linewidth=2, markersize=4, label=name)
ax1.axhline(y=5, color='k', linestyle='--', alpha=0.5, label='Viscous Sublayer')
ax1.axhline(y=30, color='k', linestyle=':', alpha=0.5, label='Buffer Layer')
ax1.axhline(y=300, color='k', linestyle='-.', alpha=0.5, label='Log Layer End')
ax1.set_xlabel('Grid Index', fontsize=11)
ax1.set_ylabel(r'$y^+$', fontsize=11)
ax1.set_title('y+ Distribution for Different Grids', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 边界层子层区域
ax2 = axes[0, 1]
y_plus_range = np.logspace(-1, 3, 100)
ax2.fill_betweenx(y_plus_range, 0, 1, where=(y_plus_range < 5), 
                  alpha=0.3, color='blue', label='Viscous Sublayer')
ax2.fill_betweenx(y_plus_range, 0, 1, where=((y_plus_range >= 5) & (y_plus_range < 30)), 
                  alpha=0.3, color='yellow', label='Buffer Layer')
ax2.fill_betweenx(y_plus_range, 0, 1, where=((y_plus_range >= 30) & (y_plus_range < 500)), 
                  alpha=0.3, color='green', label='Logarithmic Layer')
ax2.set_xlim([0, 1])
ax2.set_ylim([0.1, 1000])
ax2.set_yscale('log')
ax2.set_ylabel(r'$y^+$', fontsize=11)
ax2.set_title('Boundary Layer Sublayers', fontsize=12)
ax2.legend(fontsize=9, loc='upper right')

# 网格高度分布
ax3 = axes[1, 0]
for i, (name, data) in enumerate(results_yplus.items()):
    ax3.plot(range(len(data['y'])), data['y']*1000, 
            'o-', color=colors[i], linewidth=2, markersize=4, label=name)
ax3.set_xlabel('Grid Index', fontsize=11)
ax3.set_ylabel('y (mm)', fontsize=11)
ax3.set_title('Grid Height Distribution', fontsize=12)
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 第一层网格高度vs目标y+
ax4 = axes[1, 1]
target_yplus_range = np.array([0.1, 0.5, 1, 5, 10, 30, 50, 100, 300])
y1_heights = []
for ty in target_yplus_range:
    y, _ = create_grid_yplus(ty, U_inf, nu, N=10, growth_rate=1.2)
    y1_heights.append(y[0]*1000)

ax4.loglog(target_yplus_range, y1_heights, 'bo-', linewidth=2, markersize=8)
ax4.axvline(x=5, color='b', linestyle='--', alpha=0.5, label='Viscous Sublayer')
ax4.axvline(x=30, color='r', linestyle='--', alpha=0.5, label='Buffer Layer Start')
ax4.set_xlabel('Target y+', fontsize=11)
ax4.set_ylabel('First Layer Height (mm)', fontsize=11)
ax4.set_title('First Layer Height vs Target y+', fontsize=12)
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('case1_yplus_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ y+分布分析图已保存")

# ============================================
# 案例二：壁面函数vs低雷诺数模型
# ============================================
print("\n案例二：壁面函数vs低雷诺数模型")
print("-" * 50)

def wall_function_velocity(y_plus):
    """
    标准壁面函数速度分布
    """
    kappa = 0.41
    E = 9.8
    U_plus = np.zeros_like(y_plus)
    
    # 粘性底层
    mask_viscous = y_plus < 5
    U_plus[mask_viscous] = y_plus[mask_viscous]
    
    # 缓冲层和对数律层
    mask_log = y_plus >= 5
    U_plus[mask_log] = 1/kappa * np.log(E * y_plus[mask_log])
    
    return U_plus

def low_re_velocity(y_plus, Re_tau=1000):
    """
    低雷诺数模型速度分布（使用Spalding公式近似）
    """
    kappa = 0.41
    B = 5.0
    
    # Spalding公式
    U_plus = np.zeros_like(y_plus)
    for i, yp in enumerate(y_plus):
        if yp < 5:
            U_plus[i] = yp
        else:
            # 迭代求解
            def spalding_eq(U):
                return yp - U - 1/kappa * (np.exp(kappa*U) - 1 - kappa*U - (kappa*U)**2/2)
            
            from scipy.optimize import fsolve
            try:
                U_plus[i] = fsolve(spalding_eq, np.log(yp))[0]
            except:
                U_plus[i] = 1/kappa * np.log(9.8 * yp)
    
    return U_plus

def boundary_layer_with_wall_function(Re_theta=1000, N=50, use_wall_function=True):
    """
    求解边界层，使用壁面函数或低雷诺数模型
    """
    # 边界层厚度
    theta = 1.0
    delta = 10 * theta
    
    # 网格
    if use_wall_function:
        # 壁面函数：第一个点在对数律层
        y_wall = 30 * theta / Re_theta
        y = np.linspace(y_wall, delta, N)
    else:
        # 低雷诺数：精细网格到壁面
        eta = np.linspace(0, 1, N)
        y = delta * eta**2  # 二次分布，壁面附近加密
    
    # 物性
    nu = theta / Re_theta
    
    # 迭代求解
    U = np.ones(N) * 0.5
    if not use_wall_function:
        U[0] = 0
    U[-1] = 1.0
    
    for iteration in range(1000):
        U_old = U.copy()
        
        # 速度梯度
        dUdy = np.gradient(U, y)
        
        # 涡粘度（混合长度模型）
        l_mix = 0.41 * y * (1 - np.exp(-y * np.sqrt(0.5) / (26 * nu)))
        nu_t = l_mix**2 * np.abs(dUdy)
        
        # 求解速度方程
        A = np.zeros((N, N))
        b = np.zeros(N)
        
        for i in range(1, N-1):
            nu_eff = nu + nu_t[i]
            dy = y[i+1] - y[i-1]
            
            A[i, i-1] = nu_eff / ((y[i] - y[i-1]) * dy)
            A[i, i] = -2 * nu_eff / ((y[i+1] - y[i-1]) * 0.5 * dy)
            A[i, i+1] = nu_eff / ((y[i+1] - y[i]) * dy)
        
        # 边界条件
        if use_wall_function:
            # 壁面函数边界
            u_tau = 0.05
            y_plus_1 = y[0] * u_tau / nu
            U_plus_1 = 1/0.41 * np.log(9.8 * y_plus_1)
            A[0, 0] = 1.0
            b[0] = U_plus_1 * u_tau
        else:
            A[0, 0] = 1.0
            b[0] = 0.0
        
        A[-1, -1] = 1.0
        b[-1] = 1.0
        
        U = np.linalg.solve(A, b)
        
        if np.linalg.norm(U - U_old) < 1e-6:
            break
    
    return y, U

# 运行两种方法
y_wf, U_wf = boundary_layer_with_wall_function(Re_theta=1000, N=50, use_wall_function=True)
y_lr, U_lr = boundary_layer_with_wall_function(Re_theta=1000, N=80, use_wall_function=False)

# 转换为壁面单位
nu = 1/1000
u_tau_wf = np.sqrt(nu * (U_wf[1] - U_wf[0]) / (y_wf[1] - y_wf[0]))
u_tau_lr = np.sqrt(nu * (U_lr[1] - U_lr[0]) / (y_lr[1] - y_lr[0]))

y_plus_wf = y_wf * u_tau_wf / nu
y_plus_lr = y_lr * u_tau_lr / nu
U_plus_wf = U_wf / u_tau_wf
U_plus_lr = U_lr / u_tau_lr

# 理论对数律
y_plus_theory = np.linspace(30, 500, 50)
U_plus_theory = 1/0.41 * np.log(9.8 * y_plus_theory)

print(f"  壁面函数法: u_tau = {u_tau_wf:.4f}")
print(f"  低雷诺数模型: u_tau = {u_tau_lr:.4f}")

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

# 速度剖面对比
ax1 = axes[0, 0]
ax1.semilogx(y_plus_wf, U_plus_wf, 'ro-', linewidth=2, markersize=4, label='Wall Function')
ax1.semilogx(y_plus_lr, U_plus_lr, 'bs-', linewidth=2, markersize=4, label='Low-Re Model')
ax1.semilogx(y_plus_theory, U_plus_theory, 'k--', linewidth=2, label='Log Law')
ax1.set_xlabel(r'$y^+$', fontsize=11)
ax1.set_ylabel(r'$U^+$', fontsize=11)
ax1.set_title('Velocity Profile Comparison', fontsize=12)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([1, 1000])

# 壁面法向速度分布
ax2 = axes[0, 1]
ax2.plot(y_wf, U_wf, 'r-', linewidth=2, label='Wall Function')
ax2.plot(y_lr, U_lr, 'b-', linewidth=2, label='Low-Re Model')
ax2.set_xlabel('y', fontsize=11)
ax2.set_ylabel('U', fontsize=11)
ax2.set_title('Velocity Distribution (Physical Coordinates)', fontsize=12)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 误差分析
ax3 = axes[1, 0]
# 插值到相同y+进行比较
from scipy.interpolate import interp1d
f_wf = interp1d(y_plus_wf, U_plus_wf, bounds_error=False, fill_value='extrapolate')
f_lr = interp1d(y_plus_lr, U_plus_lr, bounds_error=False, fill_value='extrapolate')

y_plus_common = np.linspace(30, 300, 50)
U_wf_interp = f_wf(y_plus_common)
U_lr_interp = f_lr(y_plus_common)
U_theory_interp = 1/0.41 * np.log(9.8 * y_plus_common)

error_wf = np.abs(U_wf_interp - U_theory_interp) / U_theory_interp * 100
error_lr = np.abs(U_lr_interp - U_theory_interp) / U_theory_interp * 100

ax3.semilogy(y_plus_common, error_wf, 'r-o', linewidth=2, markersize=4, label='Wall Function Error')
ax3.semilogy(y_plus_common, error_lr, 'b-s', linewidth=2, markersize=4, label='Low-Re Error')
ax3.set_xlabel(r'$y^+$', fontsize=11)
ax3.set_ylabel('Relative Error (%)', fontsize=11)
ax3.set_title('Comparison with Log Law', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 方法对比总结
ax4 = axes[1, 1]
ax4.axis('off')

comparison_text = """
Wall Function vs Low-Re Model:

Wall Function:
  ✓ Coarser grid (y+ = 30-300)
  ✓ Lower computational cost
  ✓ Good for high-Re flows
  ✗ Limited accuracy near wall
  ✗ Poor for separated flows
  ✗ Sensitive to y+ value

Low-Re Model:
  ✓ High accuracy near wall
  ✓ Good for separated flows
  ✓ Direct integration to wall
  ✗ Fine grid required (y+ < 1)
  ✗ Higher computational cost
  ✗ More grid points needed

Recommendation:
  • Use wall function for:
    - High Reynolds numbers
    - Attached flows
    - Preliminary studies
  
  • Use low-Re model for:
    - Flow separation
    - Accurate wall quantities
    - Research applications
"""

ax4.text(0.05, 0.95, comparison_text, fontsize=9, verticalalignment='top',
        family='monospace', bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.5))

plt.tight_layout()
plt.savefig('case2_wall_function_vs_lowre.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 壁面函数vs低雷诺数模型对比图已保存")

# ============================================
# 案例三：圆管流动网格敏感性
# ============================================
print("\n案例三：圆管流动网格敏感性")
print("-" * 50)

def pipe_flow_grid_study(Re=10000, Pr=0.7, N_list=[20, 40, 80, 160]):
    """
    网格敏感性研究
    """
    R = 1.0
    nu = 1.0 / Re
    
    results = {}
    
    for N in N_list:
        # 创建网格（壁面附近加密）
        eta = np.linspace(0, 1, N)
        r = R * (1 - np.tanh(3*(1-eta))/np.tanh(3))
        dr = np.diff(r)
        dr = np.concatenate([[dr[0]], dr])
        
        # 初始化
        U = np.zeros(N)
        T = np.ones(N)
        
        # 迭代求解
        for iteration in range(2000):
            U_old = U.copy()
            
            # 速度方程
            A = np.zeros((N, N))
            b = np.zeros(N)
            
            for i in range(1, N-1):
                r_p = 0.5 * (r[i] + r[i+1])
                r_m = 0.5 * (r[i] + r[i-1])
                
                A[i, i-1] = r_m / (r[i] * dr[i]**2)
                A[i, i] = -(r_p + r_m) / (r[i] * dr[i]**2)
                A[i, i+1] = r_p / (r[i] * dr[i]**2)
                b[i] = -0.1  # 压力梯度
            
            # 边界条件
            A[0, 0] = 1.0
            A[0, 1] = -1.0
            b[0] = 0.0
            A[-1, -1] = 1.0
            b[-1] = 0.0
            
            U = np.linalg.solve(A, b)
            
            if np.linalg.norm(U - U_old) < 1e-8:
                break
        
        # 温度方程
        dTdx = -0.01
        A_T = np.zeros((N, N))
        b_T = np.zeros(N)
        
        for i in range(1, N-1):
            r_p = 0.5 * (r[i] + r[i+1])
            r_m = 0.5 * (r[i] + r[i-1])
            
            A_T[i, i-1] = r_m / (r[i] * dr[i]**2)
            A_T[i, i] = -(r_p + r_m) / (r[i] * dr[i]**2)
            A_T[i, i+1] = r_p / (r[i] * dr[i]**2)
            b_T[i] = U[i] * dTdx
        
        A_T[0, 0] = 1.0
        A_T[0, 1] = -1.0
        b_T[0] = 0.0
        A_T[-1, -1] = 1.0
        b_T[-1] = 0.0
        
        T = np.linalg.solve(A_T, b_T)
        
        # 计算努塞尔数
        dTdr_wall = (T[-1] - T[-2]) / (r[-1] - r[-2])
        qw = -dTdr_wall
        Tb = 2 * np.trapezoid(U * T * r, r) / np.trapezoid(U * r, r)
        Tw = T[-1]
        Nu = qw * 2 * R / abs(Tw - Tb)
        
        # 计算摩擦因子
        tau_w = nu * (U[-1] - U[-2]) / (r[-1] - r[-2])
        f = 8 * tau_w / (np.mean(U)**2)
        
        # 计算第一层y+
        u_tau = np.sqrt(abs(tau_w))
        y1 = r[-1] - r[-2]
        y1_plus = y1 * u_tau / nu
        
        results[N] = {
            'r': r, 'U': U, 'T': T, 'Nu': Nu, 'f': f,
            'y1_plus': y1_plus
        }
        
        print(f"  N={N:3d}: Nu={Nu:.3f}, f={f:.5f}, y1+={y1_plus:.3f}")
    
    return results

results_grid = pipe_flow_grid_study(Re=10000, Pr=0.7)

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

# 速度分布对比
ax1 = axes[0, 0]
colors = plt.cm.viridis(np.linspace(0, 1, len(results_grid)))
for i, (N, data) in enumerate(results_grid.items()):
    ax1.plot(data['r'], data['U'], color=colors[i], linewidth=2, label=f'N={N}')
ax1.set_xlabel('r/R', fontsize=11)
ax1.set_ylabel('U', fontsize=11)
ax1.set_title('Velocity Profile (Grid Convergence)', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 温度分布对比
ax2 = axes[0, 1]
for i, (N, data) in enumerate(results_grid.items()):
    ax2.plot(data['r'], data['T'], color=colors[i], linewidth=2, label=f'N={N}')
ax2.set_xlabel('r/R', fontsize=11)
ax2.set_ylabel('T', fontsize=11)
ax2.set_title('Temperature Profile (Grid Convergence)', fontsize=12)
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 努塞尔数收敛
ax3 = axes[1, 0]
N_values = list(results_grid.keys())
Nu_values = [results_grid[N]['Nu'] for N in N_values]
ax3.plot(N_values, Nu_values, 'bo-', linewidth=2, markersize=8)
ax3.axhline(y=Nu_values[-1], color='r', linestyle='--', label='Reference (Finest Grid)')
ax3.set_xlabel('Number of Grid Points', fontsize=11)
ax3.set_ylabel('Nu', fontsize=11)
ax3.set_title('Nusselt Number Convergence', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 摩擦因子收敛
ax4 = axes[1, 1]
f_values = [results_grid[N]['f'] for N in N_values]
f_theory = 0.3164 / 10000**0.25  # Blasius公式
ax4.plot(N_values, f_values, 'go-', linewidth=2, markersize=8, label='Computed')
ax4.axhline(y=f_theory, color='r', linestyle='--', label=f'Blasius: {f_theory:.5f}')
ax4.set_xlabel('Number of Grid Points', fontsize=11)
ax4.set_ylabel('Friction Factor', fontsize=11)
ax4.set_title('Friction Factor Convergence', fontsize=12)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('case3_grid_sensitivity.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 网格敏感性分析图已保存")

# ============================================
# 案例四：后台阶流动的近壁面处理
# ============================================
print("\n案例四：后台阶流动的近壁面处理")
print("-" * 50)

def backward_step_wall_treatment(Re_h=5000, method='standard'):
    """
    后台阶流动，不同近壁面处理
    
    方法:
        'standard': 标准壁面函数
        'non_equilibrium': 非平衡壁面函数
        'low_re': 低雷诺数模型
    """
    # 简化模型
    x = np.linspace(-5, 20, 100)
    
    # 再附长度（经验值）
    if method == 'standard':
        Xr_h = 6.0  # 标准壁面函数通常低估
    elif method == 'non_equilibrium':
        Xr_h = 6.8  # 非平衡壁面函数改进
    else:  # low_re
        Xr_h = 7.2  # 低雷诺数模型最准确
    
    # 简化的壁面剪切应力分布
    tau_w = np.zeros_like(x)
    for i, xi in enumerate(x):
        if xi < 0:
            tau_w[i] = 1.0
        elif xi < Xr_h:
            # 分离区
            tau_w[i] = -0.5 * np.exp(-xi/2)
        else:
            # 再附后恢复
            tau_w[i] = 1.0 * (1 - np.exp(-(xi - Xr_h)/3))
    
    # 速度剖面（简化）
    U_sep = np.zeros_like(x)
    for i, xi in enumerate(x):
        if 0 < xi < Xr_h:
            U_sep[i] = 0.3 * np.sin(np.pi * xi / Xr_h)
    
    return x, tau_w, U_sep, Xr_h

methods = ['standard', 'non_equilibrium', 'low_re']
results_step = {}

for method in methods:
    x, tau_w, U_sep, Xr = backward_step_wall_treatment(Re_h=5000, method=method)
    results_step[method] = {'x': x, 'tau_w': tau_w, 'U_sep': U_sep, 'Xr': Xr}
    print(f"  {method:20s}: Xr/h = {Xr:.1f}")

# 实验参考值
Xr_exp = 7.0
print(f"  {'Experimental':20s}: Xr/h = {Xr_exp:.1f}")

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

# 壁面剪切应力分布
ax1 = axes[0, 0]
colors = ['#E74C3C', '#F39C12', '#27AE60']
for i, (method, data) in enumerate(results_step.items()):
    label = method.replace('_', ' ').title()
    ax1.plot(data['x'], data['tau_w'], color=colors[i], linewidth=2, label=label)
ax1.axhline(y=0, color='k', linestyle='--', linewidth=1)
ax1.axvline(x=0, color='gray', linestyle=':', alpha=0.5)
ax1.set_xlabel('x/h', fontsize=11)
ax1.set_ylabel(r'$\tau_w$', fontsize=11)
ax1.set_title('Wall Shear Stress Distribution', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([-5, 15])

# 分离区速度
ax2 = axes[0, 1]
for i, (method, data) in enumerate(results_step.items()):
    label = method.replace('_', ' ').title()
    ax2.plot(data['x'], data['U_sep'], color=colors[i], linewidth=2, label=label)
ax2.set_xlabel('x/h', fontsize=11)
ax2.set_ylabel('Reverse Flow Velocity', fontsize=11)
ax2.set_title('Separation Bubble Velocity', fontsize=12)
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)
ax2.set_xlim([0, 10])

# 再附长度对比
ax3 = axes[1, 0]
method_names = ['Standard\nWF', 'Non-Equilibrium\nWF', 'Low-Re\nModel', 'Experimental']
Xr_values = [results_step['standard']['Xr'], 
             results_step['non_equilibrium']['Xr'],
             results_step['low_re']['Xr'], Xr_exp]
colors_bar = ['#E74C3C', '#F39C12', '#27AE60', '#3498DB']
bars = ax3.bar(method_names, Xr_values, color=colors_bar, alpha=0.8, edgecolor='black')
ax3.axhline(y=Xr_exp, color='k', linestyle='--', linewidth=2, label='Experimental')
ax3.set_ylabel('Reattachment Length (Xr/h)', fontsize=11)
ax3.set_title('Reattachment Length Comparison', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3, axis='y')
for bar, xr in zip(bars, Xr_values):
    height = bar.get_height()
    ax3.text(bar.get_x() + bar.get_width()/2., height,
            f'{xr:.1f}', ha='center', va='bottom', fontsize=10)

# 误差分析
ax4 = axes[1, 1]
errors = [abs(xr - Xr_exp)/Xr_exp * 100 for xr in Xr_values[:-1]]
method_names_short = ['Standard WF', 'Non-Equilibrium WF', 'Low-Re Model']
bars = ax4.bar(method_names_short, errors, color=colors_bar[:-1], alpha=0.8, edgecolor='black')
ax4.set_ylabel('Error (%)', fontsize=11)
ax4.set_title('Reattachment Length Error', fontsize=12)
ax4.grid(True, alpha=0.3, axis='y')
for bar, err in zip(bars, errors):
    height = bar.get_height()
    ax4.text(bar.get_x() + bar.get_width()/2., height,
            f'{err:.1f}%', ha='center', va='bottom', fontsize=10)

plt.tight_layout()
plt.savefig('case4_backward_step_wall_treatment.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 后台阶流动近壁面处理图已保存")

# ============================================
# 案例五：y+计算器与网格优化
# ============================================
print("\n案例五：y+计算器与网格优化")
print("-" * 50)

def yplus_calculator(U_inf, L, rho, mu, target_yplus, method='turbulent'):
    """
    y+计算器：估算第一层网格高度
    
    参数:
        U_inf: 自由流速度 (m/s)
        L: 特征长度 (m)
        rho: 密度 (kg/m³)
        mu: 动力粘度 (Pa·s)
        target_yplus: 目标y+值
        method: 'laminar' 或 'turbulent'
    
    返回:
        推荐的网格参数
    """
    nu = mu / rho
    Re = rho * U_inf * L / mu
    
    if method == 'turbulent':
        # 湍流摩擦系数
        Cf = 0.0592 * Re**(-0.2)
        tau_w = 0.5 * Cf * rho * U_inf**2
    else:
        # 层流
        Cf = 1.328 / np.sqrt(Re)
        tau_w = 0.5 * Cf * rho * U_inf**2
    
    u_tau = np.sqrt(tau_w / rho)
    
    # 第一层网格高度
    y1 = target_yplus * nu / u_tau
    
    # 边界层厚度估计
    if method == 'turbulent':
        delta = 0.37 * L / Re**0.2
    else:
        delta = 5 * L / np.sqrt(Re)
    
    # 推荐网格点数
    if target_yplus < 5:  # 低雷诺数模型
        N_recommended = int(20 * (delta / y1)**0.5)
    else:  # 壁面函数
        N_recommended = int(10 * np.log(delta / y1) / np.log(1.2))
    
    return {
        'Re': Re,
        'Cf': Cf,
        'u_tau': u_tau,
        'y1': y1,
        'delta': delta,
        'N_recommended': N_recommended
    }

def optimize_grid_growth(y1, delta, N, max_growth=1.3):
    """
    优化网格增长率
    
    参数:
        y1: 第一层高度
        delta: 边界层厚度
        N: 网格点数
        max_growth: 最大增长率
    
    返回:
        最优增长率和网格分布
    """
    from scipy.optimize import minimize_scalar
    
    def objective(growth):
        # 计算网格总高度
        y_total = y1 * (growth**N - 1) / (growth - 1)
        # 目标：使总高度接近边界层厚度
        return abs(y_total - delta)
    
    result = minimize_scalar(objective, bounds=(1.01, max_growth), method='bounded')
    optimal_growth = result.x
    
    # 生成网格
    y = np.zeros(N)
    y[0] = y1
    for i in range(1, N):
        y[i] = y[i-1] * optimal_growth
    
    return optimal_growth, y

# 示例计算
print("  y+计算器示例:")
print("  " + "-" * 50)

test_cases = [
    {'name': 'Airfoil (High Re)', 'U_inf': 50, 'L': 1.0, 'rho': 1.225, 'mu': 1.8e-5, 'y+': 30},
    {'name': 'Pipe Flow', 'U_inf': 2, 'L': 0.1, 'rho': 1000, 'mu': 1e-3, 'y+': 1},
    {'name': 'Low Re Channel', 'U_inf': 0.1, 'L': 0.05, 'rho': 1000, 'mu': 1e-3, 'y+': 0.5},
]

results_calculator = {}
for case in test_cases:
    result = yplus_calculator(
        case['U_inf'], case['L'], case['rho'], case['mu'], 
        case['y+'], method='turbulent'
    )
    
    # 优化网格
    growth_opt, y_grid = optimize_grid_growth(
        result['y1'], result['delta'], result['N_recommended']
    )
    
    results_calculator[case['name']] = {
        'params': result,
        'growth': growth_opt,
        'y_grid': y_grid
    }
    
    print(f"  {case['name']}:")
    print(f"    Re = {result['Re']:.2e}")
    print(f"    y1 = {result['y1']*1e6:.2f} μm")
    print(f"    δ = {result['delta']*1000:.2f} mm")
    print(f"    N_recommended = {result['N_recommended']}")
    print(f"    Optimal growth = {growth_opt:.3f}")
    print()

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

# y1 vs Re for different target y+
ax1 = axes[0, 0]
Re_range = np.logspace(4, 7, 50)
y1_targets = [0.5, 1, 5, 30, 100]
colors_y1 = plt.cm.viridis(np.linspace(0, 1, len(y1_targets)))

for i, yp in enumerate(y1_targets):
    y1_values = []
    for Re in Re_range:
        U_inf = 1.0
        L = Re * 1.8e-5 / 1.225 / U_inf
        result = yplus_calculator(U_inf, L, 1.225, 1.8e-5, yp)
        y1_values.append(result['y1']*1e6)
    ax1.loglog(Re_range, y1_values, color=colors_y1[i], linewidth=2, 
              label=f'y+ = {yp}')

ax1.set_xlabel('Reynolds Number', fontsize=11)
ax1.set_ylabel('First Layer Height (μm)', fontsize=11)
ax1.set_title('First Layer Height vs Re (Air)', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3, which='both')

# 网格分布对比
ax2 = axes[0, 1]
for i, (name, data) in enumerate(results_calculator.items()):
    y_norm = data['y_grid'] / data['y_grid'][-1]
    ax2.plot(range(len(y_norm)), y_norm, 'o-', linewidth=2, markersize=4,
            label=name)
ax2.set_xlabel('Grid Index', fontsize=11)
ax2.set_ylabel('Normalized Height', fontsize=11)
ax2.set_title('Optimized Grid Distribution', fontsize=12)
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 增长率vs网格点数
ax3 = axes[1, 0]
N_range = np.arange(10, 101, 10)
growth_values = []
for N in N_range:
    growth, _ = optimize_grid_growth(1e-5, 0.01, N)
    growth_values.append(growth)
ax3.plot(N_range, growth_values, 'bo-', linewidth=2, markersize=8)
ax3.axhline(y=1.2, color='r', linestyle='--', label='Typical Growth Rate')
ax3.set_xlabel('Number of Grid Points', fontsize=11)
ax3.set_ylabel('Optimal Growth Rate', fontsize=11)
ax3.set_title('Optimal Growth Rate vs Grid Points', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# y+选择指南
ax4 = axes[1, 1]
ax4.axis('off')

guide_text = """
y+ Selection Guide:

For Wall Functions:
  Target: 30 < y+ < 100
  Avoid: 5 < y+ < 30 (buffer layer)
  First cell: In logarithmic layer
  Growth rate: 1.2-1.3
  BL points: 10-20

For Low-Re Models:
  Target: y+ < 1
  First cell: In viscous sublayer
  Growth rate: 1.1-1.2
  BL points: 15-30
  Viscous sublayer: 5-10 points

For Enhanced Wall Treatment:
  Target: y+ ≈ 1 or y+ > 30
  Avoid: 5 < y+ < 10
  Can handle mixed y+
  Automatic blending

Grid Quality Checklist:
  □ y+ within target range
  □ Smooth growth rate
  口 Sufficient BL points
  口 Orthogonal to wall
  口 Aspect ratio < 100
"""

ax4.text(0.05, 0.95, guide_text, fontsize=9, verticalalignment='top',
        family='monospace', bbox=dict(boxstyle='round', facecolor='lightcyan', alpha=0.5))

plt.tight_layout()
plt.savefig('case5_yplus_calculator.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ y+计算器与网格优化图已保存")

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

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

# 1. 边界层结构
ax1 = axes[0, 0]
y_plus_full = np.logspace(-1, 3, 100)
U_plus_viscous = y_plus_full
U_plus_log = 1/0.41 * np.log(9.8 * y_plus_full)

ax1.fill_betweenx(y_plus_full, 0, 1, where=(y_plus_full < 5), 
                  alpha=0.3, color='blue', label='Viscous')
ax1.fill_betweenx(y_plus_full, 0, 1, where=((y_plus_full >= 5) & (y_plus_full < 30)), 
                  alpha=0.3, color='yellow', label='Buffer')
ax1.fill_betweenx(y_plus_full, 0, 1, where=((y_plus_full >= 30) & (y_plus_full < 500)), 
                  alpha=0.3, color='green', label='Logarithmic')

mask_viscous = y_plus_full < 5
mask_log = y_plus_full >= 30
ax1.plot(U_plus_viscous[mask_viscous], y_plus_full[mask_viscous], 'b-', linewidth=2)
ax1.plot(U_plus_log[mask_log], y_plus_full[mask_log], 'g-', linewidth=2)
ax1.set_xlabel(r'$U^+$', fontsize=10)
ax1.set_ylabel(r'$y^+$', fontsize=10)
ax1.set_title('Boundary Layer Structure', fontsize=11)
ax1.set_yscale('log')
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.3)

# 2. 壁面函数公式
ax2 = axes[0, 1]
y_plus_wf = np.linspace(1, 500, 100)
U_plus_standard = wall_function_velocity(y_plus_wf)
U_plus_lam = y_plus_wf

ax2.plot(y_plus_wf, U_plus_standard, 'b-', linewidth=2, label='Wall Function')
ax2.plot(y_plus_wf, U_plus_lam, 'r--', linewidth=2, label='Linear')
ax2.axvline(x=5, color='gray', linestyle=':', alpha=0.5)
ax2.axvline(x=30, color='gray', linestyle=':', alpha=0.5)
ax2.set_xlabel(r'$y^+$', fontsize=10)
ax2.set_ylabel(r'$U^+$', fontsize=10)
ax2.set_title('Wall Function Velocity Profile', fontsize=11)
ax2.set_xscale('log')
ax2.legend(fontsize=8)
ax2.grid(True, alpha=0.3)

# 3. 方法对比
ax3 = axes[0, 2]
methods_comp = ['Standard\nWF', 'Non-Equil.\nWF', 'Low-Re\nModel', 'EWT']
accuracy = [3, 4, 5, 4.5]
cost = [2, 3, 5, 4]
x_pos = np.arange(len(methods_comp))
width = 0.35
bars1 = ax3.bar(x_pos - width/2, accuracy, width, label='Accuracy', color='green', alpha=0.7)
bars2 = ax3.bar(x_pos + width/2, cost, width, label='Cost', color='red', alpha=0.7)
ax3.set_ylabel('Score (1-5)', fontsize=10)
ax3.set_title('Method Comparison', fontsize=11)
ax3.set_xticks(x_pos)
ax3.set_xticklabels(methods_comp, fontsize=8)
ax3.legend(fontsize=8)
ax3.grid(True, alpha=0.3, axis='y')

# 4. y+分布示例
ax4 = axes[1, 0]
for i, (name, data) in enumerate(results_yplus.items()):
    ax4.semilogy(data['y_plus'], 'o-', markersize=3, linewidth=1.5, label=name)
ax4.axhline(y=5, color='k', linestyle='--', alpha=0.3)
ax4.axhline(y=30, color='k', linestyle='--', alpha=0.3)
ax4.set_xlabel('Grid Index', fontsize=10)
ax4.set_ylabel(r'$y^+$', fontsize=10)
ax4.set_title('Grid y+ Distribution Examples', fontsize=11)
ax4.legend(fontsize=8)
ax4.grid(True, alpha=0.3)

# 5. 网格收敛
ax5 = axes[1, 1]
N_conv = list(results_grid.keys())
Nu_conv = [results_grid[N]['Nu'] for N in N_conv]
ax5.plot(N_conv, Nu_conv, 'bo-', linewidth=2, markersize=8)
ax5.set_xlabel('Grid Points', fontsize=10)
ax5.set_ylabel('Nu', fontsize=10)
ax5.set_title('Grid Convergence', fontsize=11)
ax5.grid(True, alpha=0.3)

# 6. 最佳实践总结
ax6 = axes[1, 2]
ax6.axis('off')

best_practice_text = """
Best Practices:

1. Grid Generation:
   • Calculate y1 based on target y+
   • Use smooth growth rate (1.1-1.3)
   • Ensure sufficient BL points
   • Check y+ after simulation

2. Method Selection:
   • High Re + attached → WF
   • Separation → Low-Re or EWT
   • Accurate wall data → Low-Re
   • Resource limited → WF

3. Validation:
   • Compare with correlations
   • Check grid independence
   • Verify y+ distribution
   • Monitor wall quantities

4. Common Pitfalls:
   • y+ in buffer layer
   • Sudden grid expansion
   • Insufficient BL resolution
   • Mixed wall treatments
"""

ax6.text(0.05, 0.95, best_practice_text, fontsize=8.5, verticalalignment='top',
        family='monospace', bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.5))

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_yplus_analysis.png - y+分布分析")
print("  2. case2_wall_function_vs_lowre.png - 壁面函数vs低雷诺数模型")
print("  3. case3_grid_sensitivity.png - 网格敏感性分析")
print("  4. case4_backward_step_wall_treatment.png - 后台阶流动近壁面处理")
print("  5. case5_yplus_calculator.png - y+计算器与网格优化")
print("  6. comprehensive_comparison.png - 综合对比")
print("\n输出目录：d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题014_近壁面处理方法")
print("=" * 60)