1D Heat Equation Explicit Scheme (fixed)

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')

Answer 1

[Torsten]

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.

Answer 2

[Alan Stevens]

With

for i=0:0.1:50

you are trying to index variables with zero. However, Matlab's indexing starts at 1.