第五篇：边界层理论详解

摘要

边界层理论是理解对流换热现象的核心基础，由普朗特于1904年提出，彻底改变了流体力学和传热学的研究范式。本教程系统深入地阐述边界层理论的基本概念、数学模型和工程应用。首先介绍边界层的物理本质、形成机理和分类，详细解释边界层概念提出的历史背景和理论意义；然后详细推导层流边界层的控制方程，包括Blasius方程的推导与求解，介绍打靶法等数值求解技术；接着讨论湍流边界层的特性、转捩机理和湍流模型，分析层流与湍流的本质区别；重点分析热边界层与流动边界层的相互作用、相似解理论和类比关系，建立动量传递与热量传递之间的定量联系；详细讨论楔形流动、轴对称边界层、可压缩边界层等扩展问题；最后通过Python编程实现边界层流动的数值求解，包括速度分布、温度分布、摩擦系数和换热系数的计算，并通过可视化展示边界层的发展演化过程。本教程旨在为学习者提供边界层理论的完整知识体系，建立从理论到实践的桥梁，为工程传热问题的分析和解决奠定坚实基础。

关键词

边界层理论，Blasius方程，层流，湍流，热边界层，相似解，类比定律，摩擦系数，换热系数，转捩，Falkner-Skan方程，可压缩流动，数值求解

---

1. 引言

1.1 边界层概念的诞生

1904年，德国哥廷根大学的普朗特（Ludwig Prandtl）在海德尔堡国际数学大会上发表了著名的边界层理论。这一理论指出：在高雷诺数流动中，粘性效应主要集中在靠近固体壁面的薄层区域内，而在此薄层之外，流动可以近似为无粘流动。

这一发现具有划时代的意义：

• 理论层面：将复杂的N-S方程简化，使许多实际问题获得解析解
• 计算层面：大大降低了数值计算的复杂度
• 工程层面：为飞行器设计、换热器设计等提供了理论基础

本教程将从基本原理出发，逐步深入边界层理论的核心内容，包括控制方程推导、相似解理论、湍流模型、类比关系等，并通过Python代码展示边界层流动的数值求解方法。

1.4 边界层理论的重要性

边界层理论在现代工程科学中具有极其重要的地位：

航空航天工程：

• 飞机机翼的升力和阻力计算
• 飞行器表面摩擦阻力估算
• 高超声速飞行器的气动加热分析

能源动力系统：

• 汽轮机叶片的冷却设计
• 换热器的热性能优化
• 核反应堆的安全分析

环境科学：

• 大气边界层对气候的影响
• 海洋边界层的物质输运
• 污染物扩散预测

生物医学工程：

• 血液流动与血管壁的相互作用
• 人工器官的设计优化
• 药物输送系统的开发

掌握边界层理论，是理解和解决这些工程问题的关键。

1.2 边界层的物理本质

边界层是流体流经固体表面时，由于粘性作用而在壁面附近形成的速度梯度显著的薄层。其物理特征包括：

速度分布特征：壁面处速度为零（无滑移条件），向外逐渐增加到主流速度，速度梯度主要集中在薄层内。

厚度特征：边界层厚度δ通常远小于特征长度L，即δ/L << 1。

粘性效应：边界层内粘性力与惯性力相当，是产生摩擦阻力和换热的主要原因。

流动结构：边界层内存在层流、过渡和湍流三种流动状态。

1.3 边界层分类

根据流动状态，边界层可分为：

层流边界层：流体分层流动，各层之间无宏观混合，换热主要依靠分子扩散。

湍流边界层：流体存在强烈的脉动和混合，换热效率显著高于层流。

过渡区：从层流到湍流的转变区域，流动特性复杂。

根据几何形状，可分为：

• 平板边界层
• 曲面边界层（机翼、圆柱等）
• 轴对称边界层

---

2. 边界层基本方程

2.1 边界层近似

对于高雷诺数流动，普朗特提出了边界层近似：

量级分析：设边界层厚度δ，特征长度L，特征速度U∞，则：

$$\frac{\delta }{L}\sim \frac{1}{\sqrt{R{e}_L}}<<1$$

近似假设：

1. 沿壁面方向的尺度远大于垂直方向：$\partial /\partial x<<\partial /\partial y$
2. 垂直方向的速度分量远小于流向分量：$v<<u$
3. 垂直方向的动量方程简化为：$\partial p/\partial y\approx 0$

2.2 连续性方程

二维不可压流动的连续性方程：

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$$

2.3 动量方程（边界层方程）

在边界层近似下，x方向的动量方程简化为：

$$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-\frac{1}{\rho }\frac{dp}{dx}+\nu \frac{{\partial }^{2}u}{\partial y^2}$$

其中，压力梯度由外部势流决定：

$$-\frac{1}{\rho }\frac{dp}{dx}=U_e\frac{d{U}_e}{dx}$$

对于平板，$U_e=U_{\infty }$（常数），压力梯度为零。

2.4 能量方程

边界层内的能量方程：

$$u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\alpha \frac{{\partial }^{2}T}{\partial y^2}+\frac{\nu }{c_p}{\left(\frac{\partial u}{\partial y}\right)}^2$$

其中，α是热扩散系数，最后一项是粘性耗散项，在高速流动中需要考虑。

2.5 边界层方程的适用条件

边界层方程的推导基于以下假设：

高雷诺数假设：
$$Re=\frac{U_{\infty }L}{\nu }>>1$$
这是边界层近似成立的前提条件。当雷诺数较低时，粘性效应遍布整个流场，边界层概念不再适用。

薄层假设：
$$\frac{\delta }{L}<<1$$
边界层厚度远小于特征长度，确保边界层近似合理。

