如何用scipy solve_bvp求解带连续条件的双层一维稳态导热系统

问题描述 #

需要求解的一维稳态扩散导热方程（补充负号后）：
k T'' = k (d²/dz²) T = -H
其中：

• k (>0 W/m·K)：热导率
• H (≥0 W/m³)：体积加热率

系统分为两层：

• 0 ≤ z ≤ z₁：k=k₁，H=H₁
• z₁ < z ≤ z₂：k=k₂，H=H₂（k₂≥k₁）

边界条件：

• T(z=0)=Tᵤ，T(z=z₂)=T_b
• z=z₁处需满足温度连续和热通量连续（k₁ dT₁/dz(z=z₁) = k₂ dT₂/dz(z=z₁)）

使用scipy.integrate.solve_bvp数值求解时，出现以下问题：

• 当k₁=k₂时，数值解与解析解一致
• 当k₁≠k₂时，数值解错误，核心原因是solve_bvp默认强制状态变量（温度T和温度梯度dT/dr）连续，但实际需要的是T和热通量k dT/dz连续。

原始代码 #

import numpy as np
import scipy
from scipy.integrate import solve_bvp

def H_func(x, z_1, z_2, H_1, H_2, R):
    #x is radius 
    R_1 = R - z_1
    R_2 = R - z_2
    if type(x) == np.ndarray:
        H_arr = np.array([])
        for i in np.arange(len(x)):
            x_val = x[i]
            if (x_val > R_1 and x_val <= R):
                H_val = H_1
            elif (x_val >= R_2 and x_val <= R_1):
                H_val = H_2
            else:
                1/0
            H_arr = np.append(H_arr, H_val)
        return H_arr
    elif type(x) == np.float64:
        if (x > R_1 and x <= R):
            H_val = H_1
        elif (x >= R_2 and x <= R_1):
            H_val = H_2
        else:
            1/0
        return H_val
    else:
        1/0 #failure
        return

def k_func(x, z_1, z_2, k_1, k_2, R):
    R_1 = R - z_1
    R_2 = R - z_2
    if type(x) == np.ndarray:
        k_arr = np.array([])
        for j in np.arange(len(x)):
            x_val = x[j]
            if (x_val > R_1 and x_val <= R):
                 k_val = k_1
            elif (x_val >= R_2 and x_val <= R_1):
                 k_val = k_2
            else:
                 1/0
            k_arr = np.append(k_arr, k_val)
        return k_arr
    elif type(x) == np.float64:
        if (x > R_1 and x <= R):
            k_val = k_1
        elif (x >= R_2 and x <= R_1):
            k_val = k_2
        else:
            1/0
        return k_val
    else:
        1/0
        return 

def fun(x, y, z_1, z_2, k_1, k_2, H_1, H_2, R):

    second_deriv_arr = np.array([])
    for i in np.arange(len(x)):
        x_Val = x[i]
        z_val = R - x_Val  # 修正原代码中的R_p为R
        T_val = y[0][i]
        dTdr_val = y[1][i]
        dTdz_val = -1.*dTdr_val
        if (0. <= z_val and z_val <= z_1):
            second_deriv_val = -1.*H_func(x_Val, z_1, z_2, H_1, H_2, R)/k_func(x_Val, z_1, z_2, k_1, k_2, R)
        elif (z_1 < z_val and z_val <= z_2):
            second_deriv_val = -1.*H_func(x_Val, z_1, z_2, H_1, H_2, R)/k_func(x_Val, z_1, z_2, k_1, k_2, R)
        else:
            print("Probing incorrect parameter space.")
            1/0            

        second_deriv_arr = np.append(second_deriv_arr, second_deriv_val)

    return np.vstack((y[1], second_deriv_arr))

def bc(ya, yb):
    return np.array([ya[0] - Tu_global, yb[0] - Tb_global])

def get_XY_arrays(z_1, z_2, T_u, T_b, k_1, k_2, H_1, H_2, R, spacing1 = 25, spacing2 = 50, tolerance = 1.e-5):

    R_1 = R - z_1
    R_2 = R - z_2

    #an array of radii
    x_in = np.linspace(R_2, R, spacing1)  # 修正原代码的z_2为R_2
    #an array of temperatures and temperature gradients
    y_in = np.zeros((2, x_in.size))

    y_temp = np.linspace(T_b, T_u, spacing1)
    y_in[0] = y_temp
    y_in[1] = -(T_b - T_u)/z_2

    global Tb_global, Tu_global
    Tb_global = T_b
    Tu_global = T_u
    
    res_func = solve_bvp(lambda x,y: fun(x, y, z_1, z_2, k_1, k_2, H_1, H_2, R), bc, x_in, y_in, tol = tolerance)

    x_arr = np.linspace(R_2, R, spacing2)
    y_arr = res_func.sol(x_arr)[0]

    return x_arr, y_arr, res_func

注：原代码中R_p应为R，x_in的起始点应为R_2，已在代码中修正。

解决方案 #

方法1：修改状态变量适配solve_bvp的连续条件 #

核心思路是将状态变量从[T, dT/dr]替换为[T, q]，其中q是热通量（q = -k dT/dz），这样solve_bvp默认的状态变量连续条件正好匹配物理上的温度连续和热通量连续要求。

推导r为自变量的ODE：

已知z = R - r，因此dr = -dz，d/dr = -d/dz。

