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

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

[Torsten]

could you please figure out how to specify the boundary condition here?

As said, this is not a question for a MATLAB forum. Maybe you should ask this in a physics forum.

Perhaps

layer1_k*dT/dy = 138000 / upper surface area  [W/m^2]

on the upper surface ? But this radiative heat flux looks much too strong. So: You must find how much of these 138 000 W is netto energy input to the upper surface of your object, and divide this amount by the upper surface area. But the surface temperature will rise during the process so that less energy will be absorbed. Here, the difference of the temperatures to the power of 4 again comes into play.

[Joydeep Bagchi]

Thankyou @Torsten

Matlab do not allow me to post comment on the answer so I have to write other question. I do not know why it is not allowing. I have values like absorbance and performance as 0.45 and 0.9 but still after that the effect will be too much.

You told me to ask in the physicsforum, could you please mention which one you are referring to?

also as you referred, i tried in the code

% Constants
T = 40;             % Total simulation time (seconds)
Ny = 50;            % Number of spatial grid points
Nt = 5000;           % Number of time steps
velocity_x = 0.2; % line speed in m/s

% Thermal properties of different layers
layer1_thickness = 550e-6;  % Thickness of layer 1 (meters)
layer1_k = 0.15;            % Thermal conductivity of layer 1 (W/m-K)
layer2_thickness = 70e-6;  % Thickness of layer 2 (meters)
layer2_k = 0.16;            % Thermal conductivity of layer 2 (W/m-K)
layer3_thickness = 3000e-6;  % Thickness of layer 3 (meters)
layer3_k = 0.35;            % Thermal conductivity of layer 3 (W/m-K)
rho_1 = 1250;           %density of layer 1 (kg/m³)
rho_2 = 1400;           %density of layer 2 (kg/m³)
rho_3 = 1800;           %density of layer 3(kg/m3)
cp_1 = 1350;            %specific heat capacity of layer 1 (J/kg-K)
cp_2 = 1100;            %specific heat capacity of layer 2 (J/kg-K)
cp_3 = 1800;            %specific heat capacity of layer 3 (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)

x0  = layer3_thickness;
x1 = layer2_thickness;
x2 = layer1_thickness;
x01 = layer3_thickness+layer2_thickness;



Ly = layer1_thickness+layer2_thickness+layer3_thickness;         % Total thickness of the system (meters)

% 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 = 200e3; %W
radiative_flux = 50e3; %W/m2
absorption = 45; % in %
performance = 90; % in %
emissivity = 0.90;

%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
dy = Ly / (Ny - 1);
dt = T / Nt;
y = linspace(0, Ly, Ny);
t = linspace(0, T, Nt);


% Initialize temperature matrix
u = zeros(Ny, Nt);
initial_temperature1 = 25;
initial_temperature2 = 25;
initial_temperature3 = 25;

%Assigning initial values to the temperature matrix
for i=1:Ny

    if y(i) <= x0
        u(i,1) = initial_temperature3;

    elseif y(i)>x0 && y(i) <= x01
        u(i,1) = initial_temperature2;

    else
        u(i,1) = initial_temperature1;
 
    end

end


distance_travelled = zeros(1,Nt);   %initial distance travelled
dia_roller = 1.0;                   %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 heating roller
y_3 = 1.28;                                              %in contact with IR
y_4 = 0.3;                                              % after IR
y_5 = pi*dia_cooling_roller*contact_angle_cool1/360;    %in contact with cooling roller1
y_6 = 0.5;                                              %after cooling roller 1
y_7 = pi*dia_cooling_roller*contact_angle_cool2/360;    %in contact with cooling roller2
y_8 = 0.5;                                              %after cooling roller 2
y_9 = pi*dia_cooling_roller*contact_angle_cool3/360;    %in contact with cooling roller3
y_10 = 1.0;                                              %after cooling roller 3

p0 = y_1;
p1 = p0+y_2;
p2 = p1 + y_3;
p3 = p2 + y_4;
p4 = p3 + y_5;
p5 = p4 + y_6;
p6 = p5 + y_7;
p7 = p6 + y_8;
p8 = p7 + y_9;
p9 = p8 + y_10;

% Forward Euler method
%r = alpha * dt / (dy^2);
%it is better to keep the parameter 'r' between 0.2 - 0.3 for better stability
%purpose  and if >0.5 then it is highly unstable

