implicit numerical method: 1-D unsteady state heat transfer using finite difference method

% 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;

%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;

%Assigning initial values to the temperature matrix
for i=1:Ny
    
    if y(i) <= layer2_thickness
        u(i,1) = initial_temperature1;
        
    else
        u(i,1) = initial_temperature2;
        
    end
    
end


distance_travelled = zeros(1,Nt);   %initial distance travelled
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 heating roller

p0 = y_1;
p1 = p0+y_2;

% 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) <= layer2_thickness
                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
            
            
            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) + (layer2_k * u(Ny - 1, n + 1))) / (layer2_k + h_roller * dy);
        u(1, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(2, n + 1))) / (layer1_k + h_air * dy);
        
        
    else
        for i = 2:Ny - 1
            % Determine the thermal conductivity based on the layer
            if y(i) <= layer2_thickness
                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
            
            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) + (layer2_k * u(2, n + 1))) / (layer2_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) >= p9
        
        break;
    end
    
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time-stepping loop (crank nicholson 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) <= layer2_thickness
                k = layer2_k;
                rho = rho_2;
                cp = cp_2;
                alpha = k/(rho*cp);
                r = alpha*dt/(2*dy^2);
                
            else
                k = layer1_k;
                rho = rho_1;
                cp = cp_1;
                alpha = k/(rho*cp);
                r = alpha*dt/(2*dy^2);
            end
            
            
            u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n) + u(i + 1, n + 1) - 2 * u(i, n + 1) + u(i - 1, n + 1));
        end
        
        
        % Apply convective boundary condition
        u(Ny, n + 1) = ((h_roller * T_roller * dy) + (layer2_k * u(Ny - 1, n + 1))) / (layer2_k + h_roller * dy);
        u(1, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(2, n + 1))) / (layer1_k + h_air * dy);
        
        
    else
        for i = 2:Ny - 1
            % Determine the thermal conductivity based on the layer
            if y(i) <= layer2_thickness
                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
            
            u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n) + u(i + 1, n + 1) - 2 * u(i, n + 1) + u(i - 1, n + 1));
        end
        
        
        
        % Apply convective boundary condition
        u(1, n + 1) = ((h_air * T_ambient * dy) + (layer2_k * u(2, n + 1))) / (layer2_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) >= p1
        
        break;
    end
    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Time-stepping loop (backward euler 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) <= layer2_thickness
                k = layer2_k;
                rho = rho_2;
                cp = cp_2;
                %alpha = k/(rho*cp);
                %r = alpha*dt/(2*dy^2);
                
            else
                k = layer1_k;
                rho = rho_1;
                cp = cp_1;
                %alpha = k/(rho*cp);
                %r = alpha*dt/(dy^2);
            end
            
            % Calculate coefficient matrix for implicit equation
            alpha = k / (rho * cp);
            r = alpha * dt / dy^2;
            A = diag(1 + 2 * r * ones(Ny-2, 1)) + diag(-r * ones(Ny-3, 1), 1) + diag(-r * ones(Ny-3, 1), -1);
            
            %u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n) + u(i + 1, n + 1) - 2 * u(i, n + 1) + u(i - 1, n + 1));
            %end
            
            
            % Apply convective boundary condition
            u(Ny, n + 1) = ((h_roller * T_roller * dy) + (layer2_k * u(Ny - 1, n + 1))) / (layer2_k + h_roller * dy);
            u(1, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(2, n + 1))) / (layer1_k + h_air * dy);
            
            % Right-hand side of the implicit equation
            b = u(2:end-1, n) + r * (u(3:end, n) - 2 * u(2:end-1, n) + u(1:end-2, n));
            
            % Solve the implicit equation
            u(2:end-1, n + 1) = A \ b;
        end
        
        
    else
        for i = 2:Ny - 1
            % Determine the thermal conductivity based on the layer
            if y(i) <= layer2_thickness
                k = layer2_k;
                rho = rho_2;
                cp = cp_2;
                %alpha = k/(rho*cp);
                %r = alpha*dt/(2*dy^2);
                
            else
                k = layer1_k;
                rho = rho_1;
                cp = cp_1;
                %alpha = k/(rho*cp);
                %r = alpha*dt/(dy^2);
            end
            
            % Calculate coefficient matrix for implicit equation
            alpha = k / (rho * cp);
            r = alpha * dt / dy^2;
            A = diag(1 + 2 * r * ones(Ny-2, 1)) + diag(-r * ones(Ny-3, 1), 1) + diag(-r * ones(Ny-3, 1), -1);
            
            %u(i, n + 1) = u(i, n) + r * (u(i + 1, n) - 2 * u(i, n) + u(i - 1, n) + u(i + 1, n + 1) - 2 * u(i, n + 1) + u(i - 1, n + 1));
            %end
            
            
            % Apply convective boundary condition
            u(Ny, n + 1) = ((h_roller * T_roller * dy) + (layer2_k * u(Ny - 1, n + 1))) / (layer2_k + h_roller * dy);
            u(1, n + 1) = ((h_air * T_ambient * dy) + (layer1_k * u(2, n + 1))) / (layer1_k + h_air * dy);
            
            % Right-hand side of the implicit equation
            b = u(2:end-1, n) + r * (u(3:end, n) - 2 * u(2:end-1, n) + u(1:end-2, n));
            
            % Solve the implicit equation
            u(2:end-1, n + 1) = A \ b;
        end
    end
    
    % Break the loop when the total distance is reached
    if distance_travelled(n + 1) >= p1
        
        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(50, 1:last_valid_time_step), ...
    'r', 'DisplayName', 'Oberseite');
