1D Heat Conduction using Eulers Explicit discretisation

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.

Answer 1

[William Rose]

@Yashraj Randad,

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.

[Yashraj Randad]

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.

[William Rose]

@Yashraj Randad, you're welcome.