对流换热仿真-主题016_非牛顿流体对流换热-非牛顿流体对流换热

kkchenjj

35人浏览 · 2026-04-17 21:52:04

kkchenjj · 2026-04-17 21:52:04 发布

主题016：非牛顿流体对流换热

非牛顿流体在化工、食品、石油和生物医学等工业领域中广泛存在，其对流换热特性与牛顿流体有显著差异。理解非牛顿流体的流变特性及其对换热的影响对于相关工程应用至关重要。本教程系统介绍非牛顿流体的分类与流变特性，深入分析幂律流体、宾汉流体和粘弹性流体的对流换热机理。详细推导非牛顿流体边界层方程，建立适用于不同流变模型的努塞尔数关联式。通过管道流动、平板边界层、自然对流等典型案例，展示非牛顿效应对换热性能的影响规律。教程包含完整的Python仿真代码，实现非牛顿流体流动与换热的数值模拟，为相关工程应用提供理论基础和计算工具。

摘要

非牛顿流体在化工、食品、石油和生物医学等工业领域中广泛存在，其对流换热特性与牛顿流体有显著差异。本教程系统介绍非牛顿流体的分类与流变特性，深入分析幂律流体、宾汉流体和粘弹性流体的对流换热机理。详细推导非牛顿流体边界层方程，建立适用于不同流变模型的努塞尔数关联式。通过管道流动、平板边界层、自然对流等典型案例，展示非牛顿效应对换热性能的影响规律。教程包含完整的Python仿真代码，实现非牛顿流体流动与换热的数值模拟，为相关工程应用提供理论基础和计算工具。

关键词

非牛顿流体，幂律模型，宾汉流体，粘弹性，对流换热，流变特性，边界层理论，数值模拟

在这里插入图片描述

1. 引言

1.1 非牛顿流体的定义与重要性

非牛顿流体是指其剪切应力与剪切速率之间不满足线性关系的流体，即粘度随剪切条件而变化。与牛顿流体（如水、空气）不同，非牛顿流体的流变行为更加复杂多样。

在工程实践中，非牛顿流体无处不在：

聚合物溶液和熔体：塑料加工、纤维制造
石油产品：原油输送、钻井液
食品工业：番茄酱、蛋黄酱、巧克力
生物流体：血液、淋巴液、粘液
化工浆料：涂料、陶瓷浆料、纸浆

1.2 非牛顿效应对换热的影响

非牛顿特性显著影响对流换热过程：

粘度变化：剪切变稀或剪切增稠改变速度分布
流动阻力：高粘度导致压降增大，影响泵送功率
边界层特性：速度剖面改变影响热边界层发展
湍流特性：粘弹性抑制或促进湍流，影响换热强度

1.3 本教程内容安排

本教程将从以下几个方面系统阐述非牛顿流体对流换热：

第2节：非牛顿流体分类与流变模型
第3节：非牛顿流体边界层理论
第4节：管道流动换热分析
第5节：外部流动换热特性
第6节：自然对流换热
第7节：Python仿真实现
第8节：工程应用与总结

2. 非牛顿流体分类与流变模型

2.1 非牛顿流体分类

根据流变特性，非牛顿流体可分为以下几类：

2.1.1 纯粘性非牛顿流体

剪切变稀流体（假塑性）：

表观粘度随剪切速率增加而降低
典型例子：聚合物溶液、血液、番茄酱
工程意义：高剪切区粘度降低，有利于流动但可能影响换热

剪切增稠流体（胀塑性）：

表观粘度随剪切速率增加而增大
典型例子：浓淀粉溶液、石英砂悬浮液
工程意义：高剪切区阻力增大，需要更高的泵送功率

宾汉塑性流体：

存在屈服应力，超过屈服应力后才流动
典型例子：牙膏、钻井泥浆、油脂
工程意义：管道中可能存在塞流核心

触变性流体：

表观粘度随剪切时间延长而降低
典型例子：油漆、印刷油墨
工程意义：需要预剪切才能达到稳定状态

震凝性流体：

表观粘度随剪切时间延长而增大
典型例子：某些胶体悬浮液
工程意义：长时间剪切可能导致流动困难

2.1.2 粘弹性流体

同时具有粘性和弹性特性的流体：

记忆效应：应力依赖于变形历史
法向应力差：产生爬杆效应、挤出胀大
典型例子：聚合物熔体、粘弹性表面活性剂溶液

2.1.3 触变性和震凝性流体

触变性：粘度随时间剪切而降低（可恢复）
震凝性：粘度随时间剪切而增加

时间依赖性流体的工程意义：

启动过程需要克服初始高粘度
流动稳定性受剪切历史影响
设计时需要考虑瞬态效应

2.2 常用流变模型

2.2.1 幂律模型（Ostwald-de Waele）

最常用的非牛顿流体模型：

$τ=K(dudy)n=Kγ˙n\tau = K \left(\frac{du}{dy}\right)^n = K \dot{\gamma}^n$

其中：

$τ\tau$ ：剪切应力（Pa）
$K$ ：稠度系数（Pa·sⁿ）
$n$ ：幂律指数（无量纲）
$γ˙=du/dy\dot{\gamma} = du/dy$ ：剪切速率（s⁻¹）

表观粘度：

$μapp=Kγ˙n−1\mu_{app} = K \dot{\gamma}^{n-1}$

流动特性：

$n < 1$ ：剪切变稀（假塑性）
$n = 1$ ：牛顿流体
$n > 1$ ：剪切增稠（胀塑性）

2.2.2 宾汉模型

适用于具有屈服应力的流体：

$τ=τy+μpdudy当τ>τy\tau = \tau_y + \mu_p \frac{du}{dy} \quad \text{当} \quad \tau > \tau_y$

$dudy=0当τ≤τy\frac{du}{dy} = 0 \quad \text{当} \quad \tau \leq \tau_y$

其中：

$τy\tau_y$ ：屈服应力（Pa）
$μp\mu_p$ ：塑性粘度（Pa·s）

2.2.3 Herschel-Bulkley模型

结合宾汉模型和幂律模型：