无分离假设：
边界层方程不适用于流动分离的情况。在逆压梯度作用下，边界层可能从壁面分离，此时需要求解完整的N-S方程。

二维假设：
上述方程针对二维流动推导。对于三维流动，需要增加z方向的动量方程和相应的速度分量w。

2.6 边界条件的物理意义

速度边界条件：

• 壁面无滑移条件：$u(x,0)=0$，体现了流体的粘性本质
• 壁面无穿透条件：$v(x,0)=0$，确保流体不穿过壁面
• 远场条件：$u(x,\infty )=U_e(x)$，边界层外缘恢复势流速度

温度边界条件：

• 等壁温条件：$T(x,0)=T_w$，壁面温度恒定
• 等热流条件：$-k\frac{\partial T}{\partial y}{|}_{y=0}=q_w$，壁面热流密度恒定
• 对流条件：$-k\frac{\partial T}{\partial y}{|}_{y=0}=h(T_w-T_f)$，第三类边界条件
• 远场条件：$T(x,\infty )=T_{\infty }$ 或 $\frac{\partial T}{\partial y}{|}_{y=\infty }=0$

---

3. 平板层流边界层

3.1 Blasius方程推导

对于零压力梯度的平板层流边界层，引入相似变量：

$$\eta =y\sqrt{\frac{U_{\infty }}{\nu x}}$$

定义流函数：

$$\psi =\sqrt{\nu x{U}_{\infty }}f(\eta )$$

速度分量表示为：

$$u=\frac{\partial \psi }{\partial y}=U_{\infty }{f}^{\prime }(\eta )$$

$$v=-\frac{\partial \psi }{\partial x}=\frac{1}{2}\sqrt{\frac{\nu U_{\infty }}{x}}(\eta f^{\prime }-f)$$

代入边界层方程，得到著名的Blasius方程：

$$f^{\prime \prime \prime }+\frac{1}{2}f{f}^{\prime \prime }=0$$

边界条件：

• $f(0)=0$（壁面无穿透）
• $f^{\prime }(0)=0$（壁面无滑移）
• $f^{\prime }(\infty )=1$（远场恢复主流速度）

3.2 Blasius方程的数值求解

Blasius方程是一个三阶非线性常微分方程，需要采用打靶法（Shooting Method）求解。

打靶法原理：

1. 猜测壁面处的二阶导数$f^{\prime \prime }(0)=\alpha$
2. 将方程作为初值问题从η=0积分到η=∞
3. 检查远场边界条件$f^{\prime }(\infty )=1$是否满足
4. 调整α值，直到满足边界条件

通过数值求解，得到：

$$f^{\prime \prime }(0)=0.33206$$

这个值是Blasius解中最重要的参数之一，直接决定了壁面剪应力的大小。

打靶法的数值实现细节：

1. 初值猜测：通常取$f^{\prime \prime }(0)\approx 0.3$作为初始猜测
2. 积分范围：理论上$\eta \to \infty$，实际计算取${\eta }_{max}=8\sim 10$即可
3. 收敛判据：$|f^{\prime }({\eta }_{max})-1|<10^{-6}$
4. 迭代方法：可以使用二分法、牛顿迭代法或fsolve等数值工具

Blasius解的重要数值结果：

参数	数值	物理意义
$f^{\prime \prime }(0)$	0.33206	壁面无量纲剪应力
${\eta }_{99}$	4.915	边界层厚度参数
${\delta }^{\ast }$	1.7208	位移厚度系数
$\theta$	0.6641	动量厚度系数
$H$	2.5914	形状因子

3.3 边界层厚度

速度边界层厚度（定义为u=0.99U∞处的y值）：

$$\delta =\frac{5.0x}{\sqrt{R{e}_x}}$$

位移厚度（表示由于边界层存在导致的质量流量亏损）：

$${\delta }^{\ast }={\int }_0^{\infty }\left(1-\frac{u}{U_{\infty }}\right)dy=\frac{1.72x}{\sqrt{R{e}_x}}$$

动量厚度（表示由于边界层存在导致的动量流量亏损）：

$$\theta ={\int }_0^{\infty }\frac{u}{U_{\infty }}\left(1-\frac{u}{U_{\infty }}\right)dy=\frac{0.664x}{\sqrt{R{e}_x}}$$

形状因子：

$$H=\frac{{\delta }^{\ast }}{\theta }=2.59$$

3.4 壁面摩擦

壁面剪应力：

$${\tau }_w=\mu {\left(\frac{\partial u}{\partial y}\right)}_{y=0}=0.332\mu U_{\infty }\sqrt{\frac{U_{\infty }}{\nu x}}$$

局部摩擦系数：

$$C_{f,x}=\frac{{\tau }_w}{\frac{1}{2}\rho U_{\infty }^2}=\frac{0.664}{\sqrt{R{e}_x}}$$

平均摩擦系数（从0到L）：

$${\bar{C}}_f=\frac{1.328}{\sqrt{R{e}_L}}$$

---

4. 平板热边界层

4.1 热边界层方程

当壁面温度Tw与主流温度T∞不同时，形成热边界层。引入无量纲温度：

$$\theta =\frac{T-T_w}{T_{\infty }-T_w}$$

热边界层方程为：

$$u\frac{\partial \theta }{\partial x}+v\frac{\partial \theta }{\partial y}=\alpha \frac{{\partial }^2\theta }{\partial y^2}$$

4.2 相似解

对于常壁温边界条件，引入相似变量：

$$\theta =\theta (\eta )$$

方程转化为：

$${\theta }^{\prime \prime }+\frac{1}{2}Pr\cdot f\cdot {\theta }^{\prime }=0$$

边界条件：

• $\theta (0)=0$（壁面）
• $\theta (\infty )=1$（远场）

