传热学仿真-主题003-导热微分方程的推导

kkchenjj

418人浏览 · 2026-03-17 21:24:00

kkchenjj · 2026-03-17 21:24:00 发布

第3篇：导热微分方程的推导

摘要

导热微分方程是描述热量传递过程的数学基础，是传热学理论体系的核心组成部分。本教程系统介绍导热微分方程的完整推导过程，从能量守恒定律和傅里叶定律出发，逐步建立描述热量传递的数学模型。教程详细讨论能量平衡分析方法，阐述温度梯度、热流密度与内能变化之间的物理关系，推导三维非稳态导热微分方程的一般形式。通过Python数值模拟，展示一维非稳态导热问题的求解过程，分析不同时间步长和空间网格对数值解精度的影响，验证显式格式和隐式格式的稳定性条件。本教程为读者深入理解传热学理论、掌握数值模拟方法提供系统性的知识框架和实践指导。

关键词

导热微分方程，能量守恒定律，傅里叶定律，非稳态导热，热扩散系数，有限差分法，显式格式，隐式格式

在这里插入图片描述

一、导热微分方程的理论基础

1.1 能量守恒定律的数学表达

能量守恒定律是自然界最基本的物理定律之一，也是推导导热微分方程的出发点。该定律指出：在一个孤立系统中，能量既不会凭空产生，也不会凭空消失，它只会从一种形式转化为另一种形式，或者从一个物体转移到另一个物体，而能量的总量保持不变。

对于传热问题，热力学第一定律可以表述为：单位时间内物体热力学能的增加等于单位时间内传入物体的热量与物体内部热源产生的热量之和。用数学语言描述，考虑一个体积为 $V$ 、表面积为 $A$ 的控制体，能量守恒方程可以写成：

$dEdt=Q˙in−Q˙out+Q˙gen\frac{dE}{dt} = \dot{Q}_{in} - \dot{Q}_{out} + \dot{Q}_{gen}$

其中， $dEdt\frac{dE}{dt}$ 表示控制体内能量的变化率， $Q˙in\dot{Q}_{in}$ 和 $Q˙out\dot{Q}_{out}$ 分别表示通过控制面传入和传出的热流量， $Q˙gen\dot{Q}_{gen}$ 表示控制体内热源产生的热流量。

在导热问题中，控制体的能量主要表现为内能，内能的变化与温度变化直接相关。对于固体材料，内能的变化可以表示为：

$dEdt=∫Vρcp∂T∂tdV\frac{dE}{dt} = \int_V \rho c_p \frac{\partial T}{\partial t} dV$

这里， $ρ\rho$ 是材料密度， $c_p$ 是比定压热容， $T$ 是温度， $t$ 是时间。这个积分表达式表明，控制体内能量的变化率等于各点温度变化率的体积分。

1.2 微元体能量平衡分析

为了建立导热微分方程，我们需要从微元体的角度进行能量平衡分析。考虑一个位于连续介质中的微小矩形六面体，其边长分别为 $d x$ 、 $d y$ 、 $d z$ ，体积为 $\cdot dy \cdot dz$ 。

根据傅里叶定律，热量通过导热方式传递，热流密度矢量 $q\mathbf{q}$ 与温度梯度成正比：

$q=−k∇T\mathbf{q} = -k \nabla T$

其中， $k$ 是材料的热导率， $∇T\nabla T$ 是温度梯度。负号表示热量从高温区域流向低温区域。

在 $x$ 方向上，通过微元体左侧面（位于 $x$ 处）传入的热流量为：

$Q˙x=qx⋅dy⋅dz\dot{Q}_x = q_x \cdot dy \cdot dz$

通过右侧面（位于 $x + d x$ 处）传出的热流量为：

$Q˙x+dx=(qx+∂qx∂xdx)⋅dy⋅dz\dot{Q}_{x+dx} = \left(q_x + \frac{\partial q_x}{\partial x}dx\right) \cdot dy \cdot dz$

