第二十一篇：平板导热：理论与数值对比

摘要

平板导热是传热学中最基础且应用广泛的问题之一。本文系统分析了一维稳态平板导热的理论解与数值解的对比验证。详细推导了单层平板、多层复合平板以及含内热源平板的温度分布解析解。采用有限差分法和有限体积法建立数值模型，通过网格收敛性研究验证数值方法的准确性。定量分析了不同边界条件（Dirichlet、Neumann、Robin）对温度分布的影响。通过Python仿真实现了理论解与数值解的对比，验证了数值方法的可靠性，为后续复杂导热问题的数值模拟奠定了基础。

关键词

平板导热，理论解，数值解，有限差分法，网格收敛性，边界条件，温度分布，误差分析

---

1. 引言

1.1 平板导热的重要性

平板导热模型在工程中无处不在：

• 建筑墙体保温设计
• 电子器件散热基板
• 换热器壁面传热
• 隔热材料性能评估

1.2 研究目标

本文旨在：

1. 建立平板导热的理论分析框架
2. 开发可靠的数值求解方法
3. 验证数值解的准确性
4. 分析误差来源和收敛特性

---

2. 理论分析

2.1 单层平板稳态导热

控制方程：
$$\frac{d^{2}T}{d{x}^2}=0$$

通解：
$$T(x)=C_{1}x+C_2$$

边界条件：

• $T(0)=T_1$
• $T(L)=T_2$

温度分布：
$$T(x)=T_1+(T_2-T_1)\frac{x}{L}$$

热流密度（常数）：
$$q=-k\frac{dT}{dx}=\frac{k(T_1-T_2)}{L}$$

2.2 多层复合平板

对于 $n$ 层材料，第 $i$ 层的温度分布：
$$T_i(x)=T_{i-1}+(T_i-T_{i-1})\frac{x-x_{i-1}}{L_i}$$

总热阻：
$$R_{total}=\sum _{i=1}^n\frac{L_i}{k_i}$$

总热流：
$$q=\frac{T_1-T_{n+1}}{R_{total}}$$

2.3 含内热源的平板

控制方程：
$$\frac{d^{2}T}{d{x}^2}+\frac{\dot{q}}{k}=0$$

通解：
$$T(x)=-\frac{\dot{q}}{2k}{x}^2+C_{1}x+C_2$$

对称边界条件（中心绝热）：

• $\frac{dT}{dx}{|}_{x=0}=0$
• $T(L)=T_w$

温度分布：
$$T(x)=T_w+\frac{\dot{q}{L}^2}{2k}\left(1-\frac{x^2}{L^2}\right)$$

最高温度（中心）：
$$T_{max}=T_w+\frac{\dot{q}{L}^2}{2k}$$

---

3. 数值方法

3.1 有限差分离散

内部节点（中心差分）：
$$\frac{T_{i+1}-2T_i+T_{i-1}}{\Delta x^2}=0$$

边界条件离散：

• Dirichlet：$T_0=T_1$（已知）
• Neumann：$\frac{T_1-T_{-1}}{2\Delta x}=0$
• Robin：$-k\frac{T_1-T_0}{\Delta x}=h(T_0-T_{\infty })$

3.2 有限体积法

对控制体积分：
$${\int }_{x_w}^{x_e}\frac{d^{2}T}{d{x}^2}dx={\frac{dT}{dx}|}_{x_e}-{\frac{dT}{dx}|}_{x_w}=0$$

界面热流连续：
$$q_e=q_w$$

3.3 网格收敛性分析

Richardson外推：
$$f_{exact}\approx f_h+\frac{f_h-f_{2h}}{2^p-1}$$

其中 $p$ 为方法精度阶数。

---

4. Python仿真实现

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
import os

output_dir = r'd:\文档\500仿真领域\工程仿真\传热学仿真\主题021'
os.makedirs(output_dir, exist_ok=True)

print("="*60)
print("仿真1：单层平板导热 - 理论vs数值")
print("="*60)

# 参数
L = 0.1  # 平板厚度，m
k = 50   # 导热系数，W/(m·K)
T1 = 100  # 左边界温度，°C
T2 = 20   # 右边界温度，°C

# 理论解
x_exact = np.linspace(0, L, 1000)
T_exact = T1 + (T2 - T1) * x_exact / L
q_exact = k * (T1 - T2) / L

print(f"理论热流密度: {q_exact:.2f} W/m²")

# 数值解 - 不同网格
grid_sizes = [11, 21, 41, 81, 161]
numerical_results = {}