% Time-stepping loop (explicit method)
for n = 1:Nt - 1

    % Update position in x-direction
    distance_travelled(n+1) = distance_travelled(n) + velocity_x * dt;

    % Check if the total distance is within the contact length
    if distance_travelled(n) <= p0

    for i = 2:Ny - 1
        % Determine the thermal conductivity based on the layer
        if y(i) <= x0
            k = layer3_k;
            rho = rho_3;
            cp = cp_3;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);

        elseif y(i)>x0 && y(i)<= x01
            k = layer2_k;
            rho = rho_2;
            cp = cp_2;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);


        else
            k = layer1_k;
            rho = rho_1;
            cp = cp_1;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);
        end
        alpha = k/(rho*cp);
        r = alpha*dt/(dy^2);                

        u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
    end
    

    % Apply convective boundary condition
    u(Ny, n + 1) = ((h_roller * T_roller * dy) + (layer1_k * u(Ny - 1, n + 1))) / (layer1_k + h_roller * dy);
    u(1, n + 1) = ((h_air * T_ambient * dy) + (layer3_k * u(2, n + 1))) / (layer3_k + h_air * dy);


    elseif distance_travelled(n) <= p1
        for i = 2:Ny - 1
        % Determine the thermal conductivity based on the layer
        if y(i) <= x0
            k = layer3_k;
            rho = rho_3;
            cp = cp_3;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);

        elseif y(i)>x0 && y(i)<= x01
            k = layer2_k;
            rho = rho_2;
            cp = cp_2;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);


        else
            k = layer1_k;
            rho = rho_1;
            cp = cp_1;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);
        end
        alpha = k/(rho*cp);
        r = alpha*dt/(dy^2);                

        u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
    end



    % Apply convective boundary condition
    u(1, n + 1) = ((h_air * T_ambient * dy) + (layer3_k * u(2, n + 1))) / (layer3_k + h_air * dy);
    u(Ny, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(Ny - 1, n + 1))) / (layer1_k + h_air * dy);

    elseif distance_travelled(n) <= p2
        for i = 2:Ny - 1
        % Determine the thermal conductivity based on the layer
        if y(i) <= x0
            k = layer3_k;
            rho = rho_3;
            cp = cp_3;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);

        elseif y(i)>x0 && y(i)<= x01
            k = layer2_k;
            rho = rho_2;
            cp = cp_2;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);


        else
            k = layer1_k;
            rho = rho_1;
            cp = cp_1;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);
        end
        alpha = k/(rho*cp);
        r = alpha*dt/(dy^2);                

        u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
    end


    % Apply convective boundary condition
    u(1, n + 1) = ((h_air * T_ambient * dy) + (layer3_k * u(2, n + 1))) / (layer3_k + h_air * dy);
    u(Ny, n + 1) = u(Ny - 1, n + 1) + ((Net_radiative_intensity * dy) / layer1_k);
    
 else
        for i = 2:Ny - 1
        % Determine the thermal conductivity based on the layer
        if y(i) <= x0
            k = layer3_k;
            rho = rho_3;
            cp = cp_3;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);

        elseif y(i)>x0 && y(i)<= x01
            k = layer2_k;
            rho = rho_2;
            cp = cp_2;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);


        else
            k = layer1_k;
            rho = rho_1;
            cp = cp_1;
            %alpha = k/(rho*cp);
            %r = alpha*dt/(dy^2);
        end
        alpha = k/(rho*cp);
        r = alpha*dt/(dy^2);                

        u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n));
    end

    % Apply convective boundary condition
    u(1, n + 1) = ((h_air * T_ambient * dy) + (layer3_k * u(2, n + 1))) / (layer3_k + h_air * dy);
    u(Ny, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(Ny - 1, n + 1))) / (layer1_k + h_air * dy);


    end

    % Break the loop when the total distance is reached
    if distance_travelled(n + 1) >= p3
       
        break;
    end

end

% Plot the temperature distribution over time
% Find the last valid time step where u(:, Nt) is non-zero
last_valid_time_step = find(u(1, :) ~= 0, 1, 'last');
figure();
plot(t(1:last_valid_time_step), u(1, 1:last_valid_time_step), ...
    'b', 'DisplayName', 'Unterseite');
hold on
plot(t(1:last_valid_time_step), u((Ny), 1:last_valid_time_step), ...
   'r', 'DisplayName', 'Oberseite');
hold on
plot(t(1:last_valid_time_step), u(((Ny)/2), 1:last_valid_time_step), ...
   'g', 'DisplayName', 'Interface');
xlabel('Time (seconds)');
ylabel('Temperature (°C)');
yticks(20:10:250);
xticks(0 :1 :t(1,last_valid_time_step));
legend('Location', 'Best');
title('Temperature Distribution in composite');
grid on;
axis([0 t(1,last_valid_time_step) 20 250]);
hold off

you can see that there is sharp increase in the temperature which is not the actual case.

The actual case somewhat looks like this

you can see the second peak which is somewhere around 190

[Torsten]

You told me to ask in the physicsforum, could you please mention which one you are referring to?

I've not yet used a physics forum. I just googled "physics forum" and found quite a few. Maybe you can search in advance if one of them is specialized on modelling and simulation.

You set the boundary condition

layer1_k*dT/dy = Net_radiative_intensity

This boundary condition can also be set in the "pdepe" version.

I thought you had problems dealing with a volumetric heat source term from the radiation.

[Joydeep Bagchi]

I've not yet used a physics forum. I just googled "physics forum" and found quite a few. Maybe you can search in advance if one of them is specialized on modelling and simulation.

Ok, I will look for it.

This boundary condition can also be set in the "pdepe" version.

Could you please tell how, because I have tried by putting pr=Net_radiative_intensity and qr=0 then it is giving error.

I thought you had problems dealing with a volumetric heat source term from the radiation.

Since the last said boundary condition is not working for me, so I go by definition and then try to put the radiation as an external source.

[Torsten]

Could you please tell how, because I have tried by putting pr=Net_radiative_intensity and qr=0 then it is giving error.

pr = -Net_radiative_intensity
qr = 1

[Joydeep Bagchi]

Thankyou @Torsten for the reply.

Could you please explain why you have selected qr=1 while value of pr is a constant?

also, 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.

[Torsten]

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.

[Joydeep Bagchi]

Thankyou @Torsten for the insight and the explanation of the boundary condition.