$τ=τy+Kγ˙n当τ>τy\tau = \tau_y + K \dot{\gamma}^n \quad \text{当} \quad \tau > \tau_y$

适用于大多数真实非牛顿流体。

2.2.4 Carreau模型

描述从牛顿平台到幂律行为的转变：

$μ−μ∞μ0−μ∞=[1+(λγ˙)2](n−1)/2\frac{\mu - \mu_\infty}{\mu_0 - \mu_\infty} = \left[1 + (\lambda \dot{\gamma})^2\right]^{(n-1)/2}$

其中：

$μ0\mu_0$ ：零剪切粘度
$μ∞\mu_\infty$ ：无限剪切粘度
$λ\lambda$ ：特征时间

2.2.5 粘弹性模型

Maxwell模型：

$τ+λdτdt=μdγdt\tau + \lambda \frac{d\tau}{dt} = \mu \frac{d\gamma}{dt}$

Oldroyd-B模型：

更复杂的粘弹性模型，适用于聚合物溶液：

$τ+λ1∇τ=μ(γ˙+λ2∇γ˙)\tau + \lambda_1 \frac{\nabla}{\tau} = \mu \left(\dot{\gamma} + \lambda_2 \frac{\nabla}{\dot{\gamma}}\right)$

其中：

$λ1\lambda_1$ ：应力松弛时间
$λ2\lambda_2$ ：变形延迟时间
$∇τ\frac{\nabla}{\tau}$ ：上随体导数

Giesekus模型：

考虑非线性效应的粘弹性模型：

$τ+λ∇τ+αλμτ⋅τ=μγ˙\tau + \lambda \frac{\nabla}{\tau} + \alpha \frac{\lambda}{\mu} \tau \cdot \tau = \mu \dot{\gamma}$

其中 $α\alpha$ 为迁移参数，控制剪切变稀程度。

Oldroyd-B模型：适用于稀聚合物溶液。

2.3 典型非牛顿流体的流变参数

流体	模型	K (Pa·sⁿ)	n	τᵧ (Pa)
血液	幂律	0.015	0.7	-
番茄酱	幂律	10.0	0.3	-
钻井泥浆	宾汉	-	-	10.0
聚合物溶液	幂律	0.5	0.6	-
原油	幂律	1.0	0.8	-
蛋黄酱	幂律	15.0	0.35	-
巧克力浆	幂律	5.0	0.5	-
纸浆悬浮液	宾汉	-	-	5.0

参数说明：

K值越大，流体越粘稠
n值越小，剪切变稀效应越明显
τᵧ越大，屈服应力效应越显著

温度影响：
大多数非牛顿流体的流变参数对温度敏感，通常温度升高导致K值降低。

3. 非牛顿流体边界层理论

3.1 边界层控制方程

对于幂律流体，二维稳态边界层方程为：

连续性方程：
$∂u∂x+∂v∂y=0\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0$

动量方程：
$u∂u∂x+v∂u∂y=U∞dU∞dx+1ρ∂τxy∂yu\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = U_\infty \frac{dU_\infty}{dx} + \frac{1}{\rho}\frac{\partial \tau_{xy}}{\partial y}$

其中剪切应力：
$τxy=K∣∂u∂y∣n−1∂u∂y\tau_{xy} = K \left|\frac{\partial u}{\partial y}\right|^{n-1} \frac{\partial u}{\partial y}$

能量方程：
$u∂T∂x+v∂T∂y=α∂2T∂y2u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y} = \alpha \frac{\partial^2 T}{\partial y^2}$

边界条件：

壁面处（ $y = 0$ ）： $u = 0$ ， $v = 0$ ， $T=T_w$
边界层外缘（ $\to \infty$ ）： $u=U∞u=U_\infty$ ， $T=T∞T=T_\infty$

非牛顿效应：

幂律指数n影响速度剖面形状
剪切变稀（n<1）使壁面速度梯度增大
剪切增稠（n>1）使壁面速度梯度减小

3.2 相似解方法

对于平板层流，引入相似变量：

$η=y(ρU∞2−nKx)1/(n+1)\eta = y \left(\frac{\rho U_\infty^{2-n}}{K x}\right)^{1/(n+1)}$

定义流函数：
$ψ=(KxU∞2n−1ρ)1/(n+1)f(η)\psi = \left(\frac{K x U_\infty^{2n-1}}{\rho}\right)^{1/(n+1)} f(\eta)$

速度分量：
$U_\infty f'(\eta)$

动量方程转化为：

