Why is there discrepancy in the result with matlab and python?

I want to solve the 1-D transient heat transfer equation.

Define a 1-D geometry (a line) in y-direction

dT/dt = (k/(rho*cp))*d²T/dy²

I.C.: @t=0, temperature of all elements is 25°C

B.C.: The geometry is divided into several elements.

1. The bottom element is exposed to hot roller -kdT/dy = h(T_roller-T)

2. The upper element is exposed to atmosphere -k*dT/dy = h(T-T_ambient)

But I am getting different results in the Matlab file and Python file. The radiation boundary condition works perfectly for both of the files, but the roller in the contact boundary condition produces different results. The results differ by 4-5°C when the temperature of the roller is low and the difference increases to about 10-11°C when the temperature of the roller increases by 100°C.

Matlab code:

function pde()
% Constants

Ny = 1000;            % Number of spatial grid points
Nt = 5000;           % Number of time steps
velocity_x = 15/60; % line speed in m/s
% Thermal properties of different layers
layer1_thickness = 0.001;  % Thickness of layer 1 (meters)
layer1_k = 2.15;            % Thermal conductivity of layer 1 (W/m-K)
layer2_thickness = 0.001;  % Thickness of layer 2 (meters)
layer2_k = layer1_k;            % Thermal conductivity of layer 2 (W/m-K)
rho_1 = 1950;           %density of layer 1 (kg/m³)
rho_2 = 1950;           %density of layer 2 (kg/m³)
cp_1 = 2100;            %specific heat capacity of layer 1 (J/kg-K)
cp_2 = 2100;            %specific heat capacity of layer 2 (J/kg-K)
T_roller = 120;         %Temperature of hot roller (°C)
h_roller = 500;         %Convective heat transfer coefficient (W/m2-K)



Ly = layer1_thickness+layer2_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)

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

radiative_flux = 150e3;
absorption = 50;
performance = 80;
Net_radiative_intensity = radiative_flux*(absorption/100)*(performance/100);

dia_roller = 0.5;                   %diameter of hot roller (m)
contact_angle = 180;                %contact angle/wrap angle for hot_roller in °

y_1 = pi*dia_roller*contact_angle/360;              %in contact with hot roller
y_2 = 1.3;                                          %after hot roller


t1 = y_1/velocity_x;
t2 = y_2/velocity_x;


T = t1+t2;

nt1 = ceil(t1/(t1+t2)*Nt);
nt2 = Nt-nt1;


t11 = linspace(0,t1,nt1);
t12 = linspace(t1,t1+t2,nt2);



%pdepe settings
m = 0;  %for 1-D cartesian coordinates with no symmetry
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_roller*(ul-T_roller); %% in contact with roller at bottom element
    pl = Net_radiative_intensity; %% radiation at bottom element
    %pl = -h_air*(ul-T_ambient);
    ql = 1;   
    pr = h_air*(ur-T_ambient);
    %pr = -Net_radiative_intensity;   
    qr = 1;

  else
    pl = -h_air*(ul-T_ambient);
    ql = 1;
    pr = h_air*(ur-T_ambient);
    %pr = -Net_radiative_intensity;
    qr = 1;
  end
end
end

Result: Matlab_radiation_bottom_element

Matlab_roller_contact_bottom_element

Python code:

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# Constants
Ny = 1000           # Number of spatial grid points
Nt = 5000         # Number of time steps
velocity_x = 15 / 60  # Line speed in m/s

# Thermal properties of different layers
layer1_thickness = 0.001  # Thickness of layer 1 (meters)
layer1_k = 2.15          # Thermal conductivity of layer 1 (W/m-K)
rho_1 = 1950              # Density of layer 1 (kg/m³)
cp_1 = 2100              # Specific heat capacity of layer 1 (J/kg-K)

layer2_thickness = 0.001  # Thickness of layer 2 (meters)
layer2_k = 2.15           # Thermal conductivity of layer 2 (W/m-K)
rho_2 = rho_1             # Density of layer 2 (kg/m³)
cp_2 = cp_1               # Specific heat capacity of layer 2 (J/kg-K)

T_roller = 120           # Temperature of hot roller (°C)
h_roller = 500           # Convective heat transfer coefficient (W/m²-K)

T_ambient = 25           # Ambient temperature (°C)
h_air = 12               # Convective heat transfer coefficient of air (W/m²-K)

Ly = layer1_thickness + layer2_thickness

# Discretization in space
ny1 = int(np.ceil(layer2_thickness / Ly * (Ny)))  
ny2 = Ny - ny1
                                        
y1 = np.linspace(0, layer2_thickness, ny1)
y2 = np.linspace(layer2_thickness, Ly, ny2 + 1)[1:]   
y = np.concatenate((y1, y2))


# Initial condition function
def icfun(y, ny1, ny2):
    u_initial = np.zeros_like(y)
    u_initial[:ny1] = 25  # Temperature for layer 2
    u_initial[ny1:] = 25  # Temperature for layer 1
    
    return u_initial


# Initial condition
initial_condition = icfun(y, ny1, ny2)