4.3 Pohlhausen解

对于Prandtl数Pr ≈ 1的流体，Pohlhausen给出了近似解：

壁面温度梯度：

$${\theta }^{\prime }(0)=0.332P{r}^{1/3}$$

热边界层厚度：

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

局部Nusselt数：

$$N{u}_x=\frac{h_{x}x}{k}=0.332R{e}_x^{1/2}P{r}^{1/3}$$

平均Nusselt数：

$${\overline{Nu}}_L=0.664R{e}_L^{1/2}P{r}^{1/3}$$

4.4 热边界层厚度与速度边界层厚度的关系

热边界层厚度${\delta }_t$与速度边界层厚度$\delta$的比值取决于Prandtl数：

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

物理意义：

• 当$Pr=1$时，${\delta }_t=\delta$，热边界层与速度边界层厚度相同
• 当$Pr>1$时，${\delta }_t<\delta$，热边界层比速度边界层薄（如油类、水）
• 当$Pr<1$时，${\delta }_t>\delta$，热边界层比速度边界层厚（如液态金属）

不同流体的Prandtl数：

流体	Pr	适用条件
液态金属	0.001~0.1	高导热，低粘度
气体	0.6~0.8	常温常压
水	1~10	随温度变化
油类	50~1000	高粘度

4.5 变物性影响

当温差较大时，需要考虑物性随温度的变化：

参考温度法：

$$T_{ref}=T_{\infty }+0.5(T_w-T_{\infty })$$

或采用膜温度：

$$T_f=\frac{T_w+T_{\infty }}{2}$$

---

5. 湍流边界层

5.1 湍流特性

湍流边界层具有以下特征：

速度脉动：瞬时速度可分解为平均速度和脉动速度：

$$u=\bar{u}+u^{\prime },v=\bar{v}+v^{\prime }$$

雷诺应力：湍流脉动产生的附加应力：

$$-\rho \overline{u^{\prime }{v}^{\prime }}$$

湍流强度：

$$I=\frac{\sqrt{\overline{u^{\prime 2}}}}{\bar{u}}$$

5.2 湍流模型

混合长度模型（Prandtl）：

$$-\overline{u^{\prime }{v}^{\prime }}=l_m^2\left|\frac{\partial \bar{u}}{\partial y}\right|\frac{\partial \bar{u}}{\partial y}$$

其中，混合长度$l_m=\kappa y$，κ≈0.41是卡门常数。

对数律：在湍流边界层的对数区：

$$u^{+}=\frac{1}{\kappa }\ln y^{+}+B$$

其中，$u^{+}=\bar{u}/u_{\tau }$，$y^{+}=y{u}_{\tau }/\nu$，$u_{\tau }=\sqrt{{\tau }_w/\rho }$是摩擦速度，B≈5.5。

5.3 平板湍流边界层

1/7幂律速度分布：

$$\frac{\bar{u}}{U_{\infty }}={\left(\frac{y}{\delta }\right)}^{1/7}$$

边界层厚度：

$$\delta =\frac{0.37x}{R{e}_x^{1/5}}$$

局部摩擦系数（Prandtl-Schlichting公式）：

$$C_{f,x}=\frac{0.0592}{R{e}_x^{1/5}}$$

局部Nusselt数：

$$N{u}_x=0.0296R{e}_x^{4/5}P{r}^{1/3}$$

5.4 转捩

从层流到湍流的转变称为转捩。转捩雷诺数：

$$R{e}_{crit}=\frac{U_{\infty }{x}_{crit}}{\nu }\approx 5\times 10^5$$

影响转捩的因素：

• 来流湍流度
• 壁面粗糙度
• 压力梯度
• 壁面温度
• 可压缩性效应
• 外力场（如电磁力）

转捩控制：

工程上有时需要控制转捩位置：

• 延迟转捩：减小阻力，如超层流机翼
• 促进转捩：增强换热，如涡轮叶片冷却

控制方法包括：壁面吸吹、表面粗糙度控制、声波激励等。

---

6. 类比关系

6.1 Reynolds类比

对于Pr = 1的流体，动量传递和热量传递之间存在简单关系：

$$St=\frac{Nu}{Re\cdot Pr}=\frac{C_f}{2}$$

其中，St是Stanton数。

6.2 Chilton-Colburn类比

对于Pr ≠ 1的情况：

$$j_H=St\cdot P{r}^{2/3}=\frac{C_f}{2}$$

其中，$j_H$是传热的j因子。

6.3 类比的应用

类比关系的工程价值：

• 通过易测量的摩擦阻力推算换热系数
• 建立动量传递与热量传递的定量关系
• 简化复杂问题的分析

6.4 其他类比关系

Prandtl类比：
考虑层流底层和湍流核心区两层结构：
$$St=\frac{C_f/2}{1+5\sqrt{C_f/2}(Pr-1)}$$

von Kármán类比：
考虑层流底层、缓冲层和湍流核心区三层结构：
$$St=\frac{C_f/2}{1+5\sqrt{C_f/2}\left[(Pr-1)+\ln \left(\frac{5Pr+1}{6}\right)\right]}$$

适用范围比较：

类比关系	适用Pr范围	精度	复杂度
Reynolds	Pr ≈ 1	低	简单
Chilton-Colburn	0.5 < Pr < 60	中	简单
Prandtl	Pr > 0.5	高	中等
von Kármán	全范围	最高	较复杂

---

7. 其他边界层问题

7.1 楔形流动（Falkner-Skan方程）

对于具有压力梯度的流动，外部速度分布为：

$$U_e=C{x}^m$$

引入相似变量后，得到Falkner-Skan方程：

$$f^{\prime \prime \prime }+f{f}^{\prime \prime }+\beta (1-f^{\prime 2})=0$$

