solving 1D transient heat equation using finite difference method
210 views (last 30 days)
Show older comments
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');
5 Comments
Torsten
on 21 Nov 2023
Edited: Torsten
on 22 Nov 2023
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');
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!