📚 课程目录第2章热传导控制方程与定解条件

2. 热传导控制方程与定解条件

数值模拟温度场的第一步，不是打开软件画网格，而是把“温度在空间和时间中如何变化”写成可求解的数学问题。这个数学问题由两部分组成：控制方程描述区域内部温度演化的普遍规律，定解条件则把普遍规律约束到具体问题上——包括初始时刻的温度分布和边界上的换热方式。缺少任何一部分，温度场都不唯一。

一、傅里叶定律：热流与温度梯度的本构关系

傅里叶定律是导热问题的本构方程，它回答“热量往哪里流、流多少”。在一维情形下，热流密度 q（W/m²）与温度梯度成正比：

q_x = -λ · ∂T/∂x

其中 λ 为导热系数（W/(m·K)），负号表示热量沿温度降低的方向传递。三维一般形式写成矢量：q = −λ∇T。这里必须约定符号：本文取“热流密度沿坐标轴正方向为正”，因此当温度沿 x 增大而降低时 ∂T/∂x < 0，q_x > 0，热量向 +x 方向流动，与物理直觉一致。后续推导和边界条件均沿用这一约定。

傅里叶定律只说明“有梯度就有热流”，它本身不保证能量守恒。要得到温度场的控制方程，还需把热流代入能量守恒。

二、能量守恒推导导热微分方程

在导热体内取边长 dx、dy、dz 的微元控制体。按“流入能量 − 流出能量 + 内热源生成 = 内能增量”列衡算。以 x 方向为例，净流入热量为 q_x − q_{x+dx} ≈ −(∂q_x/∂x)dx，三个方向相加得净导热流入量 −∇·q dV。设单位体积内热源强度为 Φ（W/m³），内能增量为 ρc_p(∂T/∂t)dV，其中 ρ 为密度、c_p 为比热容。于是：

ρc_p ∂T/∂t = −∇·q + Φ

将 q = −λ∇T 代入，得到导热微分方程的一般形式：

ρc_p ∂T/∂t = ∇·(λ∇T) + Φ

若材料各向同性且 λ 为常数，可写成 ρc_p ∂T/∂t = λ∇²T + Φ，进一步定义热扩散率 a = λ/(ρc_p)（m²/s），则 ∂T/∂t = a∇²T + Φ/(ρc_p)。a 越大，温度扰动在材料中传播并趋于均匀的速度越快，这是判断瞬态过程快慢的关键参数。稳态无内热源时方程退化为 ∇²T = 0，即拉普拉斯方程。

数值实现中，λ 随温度变化时不要提前把它提出散度算子，而应保留 ∇·(λ∇T) 的守恒形式，否则在物性突变界面附近会引入额外误差。

三、初始条件与三类边界条件

初始条件只对瞬态问题需要，给出 t = 0 时刻整个区域的温度分布：T(x,y,z,0) = T₀(x,y,z)。工程上常见两种简化：全域均匀温度 T₀，或按实测/解析分布给定。初始条件与边界条件必须相容，例如边界温度恒为 100 ℃ 而初始全场为 20 ℃，会在 t = 0⁺ 产生温度跳跃，数值上表现为初始时刻边界附近梯度极大，需要更细的时间步。

边界条件分三类，是数值模拟中最容易设错的部分：

第一类（Dirichlet）：给定边界温度，T|_Γ = T_w(t)。适用于表面与恒温源紧密接触、相变界面等情形。
第二类（Neumann）：给定边界热流密度，−λ∂T/∂n|_Γ = q_w(t)，n 为外法向。q_w = 0 即绝热边界；给定热流加热时 q_w 取负值（热量流入）。
第三类（Robin）：给定表面与流体的对流换热，−λ∂T/∂n|_Γ = h(T_w − T_f)，h 为对流换热系数（W/(m²·K)），T_f 为环境流体温度。它把边界温度与外部流体温度通过热阻联系起来。

注意第三类条件中 h 与 λ 的量纲差异：h 很大时 T_w → T_f，退化为第一类；h → 0 时退化为绝热。这一极限关系可用于检查边界设置是否合理。

四、把方程与条件落到数值模型上

以一个可复现的教学示例说明如何组装定解问题：一块 200 mm × 100 mm 的钢板（λ = 45 W/(m·K)，ρ = 7850 kg/m³，c_p = 470 J/(kg·K)，无内热源），初始温度 20 ℃。左端面与 200 ℃ 热源理想接触，右端面与 25 ℃ 空气对流换热（h = 15 W/(m²·K)），上下两面绝热。则：

控制方程：ρc_p ∂T/∂t = λ(∂²T/∂x² + ∂²T/∂y²)
初始条件：T(x,y,0) = 20 ℃
左边界 x=0：T = 200 ℃                （第一类）
右边界 x=L：−λ∂T/∂x = h(T − 25)      （第三类）
上下边界：∂T/∂y = 0                   （第二类，绝热）

由此可估算热扩散率 a = 45/(7850×470) ≈ 1.22×10⁻⁵ m²/s。若用显式格式，稳定性要求时间步满足 Δt ≤ Δx²/(4a)（二维）。取 Δx = 2 mm，则 Δt ≤ (0.002)²/(4×1.22×10⁻⁵) ≈ 0.082 s，实际取 0.05 s 留出裕量。这一估算直接决定网格与时间步的匹配，是数值模拟前必须完成的参数检查。

边界条件的离散要特别小心：第二类条件用单侧差分或引入虚拟节点实现，第三类条件通常把 h 与 λ 组合成边界节点的等效热阻。若把第三类误设为第一类恒温，会人为提高换热能力，导致右端温度偏高、整体升温偏快；反之把绝热面误设为对流面，会引入本不存在的散热通道。判断结果是否可信，可先看稳态趋势：本例中热量从左向右传递，右端温度应介于 25 ℃ 与 200 ℃ 之间，且沿 x 单调下降，若出现非单调或超出区间，应优先回查边界设置与符号约定。

三类边界条件不可混用同一套离散模板：第一类直接约束节点温度，第二、三类约束的是温度梯度或梯度与温差的组合，需在代数方程组中体现为不同的系数项。修改边界类型后，务必重新检查系数矩阵是否仍保持对角占优。

综上，控制方程给出温度场的演化规律，初始条件锁定时间起点，三类边界条件锁定空间界面的换热方式，三者共同构成适定的定解问题。数值模拟的可靠性，首先取决于这套数学表达是否与物理对象一致，其次才是离散格式与网格质量。