solving 1D transient heat equation using finite difference method

Hello,

I am trying to solve a 1D transient heat equation using the finite difference method for different radii from 1 to 5 cm, with adiabatic bounday conditions as shown in the picture. I have attached my discretization method as well to give a better insight into my problem. I do not know where I did wrong that the results are not as expected. The solution should obviously be exponential. What would you suggest?

% 1D Cylindrical Heat Equation with Heat Generation - Explicit FD
clear
close all
clc
% Parameters
R = 1*10^-2;         % Radius of the cylinder (meters)
t_final = 300;           % Total simulation time (seconds)
cp = 2100;    % specific heat
h = 200;    %conduction coefficient
k = 0.5;    %thermal conductivity
rho = 1100; %density
alpha = k/(rho*cp);
heatGeneration = 0.1;  % Heat generation term

% Spatial discretization
Nr = 100;         % Number of radial nodes
dr = R / (Nr - 1); % Radial step size
r = 0:dr:R;

% Temporal discretization
Nt = 100;        % Number of time steps
dt = t_final/(Nt-1);     % Temporal step size
time = 0:dt:t_final;

% Initialize temperature field
u = zeros(Nr, Nt);

% Initial condition
u(:, 1) = -196;  % Initial temperature distribution

% Explicit finite difference scheme
for n = 1:Nt-1
    for i = 2:Nr-1
        % Discretization in radial direction
        du_dr = (u(i+1, n) - u(i-1, n)) / (2 * dr);
        d2u_dr2 = (u(i+1, n) - 2 * u(i, n) + u(i-1, n)) / dr^2;
        
        % Update temperature
        u(i, n+1) = u(i, n) + alpha * dt * (1/r(i) * du_dr + d2u_dr2) + (heatGeneration*dt)/(rho*cp);
    end
    
    % Boundary conditions
    u(1,n+1) = u(2,n+1);
    u(Nr, n+1) = u(Nr-1, n+1);
end

% Plot the temperature distribution over time
figure;
plot(time,u);
title('1D Cylindrical Heat Equation with Heat Generation');
xlabel('Time (s)');
ylabel('Temperature');

[Torsten]

You start with a flat profile for T, you have adiabatic boundary conditions and your heat generation rate is constant and independent of r. As a consequence, you get flat profiles over time. What do you expect to be exponential ?

[Hasan Humeidat]

Hmm, but I solved this analytically and obtained an exponential function. Regardless, in case I want to change the conditions where the boundary condition at 1cm is a convective heat transfer, how can I translate this to the method I am using? (discretization wise).

[Torsten]

For comparison with your discretization method.

In order to obtain similar results with convective heat transfer at 1 cm, you will have to switch to implicit instead of explicit Euler. Else the time-stepsize to get physical senseful solutions seems to become extremely small (CFL condition says dt < 0.25*dr^2, thus dt < 0.25*(1e-2*1e-2)^2, thus dt < 2.5e-9.)

exercise()

function exercise
% Parameters
R = 1*10^-2;         % Radius of the cylinder (meters)
t_final = 3000;           % Total simulation time (seconds)
cp = 2100;    % specific heat
k = 0.5;    %thermal conductivity
rho = 1100; %density
alpha = k/(rho*cp);
heatGeneration = 10000;  % Heat generation term
kinf = 10;   %heat transfer coefficient
Tinf = -196;  % external temperature
x = linspace(0,R,100);
t = linspace(0,t_final,100);
m = 1;
sol = pdepe(m,@heatcyl,@heatic,@heatbc,x,t);
u = sol(:,:,1);
figure()
plot(t,sol(:,1))
xlabel('Time')
ylabel('Temperature u(0,t)')
title('Temperature change at center of cylinder')
figure()
plot(x,sol(end,:))
xlabel('Radius')
ylabel('Temperature')
title('Temperature at final time')

function [c,f,s] = heatcyl(x,t,u,dudx)
  c = rho*cp;
  f = k*dudx;
  s = heatGeneration;
end
function u0 = heatic(x)
  u0 = -196;
end
function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t)
  pl = 0; %ignored by solver since m=1
  ql = 0; %ignored by solver since m=1
  pr = kinf*(ur-Tinf);
  qr = 1; 
end

end

[Torsten]

You need dt = 3000/200000 = 0.015 to get a stable solution with Explicit Euler:

clear
close all
clc
% Parameters
R = 1*10^-2;         % Radius of the cylinder (meters)
t_final = 3000;           % Total simulation time (seconds)
cp = 2100;    % specific heat
k = 0.5;    %thermal conductivity
rho = 1100; %density
alpha = k/(rho*cp);
heatGeneration = 10000;  % Heat generation term
kinf = 10;  % exterior heat transfer coefficient
Tinf = -196;   % exterior temperature

% Spatial discretization
Nr = 100;         % Number of radial nodes
dr = R / (Nr - 1); % Radial step size
r = 0:dr:R;

% Temporal discretization
Nt = 200000;        % Number of time steps
dt = t_final/(Nt-1);     % Temporal step size
time = 0:dt:t_final;

% Initialize temperature field
u = zeros(Nr, Nt);

% Initial condition
u(:, 1) = -196;  % Initial temperature distribution

% Explicit finite difference scheme
for n = 1:Nt-1
    for i = 2:Nr-1
        % Discretization in radial direction
        du_dr = (u(i+1, n) - u(i-1, n)) / (2 * dr);
        d2u_dr2 = (u(i+1, n) - 2 * u(i, n) + u(i-1, n)) / dr^2;
        
        % Update temperature
        u(i, n+1) = u(i, n) + alpha * dt * (1/r(i) * du_dr + d2u_dr2) + (heatGeneration*dt)/(rho*cp);
        %u(i,n+1) = u(i,n) + alpha * dt * 1/r(i)*((r(i)+dr/2)*(u(i+1,n)-u(i,n)) - (r(i)-dr/2)*(u(i,n)-u(i-1,n)))/dr^2 + (heatGeneration*dt)/(rho*cp);
    end
    % Boundary conditions
    u(1,n+1) = u(2,n+1);
    %u(Nr,n+1) = u(Nr-1,n+1);
    u(Nr,n+1) = u(Nr,n) + alpha * dt * 1/R*(R*(-kinf/k)*(u(Nr,n)-Tinf)-(R-dr/2)*(u(Nr,n)-u(Nr-1,n))/dr)/(dr/2) + (heatGeneration*dt)/(rho*cp); 
    %u(Nr, n+1) = u(Nr-1, n+1);
end

% Plot the temperature distribution over time
figure;
plot(time,u(1,:));
title('1D Cylindrical Heat Equation with Heat Generation');
xlabel('Time (s)');
ylabel('Temperature');

figure;
plot(r,u(:,end));
title('1D Cylindrical Heat Equation with Heat Generation');
xlabel('Position (m)');
ylabel('Temperature');

[Hasan Humeidat]

Yes, now I see where I had a problem in. Thank you very much, this is very helpful.