for N in grid_sizes:
    dx = L / (N - 1)
    x_num = np.linspace(0, L, N)
    
    # 构建三对角矩阵
    main = np.ones(N) * 2
    off = np.ones(N-1) * (-1)
    A = diags([off, main, off], [-1, 0, 1], format='csr')
    
    # 右端项
    b = np.zeros(N)
    
    # 边界条件
    A[0, 0] = 1
    A[0, 1] = 0
    b[0] = T1
    
    A[-1, -1] = 1
    A[-1, -2] = 0
    b[-1] = T2
    
    # 求解
    T_num = spsolve(A, b)
    
    # 计算热流（中心差分）
    q_num = -k * (T_num[1] - T_num[0]) / dx
    
    numerical_results[N] = (x_num, T_num, q_num)
    
    error = abs(q_num - q_exact) / q_exact * 100
    print(f"网格 {N}: 数值热流 = {q_num:.2f} W/m², 误差 = {error:.4f}%")

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

# 温度分布对比
ax1 = axes[0]
ax1.plot(x_exact*1000, T_exact, 'k-', linewidth=2, label='Analytical')
for N, (x_num, T_num, _) in numerical_results.items():
    ax1.plot(x_num*1000, T_num, 'o-', markersize=3, label=f'N={N}')
ax1.set_xlabel('x (mm)', fontsize=11)
ax1.set_ylabel('Temperature (°C)', fontsize=11)
ax1.set_title('Temperature Distribution Comparison', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# 网格收敛性
ax2 = axes[1]
N_list = list(numerical_results.keys())
q_errors = [abs(numerical_results[N][2] - q_exact) / q_exact * 100 for N in N_list]
ax2.loglog(N_list, q_errors, 'bo-', linewidth=2, markersize=8)
dx_list = [L/(N-1)*1000 for N in N_list]
ax2.loglog(N_list, [100*dx**2 for dx in dx_list], 'r--', label='2nd Order Slope')
ax2.set_xlabel('Grid Size N', fontsize=11)
ax2.set_ylabel('Relative Error (%)', fontsize=11)
ax2.set_title('Grid Convergence Study', fontsize=12, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3, which='both')

plt.tight_layout()
plt.savefig(f'{output_dir}/plate_conduction_comparison.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n图1：平板导热对比已保存")

print("\n" + "="*60)
print("仿真2：含内热源平板")
print("="*60)

# 含内热源
q_dot = 1e6  # 内热源，W/m³

# 理论解
T_w = 20  # 壁面温度
x_source = np.linspace(0, L, 1000)
T_source_exact = T_w + q_dot * L**2 / (2*k) * (1 - (x_source/L)**2)
T_max_exact = T_w + q_dot * L**2 / (2*k)

print(f"理论最高温度: {T_max_exact:.2f}°C")

# 数值解
N = 81
dx = L / (N - 1)
x_num = np.linspace(0, L, N)

main = np.ones(N) * 2
off = np.ones(N-1) * (-1)
A = diags([off, main, off], [-1, 0, 1], format='csr')

# 右端项（内热源）
b = np.ones(N) * q_dot * dx**2 / k

# 边界条件（对称+Dirichlet）
A[0, 0] = 1
A[0, 1] = -1  # 绝热边界
b[0] = 0

A[-1, -1] = 1
A[-1, -2] = 0
b[-1] = T_w

T_source_num = spsolve(A, b)
T_max_num = np.max(T_source_num)

print(f"数值最高温度: {T_max_num:.2f}°C")
print(f"温度误差: {abs(T_max_num - T_max_exact):.4f}°C")

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(x_source*1000, T_source_exact, 'k-', linewidth=2, label='Analytical')
ax.plot(x_num*1000, T_source_num, 'ro', markersize=4, label='Numerical (N=81)')
ax.set_xlabel('x (mm)', fontsize=11)
ax.set_ylabel('Temperature (°C)', fontsize=11)
ax.set_title('Temperature Distribution with Internal Heat Source', fontsize=12, fontweight='bold')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/plate_with_source.png', dpi=150, bbox_inches='tight')
plt.close()

print("图2：含内热源平板已保存")

print("\n所有仿真完成！")

---

5. 误差分析与讨论

5.1 截断误差

有限差分法的截断误差：
$$\tau =\frac{\Delta x^2}{12}\frac{{\partial }^{4}T}{\partial x^4}+O(\Delta x^4)$$

对于线性温度分布，理论误差为零。

5.2 舍入误差

计算机浮点数精度限制：

• 单精度：约7位有效数字
• 双精度：约15位有效数字

5.3 收敛阶数验证

通过网格细化验证二阶收敛：
$$\frac{\Vert e_{2h}\Vert }{\Vert e_h\Vert }\approx 4$$

---

6. 边界条件影响分析

6.1 不同边界条件组合

第一类边界条件（Dirichlet）：
$$T(0)=T_1,T(L)=T_2$$
温度分布为线性，与材料属性无关。

第二类边界条件（Neumann）：
$$-k\frac{dT}{dx}{|}_{x=0}=q_1,-k\frac{dT}{dx}{|}_{x=L}=q_2$$
需要满足$q_1=q_2$才能保证稳态解存在。

第三类边界条件（Robin）：
$$-k\frac{dT}{dx}{|}_{x=0}=h_1(T_1-T_{\infty 1})$$
$$-k\frac{dT}{dx}{|}_{x=L}=h_2(T_2-T_{\infty 2})$$

混合边界条件：
一端固定温度，一端对流换热，工程中最常见。

6.2 接触热阻影响

实际工程中，界面接触不完美，存在接触热阻：
$$R_{contact}=\frac{\Delta T_{interface}}{q}$$

总热阻：
$$R_{total}=R_{cond}+R_{contact}=\frac{L}{kA}+R_{contact}$$

减小接触热阻的方法：

1. 增加接触压力
2. 使用导热脂/导热垫
3. 表面抛光处理
4. 焊接或钎焊连接

6.3 变导热系数处理

当导热系数随温度变化时：
$$k(T)=k_0(1+\beta T)$$

Kirchhoff变换：
定义新变量$\theta ={\int }_0^{T}k(T^{\prime })d{T}^{\prime }$，将方程线性化。

数值处理：
采用迭代法或局部线性化处理。

---

7. 多层平板导热

7.1 串联热阻

对于n层材料串联：
$$R_{total}=\sum _{i=1}^n\frac{L_i}{k_{i}A}$$

界面温度：
$$T_i=T_1-q\sum _{j=1}^{i-1}{R}_j$$

7.2 并联热阻

热量分流：
$$q_{total}=q_1+q_2=\frac{T_1-T_2}{R_1}+\frac{T_1-T_2}{R_2}$$

等效热阻：
$$\frac{1}{R_{eq}}=\frac{1}{R_1}+\frac{1}{R_2}$$

7.3 工程应用

复合墙体：

• 保温层（低k）
• 结构层（高k）
• 饰面层（保护）

热桥分析：
局部高导热路径导致的热流集中。

---

8. 瞬态平板导热

8.1 控制方程

$$\rho c_p\frac{\partial T}{\partial t}=k\frac{{\partial }^{2}T}{\partial x^2}$$

无量纲化：
$$\theta =\frac{T-T_{\infty }}{T_i-T_{\infty }},\tau =\frac{\alpha t}{L^2},X=\frac{x}{L}$$

8.2 解析解

分离变量法：
$$\theta (X,\tau )=\sum _{n=1}^{\infty }{C}_n\sin ({\lambda }_{n}X)e^{-{\lambda }_n^2\tau }$$

特征值：
$${\lambda }_n=(2n-1)\frac{\pi }{2}$$

8.3 半无限大近似

短时间或厚板：
$$\frac{T-T_i}{T_s-T_i}=\text{erfc}\left(\frac{x}{2\sqrt{\alpha t}}\right)$$

---

9. 数值方法深入

9.1 有限差分格式对比

格式	精度	稳定性	适用性
中心差分	二阶	无条件	稳态
向前差分	一阶	有条件	瞬态
向后差分	一阶	无条件	瞬态
Crank-Nicolson	二阶	无条件	瞬态

9.2 边界条件离散

虚拟节点法：

对于Neumann边界$\frac{dT}{dx}=0$：
$$\frac{T_1-T_{-1}}{2\Delta x}=0\Rightarrow T_{-1}=T_1$$

代入内部节点方程消去虚拟节点。

9.3 网格无关性验证

Richardson外推：
$$T_{exact}\approx \frac{4T_{h/2}-T_h}{3}$$

网格细化策略：

• 从粗网格开始
• 逐步细化
• 比较结果差异
• 确定收敛网格

---

10. 本章小结

本章系统对比了平板导热的理论解与数值解：

核心理论

1. 理论解特点：

   • 无内热源：线性分布
   • 有内热源：抛物线分布
   • 解析解提供基准验证

2. 数值方法：

   • 有限差分法
   • 二阶精度
   • 边界条件处理

3. 误差分析：

   • 截断误差$O(\Delta x^2)$
   • 舍入误差
   • 网格收敛性

扩展内容

1. 边界条件：Dirichlet、Neumann、Robin及其组合
2. 接触热阻：界面效应及减小方法
3. 多层平板：串联/并联热阻网络
4. 瞬态导热：分离变量法与半无限大近似
5. 数值深入：不同格式对比与网格验证

工程价值

1. 设计指导：平板是工程中最基本的几何形状
2. 验证基准：用于检验数值代码正确性
3. 简化分析：复杂问题可局部简化为平板模型
4. 教学价值：理解导热基本规律的最佳起点

平板导热理论是传热学的基础，掌握其解析与数值方法对于解决复杂工程问题至关重要。

• 简单问题的精确解
• 验证数值方法的标准
• 物理洞察清晰

2. 数值解优势：

   • 处理复杂几何
   • 非线性问题
   • 非均匀材料

3. 验证方法：

   • 网格收敛性研究
   • Richardson外推
   • 误差阶数验证

4. 工程应用：

   • 数值方法可靠
   • 合理选择网格
   • 关注边界条件处理

理论-数值对比是验证计算模型可靠性的标准流程，在工程实践中具有重要意义。