使用solve_ivp求解钢棒多维热传导问题的异常求助

你的代码存在几个核心错误，导致温度求解和绘图出现问题，逐一说明并修正：

1. 能量方程的核心错误：温度变化率计算缺失质量项 #

能量守恒方程中，能量变化率 = 质量 × 比热容 × 温度变化率，你直接将能量变化率赋值给温度变化率，完全漏掉了质量和比热容的除法，导致温度变化的数值完全错误。

每个分段的质量为：mass = rho * 体积，体积是径向环形面积×轴向长度，即A * z，因此温度变化率应为：

dTdt_array[i, j] = dEi_dt / (mass * Cp)

2. 轴向/径向分段长度的变量混淆 #

你将径向分段数n和轴向分段数m搞反了：

• 径向从中心到边缘（总半径D/2）分为n段，因此径向每段增量应为 r = (D/2)/n
• 轴向总长度L分为m段，因此轴向每段长度应为 z = L/m

原代码中这两个变量的计算完全颠倒，导致所有面积、换热的几何参数错误。

3. 求解结果的维度reshape错误 #

solve_ivp返回的sol.y形状为**(状态变量数, 时间点数)**，你的初始状态是(m,n)的矩阵，flatten后有m*n个状态变量，因此需要先转置再reshape，才能得到(时间点数, 轴向分段数, 径向分段数)的结构：

temperatures = sol.y.T.reshape((len(t), m, n))

原代码的reshape顺序完全错误，导致时间、轴向、径向的维度混乱，绘图时自然无法正确对应。

4. 径向侧面积的计算误差 #

原代码用径向分段的外半径计算侧面积，更准确的做法是用分段的平均半径，即该段内半径和外半径的平均值：

avg_radius = (j*r + (j+1)*r) / 2
S = 2 * np.pi * avg_radius * z

修正后的完整代码 #

import matplotlib.pyplot as plt 
import numpy as np 
from scipy.integrate import solve_ivp

# Parameters
n = 5  # 径向分段数
m = 3  # 轴向分段数
D = 0.1  # 钢棒直径 (m) 
Cp = 500  # 比热容 (J/kg·K) 
rho = 8000  # 密度 (kg/m³) 
k = 20  # 热导率 (W/m·K) 
L = 0.5  # 钢棒长度 (m) 
r = (D/2)/n  # 径向每段半径增量 (m)
z = L/m      # 轴向每段长度 (m)
U = 18  # 侧面总换热系数 (W/m²·K) 
U_top = 18  # 顶部总换热系数 (W/m²·K) 
T_amb = 290  # 环境温度 (K) 
T_plate = 500  # 加热板温度 (K) 

def Energy_eqns(t, y):
    y = y.reshape((m, n))  # 1D转2D：(轴向分段, 径向分段)
    dTdt_array = np.zeros((m, n))

    for i in range(m):  # 遍历轴向分段
        for j in range(n):  # 遍历该轴向分段内的径向分段
            # 径向环形面积（轴向端面面积）
            A = np.pi * (((j + 1) * r) ** 2 - (j * r) ** 2)
            # 分段平均半径
            avg_radius = (j*r + (j+1)*r) / 2
            # 径向侧面积
            S = 2 * np.pi * avg_radius * z

            # 轴向流入热流
            if i == 0:
                # 底部与加热板接触，热阻为z/2（分段中心到加热板的距离）
                Qi_1 = (k * A * (T_plate - y[i, j])) / (z / 2)
            else:
                Qi_1 = (k * A * (y[i - 1, j] - y[i, j])) / z

            # 轴向流出热流
            if i < m - 1:
                Qi = k * A * (y[i, j] - y[i + 1, j]) / z
            else:
                # 顶部与环境换热
                Qi = U_top * A * (y[i, j] - T_amb)

            # 径向流入热流
            if j > 0:
                Qrj_1 = k * S * (y[i, j - 1] - y[i, j]) / r
            else:
                # 中心无径向热流
                Qrj_1 = 0

            # 径向流出热流
            if j < n - 1:
                Qrj = k * S * (y[i, j] - y[i, j + 1]) / r
            else:
                # 外侧与环境换热
                Qrj = U * S * (y[i, j] - T_amb)

            # 能量变化率
            dEi_dt = Qi_1 + Qrj_1 - Qi - Qrj
            # 计算温度变化率：dT/dt = dEi_dt / (质量*Cp)
            mass = rho * A * z
            dTdt_array[i, j] = dEi_dt / (mass * Cp)

    return dTdt_array.flatten()

# 初始温度：所有分段为环境温度
initial_temperatures = np.ones((m, n)) * T_amb
# 积分时间范围
t_span = (0, 100)
# 指定输出时间点
t_eval = np.linspace(0, 100, 20)  # 减少绘图数量，避免弹窗过多
# 求解微分方程组
sol = solve_ivp(Energy_eqns, t_span, initial_temperatures.flatten(), method='LSODA', t_eval=t_eval) 

t = sol.t
# 修正reshape：(时间点, 轴向分段, 径向分段)
temperatures = sol.y.T.reshape((len(t), m, n))

# 绘图：选取部分时间点绘制（避免100个弹窗）
plot_indices = [0, 5, 10, 15, 19]  # 初始、25s、50s、75s、100s
vmin = np.min(temperatures)
vmax = np.max(temperatures) + 20  # 扩大颜色范围，增强对比

for idx in plot_indices:
    time = t[idx]
    temp_profile = temperatures[idx]  # (轴向分段, 径向分段)
    
    fig, ax = plt.subplots(figsize=(8, 4))
    # imshow输入为(径向, 轴向)，对应ylabel是径向，xlabel是轴向
    im = ax.imshow(temp_profile.T, cmap='hot', aspect='auto', origin='lower', vmin=vmin, vmax=vmax)
    
    ax.set_xlabel('轴向分段 (从底部到顶部)')
    ax.set_ylabel('径向分段 (从中心到边缘)')
    ax.set_title(f't={time:.1f}s 时的温度分布')
    plt.colorbar(im, ax=ax, label='温度 (K)')
    plt.tight_layout()
    plt.show()

额外说明 #

1. 我减少了t_eval的点数并只选取关键时间点绘图，避免弹出100个窗口影响体验，你可以根据需求调整。
2. 修正后的代码会正确计算每个分段的温度变化，热力图的轴向（从底部加热板到顶部）和径向（从中心到边缘）维度完全对应，时间关联正常。

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