如何使用fsolve求解theta2？多参数非线性方程求解问题

核心问题梳理 #

你需要求解热平衡方程的根θ₂，对应7种绝缘厚度、4种θ₁值、2种发射率，共7×4×2=56组参数组合。基于热传递原理，完整的热平衡方程应为：
绝缘层热传导速率 = 外表面对流换热速率 + 外表面与环境辐射换热速率

对应的数学表达式（假设为单位长度管道，温度需转换为绝对温标°R，辐射计算要求绝对温度）：

Q_cond = (2πk(θ₁ᴿ - θ₂ᴿ)) / ln(d₂/d₁)
Q_conv = hc·πd₂(θ₂ᴿ - θ₃ᴿ)
Q_rad = σπd₂·(1/(1/εₘ + 1/εc - 1))·(θ₂ᴿ⁴ - θ₃ᴿ⁴)
目标函数：f(θ₂) = Q_cond - Q_conv - Q_rad = 0

其中：θᴿ = θ°F + 459.67（华氏度转兰金温标）

fsolve错误原因分析 #

你遇到的Result from function call is not a proper array of floats错误，本质是：

• 直接用sympy符号变量θ₂构建函数，fsolve无法处理符号对象，要求输入函数必须返回纯浮点数（标量或数组）。
• 未将符号表达式转换为可计算数值的函数，导致返回值非数值类型。

正确实现步骤与代码示例 #

1. 导入库与定义基础参数 #

import math
from scipy.optimize import fsolve

# 给定参数
k = 0.5
d1 = 20/12  # 单位：ft
Linsulation = [2/12, 3/12, 4/12, 5/12, 6/12, 7/12, 8/12]  # 绝缘厚度，单位：ft
em_list = [0.09, 0.9]
sigma = 0.171e-8  # 单位：Btu/(hr·ft²·°R⁴)
theta1_list = [800, 900, 1000, 1100]  # 单位：°F
theta3_F = 70  # 环境温度，单位：°F

# 转换为绝对温标
theta3_R = theta3_F + 459.67

# 计算外直径d2数组
d2 = [d1 + 2*L for L in Linsulation]

2. 定义纯数值目标函数 #

目标函数需接收数值类型的θ₂ᴿ，返回浮点数形式的目标函数值：

def heat_balance(theta2_R, theta1_R, d2_val, em_val, ec_val=0.9):
    # 计算传导速率（单位长度）
    Q_cond = (2 * math.pi * k * (theta1_R - theta2_R)) / math.log(d2_val / d1)
    
    # 计算对流换热系数hc
    delta_T_conv = theta2_R - theta3_R
    hc = 0.270 * (delta_T_conv ** 0.25) * (d2_val ** -0.25)
    
    # 计算对流速率
    Q_conv = hc * math.pi * d2_val * delta_T_conv
    
    # 计算辐射等效发射率与辐射速率
    eps_eff = 1 / (1/em_val + 1/ec_val - 1)
    Q_rad = sigma * math.pi * d2_val * eps_eff * (theta2_R**4 - theta3_R**4)
    
    # 返回热平衡差值（目标为0）
    return Q_cond - Q_conv - Q_rad

3. 批量求解所有参数组合 #

通过嵌套循环遍历所有参数组合，调用fsolve求解，同时给出合理的初始猜测值：

# 存储所有结果
results = []

for theta1_F in theta1_list:
    theta1_R = theta1_F + 459.67
    for d2_val in d2:
        insulation_thickness = (d2_val - d1)/2  # 还原绝缘厚度
        for em_val in em_list:
            # 初始猜测值：取θ₁与θ₃的中间绝对温度
            initial_guess = (theta1_R + theta3_R) / 2
            # 调用fsolve，通过args传递固定参数
            theta2_R_sol, = fsolve(heat_balance, initial_guess, args=(theta1_R, d2_val, em_val))
            # 转换回华氏度
            theta2_F_sol = theta2_R_sol - 459.67
            # 记录结果
            results.append({
                "theta1(°F)": theta1_F,
                "绝缘厚度(in)": round(insulation_thickness*12, 1),
                "发射率em": em_val,
                "theta2(°F)": round(theta2_F_sol, 2)
            })

# 打印前5组结果示例
for idx, res in enumerate(results[:5], 1):
    print(f"第{idx}组: θ₁={res['theta1(°F)']}, 绝缘厚度={res['绝缘厚度(in)']}in, em={res['发射率em']}, θ₂={res['theta2(°F)']}°F")

关键注意事项 #

• 绝对温标：辐射计算必须用绝对温度，否则结果完全错误。
• 初始猜测值：fsolve依赖合理初始值，建议取θ₁与θ₃的中间值，避免收敛到不合理的根。
• 纯数值化函数：确保函数内所有运算为浮点数操作，不混入sympy符号变量，这是解决错误的核心。
• 参数传递：通过fsolve的args参数传递固定参数，避免硬编码，提升代码复用性。

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