1D Heat Conduction using Eulers Explicit discretisation

11 views (last 30 days)
I am trying to solve a problem regarding heat conduction. I have written down a code and using the method that I know I tried to solve it but when I'm making the T0 matrix it gives me error because I have to start at 0 Kelvin and it gives me an error. Attached is an image of the question I'm trying to solve here.
%% Clearing
clc
clear
close all
%% Initialisation
L=1;
t=100;
k=0.01;
nt=500;
dx=0.05;
n=L/dx;
dt=0.5;
x=0:dx:L;
alpha=k*dt/dx^2;
T0=0*ones(1,n);
T1=20*ones(1,n);
T0(1) = 0;
T0(end) = 20;
%% Solving
for j=1:nt
for i=2:n-1
T1(i)=T0(i)+alpha*(T0(i+1)-2*T0(i)+T0(i-1));
end
T0=T1;
end
%% plotting
timeset=[0 1 5 10 100]; % time instants required
for i=1:length(timeset) % 1 to number of elements in timeset
z=timeset(i); % time instant
ntset=z/dt+1; % time instant extracted for the corresponding step
figure, plot(x,T1(:,ntset)) % plotting temperature vs distance
ylabel('Temperature (°C)') % labels y axis
xlabel('Distance (m)') % label x axis
title('Distribution of Temperature') % title of graph
grid minor % detailed view
end
Index in position 2 exceeds array bounds. Index must not exceed 20.

Accepted Answer

William Rose
William Rose on 8 Aug 2022
Note that length(x)=21, but your n=20. This causes problems. I recommend that you initialize a bit differently, as shown below.
I do not understand why you have vectors for T0 and T1. I recommend that T0 and T1 be scalars. I recommend that you create a single array, called T. Each row of T corresponds to one time. Each column of T corresponds to one position along the bar.
%% Initialisation
L=1.0; dx=0.05; x=0:dx:L; nx=length(x); %space
tmax=30; dt=0.05; t=0:dt:tmax; nt=length(t); %time
alpha=0.01; %thermal diffusivity, m^2/s
T0=0; %left temperature
T1=20; %right termperature
T=zeros(nt,nx); %allocate array for temperature
T(:,1)=T0; %make the left temp = T0 at all times
T(1,:)=T0; %make the initial temp=T0 at all locations
T(:,end)=T1; %make the right temp = T1 at all times
Solve the problem.
for i=2:nt
for j=2:nx-1
T(i,j)=T(i-1,j)+alpha*(dt/dx^2)*(T(i-1,j+1)-2*T(i-1,j)+T(i-1,j-1));
end
end
Plot some results.
plot(x,T(1,:),'-r',x,T(101,:),'-y',x,T(201,:),'-g',x,T(301,:),'-c',...
x,T(401,:),'-b',x,T(501,:),'-m',x,T(601,:),'-r');
grid on; title('Temperature versus position and time')
xlabel('Position (m)'); ylabel('Temperature (K)')
legend('t=0 s','t=5 s','t=10 s','t=15 s','t=20 s','t=25 s','t=30 s');
When I ran the code with dt=0.5 s, the solution oscillated, and gave results such as T=+/-10^167 after 100 seconds. This happens if the time step is too large. So I reduced the time step by a factor of 10: dt=0.05 s. Then I got reasonable results, as seen above.
Learn how to use surf() to make a nice 3D plot with axes position, time, and temperature.
  2 Comments
Yashraj Randad
Yashraj Randad on 8 Aug 2022
Thank you William. I appreciate the answer. You're right about the T0 and T1 being scalars, initialising them as vectors added an unnecessary layer of complication that wasn't needed.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!