$n∣f′′∣n−1f′′′+1n+1ff′′=0n|f''|^{n-1}f''' + \frac{1}{n+1}ff'' = 0$

3.3 边界层厚度

对于幂律流体平板边界层：

流动边界层厚度：
$δx=C(n)Rex−1/(n+1)\frac{\delta}{x} = C(n) Re_x^{-1/(n+1)}$

其中修正雷诺数：
$Rex=ρU∞2−nxnKRe_x = \frac{\rho U_\infty^{2-n} x^n}{K}$

与牛顿流体的比较：

剪切变稀流体（n < 1）：边界层更薄，速度梯度更陡峭
剪切增稠流体（n > 1）：边界层更厚，速度分布更平缓

边界层厚度与换热的关系：

边界层厚度减小，壁面温度梯度增大
剪切变稀流体通常具有更高的努塞尔数
边界层厚度比 $δt/δ\delta_t/\delta$ 受普朗特数和幂律指数共同影响

3.4 摩擦系数

局部摩擦系数：
$Cf,x=τw12ρU∞2=2[(n+1)C(n)nRex]1/(n+1)C_{f,x} = \frac{\tau_w}{\frac{1}{2}\rho U_\infty^2} = 2\left[\frac{(n+1)C(n)^n}{Re_x}\right]^{1/(n+1)}$

4. 管道流动换热分析

4.1 充分发展流动

对于幂律流体在圆管中的层流流动，速度分布为：

$U_{max} \left[1 - \left(\frac{r}{R}\right)^{(n+1)/n}\right]$

其中最大速度：
$Umax=(ΔP2KL)1/nnR(n+1)/nn+1U_{max} = \left(\frac{\Delta P}{2KL}\right)^{1/n} \frac{nR^{(n+1)/n}}{n+1}$

平均速度与最大速度之比：
$uˉUmax=n+13n+1\frac{\bar{u}}{U_{max}} = \frac{n+1}{3n+1}$

4.2 努塞尔数关联式

恒壁温条件：

对于幂律流体，充分发展的努塞尔数：

$NuD=C(n)⋅3.66Nu_D = C(n) \cdot 3.66$

其中C(n)为幂律指数修正因子：

n = 0.5时，Nu ≈ 4.2
n = 1.0时，Nu = 3.66
n = 1.5时，Nu ≈ 3.2

入口段效应：

热入口段长度：
$Lt,fdD=0.05ReMRPreff\frac{L_{t,fd}}{D} = 0.05 Re_{MR} Pr_{eff}$

其中Metzner-Reed雷诺数：
$ReMR=ρuˉ2−nDnK′Re_{MR} = \frac{\rho \bar{u}^{2-n} D^n}{K'}$

4.3 宾汉流体管道流动

对于宾汉流体，存在塞流核心（plug flow region）：

塞流核心半径：
$rp=2τyLΔPr_p = \frac{2\tau_y L}{\Delta P}$

速度分布：

在塞流区（r < rₚ）：
$u_p = \text{常数}$

在剪切区（rₚ < r < R）：
$\frac{1}{\mu_p}\left(\frac{\Delta P}{4L}\right)(R - r_p)^2 - \frac{\tau_y}{\mu_p}(R - r)$

换热特性：

塞流核心区温度分布均匀
换热主要发生在剪切区
有效努塞尔数低于牛顿流体

宾汉数的影响：
宾汉数 $\tau_y D / (\mu_p \bar{u})$ 表征屈服应力效应的相对重要性：

Bi = 0：牛顿流体
Bi增大：塞流区扩大，换热减弱
Bi > 10：塞流效应显著，需要特殊设计

工程应用建议：

对于高宾汉数流体，考虑使用脉动流或插入物破坏塞流
增加管长以提高换热效率
采用小直径管道增强换热

5. 外部流动换热特性

5.1 平板对流换热

对于幂律流体平板层流，局部努塞尔数：

$Nu_x = C(n) Re_x^{1/(n+1)} Pr_{eff}^{1/3}$

其中有效普朗特数：
$Preff=Kcpk(U∞x)n−1Pr_{eff} = \frac{K c_p}{k} \left(\frac{U_\infty}{x}\right)^{n-1}$

幂律指数的影响：

n < 1（剪切变稀）：Nu高于牛顿流体
n > 1（剪切增稠）：Nu低于牛顿流体

5.2 绕圆柱换热

非牛顿流体绕圆柱的换热关联式：

$NuD=C⋅ReDm⋅Preff1/3⋅f(n)Nu_D = C \cdot Re_D^m \cdot Pr_{eff}^{1/3} \cdot f(n)$

其中f(n)为幂律修正函数。

流动分离特性：

剪切变稀流体：分离点延后，尾流区减小
剪切增稠流体：分离点提前，尾流区增大
粘弹性流体：可能抑制涡脱落

工程应用：

换热器管束设计需考虑非牛顿效应
聚合物加工中的冷却段设计
生物反应器中的传热优化

—### 5.3 射流冲击换热

非牛顿流体射流冲击的换热特性：

停滞点换热：
$Nustag∝Rej2/(n+1)Preff1/3Nu_{stag} \propto Re_j^{2/(n+1)} Pr_{eff}^{1/3}$

幂律效应：

剪切变稀增强停滞点换热
粘弹性可能降低湍流射流的换热强度

射流冲击的工程应用：

聚合物薄膜的冷却和固化
食品加工中的快速冷却
电子器件的冷却

6. 自然对流换热

6.1 控制方程

非牛顿流体自然对流的边界层方程：

动量方程：
$u∂u∂x+v∂u∂y=gβ(T−T∞)+1ρ∂τxy∂yu\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = g\beta(T - T_\infty) + \frac{1}{\rho}\frac{\partial \tau_{xy}}{\partial y}$

能量方程：
$u∂T∂x+v∂T∂y=α∂2T∂y2u\frac{\partial T}{\partial x} + v\frac{\partial T}{\partial y} = \alpha \frac{\partial^2 T}{\partial y^2}$

6.2 相似解

引入广义Grashof数：
$Grx∗=ρ2gβ(Tw−T∞)x2n+1K2Gr_x^* = \frac{\rho^2 g \beta (T_w - T_\infty) x^{2n+1}}{K^2}$

局部努塞尔数关联式：
$Nu_x = C(n) (Gr_x^* Pr_{eff})^{1/(2n+2)}$

6.3 非牛顿效应

剪切变稀流体：

近壁区粘度降低，速度增加
换热增强，Nu高于牛顿流体

粘弹性流体：

弹性应力影响流动结构
可能抑制或增强自然对流

自然对流的工程应用：

聚合物熔体的冷却过程
地热系统中高粘度流体的流动
食品储存中的热传递

7. Python仿真实现

7.1 案例一：幂律流体管道流动

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

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

print("=" * 60)
print("主题016：非牛顿流体对流换热")
print("=" * 60)

# ============================================
# 案例一：幂律流体管道流动
# ============================================
print("\n案例一：幂律流体管道流动")
print("-" * 50)

def power_law_pipe_velocity(r, R, n, K, dp_dx):
    """
    计算幂律流体在圆管中的速度分布
    
    参数:
        r: 径向位置
        R: 管半径
        n: 幂律指数
        K: 稠度系数
        dp_dx: 压力梯度
    """
    # 最大速度位置
    r_max = R
    
    # 速度分布
    u = np.zeros_like(r)
    for i, ri in enumerate(r):
        if ri <= R:
            # 幂律流体速度分布
            u[i] = (dp_dx / (2 * K))**(1/n) * n / (n + 1) * \
                   (R**((n+1)/n) - ri**((n+1)/n))
    
    return u

def pipe_flow_nusselt(n, Re_MR, Pr_eff):
    """
    计算幂律流体管道流动的努塞尔数
    """
    # 牛顿流体基准值
    Nu_newtonian = 3.66
    
    # 幂律修正因子（经验关联）
    if n < 1:  # 剪切变稀
        correction = 1.0 + 0.2 * (1 - n)
    else:  # 剪切增稠
        correction = 1.0 - 0.15 * (n - 1)
    
    Nu = Nu_newtonian * correction
    
    return Nu

# 参数设置
R = 1.0  # 管半径
dp_dx = -100  # 压力梯度
K = 1.0  # 稠度系数

# 不同幂律指数
n_values = [0.3, 0.5, 0.7, 1.0, 1.3, 1.5]

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

# 速度分布对比
ax1 = axes[0, 0]
r = np.linspace(0, R, 100)
colors = plt.cm.viridis(np.linspace(0, 1, len(n_values)))

for i, n in enumerate(n_values):
    u = power_law_pipe_velocity(r, R, n, K, dp_dx)
    u_max = np.max(u)
    ax1.plot(r/R, u/u_max, color=colors[i], linewidth=2, label=f'n = {n}')

ax1.set_xlabel('r/R', fontsize=11)
ax1.set_ylabel('u/u_max', fontsize=11)
ax1.set_title('Velocity Profiles for Different n', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 流量与压降关系
ax2 = axes[0, 1]
n_range = np.linspace(0.2, 1.8, 50)
Q_ratio = []
for n in n_range:
    # 流量（相对牛顿流体）
    q_ratio = (n + 1) / (3*n + 1) * (2 / (1/n + 1))**n
    Q_ratio.append(q_ratio)

ax2.plot(n_range, Q_ratio, 'b-', linewidth=2)
ax2.axhline(y=1.0, color='r', linestyle='--', alpha=0.5, label='Newtonian')
ax2.axvline(x=1.0, color='gray', linestyle=':', alpha=0.5)
ax2.set_xlabel('Power-law index n', fontsize=11)
ax2.set_ylabel('Relative Flow Rate', fontsize=11)
ax2.set_title('Flow Rate vs n (constant pressure drop)', fontsize=12)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 努塞尔数随n变化
ax3 = axes[1, 0]
Nu_values = []
for n in n_values:
    Nu = pipe_flow_nusselt(n, 1000, 7.0)
    Nu_values.append(Nu)
    print(f"  n = {n:.1f}: Nu = {Nu:.3f}")

bars = ax3.bar(range(len(n_values)), Nu_values, color=colors, alpha=0.7, edgecolor='black')
ax3.axhline(y=3.66, color='r', linestyle='--', linewidth=2, label='Newtonian (n=1)')
ax3.set_xticks(range(len(n_values)))
ax3.set_xticklabels([f'n={n}' for n in n_values], fontsize=9)
ax3.set_ylabel('Nu', fontsize=11)
ax3.set_title('Nusselt Number vs Power-law Index', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3, axis='y')

# 表观粘度分布
ax4 = axes[1, 1]
r_visc = np.linspace(0.1, R, 100)
gamma_dot = np.gradient(power_law_pipe_velocity(r_visc, R, 0.5, K, dp_dx), r_visc)
mu_app_03 = K * np.abs(gamma_dot)**(0.3 - 1)
mu_app_05 = K * np.abs(gamma_dot)**(0.5 - 1)
mu_app_10 = K * np.abs(gamma_dot)**(1.0 - 1)
mu_app_15 = K * np.abs(gamma_dot)**(1.5 - 1)

ax4.semilogy(r_visc/R, mu_app_03, 'b-', linewidth=2, label='n = 0.3')
ax4.semilogy(r_visc/R, mu_app_05, 'g-', linewidth=2, label='n = 0.5')
ax4.semilogy(r_visc/R, mu_app_10, 'r-', linewidth=2, label='n = 1.0')
ax4.semilogy(r_visc/R, mu_app_15, 'm-', linewidth=2, label='n = 1.5')
ax4.set_xlabel('r/R', fontsize=11)
ax4.set_ylabel('Apparent Viscosity (Pa·s)', fontsize=11)
ax4.set_title('Apparent Viscosity Distribution', fontsize=12)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig('case1_power_law_pipe.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 幂律流体管道流动图已保存")

7.2 案例二：宾汉流体流动与换热

# ============================================
# 案例二：宾汉流体流动与换热
# ============================================
print("\n案例二：宾汉流体流动与换热")
print("-" * 50)

def bingham_pipe_velocity(r, R, tau_y, mu_p, dp_dx):
    """
    计算宾汉流体在圆管中的速度分布
    """
    # 塞流核心半径
    r_p = 2 * tau_y / abs(dp_dx)
    r_p = min(r_p, R)
    
    u = np.zeros_like(r)
    
    for i, ri in enumerate(r):
        if ri <= r_p:
            # 塞流区 - 均匀速度
            u[i] = (abs(dp_dx) / (4 * mu_p)) * (R - r_p)**2
        elif ri <= R:
            # 剪切区
            u[i] = (abs(dp_dx) / (4 * mu_p)) * ((R - r_p)**2 - (ri - r_p)**2)
    
    return u, r_p

def bingham_nusselt(tau_y, mu_p, R, dp_dx):
    """
    计算宾汉流体的努塞尔数（简化模型）
    """
    # 塞流核心半径比
    r_p = 2 * tau_y / abs(dp_dx)
    r_p_ratio = min(r_p / R, 1.0)
    
    # 塞流效应降低换热效率
    Nu = 3.66 * (1 - 0.3 * r_p_ratio)
    
    return Nu, r_p_ratio

# 参数设置
R = 1.0
dp_dx = -100
mu_p = 1.0
tau_y_values = [0, 5, 10, 20, 50]

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

# 速度分布
ax1 = axes[0, 0]
r = np.linspace(0, R, 100)
colors = plt.cm.plasma(np.linspace(0, 1, len(tau_y_values)))

plug_ratios = []
for i, tau_y in enumerate(tau_y_values):
    u, r_p = bingham_pipe_velocity(r, R, tau_y, mu_p, dp_dx)
    u_max = np.max(u) if np.max(u) > 0 else 1
    ax1.plot(r/R, u/u_max, color=colors[i], linewidth=2, label=f'τᵧ = {tau_y} Pa')
    if tau_y > 0:
        ax1.axvline(x=r_p/R, color=colors[i], linestyle='--', alpha=0.3)
    plug_ratios.append(r_p/R)

ax1.set_xlabel('r/R', fontsize=11)
ax1.set_ylabel('u/u_max', fontsize=11)
ax1.set_title('Bingham Fluid Velocity Profiles', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 塞流核心大小
ax2 = axes[0, 1]
ax2.plot(tau_y_values, plug_ratios, 'bo-', linewidth=2, markersize=8)
ax2.set_xlabel('Yield Stress τᵧ (Pa)', fontsize=11)
ax2.set_ylabel('Plug Core Radius Ratio rₚ/R', fontsize=11)
ax2.set_title('Plug Flow Region Size', fontsize=12)
ax2.grid(True, alpha=0.3)

# 努塞尔数随屈服应力变化
ax3 = axes[1, 0]
Nu_bingham = []
for tau_y in tau_y_values:
    Nu, r_p_ratio = bingham_nusselt(tau_y, mu_p, R, dp_dx)
    Nu_bingham.append(Nu)
    print(f"  τᵧ = {tau_y:4.1f} Pa: Nu = {Nu:.3f}, rₚ/R = {r_p_ratio:.3f}")

bars = ax3.bar(range(len(tau_y_values)), Nu_bingham, color=colors, alpha=0.7, edgecolor='black')
ax3.axhline(y=3.66, color='r', linestyle='--', linewidth=2, label='Newtonian')
ax3.set_xticks(range(len(tau_y_values)))
ax3.set_xticklabels([f'{ty}' for ty in tau_y_values], fontsize=9)
ax3.set_xlabel('Yield Stress (Pa)', fontsize=11)
ax3.set_ylabel('Nu', fontsize=11)
ax3.set_title('Nusselt Number vs Yield Stress', fontsize=12)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3, axis='y')

# 流动状态图
ax4 = axes[1, 1]
ax4.axis('off')

flow_regime_text = """
Bingham Fluid Flow Regimes:

τ_w > τᵧ: Fluid flows
  • Plug core region (r < rₚ)
    - Uniform velocity
    - No shear deformation
  • Shear region (rₚ < r < R)
    - Velocity gradient
    - Viscous dissipation

τ_w ≤ τᵧ: No flow
  • Static fluid
  • No heat transfer by convection

Heat Transfer Characteristics:
• Plug core reduces mixing
• Lower Nu than Newtonian
• Higher τᵧ → Lower Nu
• Engineering challenge:
  - Avoid high τᵧ in heat exchangers
  - Use pulsating flow to break plug
"""

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

plt.tight_layout()
plt.savefig('case2_bingham_fluid.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 宾汉流体流动与换热图已保存")

7.3 案例三：幂律流体平板边界层

# ============================================
# 案例三：幂律流体平板边界层
# ============================================
print("\n案例三：幂律流体平板边界层")
print("-" * 50)

def power_law_boundary_layer(eta, n):
    """
    求解幂律流体平板边界层的相似解
    简化模型：使用近似速度剖面
    """
    # 近似速度剖面
    if n < 1:  # 剪切变稀
        shape_factor = 1.0 + 0.5 * (1 - n)
    else:  # 剪切增稠或牛顿
        shape_factor = 1.0 - 0.2 * (n - 1)
    
    # 无量纲速度剖面
    f_prime = 1 - np.exp(-eta * shape_factor)
    
    return f_prime

def boundary_layer_thickness(x, n, Re_x):
    """
    计算幂律流体边界层厚度
    """
    # 边界层厚度随幂律指数变化
    if n < 1:
        thickness_factor = 1.0 - 0.3 * (1 - n)
    else:
        thickness_factor = 1.0 + 0.2 * (n - 1)
    
    # 边界层厚度
    delta = 5.0 * x * thickness_factor / (Re_x**(1/(n+1)))
    
    return delta

def local_nusselt_power_law(x, n, Re_x, Pr_eff):
    """
    计算幂律流体平板的局部努塞尔数
    """
    # 幂律修正
    if n < 1:
        n_correction = 1.0 + 0.15 * (1 - n)
    else:
        n_correction = 1.0 - 0.1 * (n - 1)
    
    # 局部努塞尔数
    Nu_x = n_correction * 0.332 * Re_x**(1/(n+1)) * Pr_eff**(1/3)
    
    return Nu_x

# 参数设置
x = np.linspace(0.01, 1.0, 100)
U_inf = 1.0
rho = 1000
K = 0.1
Pr_eff = 7.0

n_values = [0.3, 0.5, 0.7, 1.0, 1.3, 1.5]

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

# 速度剖面对比
ax1 = axes[0, 0]
eta = np.linspace(0, 8, 100)
colors = plt.cm.viridis(np.linspace(0, 1, len(n_values)))

for i, n in enumerate(n_values):
    u_profile = power_law_boundary_layer(eta, n)
    ax1.plot(eta, u_profile, color=colors[i], linewidth=2, label=f'n = {n}')

ax1.set_xlabel(r'$\eta$', fontsize=11)
ax1.set_ylabel(r'$u/U_\infty$', fontsize=11)
ax1.set_title('Boundary Layer Velocity Profiles', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 边界层厚度发展
ax2 = axes[0, 1]
for i, n in enumerate(n_values):
    # 修正雷诺数
    Re_x = rho * U_inf**(2-n) * x**n / K
    delta = boundary_layer_thickness(x, n, Re_x)
    ax2.plot(x, delta/x, color=colors[i], linewidth=2, label=f'n = {n}')

ax2.set_xlabel('x (m)', fontsize=11)
ax2.set_ylabel('δ/x', fontsize=11)
ax2.set_title('Boundary Layer Thickness Development', fontsize=12)
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 局部努塞尔数
ax3 = axes[1, 0]
for i, n in enumerate(n_values):
    Re_x = rho * U_inf**(2-n) * x**n / K
    Nu_x = local_nusselt_power_law(x, n, Re_x, Pr_eff)
    ax1.plot(x, Nu_x, color=colors[i], linewidth=2, label=f'n = {n}')

ax3.set_xlabel('x (m)', fontsize=11)
ax3.set_ylabel('Nu_x', fontsize=11)
ax3.set_title('Local Nusselt Number Distribution', fontsize=12)
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 平均努塞尔数对比
ax4 = axes[1, 1]
Nu_avg_values = []
for n in n_values:
    Re_L = rho * U_inf**(2-n) * 1.0**n / K
    if n < 1:
        correction = 1.0 + 0.2 * (1 - n)
    else:
        correction = 1.0 - 0.15 * (n - 1)
    Nu_avg = correction * 0.664 * Re_L**(1/(n+1)) * Pr_eff**(1/3)
    Nu_avg_values.append(Nu_avg)
    print(f"  n = {n:.1f}: Nu_avg = {Nu_avg:.3f}")

bars = ax4.bar(range(len(n_values)), Nu_avg_values, color=colors, alpha=0.7, edgecolor='black')
ax4.axhline(y=Nu_avg_values[3], color='r', linestyle='--', linewidth=2, label='Newtonian')
ax4.set_xticks(range(len(n_values)))
ax4.set_xticklabels([f'n={n}' for n in n_values], fontsize=9)
ax4.set_ylabel('Nu_L', fontsize=11)
ax4.set_title('Average Nusselt Number vs n', fontsize=12)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3, axis='y')

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

7.4 案例四：非牛顿流体自然对流

# ============================================
# 案例四：非牛顿流体自然对流
# ============================================
print("\n案例四：非牛顿流体自然对流")
print("-" * 50)

def natural_convection_nusselt(n, Gr_star, Pr_eff):
    """
    计算幂律流体自然对流的努塞尔数
    
    参数:
        n: 幂律指数
        Gr_star: 广义Grashof数
        Pr_eff: 有效普朗特数
    """
    # 层流自然对流
    if n < 1:  # 剪切变稀增强换热
        n_factor = 1.0 + 0.25 * (1 - n)
    else:  # 剪切增稠减弱换热
        n_factor = 1.0 - 0.15 * (n - 1)
    
    # 努塞尔数关联式
    Ra_star = Gr_star * Pr_eff
    Nu = n_factor * 0.59 * Ra_star**(1/(2*n+2))
    
    return Nu

# 参数范围
n_values = [0.3, 0.5, 0.7, 1.0, 1.3, 1.5]
Gr_range = np.logspace(4, 8, 50)
Pr_eff = 7.0

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

# 不同n的Nu-Gr关系
ax1 = axes[0, 0]
colors = plt.cm.viridis(np.linspace(0, 1, len(n_values)))

for i, n in enumerate(n_values):
    Nu_values = []
    for Gr in Gr_range:
        Nu = natural_convection_nusselt(n, Gr, Pr_eff)
        Nu_values.append(Nu)
    ax1.loglog(Gr_range, Nu_values, color=colors[i], linewidth=2, label=f'n = {n}')

ax1.set_xlabel(r'$Gr^*$', fontsize=11)
ax1.set_ylabel('Nu', fontsize=11)
ax1.set_title('Natural Convection: Nu vs Gr*', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3, which='both')

# 固定Gr下的Nu随n变化
ax2 = axes[0, 1]
Gr_fixed = 1e6
Nu_vs_n = []
for n in n_values:
    Nu = natural_convection_nusselt(n, Gr_fixed, Pr_eff)
    Nu_vs_n.append(Nu)
    print(f"  n = {n:.1f}: Nu = {Nu:.3f}")

ax2.plot(n_values, Nu_vs_n, 'bo-', linewidth=2, markersize=8)
ax2.axvline(x=1.0, color='r', linestyle='--', alpha=0.5, label='Newtonian')
ax2.set_xlabel('Power-law index n', fontsize=11)
ax2.set_ylabel('Nu', fontsize=11)
ax2.set_title(f'Nu vs n (Gr* = {Gr_fixed:.0e})', fontsize=12)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 边界层温度分布
ax3 = axes[1, 0]
y_eta = np.linspace(0, 5, 100)
for i, n in enumerate([0.5, 0.7, 1.0, 1.3]):
    # 近似温度剖面
    if n < 1:
        delta_t_factor = 1.0 + 0.2 * (1 - n)
    else:
        delta_t_factor = 1.0 - 0.1 * (n - 1)
    
    theta = np.exp(-y_eta / delta_t_factor)
    ax3.plot(y_eta, theta, color=colors[i*2], linewidth=2, label=f'n = {n}')

ax3.set_xlabel(r'$\eta$', fontsize=11)
ax3.set_ylabel(r'$\theta$', fontsize=11)
ax3.set_title('Temperature Profiles in Boundary Layer', fontsize=12)
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3)

# 应用实例对比
ax4 = axes[1, 1]
applications = ['Blood\n(n=0.7)', 'Polymer\n(n=0.6)', 'Ketchup\n(n=0.3)', 
                'Water\n(n=1.0)', 'Clay\n(n=1.2)']
n_app = [0.7, 0.6, 0.3, 1.0, 1.2]
Nu_app = []

for n in n_app:
    Nu = natural_convection_nusselt(n, 1e6, Pr_eff)
    Nu_app.append(Nu)

colors_app = ['red', 'orange', 'darkred', 'blue', 'green']
bars = ax4.bar(range(len(applications)), Nu_app, color=colors_app, alpha=0.7, edgecolor='black')
ax4.axhline(y=Nu_app[3], color='b', linestyle='--', linewidth=2, label='Newtonian (Water)')
ax4.set_xticks(range(len(applications)))
ax4.set_xticklabels(applications, fontsize=8)
ax4.set_ylabel('Nu', fontsize=11)
ax4.set_title('Natural Convection for Different Fluids', fontsize=12)
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3, axis='y')

plt.tight_layout()
plt.savefig('case4_natural_convection.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 非牛顿流体自然对流图已保存")

7.5 案例五：粘弹性流体换热特性

# ============================================
# 案例五：粘弹性流体换热特性
# ============================================
print("\n案例五：粘弹性流体换热特性")
print("-" * 50)

def viscoelastic_nusselt(We, Re, Pr, elasticity_effect='moderate'):
    """
    计算粘弹性流体的努塞尔数
    
    参数:
        We: Weissenberg数（弹性与粘性之比）
        Re: Reynolds数
        Pr: Prandtl数
        elasticity_effect: 弹性效应强度
    """
    # 牛顿流体基准
    Nu_newtonian = 0.023 * Re**0.8 * Pr**0.4
    
    # 弹性效应修正
    if elasticity_effect == 'weak':
        elastic_correction = 1.0 - 0.1 * We
    elif elasticity_effect == 'moderate':
        elastic_correction = 1.0 - 0.2 * We
    else:  # strong
        elastic_correction = 1.0 - 0.35 * We
    
    # 弹性抑制湍流，降低换热
    Nu = Nu_newtonian * max(0.6, elastic_correction)
    
    return Nu

def drag_reduction(We, concentration):
    """
    计算粘弹性流体的减阻率
    """
    # 减阻率随Weissenberg数和浓度增加
    DR = 1 - np.exp(-0.5 * We * concentration)
    return DR

# 参数范围
We_range = np.linspace(0, 5, 50)
Re = 10000
Pr = 7.0

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

# Nu随Weissenberg数变化
ax1 = axes[0, 0]
elasticity_levels = ['weak', 'moderate', 'strong']
colors_elastic = ['green', 'orange', 'red']

for level, color in zip(elasticity_levels, colors_elastic):
    Nu_values = []
    for We in We_range:
        Nu = viscoelastic_nusselt(We, Re, Pr, level)
        Nu_values.append(Nu)
    ax1.plot(We_range, Nu_values, color=color, linewidth=2, label=f'{level.title()} elasticity')

Nu_newt = 0.023 * Re**0.8 * Pr**0.4
ax1.axhline(y=Nu_newt, color='b', linestyle='--', linewidth=2, label='Newtonian')
ax1.set_xlabel('Weissenberg Number', fontsize=11)
ax1.set_ylabel('Nu', fontsize=11)
ax1.set_title('Heat Transfer vs Elasticity', fontsize=12)
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 减阻效应
ax2 = axes[0, 1]
concentrations = [0.1, 0.3, 0.5, 1.0]
colors_conc = plt.cm.viridis(np.linspace(0, 1, len(concentrations)))

for i, conc in enumerate(concentrations):
    DR_values = []
    for We in We_range:
        DR = drag_reduction(We, conc)
        DR_values.append(DR * 100)
    ax2.plot(We_range, DR_values, color=colors_conc[i], linewidth=2, label=f'c = {conc}')

ax2.set_xlabel('Weissenberg Number', fontsize=11)
ax2.set_ylabel('Drag Reduction (%)', fontsize=11)
ax2.set_title('Drag Reduction in Viscoelastic Fluids', fontsize=12)
ax2.legend(fontsize=9)
ax2.grid(True, alpha=0.3)

# 粘弹性流体的速度分布
ax3 = axes[1, 0]
y_wall = np.linspace(0, 1, 100)

# 牛顿流体
u_newtonian = y_wall**0.1

# 粘弹性流体（弹性使速度分布更平坦）
u_viscoelastic_weak = y_wall**0.15
u_viscoelastic_strong = y_wall**0.25

ax3.semilogx(y_wall, u_newtonian, 'b-', linewidth=2, label='Newtonian')
ax3.semilogx(y_wall, u_viscoelastic_weak, 'g-', linewidth=2, label='Weak elastic')
ax3.semilogx(y_wall, u_viscoelastic_strong, 'r-', linewidth=2, label='Strong elastic')
ax3.set_xlabel('Distance from wall', fontsize=11)
ax3.set_ylabel('Normalized velocity', fontsize=11)
ax3.set_title('Velocity Profiles: Newtonian vs Viscoelastic', fontsize=12)
ax3.legend(fontsize=9)
ax3.grid(True, alpha=0.3, which='both')

# 工程应用总结
ax4 = axes[1, 1]
ax4.axis('off')

viscoelastic_text = """
Viscoelastic Fluid Effects:

Heat Transfer:
• Elasticity suppresses turbulence
• Reduces mixing near wall
• Lower Nu than Newtonian
• Effect increases with We

Drag Reduction:
• Toms effect (1948)
• Up to 80% drag reduction
• Polymer additives (ppm level)
• Surfactant solutions

Applications:
• Oil pipeline transport
• District heating/cooling
• Firefighting systems
• Hydraulic fracturing

Design Considerations:
• Balance drag reduction vs heat transfer
• Consider polymer degradation
• Temperature effects on elasticity
• Concentration optimization
"""

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

plt.tight_layout()
plt.savefig('case5_viscoelastic.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 粘弹性流体换热特性图已保存")

7.6 综合对比图

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

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

# 1. 不同非牛顿模型的粘度-剪切速率关系
ax1 = axes[0, 0]
gamma_dot = np.logspace(0, 3, 100)
K = 1.0

# 幂律流体
mu_power_03 = K * gamma_dot**(0.3 - 1)
mu_power_10 = K * gamma_dot**(1.0 - 1)
mu_power_15 = K * gamma_dot**(1.5 - 1)

# Carreau模型
mu_carreau = 1.0 + (10 - 1) * (1 + (0.1 * gamma_dot)**2)**(-0.35)

ax1.loglog(gamma_dot, mu_power_03, 'b-', linewidth=2, label='Power-law n=0.3')
ax1.loglog(gamma_dot, mu_power_10, 'r-', linewidth=2, label='Newtonian')
ax1.loglog(gamma_dot, mu_power_15, 'g-', linewidth=2, label='Power-law n=1.5')
ax1.loglog(gamma_dot, mu_carreau, 'm--', linewidth=2, label='Carreau')
ax1.set_xlabel('Shear Rate (1/s)', fontsize=10)
ax1.set_ylabel('Viscosity (Pa·s)', fontsize=10)
ax1.set_title('Rheological Models Comparison', fontsize=11)
ax1.legend(fontsize=8)
ax1.grid(True, alpha=0.3, which='both')

# 2. Nu随幂律指数变化（综合）
ax2 = axes[0, 1]
n_range = np.linspace(0.2, 1.8, 50)

# 管道流动
Nu_pipe = [3.66 * (1 + 0.2 * (1 - n) if n < 1 else 1 - 0.15 * (n - 1)) for n in n_range]
# 平板流动
Nu_plate = [0.5 * (1 + 0.25 * (1 - n) if n < 1 else 1 - 0.1 * (n - 1)) for n in n_range]
# 自然对流
Nu_natural = [0.6 * (1 + 0.3 * (1 - n) if n < 1 else 1 - 0.2 * (n - 1)) for n in n_range]

ax2.plot(n_range, Nu_pipe, 'b-', linewidth=2, label='Pipe flow')
ax2.plot(n_range, Nu_plate, 'r--', linewidth=2, label='Flat plate')
ax2.plot(n_range, Nu_natural, 'g:', linewidth=2, label='Natural convection')
ax2.axvline(x=1.0, color='k', linestyle='-', alpha=0.3)
ax2.set_xlabel('Power-law index n', fontsize=10)
ax2.set_ylabel('Normalized Nu', fontsize=10)
ax2.set_title('Nu vs n for Different Flows', fontsize=11)
ax2.legend(fontsize=8)
ax2.grid(True, alpha=0.3)

# 3. 不同流体的流变分类
ax3 = axes[0, 2]
fluids_data = {
    'Blood': {'n': 0.7, 'type': 'Shear-thinning'},
    'Ketchup': {'n': 0.3, 'type': 'Shear-thinning'},
    'Polymer': {'n': 0.6, 'type': 'Shear-thinning'},
    'Water': {'n': 1.0, 'type': 'Newtonian'},
    'Clay': {'n': 1.2, 'type': 'Shear-thickening'},
    'Starch': {'n': 1.5, 'type': 'Shear-thickening'}
}

x_pos = np.arange(len(fluids_data))
n_values_plot = [fluids_data[f]['n'] for f in fluids_data]
colors_fluid = ['red' if n < 1 else 'blue' if n == 1 else 'green' for n in n_values_plot]

bars = ax3.bar(x_pos, n_values_plot, color=colors_fluid, alpha=0.7, edgecolor='black')
ax3.axhline(y=1.0, color='k', linestyle='--', linewidth=2, label='Newtonian')
ax3.set_xticks(x_pos)
ax3.set_xticklabels(fluids_data.keys(), fontsize=8, rotation=45)
ax3.set_ylabel('Power-law index n', fontsize=10)
ax3.set_title('Rheological Classification', fontsize=11)
ax3.legend(fontsize=8)
ax3.grid(True, alpha=0.3, axis='y')

# 4. 工程应用分布
ax4 = axes[1, 0]
applications = ['Polymer\nProcessing', 'Oil\nTransport', 'Food\nIndustry', 
                'Biomedical', 'Chemical\nProcesses']
n_app_range = [0.6, 0.8, 0.4, 0.7, 0.5]
colors_app = plt.cm.plasma(np.linspace(0, 1, len(applications)))

bars = ax4.barh(applications, n_app_range, color=colors_app, alpha=0.7, edgecolor='black')
ax4.axvline(x=1.0, color='r', linestyle='--', linewidth=2, label='Newtonian')
ax4.set_xlabel('Typical n value', fontsize=10)
ax4.set_title('Non-Newtonian Fluids in Industry', fontsize=11)
ax4.legend(fontsize=8)
ax4.grid(True, alpha=0.3, axis='x')

# 5. 换热增强/减弱效果
ax5 = axes[1, 1]
n_effect = np.linspace(0.2, 1.8, 50)
enhancement = [(1 + 0.25 * (1 - n) if n < 1 else 1 - 0.15 * (n - 1)) for n in n_effect]

ax5.plot(n_effect, enhancement, 'b-', linewidth=2)
ax5.axhline(y=1.0, color='k', linestyle='--', alpha=0.5, label='Newtonian')
ax5.fill_between(n_effect, 1.0, enhancement, where=[e > 1 for e in enhancement], 
                  alpha=0.3, color='green', label='Enhanced')
ax5.fill_between(n_effect, 1.0, enhancement, where=[e < 1 for e in enhancement], 
                  alpha=0.3, color='red', label='Reduced')
ax5.set_xlabel('Power-law index n', fontsize=10)
ax5.set_ylabel('Heat Transfer Ratio', fontsize=10)
ax5.set_title('Non-Newtonian Effect on Heat Transfer', fontsize=11)
ax5.legend(fontsize=8)
ax5.grid(True, alpha=0.3)

# 6. 关键公式总结
ax6 = axes[1, 2]
ax6.axis('off')

formula_text = r"""
Key Equations:

Power-Law Model:
τ = K(du/dy)^n
μ_app = K(du/dy)^(n-1)

Bingham Model:
τ = τᵧ + μ_p(du/dy) for τ > τᵧ

Pipe Flow Velocity:
u = (ΔP/2KL)^(1/n) · n/(n+1) ·
    [R^((n+1)/n) - r^((n+1)/n)]

Nusselt Number:
Nu = C(n) · Nu_N
where C(n) = 1 + 0.2(1-n) for n<1

Generalized Reynolds:
Re_MR = ρu^(2-n)D^n/K'

Natural Convection:
Nu = C(n)(Gr*·Pr)^(1/(2n+2))
"""

ax6.text(0.05, 0.95, formula_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_power_law_pipe.png - 幂律流体管道流动")
print("  2. case2_bingham_fluid.png - 宾汉流体流动与换热")
print("  3. case3_power_law_bl.png - 幂律流体平板边界层")
print("  4. case4_natural_convection.png - 非牛顿流体自然对流")
print("  5. case5_viscoelastic.png - 粘弹性流体换热特性")
print("  6. comprehensive_comparison.png - 综合对比")
print("\n输出目录：d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题016_非牛顿流体对流换热")
print("=" * 60)