def pdefun(t, u):
    dudt = np.zeros_like(u)
    
    # Internal nodes for layer 2
    for i in range(1, Ny-1):
       
        if y[i] <= layer2_thickness:  # Layer 2
            dudt[i] = layer2_k / (rho_2 * cp_2) * (u[i+1] - 2*u[i] + u[i-1]) / (y[i-1] - y[i])**2
           
        else:  # Layer 1
            
            dudt[i] = layer1_k / (rho_1 * cp_1) * (u[i+1] - 2*u[i] + u[i-1]) / (y[i-1] - y[i])**2
    
    # Boundary conditions
    if t <= t1:
        #dudt[0] = layer2_k / (rho_2 * cp_2) * (u[1] - u[0]) / (y[1] - y[0])**2 - h_air / (rho_2 * cp_2) * (u[0] - T_ambient) / (y[1] - y[0])
        #dudt[0] = layer2_k / (rho_2 * cp_2) * (u[1] - u[0]) / (y[1] - y[0])**2 - h_roller / (rho_2 * cp_2) * (u[0] - T_roller) / (y[1] - y[0])     ##in contact with roller at bottom element
        dudt[0] = layer2_k / (rho_2 * cp_2) * (u[1] - u[0]) / (y[1] - y[0])**2 + Net_radiative_intensity / (rho_2 * cp_2 * (y[1] - y[0]))           ##radiation at bottom element
        dudt[-1] = layer1_k / (rho_1 * cp_1) * (u[-2] - u[-1]) / (y[-1] - y[-2])**2 - h_air / (rho_1 * cp_1) * (u[-1] - T_ambient) / (y[-1] - y[-2])
        
        
    else:
        dudt[0] = layer2_k / (rho_2 * cp_2) * (u[1] - u[0]) / (y[1] - y[0])**2 - h_air / (rho_2 * cp_2) * (u[0] - T_ambient) / (y[1] - y[0])
        dudt[-1] = layer1_k / (rho_1 * cp_1) * (u[-2] - u[-1]) / (y[-1] - y[-2])**2 - h_air / (rho_1 * cp_1) * (u[-1] - T_ambient) / (y[-1] - y[-2])
        
        
        
    return dudt



# Solve the PDE
dia_roller = 0.5
contact_angle = 180

radiative_flux = 150e3
absorption = 50
performance = 80
Net_radiative_intensity = radiative_flux*(absorption/100)*(performance/100)

y_1 = np.pi * dia_roller * contact_angle / 360
y_2 = 1.3

# Discretization in time
t1 = y_1 / velocity_x
t2 = y_2 / velocity_x
T = t1 + t2
t_eval = np.linspace(0, T, Nt)
sol = solve_ivp(pdefun, (0, T), initial_condition, method='BDF', atol=1e-4, rtol=1e-4, t_eval=t_eval)

# Extract the solution
u = sol.y.T

# Plot results
plt.plot(t_eval, u[:, 0], label='Temperature at bottom surface')
plt.plot(t_eval, u[:, ny1], label='Temperature at interface')
plt.plot(t_eval, u[:, -1], label='Temperature at upper surface')
plt.xlabel('Time (s)')
plt.ylabel('Temperature (°C)')
plt.legend()

# Set x-ticks at a step of 1 second
x_ticks = np.arange(0, int(T) + 1, 1)  # Range from 0 to T (rounded up) with a step of 1 second
plt.xticks(x_ticks)

# Set y-axis limits dynamically based on the minimum and maximum values of u[:, 0] and u[:, -1]
y_min = min(np.min(u[:, 0]), np.min(u[:, -1]))
y_max = max(np.max(u[:, 0]), np.max(u[:, -1]))

# Round down y_min to the nearest multiple of 5
y_min = np.floor(y_min / 5) * 5 - 5

# Round up y_max to the nearest multiple of 5 plus 5
y_max = np.ceil(y_max / 5) * 5 + 5

# Set y-axis limits
plt.ylim(y_min, y_max)

# Define y-ticks with a step of 5°C within the y-axis limits
y_ticks = np.arange(y_min, y_max, 5)
plt.yticks(y_ticks)
plt.grid(True)

plt.show()

Results:

Python_radiation_bottom_element

Python_roller_contact_bottom_element

I tried to check whether the time points are defined exactly or not and checked with radiation boundary condition which works alright. But the strange behaviour is that the curve shape or behaviour is the same.

The area of interest is the inclusion of roller boundary condition at bottom element in python code (dudut[0]) gives different result as compared to the boundary condition at left element/bottom element(pl).

Below, you can see the comparison of results in matlab and python when the total time steps is 100 and number of divided elements across the total thickness is 1000.

difference in results of python and matlab

Could anyone please help or suggest how can I resolve this discrepancy?

Your Matlab function is at fault. You have muddled up the ur and ul values at the left (l) and right (r) boundaries. I can reproduce this error in my Fortran code (given below) if I fake the source terms at the boundaries.

Where you have

pl = -h_roller*(ur-T_roller); %% in contact with roller at bottom element

it should be