因此，在 $x$ 方向上的净热流量为：

$Q˙x−Q˙x+dx=−∂qx∂xdx⋅dy⋅dz\dot{Q}_x - \dot{Q}_{x+dx} = -\frac{\partial q_x}{\partial x}dx \cdot dy \cdot dz$

同理，在 $y$ 方向和 $z$ 方向上的净热流量分别为：

$Q˙y−Q˙y+dy=−∂qy∂ydx⋅dy⋅dz\dot{Q}_y - \dot{Q}_{y+dy} = -\frac{\partial q_y}{\partial y}dx \cdot dy \cdot dz$

$Q˙z−Q˙z+dz=−∂qz∂zdx⋅dy⋅dz\dot{Q}_z - \dot{Q}_{z+dz} = -\frac{\partial q_z}{\partial z}dx \cdot dy \cdot dz$

三个方向上的净热流量之和就是通过导热传入微元体的总热量：

$Q˙net=−(∂qx∂x+∂qy∂y+∂qz∂z)dx⋅dy⋅dz\dot{Q}_{net} = -\left(\frac{\partial q_x}{\partial x} + \frac{\partial q_y}{\partial y} + \frac{\partial q_z}{\partial z}\right)dx \cdot dy \cdot dz$

1.3 导热微分方程的推导

将傅里叶定律代入能量平衡方程。由于：

$qx=−k∂T∂x,qy=−k∂T∂y,qz=−k∂T∂zq_x = -k\frac{\partial T}{\partial x}, \quad q_y = -k\frac{\partial T}{\partial y}, \quad q_z = -k\frac{\partial T}{\partial z}$

因此：

$∂qx∂x=−∂∂x(k∂T∂x)\frac{\partial q_x}{\partial x} = -\frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial x}\right)$

$∂qy∂y=−∂∂y(k∂T∂y)\frac{\partial q_y}{\partial y} = -\frac{\partial}{\partial y}\left(k\frac{\partial T}{\partial y}\right)$

$∂qz∂z=−∂∂z(k∂T∂z)\frac{\partial q_z}{\partial z} = -\frac{\partial}{\partial z}\left(k\frac{\partial T}{\partial z}\right)$

如果考虑微元体内存在内热源（如电阻加热、核反应等），单位体积的热源功率为 $q˙\dot{q}$ ，则热源产生的热量为 $q˙⋅dx⋅dy⋅dz\dot{q} \cdot dx \cdot dy \cdot dz$ 。

根据能量守恒，微元体内能的增加等于传入的净热量加上内热源产生的热量：

$ρcp∂T∂tdx⋅dy⋅dz=[∂∂x(k∂T∂x)+∂∂y(k∂T∂y)+∂∂z(k∂T∂z)+q˙]dx⋅dy⋅dz\rho c_p \frac{\partial T}{\partial t}dx \cdot dy \cdot dz = \left[\frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial x}\right) + \frac{\partial}{\partial y}\left(k\frac{\partial T}{\partial y}\right) + \frac{\partial}{\partial z}\left(k\frac{\partial T}{\partial z}\right) + \dot{q}\right]dx \cdot dy \cdot dz$

两边除以 $dx \cdot dy \cdot dz$ ，得到导热微分方程的一般形式：

$ρcp∂T∂t=∂∂x(k∂T∂x)+∂∂y(k∂T∂y)+∂∂z(k∂T∂z)+q˙\rho c_p \frac{\partial T}{\partial t} = \frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial x}\right) + \frac{\partial}{\partial y}\left(k\frac{\partial T}{\partial y}\right) + \frac{\partial}{\partial z}\left(k\frac{\partial T}{\partial z}\right) + \dot{q}$

这就是描述热量传递过程的基本方程，也称为热扩散方程或傅里叶-基尔霍夫方程。

1.4 特殊形式的导热方程

