1D Heat Conduction using Eulers Explicit discretisation

6 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

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!