hold on
plot(t(1:last_valid_time_step), u(25, 1:last_valid_time_step), ...
    'g', 'DisplayName', 'Interface');
xlabel('Time (seconds)');
ylabel('Temperature (°C)');
yticks(20:10:200);
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 200]);
hold off

I want to apply implicit method to the 1-D unsteady state heat transfer problem to diminsh the effect of large thermal conductivity or very small densities or specific heat capacities. I was writing the code according to the classic example given in book and youtube videos but I am unable to match upto the solution provided by Explicit numerical method. Also there is no example of using crank nicholson or backward euler using convective boundary conditions or conductive boundary conditions.

Could you explain why is this so?

[Torsten]

You didn't supply your equations in a mathematical notation, but you should test whether it's possible to use MATLAB's "pdepe" to solve.

[Joydeep Bagchi]

@Torsten - thankyou for the reply.

ok i will try pdepe. apart from trying pdepe, is there any other way to apply the implicit method?

[Torsten]

Further, you have a two-layer problem. Where do you discretize the transmission condition

T1 = T2
lambda1*dT1/dx = lambda2*dT2/dx

at the interface ?

I can only advise you to use "pdepe" (at least at the beginning) to have a valid reference solution.

The point at the interface should be part of the mesh for this solver.

I'm confused about the thermal diffusivity settings in your code. First you write

if y(i) <= layer2_thickness
    k = layer2_k;  
    ...
else
    k = layer1_k;
    ...
end

but when you incorporate the boundary conditions, you work with layer1_k at y(1) and layer2_k at y(end).

This doesn't make sense in my opinion.

[Joydeep Bagchi]

i do not need to define the interface temperature because i considered the product as a whole and then dividing the properties accoridng to the different thickness

if you consider a paper and across its thickness of length say 1mm so 0.5 mm is of maybe paper 1 type and 0.5 mm may be of paper 2 type so the properties in those 2 layers will be different

[Torsten]

if you consider a paper and across its thickness of length say 1mm so 0.5 mm is of maybe paper 1 type and 0.5 mm may be of paper 2 type so the properties in those 2 layers will be different

I know. That's why you have to establish that temperature and heat flux at the interface point are continuous:

T1 = T2
lambda1*dT1/dx = lambda2*dT2/dx

where T1 and T2 are the limiting temperatures from the left and from the right, respectively.