在实际应用中，根据具体条件的不同，导热方程可以简化为不同的形式。

常物性、无内热源情况：

如果材料的热导率 $k$ 、密度 $ρ\rho$ 和比热容 $c_p$ 都是常数，且没有内热源（ $q˙=0\dot{q}=0$ ），则方程简化为：

$∂T∂t=α(∂2T∂x2+∂2T∂y2+∂2T∂z2)\frac{\partial T}{\partial t} = \alpha \left(\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2}\right)$

其中， $α=kρcp\alpha = \frac{k}{\rho c_p}$ 称为热扩散系数（Thermal Diffusivity），单位为 $m^2/s$ 。热扩散系数反映了材料传递热量的能力， $α\alpha$ 越大，热量在材料中传播得越快。

稳态导热：

在稳态条件下，温度不随时间变化， $∂T∂t=0\frac{\partial T}{\partial t} = 0$ ，方程简化为泊松方程（有内热源）或拉普拉斯方程（无内热源）：

$∇2T+q˙k=0或∇2T=0\nabla^2 T + \frac{\dot{q}}{k} = 0 \quad \text{或} \quad \nabla^2 T = 0$

一维非稳态导热：

对于一维问题（如大平壁、长圆柱、球体等），方程简化为：

$∂T∂t=α∂2T∂x2\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}$

这是最简单也是最重要的导热方程形式，许多实际问题都可以用一维模型来近似。

二、一维非稳态导热问题的数值求解

2.1 问题描述

考虑一个半无限大固体，初始时刻温度均匀分布为 $T_0$ 。在 $t = 0$ 时刻，固体表面（ $x = 0$ 处）的温度突然升高到 $T_s$ 并保持不变。求在不同时刻、不同位置的温度分布。

这个问题的数学描述为：

控制方程： $∂T∂t=α∂2T∂x2,x>0,t>0\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}, \quad x > 0, \quad t > 0$

初始条件： $T_0, \quad x > 0$

边界条件： $T_s, \quad T(\infty, t) = T_0, \quad t > 0$

2.2 解析解

该问题存在解析解，可以通过相似变换法或拉普拉斯变换法求得。解析解的形式为：

$T(x,t)−T0Ts−T0=erfc(x2αt)\frac{T(x,t) - T_0}{T_s - T_0} = \text{erfc}\left(\frac{x}{2\sqrt{\alpha t}}\right)$

其中， $erfc(η)=1−erf(η)=2π∫η∞e−ξ2dξ\text{erfc}(\eta) = 1 - \text{erf}(\eta) = \frac{2}{\sqrt{\pi}}\int_{\eta}^{\infty}e^{-\xi^2}d\xi$ 是余误差函数（Complementary Error Function）。

解析解表明，温度分布只与无量纲变量 $η=x2αt\eta = \frac{x}{2\sqrt{\alpha t}}$ 有关，这种特性称为相似性。这意味着温度剖面的形状在不同时间是相似的，只是沿 $x$ 方向按比例缩放。

2.3 有限差分法数值求解

虽然该问题有解析解，但对于更复杂的问题（如变物性、复杂几何形状、非线性边界条件等），通常需要采用数值方法求解。有限差分法是求解导热方程最常用的数值方法之一。

空间离散：将求解区域 $[0, L]$ 划分为 $N$ 个等间距网格，网格间距 $Δx=L/N\Delta x = L/N$ 。节点位置为 $xi=iΔxx_i = i\Delta x$ ， $i = 0, 1, 2, ..., N$ 。

时间离散：将时间区间 $0, t_{max}]$ 划分为 $M$ 个时间步，时间步长 $Δt\Delta t$ 。时间层记为 $tn=nΔtt^n = n\Delta t$ ， $n = 0, 1, 2, ..., M$ 。

显式格式：用向前差分近似时间导数，用中心差分近似空间二阶导数：

