Matlab与Python一维瞬态传热计算结果差异排查求助

问题描述 #

求解y方向一维瞬态传热方程：
dT/dt = (k/(rho*cp))*d²T/dy²

• 初始条件：t=0时所有单元温度为25℃
• 边界条件：
  • 底部：热辊对流边界 k*dT/dy = h(T_roller - T)
  • 上部：大气对流边界 k*dT/dy = h(T - T_ambient)

辐射边界下两种语言求解结果完全一致，但热辊边界场景下存在显著数值差异：辊温较低时差异4-5℃，辊温升高100℃时差异扩大至10-11℃；已确认时间节点设置一致，曲线趋势匹配。

核心排查方向 #

1. 边界条件离散格式差异 #

这是最可能的原因，重点检查热辊边界的离散实现：

• 确认两侧对dT/dy的近似方式一致：比如底部边界节点（设为T₁）与相邻内部节点（T₂）的温差近似，是否都采用(T₂-T₁)/Δy（而非单侧差分的其他形式）
• 代入边界方程后的推导结果是否一致：
  把k*dT/dy = h(T_roller - T)替换为离散形式后，整理得到的T₁更新公式应为：
  T₁ = (k*T₂/Δy + h*T_roller) / (k/Δy + h)
  检查两边代码中是否存在符号错误、系数计算错误（比如把T_roller - T写成T - T_roller）

2. 数值求解器参数不匹配 #

• 时间积分方法是否一致：比如Matlab用ode15s（隐式刚性求解器），Python是否用了scipy.integrate.solve_ivp的'Radau'或'BDF'求解器？显式/隐式方法的选择会直接影响结果精度
• 收敛阈值（rtol/atol）、时间步长限制是否完全相同：隐式求解器的迭代终止条件差异会导致数值偏差
• 空间网格划分是否一致：节点数量、网格间距Δy的取值必须完全对应，不能出现Matlab用100个节点、Python用99个节点的情况

3. 物理参数与索引细节 #

• 物理参数的精度与单位：确认k、rho、cp、h等参数在两种语言中用的是相同精度（比如都是双精度），单位统一（比如不要混淆W/m·K和W/cm·K）
• 数组索引差异：Matlab数组索引从1开始，Python从0开始，底部/顶部边界节点的索引是否对应正确（比如Matlab的T(1)对应Python的T[0]），避免边界条件应用到错误节点

4. 代码逻辑细节 #

• 变量更新顺序：检查是否存在Matlab原地更新数组、Python用副本更新的差异，导致时间步迭代中变量依赖关系错误
• 边界条件的应用时机：确认每一步时间迭代中都正确应用了边界条件，没有遗漏或重复计算

验证步骤 #

1. 简化测试用例：用2-3个空间节点、极小时间步长，手动计算前2步的结果，对比Matlab和Python的输出，定位差异出现的第一步
2. 输出中间变量：在两边代码中添加打印语句，输出每一步的边界节点温度、内部节点温度、时间步长，直接对比数值差异的来源
3. 统一求解器：手动实现相同的时间积分方法（比如显式欧拉法），替换掉各自的内置求解器，看差异是否消失

内容的提问来源于stack exchange，提问作者Adya