Writing T1 = T2 and lambda1*dT1/dx = lambda2*dT2/dx does not mean that you impose temperature at the interface, but that you need these two conditions to compute the "correct" temperature that satisfies these two conditions.

[Joydeep Bagchi]

ok so I have to define for the interface and if i change the number of gridpoints everytime then I have to define it each time

for example u(Ny/2,:)=u(Ny+1/2,:)

then similarly the gradient also, right?

[Torsten]

I usually mesh from a to the interface point and from the interface point to b and concatenate the grids (like in the pdepe code below).

I don't understand what you mean by

u(Ny/2,:)=u(Ny+1/2,:)

[Joydeep Bagchi]

I got your point. Sorry for the misconception above.

Answer 1

[Torsten]

pde()

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;
%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 heating roller
p0 = y_1;
p1 = p0+y_2;
t1 = p0/velocity_x;
t2 = p1/velocity_x-t1;
nt1 = ceil(t1/(t1+t2)*Nt);
nt2 = Nt-nt1;
t11 = linspace(0,t1,nt1);
t12 = linspace(t1,t1+t2,nt2);
%pdepe settings
m = 0;

phase = 1;
sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t11);
u1 = sol(:,:,1);

phase = 2;
sol = pdepe(m,@pdefun,@icfun,@bcfun,y,t12);
u2 = sol(:,:,1);

plot([t11,t12],[[u1(:,1);u2(:,1)],[u1(:,25);u2(:,25)],[u1(:,50);u2(:,50)]])
grid on

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 = 0;
  end
end

function u = icfun(yq)
  if phase==1
    if yq <= layer2_thickness
      u = initial_temperature1;
    else
      u = initial_temperature2;
    end
  else
    u = interp1(y,u1(end,:),yq);
  end
end

function [pl ql pr qr] = bcfun(yl,ul,yr,ur,t)
  if phase==1
    pl = -h_air*(ul-T_ambient);
    ql = 1;
    pr = h_roller*(ur-T_roller);
    qr = 1;
  else
    pl = -h_air*(ul-T_ambient);
    ql = 1;
    pr = h_air*(ur-T_ambient);
    qr = 1;
  end
end

end

[Joydeep Bagchi]

Thankyou @Torsten for the input.

I will apply the full problem statement to this now as it was only the part of the code, let's see if it can handle the big change in thermal conductivity value accurately or not

Here, the condition at the interface is not defined as you have mentioned earlier, right?

[Torsten]

I will apply the full problem statement to this now as it was only the part of the code, let's see if it can handle the big change in thermal conductivity value accurately or not

Note that pdepe should only be used for parabolic-elliptic pdes. Thus if your problem statement has equations without a second-order term, you should be careful.

Here, the condition at the interface is not defined as you have mentioned earlier, right?

According to the documentation of "pdepe", continuity of temperature and heat flux are established by the "pdepe" discretization if the interface point is part of the xmesh.

[Joydeep Bagchi]

ok the heat transfer equation is itself parabolic pde

all conditions will remain same, just the boundary condition will change. I will update here about the result.

[Joydeep Bagchi]

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

Hi @Torsten, I am sorry to bother you again.

So now the problem lies in the boundary condition for the phase 3 where an IR field is applied. As far as I have understood, this is the external source term applied on the top surface for layer_1.

Have I defined the source term wrongly or defined the boundary condition for phase=3 wrongly because if I run the program then it is giving error that spatial discretization is not possible.

Also, by the literature I found out that the external source term should be in the units of W/m³ so I have tried dividing the net_radiative_intensity by layer1_thickness but it is also giving the same error.

Could you please tell me how to resolve that?

[Torsten]

Have I defined the source term wrongly or defined the boundary condition for phase=3 wrongly because if I run the program then it is giving error that spatial discretization is not possible.

At the moment, you don't define a boundary condition at yr for phase 3:

pr = 0;
%pr = h_air*(ur-T_ambient)
qr = 0;

If radiation is present, isn't it a comination of convective and radiative heat transfer at the surface ?

lambda*dT/dy = h*(T-Tambient) + e*(T^4-Tambient^4)