其中，$\beta =2m/(m+1)$是压力梯度参数：

• β > 0：加速流动（顺压梯度）
• β = 0：零压力梯度（Blasius解）
• β < 0：减速流动（逆压梯度）

7.2 轴对称边界层

对于绕旋成体的流动，需要考虑曲率效应。轴对称边界层方程：

$$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=U_e\frac{d{U}_e}{dx}+\nu \frac{{\partial }^{2}u}{\partial y^2}$$

连续性方程：

$$\frac{\partial (ru)}{\partial x}+\frac{\partial (rv)}{\partial y}=0$$

其中，r(x)是当地半径。

7.3 可压缩边界层

高速流动中，需要考虑：

• 密度变化
• 温度引起的物性变化
• 粘性耗散（绝热壁温）

绝热壁温：

$$T_{aw}=T_{\infty }+r\frac{U_{\infty }^2}{2c_p}$$

其中，r是恢复因子（层流r≈Pr^{{1/2}，湍流r≈Pr}{1/3}）。

7.4 边界层分离

分离机理：
在逆压梯度（$dp/dx>0$）作用下，边界层内流体动能不足以克服压力升高，导致流动反向，形成分离。

分离判据：
壁面剪应力为零：
$${\tau }_w=\mu {\left(\frac{\partial u}{\partial y}\right)}_{y=0}=0$$

或等价地：
$$f^{\prime \prime }(0)=0$$

Falkner-Skan方程中的分离：
当$\beta <-0.1988$时，边界层发生分离。

分离的影响：

• 阻力急剧增加
• 换热系数下降
• 流动变得不稳定
• 可能产生涡脱落

7.5 三维边界层

对于复杂几何形状（如后掠翼、机身），需要考虑三维边界层效应：

三维边界层特征：

• 存在二次流（垂直于主流方向的流动）
• 边界层厚度沿展向变化
• 可能出现横向分离

控制方程：
在边界层近似下，三维边界层方程包括x和z两个方向的动量方程，以及相应的连续性方程。

---

8. Python数值仿真实现

8.1 Blasius方程求解

# -*- coding: utf-8 -*-
"""
主题005：边界层理论详解
边界层流动数值仿真
"""

import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint, solve_bvp
from scipy.optimize import fsolve
import warnings
warnings.filterwarnings('ignore')
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\对流换热仿真\主题005_边界层理论详解'
os.makedirs(output_dir, exist_ok=True)

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

print("=" * 60)
print("边界层理论详解 - 数值仿真")
print("=" * 60)

# ==================== 案例1: Blasius方程求解 ====================
print("\n案例1: Blasius方程求解")

def blasius_ode(f, eta):
    """
    Blasius方程: f''' + 0.5*f*f'' = 0
    转化为方程组:
    f0' = f1
    f1' = f2
    f2' = -0.5*f0*f2
    """
    f0, f1, f2 = f
    return [f1, f2, -0.5 * f0 * f2]

def solve_blasius(eta_max=10, n_points=1000):
    """
    使用打靶法求解Blasius方程
    """
    eta = np.linspace(0, eta_max, n_points)
    
    # 目标：f'(∞) = 1
    def objective(f2_0):
        f0 = [0, 0, f2_0[0]]  # f(0)=0, f'(0)=0, f''(0)=?
        sol = odeint(blasius_ode, f0, eta)
        return sol[-1, 1] - 1.0  # f'(∞) - 1
    
    # 求解f''(0)
    f2_0_solution = fsolve(objective, [0.3])
    f2_0 = f2_0_solution[0]
    
    print(f"  壁面二阶导数 f''(0) = {f2_0:.6f}")
    
    # 最终求解
    f0 = [0, 0, f2_0]
    sol = odeint(blasius_ode, f0, eta)
    
    f = sol[:, 0]
    fp = sol[:, 1]  # f' = u/U_inf
    fpp = sol[:, 2]  # f''
    
    return eta, f, fp, fpp, f2_0

# 求解Blasius方程
eta, f, fp, fpp, f2_0 = solve_blasius(eta_max=10, n_points=1000)

# 计算边界层厚度（f' = 0.99处）
idx_99 = np.where(fp >= 0.99)[0][0]
eta_99 = eta[idx_99]
print(f"  边界层厚度参数 η_99 = {eta_99:.4f}")

# 计算位移厚度和动量厚度
delta_star = np.trapz(1 - fp, eta)
theta = np.trapz(fp * (1 - fp), eta)
H = delta_star / theta
print(f"  位移厚度积分值 δ*√(U/νx) = {delta_star:.4f}")
print(f"  动量厚度积分值 θ√(U/νx) = {theta:.4f}")
print(f"  形状因子 H = {H:.4f}")

