1D Heat Equation Explicit Scheme (fixed)

6 views (last 30 days)
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')

Accepted Answer

Torsten
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.

More Answers (1)

Alan Stevens
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.

Community Treasure Hunt

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

Start Hunting!