How to define the boundary condition for radiative heat transfer on a 1-D geometry?

10 views (last 30 days)
One of the users in MATLAB suggested me to use pdepe to handle the unsteady state heat transfer equation if there is large variation in thermal conductivity. However, I could not apply the boundary condition for radiation heat transfer. At that state, the product is passing through the IR lamp and the top surface is exposed to the IR radiation and the bottom surface is exposed to air so for the top surface, boundary condition could not be figured out. As the radiation source is external, I have added it as a source term but when I put pr=0 and qr = 0 then it will show this error
Error using pdepe
Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.
Error in pde_test (line 93)
sol = pdepe(m, @pdefun, @icfun, @bcfun, y, t13);
and here is the code for that, could you please figure out how to specify the boundary condition here?
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 = 1.0; % 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 = 0;
%pr = h_air*(ur-T_ambient)
qr = 0;
else
pl = -h_air*(ul-T_ambient);
ql = 1;
pr = h_air*(ur-T_ambient);
qr = 1;
end
end
end
  8 Comments
Torsten
Torsten on 29 Feb 2024
Edited: Torsten on 29 Feb 2024
Could you please explain why you have selected qr=1 while value of pr is a constant?
The boundary conditions are set in the form
p + q*f = 0
where f is set in "pdefun".
In your case,
f = layer1_k * dT/dy
Thus setting
p = -Net_radiative_intensity , q = 1
gives
-Net_radiative_intensity + 1*layer1_k*dT/dy = 0
or
layer1_k*dT/dy = Net_radiative_intensity
thus the same boundary condition that you set in your code:
u(Ny, n + 1) = u(Ny - 1, n + 1) + ((Net_radiative_intensity * dy) / layer1_k);
do you know which type of solver pdepe is using because it is nice and kinda robust. I try to read the documentation but couldn't get a clue.
It discretizes the partial differential equation in space and leaves the temporal derivatives in their continuous form. This gives a big system of ordinary differential equations for the temperatures in the grid points. To solve it, a stiff implicit integrator (in this case: ode15s) is used. The advantage over your explicit Euler method is that the method is stable despite large time steps and that the error of the solution is controlled (at least in time direction). Look up "method of lines" for more details.

Sign in to comment.

Answers (0)

Products


Release

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!