# 1D Heat Equation Explicit Scheme (fixed)

11 views (last 30 days)
Gabriel Quattromani on 19 Apr 2022
Edited: Gabriel Quattromani on 30 Apr 2022
clear all
p=2050; % density(kg/m^3)
k=0.52; % conductivity(W/m-K)
c=1840; % specific heat capacity(J/kg-K)
T_ini=20; % initial temperature(K)
T_inf=-15; % external temperature(K)
alpha=k/(p*c); % Let thermal diffusivity be alpha
t=60*60*60*24; %time in seconds
DELTA_t=100; % change in time per timestep
M=t/DELTA_t; % number of timesteps
N=300; % number of steps of depth
L=3; % maximum depth being analyzed
DELTA_x=L/(N); % change in depth per step
t_diff=t/M;
time=0:t_diff:t;
%explicit scheme
% Boundary Conditions
T(1:N,1:M+1)=T_ini;
T(1,1:M+1)=T_inf;
CFL=(alpha*DELTA_t)/(DELTA_x^2); % Courant number, which should be less than 0.5
% Finding temperature for next timestep for each depth
for j=1:M
for i=2:N-1
T(i,j+1)=T(i,j)+((alpha*DELTA_t)/((DELTA_x)^2))*(T(i+1,j)-2*T(i,j)+T(i-1,j));
end
end
Z=N*(0.1/L); % t=0.1 seconds in terms of the time step
xm=(69/300)*L;
plot(time,T(Z,:), 'b',time,T(69,:), 'r')
xlabel('Time')
ylabel('Temperature')
grid on
set(gca,'fontsize',16)
legend('0.1','xm')

Torsten on 19 Apr 2022
You use
for i=0:0.1:50
...
T(i,j) = ...
But only positive integer values are possible as array indices. Thus something like T(0,j), T(0.1,j) etc is not allowed.
Additionally, you need a boundary condition at x=xstart and x=xend. What are these conditions ?
Delta_x and Delta_t are scalar values. Thus addressing them in the if-statement
if (alpha*DELTA_t(j))/(DELTA_x(i)^2)<=0.5
as vector is wrong.

Alan Stevens on 19 Apr 2022
With
for i=0:0.1:50
you are trying to index variables with zero. However, Matlab's indexing starts at 1.

R2021b

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!