• 热通量定义：q = -k dT/dz = -k*(-dT/dr) = k dT/dr → 温度梯度dT/dr = q/k
• 原方程k T'' = -H转换为热通量的导数：d/dz(k dT/dz) = H → -dq/dz = H → dq/dr = -H（因为dq/dr = dq/dz * dz/dr = dq/dz*(-1)）

最终，r为自变量的ODE系统为：

dT/dr = y[1] / k(r)
dq/dr = -H(r)

其中y = [T, q]，边界条件：

• r=R（对应z=0）：T = Tᵤ
• r=R-z₂（对应z=z₂）：T = T_b

修改后的代码

import numpy as np
from scipy.integrate import solve_bvp

def H_func(x, z_1, z_2, H_1, H_2, R):
    R_1 = R - z_1
    R_2 = R - z_2
    if isinstance(x, np.ndarray):
        H_arr = np.where((x > R_1) & (x <= R), H_1, H_2)
        if np.any((x < R_2) | (x > R)):
            raise ValueError("x out of valid range")
        return H_arr
    elif isinstance(x, (float, np.float64)):
        if R_1 < x <= R:
            return H_1
        elif R_2 <= x <= R_1:
            return H_2
        else:
            raise ValueError("x out of valid range")
    else:
        raise TypeError("x must be float or numpy array")

def k_func(x, z_1, z_2, k_1, k_2, R):
    R_1 = R - z_1
    R_2 = R - z_2
    if isinstance(x, np.ndarray):
        k_arr = np.where((x > R_1) & (x <= R), k_1, k_2)
        if np.any((x < R_2) | (x > R)):
            raise ValueError("x out of valid range")
        return k_arr
    elif isinstance(x, (float, np.float64)):
        if R_1 < x <= R:
            return k_1
        elif R_2 <= x <= R_1:
            return k_2
        else:
            raise ValueError("x out of valid range")
    else:
        raise TypeError("x must be float or numpy array")

def fun(x, y, z_1, z_2, k_1, k_2, H_1, H_2, R):
    T = y[0]
    q = y[1]
    k_vals = k_func(x, z_1, z_2, k_1, k_2, R)
    H_vals = H_func(x, z_1, z_2, H_1, H_2, R)
    # ODE系统
    dTdr = q / k_vals
    dqdr = -H_vals
    return np.vstack((dTdr, dqdr))

def bc(ya, yb, Tu, Tb):
    # ya: r=R处的状态（z=0），yb: r=R-z2处的状态（z=z2）
    return np.array([ya[0] - Tu, yb[0] - Tb])

def get_XY_arrays(z_1, z_2, T_u, T_b, k_1, k_2, H_1, H_2, R, spacing1 = 25, spacing2 = 50, tolerance = 1.e-5):
    R_1 = R - z_1
    R_2 = R - z_2

    # 初始猜测点
    x_in = np.linspace(R_2, R, spacing1)
    # 状态向量初始猜测：T从Tb到Tu线性变化，q初始设为常数（无热源时的热通量）
    y_in = np.zeros((2, x_in.size))
    y_in[0] = np.linspace(T_b, T_u, spacing1)
    # 无热源时q=(T_u - T_b)*k_avg/z2，这里用k1做初始猜测
    q_init = k_1 * (T_u - T_b) / z_2
    y_in[1] = np.full(x_in.size, q_init)
    
    # 求解BVP
    res_func = solve_bvp(
        lambda x,y: fun(x, y, z_1, z_2, k_1, k_2, H_1, H_2, R),
        lambda ya,yb: bc(ya, yb, T_u, T_b),
        x_in, y_in, tol=tolerance
    )

    # 生成结果数组
    x_arr = np.linspace(R_2, R, spacing2)
    y_arr = res_func.sol(x_arr)[0]
    q_arr = res_func.sol(x_arr)[1]

    return x_arr, y_arr, q_arr, res_func

代码改进点：

• 用np.where替代循环，提升效率
• 去掉全局变量，将边界条件参数传入bc函数
• 状态变量改为[T, q]，适配物理连续条件
• 初始猜测更合理，热通量设为常数（符合无热源场景的物理意义）

方法2：分段求解+匹配连续条件 #

对于分段常数的介质，可以分别求解两层的ODE，再利用连续条件和边界条件联立解出积分常数：

1. 第一层（0≤z≤z₁）：解k₁ T₁'' = -H₁，得到含两个常数的解T₁(z)
2. 第二层（z₁<z≤z₂）：解k₂ T₂'' = -H₂，得到含两个常数的解T₂(z)
3. 利用条件：
   • T₁(0)=Tᵤ
   • T₂(z₂)=T_b
   • T₁(z₁)=T₂(z₁)
   • k₁ T₁’(z₁)=k₂ T₂’(z₁)
4. 联立4个方程解出4个常数，得到解析形式的解（适合当前简单场景）

该方法可扩展到更复杂的分段系统，但代码量较大，且后续含平流项的复杂场景不易处理。

方法3：有限差分/有限元法 #

对于后续含平流项的复杂系统，可采用：

• 有限差分法（FDM）：手动离散方程，构建线性方程组求解
• 有限元法（FEM）：使用FEniCS、PyVista等库，天然支持复杂介质和边界条件

验证建议 #

1. 无热源测试：设置H₁=H₂=0，此时热通量q为常数，温度分布为两段线性，数值解应与解析解完全一致
2. 有热源测试：对比解析解（可通过分段求解推导），验证温度分布和热通量的正确性
3. 无量纲化验证：保持无量纲化处理，可减少数值误差，提升求解稳定性

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