$Tin+1−TinΔt=αTi+1n−2Tin+Ti−1nΔx2\frac{T_i^{n+1} - T_i^n}{\Delta t} = \alpha \frac{T_{i+1}^n - 2T_i^n + T_{i-1}^n}{\Delta x^2}$

整理得：

$Tin+1=Tin+αΔtΔx2(Ti+1n−2Tin+Ti−1n)T_i^{n+1} = T_i^n + \frac{\alpha \Delta t}{\Delta x^2}(T_{i+1}^n - 2T_i^n + T_{i-1}^n)$

令 $\frac{\alpha \Delta t}{\Delta x^2}$ ，则：

$T_i^{n+1} = rT_{i-1}^n + (1-2r)T_i^n + rT_{i+1}^n$

显式格式的优点是计算简单，每个时间步只需要直接计算，不需要求解方程组。但显式格式有条件稳定性，要求 $\leq 0.5$ ，即 $Δt≤Δx22α\Delta t \leq \frac{\Delta x^2}{2\alpha}$ 。如果时间步长过大，数值解会出现振荡甚至发散。

隐式格式：用向后差分近似时间导数：

$Tin+1−TinΔt=αTi+1n+1−2Tin+1+Ti−1n+1Δx2\frac{T_i^{n+1} - T_i^n}{\Delta t} = \alpha \frac{T_{i+1}^{n+1} - 2T_i^{n+1} + T_{i-1}^{n+1}}{\Delta x^2}$

整理得：

$rT_{i-1}^{n+1} + (1+2r)T_i^{n+1} - rT_{i+1}^{n+1} = T_i^n$

隐式格式在每个时间步需要求解一个三对角线性方程组，但它是无条件稳定的，可以采用较大的时间步长。

三、Python仿真实现

3.1 仿真代码

# -*- coding: utf-8 -*-
"""
主题003：导热微分方程的推导
一维非稳态导热数值仿真
"""
import matplotlib
matplotlib.use('Agg')  # 非交互式后端
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc
from scipy.linalg import solve_banded
import warnings
warnings.filterwarnings('ignore')
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\传热学仿真\主题003'
os.makedirs(output_dir, exist_ok=True)

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

# 物理参数设置
L = 0.1  # 计算域长度，单位：m
alpha = 1e-5  # 热扩散系数，单位：m^2/s
T0 = 25  # 初始温度，单位：摄氏度
Ts = 100  # 表面温度，单位：摄氏度

# 数值参数设置
N = 100  # 空间网格数
dx = L / N  # 空间步长
x = np.linspace(0, L, N + 1)  # 网格点坐标

# 时间参数
dt_explicit = 0.5 * dx**2 / alpha  # 显式格式最大稳定时间步长
dt = 0.8 * dt_explicit  # 采用稳定的时间步长
nt = 500  # 时间步数

print(f"空间步长 dx = {dx*1000:.4f} mm")
print(f"时间步长 dt = {dt:.6f} s")
print(f"稳定性参数 r = {alpha*dt/dx**2:.4f}")
print("正在计算显式格式...")

# 显式格式求解
T_explicit = np.ones(N + 1) * T0  # 初始温度分布
T_explicit[0] = Ts  # 边界条件

# 存储不同时刻的温度分布
T_history_explicit = []
time_points = [0, 100, 200, 300, 400, 500]

for n in range(nt + 1):
    if n in time_points:
        T_history_explicit.append(T_explicit.copy())
    
    if n < nt:
        T_new = T_explicit.copy()
        for i in range(1, N):
            T_new[i] = T_explicit[i] + alpha * dt / dx**2 * \
                       (T_explicit[i+1] - 2*T_explicit[i] + T_explicit[i-1])
        T_new[0] = Ts  # 应用边界条件
        T_new[N] = T_new[N-1]  # 绝热边界
        T_explicit = T_new

print("显式格式计算完成")
print("正在计算解析解...")

