Heat Diffusion on a Rod over the time

In class we learned analytical solution of 1-D heat equation

Inline image 1

In this homework we will solve the above 1-D heat equation numerically. Below is the Matlab code which simulates finite difference method to solve the above 1-D heat equation.

%1-D Heat equation

%example 1 at page 782

%dT/dt=c.d^2T/dx^2

%T(x,t)=temperature along the rod

%by finite difference method

%(T(x,t+dt)-T(x,t))/dt = (T(x+dx,t)-2T(x,t)+T(x-dx,t))/dx^2

%solve for T(x,t+dt) by iteration

%heat constant

close all

clear all

c=1;

L=1 %length of the rod

N=5 %# of elements

%dicretize Xspace

x_vec=0:0.2:1;

dx=x_vec(2)-x_vec(1);

%discretize time

dt=0.5/50;

t_vec=0:dt:0.5;

%temperature matrix

T_mat=zeros(length(x_vec),length(t_vec));

%boundary conditions

T_mat(1,:)=0;

T_mat(end,:)=0;

%initial conditions

T_mat(:,1)= sin(pi.*x_vec);

[tt,xx]=meshgrid(t_vec,x_vec);

subplot(2,1,1);

mesh(xx,tt,T_mat);

lamda=(c*dt/(dx^2));

for tdx=1:length(t_vec)-1

for idx=2:length(x_vec)-1

T_mat(idx,tdx+1)=T_mat(idx,tdx)+lamda*((T_mat(idx+1,tdx)-2*T_mat(idx,tdx)+T_mat(idx-1,tdx)));

end

end

%plot

subplot(2,1,2)

[tt,xx]=meshgrid(t_vec,x_vec);

mesh(xx,tt,T_mat);

figure

plot(x_vec,T_mat(:,1),x_vec,T_mat(:,11),x_vec,T_mat(:,21),x_vec,T_mat(:,31),x_vec,T_mat(:,41),x_vec,T_mat(:,51));

xlabel('Rod length (m)'); ylabel('Temperature (C)');

legend('Initially','At t=0.1sec','At t=0.2 secs','At t=0.3 secs',' At t=0.4 secs',' At t=0.5 secs');

Part2: In this part the end points of the rod is insulated. Modify the Matlab code so that the temperature distribution is 100+ 100*cos(x) initially.

Hand calculate the analytical solution for the temperature in this case.

Plot analytical and numerical solutions.

Do you expect to see the heat diffusing the way you see in the graph? Justify your answer by comparing this graph by comparing your analytical and numerical solutions.

Submit your typed or printed answers to the above questions and plot and submit analytical solution and numerical solutions.

You can add in your code, part below to implement the insulation at the ends of this rod..

for tdx=1:length(t_vec)-1

for idx=1:length(x_vec)

if idx==1

T_mat(idx,tdx+1)=T_mat(idx+1,tdx);

elseif idx==length(x_vec)

T_mat(idx,tdx+1)=T_mat(idx-1,tdx);

else

T_mat(idx,tdx+1)=T_mat(idx,tdx)+lamda*((T_mat(idx+1,tdx)-2*T_mat(idx,tdx)+T_mat(idx-1,tdx)));

end

end

end