Heat Diffusion on a Rod over the time In class we learned an

Heat Diffusion on a Rod over the time

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

in this homework we will solve the above 1-D heat equation numerically. There is a 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\');

QUESTION...Modify this Matlab code so that the initial temperature at the rod is 200^o half length of the rod and 0^o at the rest of the rod. This rod is set to 0 temperature at the ends for all times. Plot how the temperature diffuses along the rod, over the time and label the different times

Do you expect to see the heat diffusing the way you see in the graph? Justify your answer by comparing this graph by the analytical solution that we found in the class. (you can choose a certain number of terms in your series of solution depending on how much precise you look for. Please state how many number of terms you included.)

You can plot the analytical solution on the same graph with numerical solution.

aT a-T

Solution

%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\');