# 计算解析解
t_values = [n * dt for n in time_points]
T_analytical_history = []
for t in t_values:
    if t > 0:
        T_ana = T0 + (Ts - T0) * erfc(x / (2 * np.sqrt(alpha * t)))
    else:
        T_ana = np.ones(N + 1) * T0
        T_ana[0] = Ts
    T_analytical_history.append(T_ana)

print("解析解计算完成")
print("正在生成可视化结果...")

# 创建综合可视化图
fig = plt.figure(figsize=(18, 12))

# 图1：不同时刻温度分布（显式格式）
ax1 = fig.add_subplot(2, 3, 1)
colors = plt.cm.jet(np.linspace(0, 1, len(time_points)))
for i, (T_hist, t_val) in enumerate(zip(T_history_explicit, t_values)):
    ax1.plot(x*100, T_hist, color=colors[i], linewidth=2, 
             label=f't = {t_val:.3f}s')
ax1.set_xlabel('Position x (cm)', fontsize=11)
ax1.set_ylabel('Temperature (C)', fontsize=11)
ax1.set_title('Explicit FDM: Temperature Evolution', fontsize=12, fontweight='bold')
ax1.legend(fontsize=9, loc='upper right')
ax1.grid(True, alpha=0.3)

# 图2：解析解温度分布
ax2 = fig.add_subplot(2, 3, 2)
for i, (T_ana, t_val) in enumerate(zip(T_analytical_history, t_values)):
    ax2.plot(x*100, T_ana, color=colors[i], linewidth=2, 
             label=f't = {t_val:.3f}s')
ax2.set_xlabel('Position x (cm)', fontsize=11)
ax2.set_ylabel('Temperature (C)', fontsize=11)
ax2.set_title('Analytical Solution: Temperature Evolution', fontsize=12, fontweight='bold')
ax2.legend(fontsize=9, loc='upper right')
ax2.grid(True, alpha=0.3)

# 图3：数值解与解析解对比（最终时刻）
ax3 = fig.add_subplot(2, 3, 3)
ax3.plot(x*100, T_analytical_history[-1], 'b-', linewidth=2.5, label='Analytical')
ax3.plot(x*100, T_history_explicit[-1], 'ro', markersize=4, label='Explicit FDM')
ax3.set_xlabel('Position x (cm)', fontsize=11)
ax3.set_ylabel('Temperature (C)', fontsize=11)
ax3.set_title('Numerical vs Analytical (Final Time)', fontsize=12, fontweight='bold')
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 图4：误差分布
ax4 = fig.add_subplot(2, 3, 4)
error = np.abs(T_history_explicit[-1] - T_analytical_history[-1])
ax4.semilogy(x*100, error + 1e-10, 'g-', linewidth=2)
ax4.fill_between(x*100, 1e-10, error, alpha=0.3, color='green')
ax4.set_xlabel('Position x (cm)', fontsize=11)
ax4.set_ylabel('Absolute Error (log scale)', fontsize=11)
ax4.set_title('Numerical Error Distribution', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)

# 图5：温度场时空演化热图
ax5 = fig.add_subplot(2, 3, 5)
T_matrix = np.array(T_history_explicit)
im = ax5.imshow(T_matrix, aspect='auto', cmap='hot', 
                extent=[0, L*100, t_values[-1], 0])
ax5.set_xlabel('Position x (cm)', fontsize=11)
ax5.set_ylabel('Time t (s)', fontsize=11)
ax5.set_title('Temperature Field Evolution', fontsize=12, fontweight='bold')
cbar = plt.colorbar(im, ax=ax5)
cbar.set_label('Temperature (C)', fontsize=10)

# 图6：热流密度分布
ax6 = fig.add_subplot(2, 3, 6)
# 计算热流密度 q = -k * dT/dx
k = 50  # 假设热导率 W/(m·K)
q_history = []
for T_hist in T_history_explicit:
    q = -k * np.gradient(T_hist, dx)
    q_history.append(q)

