I am developing a model for the 1-D geometry product which is exposed to different conditions on the top and bottom surface as you can see in the code. Unfortunately, for the boundary condition, it is always giving error as bc takes multiple argument which should take single argument.
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]
here i am solving transient heat transfer.I have tried different approach but every time same error is coming up
and the error is
sol = solve_bvp(pdefun, icfun, y, t_all, bc=bc)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ TypeError: solve_bvp() got multiple values for argument 'bc'
Could someone please help or tell how to avoid the multiple argument for bc?
i have a working model in matlab and i am trying to develop the same in python and as the problem type is boundary value problem, i have used solve_bvp module but as one guy mentioned in the comment that there should be only bc not ic, then is there any other method to develop this code maybe using daetools or somehting else?
function pde()
% 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)
% layer3_thickness = 0.002; % Thickness of layer 3 (meters)
% layer3_k = 0.3; % Thermal conductivity of layer 3 (W/m-K)
% 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 = ceil(layer2_thickness/Ly*Ny);
ny2 = Ny-ny1;
y1 = linspace(0,layer2_thickness,ny1);
y2 = linspace(layer2_thickness,Ly,ny2);
y = [y1,y2(2:end)];
% 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 = 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 = ceil(t1/(t1+t2+t3+t4)*Nt);
nt2 = ceil(t2/(t1+t2+t3+t4)*Nt);
nt3 = ceil(t3/(t1+t2+t3+t4)*Nt);
nt4 = Nt-(nt1+nt2+nt3);
t11 = linspace(0,t1,nt1);
t12 = linspace(t1,t1+t2,nt2);
t13 = linspace(t1+t2,t1+t2+t3,nt3);
t14 = linspace(t1+t2+t3,t1+t2+t3+t4,nt4);
%pdepe settings
m = 0; %for 1-D cartesian coordinates with no symmetry
phase = 1;
sol = pdepe(m, @pdefun, @icfun, @bcfun, y, t11);
%sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t11);
u1 = sol(:,:,1);
phase = 2;
sol = pdepe(m, @pdefun, @icfun, @bcfun, y, t12);
%sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t12);
u2 = sol(:,:,1);
phase= 3;
sol = pdepe(m, @pdefun, @icfun, @bcfun, y, t13);
%sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t13);
u3 = sol(:,:,1);
phase=4;
sol = pdepe(m, @pdefun, @icfun, @bcfun, y, t14);
%sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t14);
u4 = sol(:,:,1);
plot([t11,t12,t13,t14],[[u1(:,1);u2(:,1);u3(:,1);u4(:,1)],[u1(:,25);u2(:,25);u3(:,25);u4(:,25)],[u1(:,50);u2(:,50);u3(:,50);u4(:,50)]])
grid on
%assigning the coefficients for pde
function [c f s] = pdefun(y,t,u,dudy)
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 radiaiton
end
end
%defining the initial conditions for pde
function u = icfun(yq)
if phase==1
if yq <= layer2_thickness
u = initial_temperature1;
else
u = initial_temperature2;
end
elseif phase==2
u = interp1(y,u1(end,:),yq);
elseif phase==3
u = interp1(y,u2(end,:),yq);
else
u = interp1(y,u3(end,:),yq);
end
end
%defining the boundary conditions for pde
function [pl ql pr qr] = bcfun(yl,ul,yr,ur,t)%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;
elseif phase==2
pl = -h_air*(ul-T_ambient);
ql = 1;
pr = h_air*(ur-T_ambient);
qr = 1;
elseif phase==3
pl = -h_air*(ul-T_ambient);
ql = 1;
pr = -Net_radiative_intensity;
%pr = h_air*(ur-T_ambient)
qr = 1;
else
pl = -h_air*(ul-T_ambient);
ql = 1;
pr = h_air*(ur-T_ambient);
qr = 1;
end
end
end
According to the documentation, the second argument of solve_bvp is for the boundary condition function. So, you give it icfun as the second argument and then you also give it a keyword argument for bc. You've doubled up on what bc is so solve_bvp is "confused".
I'm not sure what icfun is actually supposed to be used for nor why you're passing it to solve_bvp. Also note that you have other syntax errors in your code. You should review the documentation and verify how the ODE function (why do you have it called PDE?) should be defined and used (y will be a numpy array, not a scalar) as well as how the bc function should be defined and used.