pl = -h_roller*(ul-T_roller); %% in contact with roller at bottom element

(i.e. ur should be replaced by ul).

Where you have

pr = h_air*(ul-T_ambient);

it should be

pr = h_air*(ur-T_ambient);

(i.e. ul should be replaced by ur).

I wrote a Fortran code to do the simple diffusion equation in the rise phase (up to t1=PI). It agrees with your Python solution in this phase (but note my earlier comments about applying a factor 2 to both your dudt[0] and dudt[-1] to accommodate the fact that the end node only controls a width (1/2)dy. If you can run Fortran then you can simulate your Matlab error by uncommenting the lines indicated.

!============================================================================================================
!                du       d2u
! Solve   rho cp --  -  k ---   =  source   (where source comes from boundary fluxes replacing diffusive flux)
!                dt       dy2
!
! The source term is written as    S + Sp.u     where Sp should be negative for stability
!
! Use Crank-Nicolson timestepping (average flux and source over "old" and "new" times)
!
! Discretise in the form
!         A_i u_i-1 + B_i u_i + C_i u_i+1 = D_i
! and solve by the tri-diagonal matrix algorithm at every timestep.
!
!============================================================================================================
program diffusion
   use iso_fortran_env, only: real64
   implicit none

   integer, parameter :: N = 1000
   integer, parameter :: Nt = 5000

   real(real64), dimension(N) :: A, B, C, D, u

   real(real64) :: dy =  0.002 / ( N - 1 )
   real(real64) :: rho = 1950.0
   real(real64) :: cp = 2100.0
   real(real64) :: k = 2.15
   real(real64) :: T_roller = 120.0
   real(real64) :: T_air    = 25.0
   real(real64) :: h_roller = 500.0
   real(real64) :: h_air = 12.0
   real(real64) :: dt = 3.14159 / Nt

   integer      j
   real(real64) t
   real(real64) kfac, tfac
   real(real64) SL, SR, SPL, SPR

   ! The correct boundary conditions (in the rise phase)
   SL = h_roller * T_roller;   SPL = -h_roller
   SR = h_air    * T_air   ;   SPR = -h_air

   kfac = k * dt / ( rho * cp * dy ** 2 )
   tfac = dt / ( rho * cp * dy )

   A = -0.5 * kfac;   A(1) = 0;   A(N) = -kfac
   C = -0.5 * kfac;   C(N) = 0;   C(1) = -kfac

   open( 10, file="u.out" )
   u = T_air
   t = 0.0

   do j = 1, Nt
      t = t + dt

!         ! If you want to fake the error in your Matlab solution, uncomment the lines below
!         SL = h_roller * ( T_roller - u(N) );   SPL = 0
!         SR = h_air    * ( T_air    - u(1) );   SPR = 0

      B = 1 + kfac;   B(1) = B(1) - tfac * SPL;   B(N) = B(N) - tfac * SPR
      D(2:N-1) = u(2:N-1) - 0.5 * kfac * ( -u(1:N-2) + 2 * u(2:N-1) - u(3:N) )
      D(1) = u(1) - kfac * ( u(1) - u(2  ) ) + tfac * ( SL + SL + SPL * u(1) )
      D(N) = u(N) - kfac * ( u(N) - u(N-1) ) + tfac * ( SR + SR + SPR * u(N) )
      call tdma( A, B, C, D, u )
      write( 10, "(3(1x,es12.4))" ) t, u(1), u(N)

   end do

contains

   subroutine tdma( A, B, C, D, X )
      !******************************************************************************
      ! Solve, using the Thomas Algorithm (TDMA), the tri-diagonal system
      !     A_i X_i-1 + B_i X_i + C_i X_i+1 = D_i,     i = 1, n
      !
      ! Effectively, this is the n x n matrix equation.
      ! a(i), b(i), c(i) are the non-zero diagonals of the matrix and d(i) is the rhs.
      ! a(1) and c(n) are not used.
      !******************************************************************************
      real(real64), intent(in ) :: A(:), B(:), C(:), D(:)
      real(real64), intent(out) :: X(:)

      integer n
      integer i
      real(real64), allocatable :: P(:), Q(:)
      real(real64) denominator

      n = size( D )
      allocate( P(n), Q(n) )
      X = 0

      ! Forward pass
      i = 1
      denominator = b(i)
      P(i) = -c(i) / denominator
      Q(i) =  d(i) / denominator
      do i = 2, n
         denominator = b(i) + a(i) * P(i-1)
         if ( abs( denominator ) < 1.0e-10 ) return
         P(i) =  -c(i)                   / denominator
         Q(i) = ( d(i) - a(i) * Q(i-1) ) / denominator
      end do

      ! Backward pass
      X(n) = Q(n)
      do i = n - 1, 1, -1
         X(i) = P(i) * X(i+1) + Q(i)
      end do

      deallocate( P, Q )

   end subroutine tdma

end program diffusion

Output (corresponds to your 4th graph; i.e. the second python graph). Note that, in the interest of simplicity and to find your original error, I have only integrated as far as t=t1. The latter phase can be continued simply by amending the source terms after this point.