for i, (q, t_val) in enumerate(zip(q_history[::2], t_values[::2])):
    ax6.plot(x*100, q/1000, linewidth=2, label=f't = {t_val:.3f}s')
ax6.set_xlabel('Position x (cm)', fontsize=11)
ax6.set_ylabel('Heat Flux (kW/m^2)', fontsize=11)
ax6.set_title('Heat Flux Distribution', fontsize=12, fontweight='bold')
ax6.legend(fontsize=9)
ax6.grid(True, alpha=0.3)

plt.suptitle('1D Transient Heat Conduction Simulation', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig(f'{output_dir}/transient_conduction_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 综合仿真图已保存")

# 收敛性分析
print("正在进行网格收敛性分析...")
grid_sizes = [20, 40, 80, 160, 320]
errors_convergence = []

for N_test in grid_sizes:
    dx_test = L / N_test
    x_test = np.linspace(0, L, N_test + 1)
    dt_test = 0.4 * dx_test**2 / alpha
    nt_test = int(t_values[-1] / dt_test)
    
    # 显式格式计算
    T_num = np.ones(N_test + 1) * T0
    T_num[0] = Ts
    
    for n in range(nt_test):
        T_new = T_num.copy()
        for i in range(1, N_test):
            T_new[i] = T_num[i] + alpha * dt_test / dx_test**2 * \
                       (T_num[i+1] - 2*T_num[i] + T_num[i-1])
        T_new[0] = Ts
        T_num = T_new
    
    # 解析解
    T_ana = T0 + (Ts - T0) * erfc(x_test / (2 * np.sqrt(alpha * t_values[-1])))
    
    # 计算误差
    max_error = np.max(np.abs(T_num - T_ana))
    errors_convergence.append(max_error)

# 绘制收敛曲线
fig, ax = plt.subplots(figsize=(10, 6))
dx_values = [L/N for N in grid_sizes]
ax.loglog(dx_values, errors_convergence, 'bo-', linewidth=2, markersize=8, label='Numerical Error')
ax.loglog(dx_values, [errors_convergence[0] * (dx_values[0]/dx)**2 for dx in dx_values], 
          'r--', linewidth=2, label='O(dx^2) Reference')
ax.set_xlabel('Grid Spacing dx (m)', fontsize=12)
ax.set_ylabel('Maximum Absolute Error (C)', fontsize=12)
ax.set_title('Grid Convergence Analysis', fontsize=13, fontweight='bold')
ax.legend(fontsize=11)
ax.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/convergence_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 收敛性分析图已保存")

# 稳定性分析
print("正在进行稳定性分析...")
r_values = [0.1, 0.3, 0.5, 0.51, 0.6]
fig, axes = plt.subplots(1, len(r_values), figsize=(15, 4))

for idx, r in enumerate(r_values):
    dt_stab = r * dx**2 / alpha
    nt_stab = int(0.5 / dt_stab)
    
    T_stab = np.ones(N + 1) * T0
    T_stab[0] = Ts
    
    for n in range(min(nt_stab, 100)):
        T_new = T_stab.copy()
        for i in range(1, N):
            T_new[i] = T_stab[i] + r * (T_stab[i+1] - 2*T_stab[i] + T_stab[i-1])
        T_new[0] = Ts
        T_stab = T_new
    
    axes[idx].plot(x*100, T_stab, linewidth=2)
    axes[idx].set_title(f'r = {r}', fontsize=11)
    axes[idx].set_xlabel('x (cm)', fontsize=10)
    axes[idx].set_ylabel('T (C)', fontsize=10)
    axes[idx].grid(True, alpha=0.3)
    axes[idx].set_ylim([0, 150])

plt.suptitle('Stability Analysis: Effect of r Value', fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig(f'{output_dir}/stability_analysis.png', dpi=150, bbox_inches='tight')
plt.close()
print("✓ 稳定性分析图已保存")

print("\n所有仿真完成！")
print(f"结果保存在: {output_dir}")