我正在开发一维几何产品的模型，该模型的顶部和底部表面暴露在不同的条件下，正如您在代码中看到的那样。不幸的是，对于边界条件，它总是给出错误，因为 bc 采用多个参数，而本应采用单个参数。
import numpy as np
from scipy.integrate import solve_bvp
# Constants
Ly = 0.0015 # Total thickness of the system (meters)
T = 40 # Total simulation time (seconds)
Ny = 51 # Number of spatial grid points
Nt = 50000 # Number of time steps
velocity_x = 0.2 # line speed in m/s
# Thermal properties of different layers
layer1_thickness = 0.001 # Thickness of layer 1 (meters)
layer1_k = 0.15 # Thermal conductivity of layer 1 (W/m-K)
layer2_thickness = 0.5e-3 # Thickness of layer 2 (meters)
layer2_k = 0.16 # Thermal conductivity of layer 2 (W/m-K)
rho_1 = 1250 # density of layer 1 (kg/m³)
rho_2 = 1520 # density of layer 2 (kg/m³)
cp_1 = 1350 # specific heat capacity of layer 1 (J/kg-K)
cp_2 = 1460 # specific heat capacity of layer 2 (J/kg-K)
T_roller = 150 # Temperature of hot roller (°C)
h_roller = 500 # Convective heat transfer coefficient (W/m2-K)
T_cooling_roller = 22 # Temperature of cooling temperature(°C)
# Ambient temperature and convection properties
T_ambient = 25 # Ambient temperature (°C)
h_air = 12 # Convective heat transfer coefficient of air (W/m²-K)
radiation_length = 1.280 # radiation length in m
radiation_width = 2.300 # radiation width in m
radiation_power = 138e3 # W
radiative_flux = 50e3 # W/m2
absorption = 45 # in %
performance = 95 # in %
emissivity = 0.91 # emissivity values range from 0.90 to 0.93
# if radiative flux is given
# Net_radiative_intensity = radiative_flux*(absorption/100)*(performance/100);
# if radiation power is given
Net_radiative_intensity = (radiation_power / (radiation_length * radiation_width)) * (absorption / 100) * (
performance / 100) * emissivity
# Discretization in space
ny1 = np.ceil(layer2_thickness / Ly * Ny).astype(int)
ny2 = Ny - ny1
y1 = np.linspace(0, layer2_thickness, ny1)
y2 = np.linspace(layer2_thickness, Ly, ny2)
y = np.concatenate((y1, y2[1:]))
# Initialize temperature matrix
initial_temperature1 = 25
initial_temperature2 = 25
dia_roller = 0.8 # diameter of hot roller (m)
dia_cooling_roller = 0.57 # diameter of cooling roller(m)
contact_angle = 180 # contact angle/wrap angle for hot_roller in °
contact_angle_cool1 = 120 # contact angle/wrap angle for cooling_roller1 in °
contact_angle_cool2 = 220 # contact angle/wrap angle for cooling_roller2 in °
contact_angle_cool3 = 190 # contact angle/wrap angle for cooling_roller3 in °
y_1 = np.pi * dia_roller * contact_angle / 360 # in contact with hot roller
y_2 = 0.3 # after hot roller
y_3 = radiation_length # in contact with IR
y_4 = 0.3 # after IR
t1 = y_1 / velocity_x
t2 = y_2 / velocity_x
t3 = y_3 / velocity_x
t4 = y_4 / velocity_x
nt1 = np.ceil(t1 / (t1 + t2 + t3 + t4) * Nt).astype(int)
nt2 = np.ceil(t2 / (t1 + t2 + t3 + t4) * Nt).astype(int)
nt3 = np.ceil(t3 / (t1 + t2 + t3 + t4) * Nt).astype(int)
nt4 = Nt - (nt1 + nt2 + nt3)
t11 = np.linspace(0, t1, nt1)
t12 = np.linspace(t1, t1 + t2, nt2)
t13 = np.linspace(t1 + t2, t1 + t2 + t3, nt3)
t14 = np.linspace(t1 + t2 + t3, t1 + t2 + t3 + t4, nt4)
# pdepe settings
def pdefun(y, t, u, dudy,phase):
if y <= layer2_thickness:
c = rho_2 * cp_2
f = layer2_k * dudy
s = 0
else:
c = rho_1 * cp_1
f = layer1_k * dudy
s = Net_radiative_intensity # source term for IR radiation
return c, f, s
# defining the initial conditions for pde
def icfun(yq,phase):
if phase == 1:
return np.where(yq <= layer2_thickness, initial_temperature1, initial_temperature2)
elif phase == 2:
return np.interp(yq, y, u1[-1, :])
elif phase == 3:
return np.interp(yq, y, u2[-1, :])
else:
return np.interp(yq, y, u3[-1, :])
def get_bc_values(phase):
if phase == 1:
pl = -h_air * (sol.y[0, 0] - T_ambient)
ql = 1
pr = h_roller * (sol.y[0, -1] - T_roller)
qr = 1
elif phase == 2:
pl = -h_air * (sol.y[1, 0] - T_ambient)
ql = 1
pr = h_air * (sol.y[1, -1] - T_ambient)
qr = 1
elif phase == 3:
pl = -h_air * (sol.y[2, 0] - T_ambient)
ql = 1
pr = -Net_radiative_intensity
qr = 1
else:
pl = -h_air * (sol.y[3, 0] - T_ambient)
ql = 1
pr = h_air * (sol.y[3, -1] - T_ambient)
qr = 1
return pl, ql, pr, qr
"""
# defining the boundary conditions for pde
def bcfun(yl, ul, yr, ur, t,phase): # right is for top and left is for bottom
if phase == 1:
pl = -h_air * (ul - T_ambient)
ql = 1
pr = h_roller * (ur - T_roller)
qr = 1
elif phase == 2:
pl = -h_air * (ul - T_ambient)
ql = 1
pr = h_air * (ur - T_ambient)
qr = 1
elif phase == 3:
pl = -h_air * (ul - T_ambient)
ql = 1
pr = -Net_radiative_intensity
qr = 1
else:
pl = -h_air * (ul - T_ambient)
ql = 1
pr = h_air * (ur - T_ambient)
qr = 1
return pl, ql, pr, qr
"""
def bc(ya, yb, phase):
return get_bc_values(phase)
# Concatenate time arrays
t_all = np.concatenate((t11, t12, t13, t14))
# Solve the system using solve_bvp
sol = solve_bvp(pdefun, icfun, y, t_all, bc=bc)
# sol = solve_bvp(pdefun, icfun, y, t_all, bc=bcfun)
# Extract solution
u1 = sol.sol(t11)[0]
u2 = sol.sol(t12)[0]
u3 = sol.sol(t13)[0]
u4 = sol.sol(t14)[0]
这里我正在解决瞬态传热问题。我尝试了不同的方法，但每次都会出现相同的错误
错误是
sol = solve_bvp(pdefun, icfun, y, t_all, bc=bc)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ TypeError: solve_bvp() got multiple values for argument 'bc'
有人可以帮忙或告诉我如何避免 bc 的多重参数吗？ 预先感谢！
根据文档，
solve_bvpicfunbcbcsolve_bvp我不确定
icfunsolve_bvpybc