or something of that kind ? Are you sure a source term should be set within the layers ?

[Joydeep Bagchi]

yes, on the bottom surface it is convective boundary condition so it is defined in pl and ql is 1 for neumann. I have read somewhere that in case of external source, both the value of the boundary condition and derivative of the boundary condition do not play any role so i defined it as 0.

for the boltzmann law equation i think you are talking about the radiaiton loss but i do not want to measure the loss here. I want to measure the temperature distribution in the upper layer when the IR field radiation kicks in.

[Torsten]

Most probably you read about a model that is not discretized in spatial direction, but determines a single temperature only. Here, the energy input is directly set in the equation. You have a partial differential equation where usually (apart from chemical reaction heat) energy input is set across the boundaries.

Anyhow: Mathematically, a partial differential equation of second order needs two boundary conditions. And you have set only one.

[Joydeep Bagchi]

ok, thankyou for the information. But I only have the data of intensity of IR lamp which is 138 kW, radiation length and width, no other data is available. What else could I set up?

[Torsten]

As I said: I'd compute the radiation heat flux at the upper surface, add it to the convective heat flux and remove the source term. But I'm not really a process engineer.

See the chapter "Heat Transfer between surfaces" under

https://en.wikipedia.org/wiki/Thermal_radiation

how the radiative heat flux between surfaces is usually modelled.

Answer 2

[William Rose]

@Joydeep Bagchi,

I think you need to divide the long term (u(i+1,n)-2*u(i,n)+...) by 2, on the right side of your Crank-Nicolson update equation. But still, I think it is wrong, because Crank-Nicolson is an implicit method which one solves with a triadiagonal matrix. It is just a different tridiagonal matrix than the one in the backward Euler implicit method. I have not evvaluatred the handling of the boundary conditions.

This is a complicated model. I would test the 3 algorithms on a simpler model first, to see if they give the same results - or close-enough-to-the-same.

How do you know the algoritms give different results? The Crank-Nicolson results appear to overwrite the explicit method results. Then the backward Euler appears to overwrite the Crank-Nicolson. Am I missing something?

You can compute r1 and r2 (for layers 1 and 2), and matrices A1 (using r1) and A2 (using r2) at the start of your code, before all the loops.

r1 = layer1_k*dt/(rho_1*cp_1*2*dy^2); % don't need alpha
r2 = layer2_k*dt/(rho_2*cp_2*2*dy^2);
% equation for A is different for backward Euler and C-N
A1be=[triagonal matix using r1];  
A2be=[triagonal matix using r2];
A1cn=[triagonal matix using r1];  
A2cn=[triagonal matix using r2];

Then, when looping over over position, you do

for i = 2:Ny - 1
    % Determine the thermal conductivity based on the layer
    if y(i)<=layer2_thickness
        r=r2; A=A2be;   % or A2cn 
    else 
        r=r1; A=A1be;   % or A1cn
    end
    u(i,n+1) = u(i,n) + r*(stuff);
end

and so on. This will make your code shorter and simpler.

Explicit method update equation (with some spaces removed)

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

Your Crank-Nicholson update equation

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

Yor Implicit method (backward Euler) update equation

% Coefficient matrix for implicit equation
A = diag(1+2*r*ones(Ny-2,1)) + diag(-r*ones(Ny-3,1),1) + diag(-r*ones(Ny-3,1),-1);
% Right-hand side of the implicit equation
b = u(2:end-1,n) + r * (u(3:end,n)-2*u(2:end-1,n)+u(1:end-2,n));
u(2:end-1,n+1) = A\b;  % solve

[Joydeep Bagchi]

Thankyou @William Rose for the reply.

I am running three different matlab files so the constants are same at the beginning, just the time stepping loop is different. all three methods should give about same results and implicit methods should be more robust and unconditionally stable. because with explicit method, i am getting the solution but it heavily depends on parameter 'r' and it depends on density,thermal conductivity and specific heat capacity. So if I change any one of the parameters by a large magnitude then my solution differs a lot which should not be the case.