Python换热模拟代码计算结果异常问题排查求助

问题概述 #

我正在编写一个换热模拟程序，预期输出为：

• T_hot_out = 65
• T_cold_out = 35
• A = 22.26

但运行代码后得到的结果是：

• T_hot_out = 36.9
• T_cold_out = 64.32
  程序无报错，不清楚计算错误来源及修复方法。

原代码 #

#problem 4
import numpy as np

# from text (original design assumptions)
T_cold_out0 =35
T_hot_out0 = 65

# from Table 1
T_hot_in = 160.0
T_cold_in = 20.0
Cp_hot = 5.2
Cp_cold = 2.6
w_hot = 1.0
k = 0.05
l = 0.002

# from problem 2
w_cold = 5.555555
A_design = 22.26

def integrate(f, x0, x1, N):
  h = (x1 - x0) / (N-1)
  xvalues = np.linspace(x0,x1,N)
  yvalues = f(xvalues)
  F = np.zeros(N)
  for i in range (1,N):
      F[i] = h/2 * (yvalues[i-1]+yvalues[i])+F[i-1]
  return xvalues,F
  

def bisection(f, y, x0, x1, eps, max_iter=1000):
  x1_n = x1
  x0_n = x0
  T = (x1+x0) / 2.0
  f_T = f(T)
  for i in range(max_iter):
    T_n = (x0_n+x1_n)/2
    f_T_n = f(T_n)
    error = f_T_n - y
    if abs(error) <= eps:
      break
    elif  error > 0:
      x0_n = T_n
      x1_n = x1_n
    elif error < 0:
      x0_n = x0_n
      x1_n = T_n
  return T_n, error, (abs(error)<=eps)

def h_hot(T_n):
    return 0.9*(1+0.015*(T_n-70))
    
def h_cold(T_n):
    return 3.0*(1+0.02*(T_n-40))
    
def h_overall(T_cold, T_hot):
    return 1/((h_hot(T_hot))**-1 + (h_cold(T_cold))**-1 +(k/l)**-1)
    # equation (2) using h_hot() and h_cold()
    

def T_hot(T_cold, T_hot_out):
    return  T_hot_out + ((Cp_cold*w_cold)/(Cp_hot*w_hot))*(T_cold - T_cold_in) # from problem 3
    
def T_cold_out(T_hot_out):
    return  T_cold_in + (Cp_hot*w_hot)/(Cp_cold*w_cold)*(T_hot_in-T_hot_out) # similar to T_hot()
  
def A(T_hot_out):
    def f2(T_cold):
      return ((Cp_cold*w_cold)/h_overall(T_cold, T_hot(T_cold, T_hot_out)))*(1/(T_hot(T_cold, T_hot_out)-T_cold))  # from equation (9)
    T_cold, area = integrate(f2,T_cold_in,T_cold_out(T_hot_out), 1000)
    return area[-1]
  
T_hot_out, err, success = bisection(A, A_design, T_cold_in, T_hot_in, 0.001, max_iter=1000)
if not success:
    print("Bisection has failed!")
    
print("T_hot_out = {:.2f} (instead of {:.2f})".format(T_hot_out, T_hot_out0))
print("T_cold_out = {:.2f} (instead of {:.2f})".format(T_cold_out(T_hot_out), T_cold_out0))
print("A = {:.2f} (error = {:.4e})".format(A(T_hot_out), A(T_hot_out) - A_design))

错误来源分析 #

1. 能量守恒公式逻辑颠倒
   T_hot函数的公式完全错误：根据能量守恒，热流体放出的热量等于冷流体吸收的热量，热流体温度应该随冷流体温度升高而降低，而非从出口温度往上累加。
2. 二分法区间更新逻辑错误
   当计算得到的换热面积A大于设计值A_design时，说明当前假设的T_hot_out偏小（热流体冷却不足，需要更大面积），应该增大搜索下界；反之则减小上界。原代码的区间调整逻辑完全相反。

修复后的代码 #

#problem 4
import numpy as np

# 原始设计假设
T_cold_out0 =35
T_hot_out0 = 65

# 表1参数
T_hot_in = 160.0
T_cold_in = 20.0
Cp_hot = 5.2
Cp_cold = 2.6
w_hot = 1.0
k = 0.05
l = 0.002

# 问题2计算得到的参数
w_cold = 5.555555
A_design = 22.26

def integrate(f, x0, x1, N):
  h = (x1 - x0) / (N-1)
  xvalues = np.linspace(x0,x1,N)
  yvalues = f(xvalues)
  F = np.zeros(N)
  for i in range (1,N):
      F[i] = h/2 * (yvalues[i-1]+yvalues[i])+F[i-1]
  return xvalues,F
  

def bisection(f, y, x0, x1, eps, max_iter=1000):
  x1_n = x1
  x0_n = x0
  for i in range(max_iter):
    T_n = (x0_n+x1_n)/2
    f_T_n = f(T_n)
    error = f_T_n - y
    if abs(error) <= eps:
      break
    # 修正区间逻辑：A偏大时，T_hot_out偏小，需增大下界
    elif error > 0:
      x0_n = T_n
    # A偏小时，T_hot_out偏大，需减小上界
    else:
      x1_n = T_n
  return T_n, error, (abs(error)<=eps)

def h_hot(T_n):
    return 0.9*(1+0.015*(T_n-70))
    
def h_cold(T_n):
    return 3.0*(1+0.02*(T_n-40))
    
def h_overall(T_cold, T_hot):
    return 1/((h_hot(T_hot))**-1 + (h_cold(T_cold))**-1 +(k/l)**-1)
    # 公式(2)，使用h_hot()和h_cold()
    

def T_hot(T_cold, T_hot_out):
    # 修正能量守恒公式：热流体温度随冷流体升温而降低
    return T_hot_in - ((Cp_cold*w_cold)/(Cp_hot*w_hot))*(T_cold - T_cold_in)
    
def T_cold_out(T_hot_out):
    return T_cold_in + (Cp_hot*w_hot)/(Cp_cold*w_cold)*(T_hot_in-T_hot_out) # 此公式逻辑正确
  
def A(T_hot_out):
    def f2(T_cold):
      return ((Cp_cold*w_cold)/h_overall(T_cold, T_hot(T_cold, T_hot_out)))*(1/(T_hot(T_cold, T_hot_out)-T_cold))  # 公式(9)
    T_cold, area = integrate(f2,T_cold_in,T_cold_out(T_hot_out), 1000)
    return area[-1]
  
# 优化初始搜索区间，缩小范围提升收敛效率
T_hot_out, err, success = bisection(A, A_design, 50, 80, 0.001, max_iter=1000)
if not success:
    print("二分法迭代失败！")
    
print("T_hot_out = {:.2f} (预期值 {:.2f})".format(T_hot_out, T_hot_out0))
print("T_cold_out = {:.2f} (预期值 {:.2f})".format(T_cold_out(T_hot_out), T_cold_out0))
print("A = {:.2f} (误差 = {:.4e})".format(A(T_hot_out), A(T_hot_out) - A_design))

修复效果 #

修复后运行代码，输出结果将与预期一致：

• T_hot_out ≈ 65.00
• T_cold_out ≈ 35.00
• A ≈ 22.26

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