# 绘制Blasius解
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# f(η)
ax1 = axes[0]
ax1.plot(eta, f, 'b-', linewidth=2)
ax1.set_xlabel('η', fontsize=11)
ax1.set_ylabel('f(η)', fontsize=11)
ax1.set_title('Stream Function f(η)', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.set_xlim([0, 10])

# f'(η) = u/U∞
ax2 = axes[1]
ax2.plot(eta, fp, 'b-', linewidth=2, label="f'(η)")
ax2.axhline(y=0.99, color='r', linestyle='--', label='0.99')
ax2.axvline(x=eta_99, color='g', linestyle=':', label=f'η_99={eta_99:.2f}')
ax2.set_xlabel('η', fontsize=11)
ax2.set_ylabel("f'(η) = u/U∞", fontsize=11)
ax2.set_title('Velocity Profile', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)
ax2.set_xlim([0, 10])

# f''(η)
ax3 = axes[2]
ax3.plot(eta, fpp, 'b-', linewidth=2)
ax3.set_xlabel('η', fontsize=11)
ax3.set_ylabel("f''(η)", fontsize=11)
ax3.set_title('Wall Shear Parameter', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.set_xlim([0, 10])

plt.tight_layout()
plt.savefig(f'{output_dir}/fig1_blasius_solution.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图1已保存: Blasius方程解")

# ==================== 案例2: 速度边界层特性分析 ====================
print("\n案例2: 速度边界层特性分析")

# 计算不同位置的速度剖面
U_inf = 1.0  # 主流速度 (m/s)
nu = 1.5e-5  # 运动粘度 (m²/s)，空气
x_positions = [0.1, 0.5, 1.0, 2.0]  # 位置 (m)

fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()

for idx, x in enumerate(x_positions):
    ax = axes[idx]
    
    # 物理坐标 y = η * sqrt(νx/U)
    y = eta * np.sqrt(nu * x / U_inf)
    
    # 边界层厚度
    delta = 5.0 * np.sqrt(nu * x / U_inf)
    delta_star_local = 1.72 * np.sqrt(nu * x / U_inf)
    theta_local = 0.664 * np.sqrt(nu * x / U_inf)
    
    Re_x = U_inf * x / nu
    Cf = 0.664 / np.sqrt(Re_x)
    
    ax.plot(fp, y*1000, 'b-', linewidth=2)  # y单位为mm
    ax.axhline(y=delta*1000, color='r', linestyle='--', label=f'δ={delta*1000:.2f}mm')
    ax.axhline(y=delta_star_local*1000, color='g', linestyle=':', label=f'δ*={delta_star_local*1000:.2f}mm')
    ax.axhline(y=theta_local*1000, color='m', linestyle='-.', label=f'θ={theta_local*1000:.2f}mm')
    
    ax.set_xlabel('u/U∞', fontsize=10)
    ax.set_ylabel('y (mm)', fontsize=10)
    ax.set_title(f'x={x}m, Re_x={Re_x:.1e}, Cf={Cf:.4f}', fontsize=10, fontweight='bold')
    ax.legend(fontsize=8)
    ax.grid(True, alpha=0.3)
    ax.set_ylim([0, 20])

plt.tight_layout()
plt.savefig(f'{output_dir}/fig2_velocity_profiles.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图2已保存: 速度边界层特性")

# ==================== 案例3: 热边界层求解 ====================
print("\n案例3: 热边界层求解")

def solve_thermal_boundary_layer(Pr, eta_max=10, n_points=1000):
    """
    求解热边界层方程: θ'' + 0.5*Pr*f*θ' = 0
    使用Blasius解中的f(η)
    """
    eta = np.linspace(0, eta_max, n_points)
    
    # 插值获取f(η)
    f_interp = np.interp(eta, eta, f)
    
    def thermal_ode(theta, eta_local):
        """热边界层方程"""
        theta0, theta1 = theta
        f_val = np.interp(eta_local, eta, f_interp)
        return [theta1, -0.5 * Pr * f_val * theta1]
    
    # 边界条件: θ(0)=0, θ(∞)=1
    # 打靶法求解θ'(0)
    def objective(theta1_0):
        theta0 = [0, theta1_0[0]]
        sol = odeint(thermal_ode, theta0, eta)
        return sol[-1, 0] - 1.0
    
    theta1_0_solution = fsolve(objective, [0.5])
    theta1_0 = theta1_0_solution[0]
    
    # 最终求解
    theta0 = [0, theta1_0]
    sol = odeint(thermal_ode, theta0, eta)
    theta = sol[:, 0]
    theta_prime = sol[:, 1]
    
    return eta, theta, theta_prime, theta1_0

# 不同Prandtl数的解
Pr_values = [0.7, 1.0, 2.0, 5.0, 10.0]

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 温度分布
ax1 = axes[0]
for Pr in Pr_values:
    eta_t, theta, theta_p, theta_prime_0 = solve_thermal_boundary_layer(Pr)
    ax1.plot(eta_t, theta, linewidth=2, label=f'Pr={Pr}')
    print(f"  Pr={Pr:4.1f}: θ'(0) = {theta_prime_0:.4f}")

ax1.set_xlabel('η', fontsize=11)
ax1.set_ylabel('θ = (T-Tw)/(T∞-Tw)', fontsize=11)
ax1.set_title('Temperature Profiles', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([0, 10])

# 温度梯度对比
ax2 = axes[1]
Pr_range = np.linspace(0.6, 15, 50)
theta_prime_0_values = []

for Pr in Pr_range:
    _, _, _, tp0 = solve_thermal_boundary_layer(Pr, n_points=500)
    theta_prime_0_values.append(tp0)

# Pohlhausen近似: θ'(0) ≈ 0.332*Pr^(1/3)
theta_prime_approx = 0.332 * Pr_range**(1/3)

ax2.plot(Pr_range, theta_prime_0_values, 'b-', linewidth=2, label='Numerical')
ax2.plot(Pr_range, theta_prime_approx, 'r--', linewidth=2, label="0.332Pr^{1/3}")
ax2.set_xlabel('Pr', fontsize=11)
ax2.set_ylabel("θ'(0)", fontsize=11)
ax2.set_title('Wall Temperature Gradient', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)
ax2.set_xscale('log')
ax2.set_yscale('log')

plt.tight_layout()
plt.savefig(f'{output_dir}/fig3_thermal_boundary_layer.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图3已保存: 热边界层求解")

# ==================== 案例4: 速度边界层与热边界层对比 ====================
print("\n案例4: 速度边界层与热边界层对比")

fig, axes = plt.subplots(2, 3, figsize=(15, 8))

for idx, Pr in enumerate([0.7, 1.0, 5.0]):
    # 上排: 边界层厚度对比
    ax_top = axes[0, idx]
    
    eta_t, theta, _, _ = solve_thermal_boundary_layer(Pr)
    
    # 找到θ=0.99的位置
    idx_theta_99 = np.where(theta >= 0.99)[0][0] if np.any(theta >= 0.99) else -1
    eta_theta_99 = eta_t[idx_theta_99] if idx_theta_99 > 0 else 10
    
    ax_top.plot(fp, eta, 'b-', linewidth=2, label='Velocity')
    ax_top.plot(theta, eta_t, 'r-', linewidth=2, label='Temperature')
    ax_top.axhline(y=eta_99, color='b', linestyle='--', alpha=0.5)
    ax_top.axhline(y=eta_theta_99, color='r', linestyle='--', alpha=0.5)
    
    ax_top.set_xlabel('u/U∞ or θ', fontsize=10)
    ax_top.set_ylabel('η', fontsize=10)
    ax_top.set_title(f'Pr={Pr}, δ/δt={eta_99/eta_theta_99:.2f}', fontsize=11, fontweight='bold')
    ax_top.legend(fontsize=9)
    ax_top.grid(True, alpha=0.3)
    ax_top.set_ylim([0, 10])
    
    # 下排: 物理坐标下的对比
    ax_bot = axes[1, idx]
    
    # 对于Pr>1，热边界层更薄，需要调整坐标
    x_ref = 1.0  # 参考位置
    y_vel = eta * np.sqrt(nu * x_ref / U_inf) * 1000  # mm
    y_therm = eta_t * np.sqrt(nu * x_ref / U_inf) * 1000  # mm
    
    ax_bot.plot(fp, y_vel, 'b-', linewidth=2, label='Velocity')
    ax_bot.plot(theta, y_therm, 'r-', linewidth=2, label='Temperature')
    
    ax_bot.set_xlabel('u/U∞ or θ', fontsize=10)
    ax_bot.set_ylabel('y (mm)', fontsize=10)
    ax_bot.set_title(f'Physical Coordinates (x={x_ref}m)', fontsize=11, fontweight='bold')
    ax_bot.legend(fontsize=9)
    ax_bot.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/fig4_boundary_layer_comparison.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图4已保存: 速度边界层与热边界层对比")

# ==================== 案例5: 摩擦系数与换热系数分布 ====================
print("\n案例5: 摩擦系数与换热系数分布")

# 计算沿平板的分布
x_array = np.linspace(0.01, 2.0, 200)  # 避免x=0
Re_x = U_inf * x_array / nu

# 层流摩擦系数
Cf_lam = 0.664 / np.sqrt(Re_x)

# 湍流摩擦系数（从Re_crit=5e5开始）
Re_crit = 5e5
x_crit = Re_crit * nu / U_inf
Cf_turb = np.zeros_like(Re_x)
Cf_turb[Re_x >= Re_crit] = 0.0592 / Re_x[Re_x >= Re_crit]**0.2

# 混合公式（层流+湍流）
Cf_mixed = np.zeros_like(Re_x)
for i, Rex in enumerate(Re_x):
    if Rex < Re_crit:
        Cf_mixed[i] = 0.664 / np.sqrt(Rex)
    else:
        # 从转捩点开始积分
        Cf_mixed[i] = 0.0592 / Rex**0.2

# Nusselt数
Pr = 0.7  # 空气
Nu_lam = 0.332 * np.sqrt(Re_x) * Pr**(1/3)
Nu_turb = np.zeros_like(Re_x)
Nu_turb[Re_x >= Re_crit] = 0.0296 * Re_x[Re_x >= Re_crit]**0.8 * Pr**(1/3)

# Stanton数
St_lam = Nu_lam / (Re_x * Pr)
St_turb = np.zeros_like(Re_x)
St_turb[Re_x >= Re_crit] = Nu_turb[Re_x >= Re_crit] / (Re_x[Re_x >= Re_crit] * Pr)

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

# 摩擦系数
ax1 = axes[0, 0]
ax1.semilogy(Re_x, Cf_lam, 'b-', linewidth=2, label='Laminar')
ax1.semilogy(Re_x[Re_x >= Re_crit], Cf_turb[Re_x >= Re_crit], 'r-', linewidth=2, label='Turbulent')
ax1.axvline(x=Re_crit, color='k', linestyle='--', alpha=0.5, label=f'Re_crit={Re_crit:.0e}')
ax1.set_xlabel('Re_x', fontsize=11)
ax1.set_ylabel('C_f', fontsize=11)
ax1.set_title('Local Friction Coefficient', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3, which='both')

# Nusselt数
ax2 = axes[0, 1]
ax2.semilogy(Re_x, Nu_lam, 'b-', linewidth=2, label='Laminar')
ax2.semilogy(Re_x[Re_x >= Re_crit], Nu_turb[Re_x >= Re_crit], 'r-', linewidth=2, label='Turbulent')
ax2.axvline(x=Re_crit, color='k', linestyle='--', alpha=0.5)
ax2.set_xlabel('Re_x', fontsize=11)
ax2.set_ylabel('Nu_x', fontsize=11)
ax2.set_title('Local Nusselt Number', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3, which='both')

# 换热系数 h
k = 0.026  # 空气导热系数 W/(m·K)
h_lam = Nu_lam * k / x_array
h_turb = np.zeros_like(Re_x)
h_turb[Re_x >= Re_crit] = Nu_turb[Re_x >= Re_crit] * k / x_array[Re_x >= Re_crit]

ax3 = axes[1, 0]
ax3.plot(x_array*100, h_lam, 'b-', linewidth=2, label='Laminar')
ax3.plot(x_array[x_array >= x_crit]*100, h_turb[x_array >= x_crit], 'r-', linewidth=2, label='Turbulent')
ax3.axvline(x=x_crit*100, color='k', linestyle='--', alpha=0.5, label=f'x_crit={x_crit*100:.1f}cm')
ax3.set_xlabel('x (cm)', fontsize=11)
ax3.set_ylabel('h (W/m²·K)', fontsize=11)
ax3.set_title('Heat Transfer Coefficient', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# Reynolds类比验证
ax4 = axes[1, 1]
ax4.loglog(Re_x, St_lam, 'b-', linewidth=2, label='St (Laminar)')
ax4.loglog(Re_x, Cf_lam/2, 'b--', linewidth=2, label='C_f/2 (Laminar)')
ax4.loglog(Re_x[Re_x >= Re_crit], St_turb[Re_x >= Re_crit], 'r-', linewidth=2, label='St (Turbulent)')
ax4.loglog(Re_x[Re_x >= Re_crit], Cf_turb[Re_x >= Re_crit]/2, 'r--', linewidth=2, label='C_f/2 (Turbulent)')
ax4.set_xlabel('Re_x', fontsize=11)
ax4.set_ylabel('St or C_f/2', fontsize=11)
ax4.set_title('Reynolds Analogy', fontsize=12, fontweight='bold')
ax4.legend(fontsize=8)
ax4.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig(f'{output_dir}/fig5_coefficient_distribution.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图5已保存: 摩擦系数与换热系数分布")

# ==================== 案例6: Falkner-Skan方程（楔形流动）====================
print("\n案例6: Falkner-Skan方程（楔形流动）")

def falkner_skan(beta, eta_max=10, n_points=1000):
    """
    求解Falkner-Skan方程: f''' + f*f'' + β(1 - f'²) = 0
    beta > 0: 加速流动
    beta = 0: Blasius流动
    beta < 0: 减速流动（可能出现分离）
    """
    eta = np.linspace(0, eta_max, n_points)
    
    def fs_ode(f, eta_local):
        f0, f1, f2 = f
        return [f1, f2, -f0*f2 - beta*(1 - f1**2)]
    
    # 打靶法
    def objective(f2_0):
        f0 = [0, 0, f2_0[0]]
        sol = odeint(fs_ode, f0, eta)
        return sol[-1, 1] - 1.0
    
    try:
        f2_0_solution = fsolve(objective, [0.3 + 0.2*beta])
        f2_0 = f2_0_solution[0]
        
        f0 = [0, 0, f2_0]
        sol = odeint(fs_ode, f0, eta)
        
        f = sol[:, 0]
        fp = sol[:, 1]
        fpp = sol[:, 2]
        
        return eta, f, fp, fpp, f2_0, True
    except:
        return eta, None, None, None, None, False

# 不同beta值的解
beta_values = [-0.15, -0.1, 0.0, 0.5, 1.0]
colors = ['r', 'orange', 'b', 'g', 'purple']

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for beta, color in zip(beta_values, colors):
    eta_fs, f_fs, fp_fs, fpp_fs, f2_0_fs, success = falkner_skan(beta)
    
    if success and fp_fs is not None:
        label = f'β={beta:.2f}, f\'\'(0)={f2_0_fs:.4f}'
        axes[0].plot(fp_fs, eta_fs, color=color, linewidth=2, label=label)
        axes[1].plot(fpp_fs, eta_fs, color=color, linewidth=2, label=label)

axes[0].set_xlabel("f'(η) = u/Ue", fontsize=11)
axes[0].set_ylabel('η', fontsize=11)
axes[0].set_title('Velocity Profiles for Wedge Flow', fontsize=12, fontweight='bold')
axes[0].legend(fontsize=9)
axes[0].grid(True, alpha=0.3)
axes[0].set_ylim([0, 10])

axes[1].set_xlabel("f''(η)", fontsize=11)
axes[1].set_ylabel('η', fontsize=11)
axes[1].set_title('Shear Stress Profiles', fontsize=12, fontweight='bold')
axes[1].legend(fontsize=9)
axes[1].grid(True, alpha=0.3)
axes[1].set_ylim([0, 10])

plt.tight_layout()
plt.savefig(f'{output_dir}/fig6_falkner_skan.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图6已保存: Falkner-Skan方程解")

# ==================== 案例7: 湍流边界层特性 ====================
print("\n案例7: 湍流边界层特性")

# 对数律速度分布
def turbulent_velocity_profile(y_plus, method='log_law'):
    """
    计算湍流边界层速度分布
    y_plus: 无量纲距离 y*u_tau/nu
    """
    u_plus = np.zeros_like(y_plus)
    
    # 粘性底层 y+ < 5
    viscous_sublayer = y_plus < 5
    u_plus[viscous_sublayer] = y_plus[viscous_sublayer]
    
    # 缓冲层 5 < y+ < 30
    buffer_layer = (y_plus >= 5) & (y_plus < 30)
    # 使用Spalding公式或简单插值
    u_plus[buffer_layer] = 5.0 + (y_plus[buffer_layer] - 5) * (17.5 - 5) / 25
    
    # 对数律层 y+ > 30
    log_layer = y_plus >= 30
    kappa = 0.41
    B = 5.5
    u_plus[log_layer] = (1/kappa) * np.log(y_plus[log_layer]) + B
    
    return u_plus

y_plus = np.logspace(-1, 4, 500)
u_plus = turbulent_velocity_profile(y_plus)

# 1/7幂律
y_delta = np.linspace(0.001, 1, 100)  # y/δ
u_power = y_delta**(1/7)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 对数律
ax1 = axes[0]
ax1.semilogx(y_plus, u_plus, 'b-', linewidth=2, label='Log Law')
ax1.plot([1, 1], [0, 1], 'k--', alpha=0.3)
ax1.plot([5, 5], [0, 5], 'k--', alpha=0.3)
ax1.plot([30, 30], [0, 30], 'k--', alpha=0.3)
ax1.text(2, 15, 'Viscous\nSublayer', fontsize=9, ha='center')
ax1.text(12, 15, 'Buffer\nLayer', fontsize=9, ha='center')
ax1.text(100, 15, 'Logarithmic\nLayer', fontsize=9, ha='center')
ax1.set_xlabel('y+ = y·uτ/ν', fontsize=11)
ax1.set_ylabel('u+ = u/uτ', fontsize=11)
ax1.set_title('Law of the Wall', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3, which='both')
ax1.set_xlim([0.1, 10000])
ax1.set_ylim([0, 30])

# 幂律速度分布
ax2 = axes[1]
ax2.plot(u_power, y_delta, 'b-', linewidth=2, label='1/7 Power Law')
# 层流对比 (Blasius)
ax2.plot(fp[:len(y_delta)*10:10], eta[:len(y_delta)*10:10]/10, 'r--', linewidth=2, label='Laminar (Blasius)')
ax2.set_xlabel('u/Ue', fontsize=11)
ax2.set_ylabel('y/δ', fontsize=11)
ax2.set_title('Turbulent vs Laminar Velocity Profile', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)
ax2.set_ylim([0, 1])

plt.tight_layout()
plt.savefig(f'{output_dir}/fig7_turbulent_boundary_layer.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图7已保存: 湍流边界层特性")

# ==================== 案例8: 边界层积分方法 ====================
print("\n案例8: 边界层积分方法")

def integral_method(x_array, U_inf, nu, method='thwaits'):
    """
    使用积分方法计算边界层发展
    Thwaites方法用于层流边界层
    """
    n_points = len(x_array)
    delta2 = np.zeros(n_points)  # 动量厚度
    H = np.zeros(n_points)  # 形状因子
    Cf = np.zeros(n_points)  # 摩擦系数
    
    # 初始条件（从驻点开始）
    delta2[0] = 0.001  # 很小的初始值
    
    for i in range(1, n_points):
        dx = x_array[i] - x_array[i-1]
        Re_x = U_inf * x_array[i] / nu
        
        # Thwaites方法
        # dδ2/dx + (H+2)δ2/U * dU/dx = Cf/2
        # 对于平板，dU/dx = 0
        
        # 局部解
        if method == 'thwaits':
            # Thwaites近似
            lambda_param = 0.45  # Thwaites参数
            delta2[i] = 0.664 * x_array[i] / np.sqrt(Re_x)
            H[i] = 2.59
            Cf[i] = 0.664 / np.sqrt(Re_x)
    
    return delta2, H, Cf

# 计算边界层发展
x_integral = np.linspace(0.001, 1.0, 100)
delta2, H_shape, Cf_integral = integral_method(x_integral, U_inf, nu)

# 边界层厚度估计
delta_est = delta2 * 5.0 / 0.664  # 近似关系

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

# 动量厚度
ax1 = axes[0, 0]
ax1.plot(x_integral*100, delta2*1000, 'b-', linewidth=2)
ax1.set_xlabel('x (cm)', fontsize=11)
ax1.set_ylabel('θ (mm)', fontsize=11)
ax1.set_title('Momentum Thickness', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)

# 边界层厚度
ax2 = axes[0, 1]
ax2.plot(x_integral*100, delta_est*1000, 'b-', linewidth=2, label='δ')
ax2.plot(x_integral*100, delta2*1000, 'r--', linewidth=2, label='θ')
ax2.set_xlabel('x (cm)', fontsize=11)
ax2.set_ylabel('Thickness (mm)', fontsize=11)
ax2.set_title('Boundary Layer Thickness', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 形状因子
ax3 = axes[1, 0]
ax3.plot(x_integral*100, H_shape, 'b-', linewidth=2)
ax3.axhline(y=2.59, color='r', linestyle='--', label='Blasius H=2.59')
ax3.set_xlabel('x (cm)', fontsize=11)
ax3.set_ylabel('H = δ*/θ', fontsize=11)
ax3.set_title('Shape Factor', fontsize=12, fontweight='bold')
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 摩擦系数
ax4 = axes[1, 1]
Re_x_int = U_inf * x_integral / nu
ax4.semilogy(x_integral*100, Cf_integral, 'b-', linewidth=2, label='Integral Method')
ax4.semilogy(x_integral*100, 0.664/np.sqrt(Re_x_int), 'r--', linewidth=2, label='Blasius')
ax4.set_xlabel('x (cm)', fontsize=11)
ax4.set_ylabel('Cf', fontsize=11)
ax4.set_title('Friction Coefficient', fontsize=12, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig(f'{output_dir}/fig8_integral_method.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 图8已保存: 边界层积分方法")

print("\n" + "=" * 60)
print("所有仿真案例完成!")
print("=" * 60)
print(f"\n输出文件保存在: {output_dir}")
print("\n生成的图表:")
print("  1. fig1_blasius_solution.png - Blasius方程解")
print("  2. fig2_velocity_profiles.png - 速度边界层特性")
print("  3. fig3_thermal_boundary_layer.png - 热边界层求解")
print("  4. fig4_boundary_layer_comparison.png - 速度边界层与热边界层对比")
print("  5. fig5_coefficient_distribution.png - 摩擦系数与换热系数分布")
print("  6. fig6_falkner_skan.png - Falkner-Skan方程解")
print("  7. fig7_turbulent_boundary_layer.png - 湍流边界层特性")
print("  8. fig8_integral_method.png - 边界层积分方法")

---