Matlab in applying Finite Difference for Temp. distribution/ rate of heat flow

[mejia198021]

Hi:

I need some assistance or feedback on a Matlab program that I working on by applying the method of finite difference to calculate and output the steady state space distribution of temperature and resulting rate of hear flow in an I Beam made of Duralumin. I have written a partial code but I am not sure if I am on the right path in solving this problem.

First: I divided the I Beam into a nodal grid, which would include 289 interior nodes (not including boundaries). The beam measures 2.4cm(.024m) x 3cm (.03m) and dx=dy=.001m.

Basically in my code I inputted all the variables given, which includes the upper/lower surface temperatures, thermoconductivity (K). Next, I created a outline of the beam by inputing the boundaries/ initial conditions. Then I applied an interior node equation that I looked up in my Heat Transfer Text:

T(m,n+1) + T(m,n-1) + T(m+1,n) +T(m-1,n) -4*T(m,n)=0

I applied this equation in 3 FOR loops for the interior nodes: Upper, center and lower sections of the I Beam. Next, I used a heat flux equation to calculate the rate of heat flow by applying: qx=k*((T1-T2)/L)

I used the same method as above to calculate this equation by stating the boundaries/ initial conditions and then apply 3 FOR loops for the three sections of the I Beam.

At this point I am not sure if I am applying the correct equations or methods. I get Temp. outputs for each interior nodes but at some pts. it seems that it is not correct. For the Heat Flux I obtain section that are 0 and neg., so that it wrong.

Any help or tips would be appreciated. Below is my code and I have attached the problem (NMProb.doc) if you would like to take a look at it. Thank you for your time.


clc
clear all;

T1 = input('Fixed Temperature at the upper surface Tt[60 deg C] : '); % deg C fixed Temperature at the upper surface.
T2 = input('Fixed Temperature at the lower surface. Tl[15 deg C] : '); % deg C fixed Temperature at the lower surface.
k = input('Enter the thermal conductivity k[Duralumin 164 W/m deg C] : '); % W/m deg C (thermal conductivity).



%Intilializing Constant Variables
k=164; %Thermoconductivity
Dn=2787; %Dn is defined as Density
Qi=0;
Qw=0;

%Outline dimensions of I beam
Wt=0.024; %Width
Ht=0.03; %Height

%Temperature Conditions
Tt=60; %Upper face Temperature (333.15K)
Tb=15; %Lower face Temperature (288.15K)
Ti=20; %Initial Temperature of interior nodes (293.15K)

%Conditions of nodal points
dx=.001; %Delta x - distance between each node
dy=.001; %Delta y - distance between each node

M=31; %Number of nodes in the Y direction -i
N=25; %Number of nodes in the X direction -j

%Applying Boundaries/ Initial Conditions to I Beam(Temperature)
T=zeros(31,25);
T(1,:)=T1; %Top Surface
T(31,:)=T2; %Bottom Surface
%Left
T(1:6,1)=0; %Upper Top Left
T(6,1:10)=0; %Bottom Top Left
T(6:26,10)=0; %Inner Center Left
T(26,1:10)=0; %Lower Bottom Left
T(26:31,1)=0; %Lower Left
%Right
T(1:6,25)=0; %Upper Top Right
T(6,16:25)=0; %Bottom Top Right
T(6:26,16)=0; %Inner Center Right
T(26,16:25)=0; %Lower Bottom Right
T(26:31,25)=0; %Lower Right

T(2:5,2:24)=Ti; %Upper Section
T(6:26,11:15)=Ti; %Center Section
T(27:30,2:24)=Ti; %Lower Section

T;

%Temperature Distribution of Interior Nodes (289 Nodal Points)
for(i=2:5)
for(j=2:24)
%Temperature of Nodal Points for Upper Portion
T(i,j)=(T(i,j+1)+T(i,j-1)+T(i+1,j)+T(i-1,j))/4;
end
end

for(i=6:26)
for(j=11:15)
%Temperature of Nodal Points for Center Portion
T(i,j)=(T(i,j+1)+T(i,j-1)+T(i+1,j)+T(i-1,j))/4;
end
end

for(i=27:30)
for(j=2:24)
%Temperature of Nodal Points for Bottom Portion
T(i,j)=(T(i,j+1)+T(i,j-1)+T(i+1,j)+T(i-1,j))/4;
end
end





Q=zeros(31,25);
Q(1,:)=Qw; %Top Surface
Q(31,:)=Qw; %Bottom Surface
%Left
Q(1:6,1)=Qw; %Upper Top Left
Q(6,1:10)=Qw; %Bottom Top Left
Q(6:26,10)=Qw; %Inner Center Left
Q(26,1:10)=Qw; %Lower Bottom Left
Q(26:31,1)=Qw; %Lower Left
%Right
Q(1:6,25)=Qw; %Upper Top Right
Q(6,16:25)=Qw; %Bottom Top Right
Q(6:26,16)=Qw; %Inner Center Right
Q(26,16:25)=Qw; %Lower Bottom Right
Q(26:31,25)=Qw; %Lower Right

Q(2:5,2:24)=Qi; %Upper Section
Q(6:26,11:15)=Qi; %Center Section
Q(27:30,2:24)=Qi; %Lower Section


%Rate of Heat Flow of Interior Nodes (289 Nodal Points
% for(i=2:5)
% for(j=2:24)
% %Temperature of Nodal Points for Upper Portion
% Q(i,j)=k*((T(i,j+1)-T(i,j))/dx);
% end
% end

for(i=6:26)
for(j=11:15)
%Temperature of Nodal Points for Center Portion
T(i,j)=(T(i,j+1)+T(i,j-1)+T(i+1,j)+T(i-1,j))/4;
end
end
%
% for(i=27:30)
% for(j=2:24)
% %Temperature of Nodal Points for Bottom Portion
% T(i,j)=(T(i,j+1)+T(i,j-1)+T(i+1,j)+T(i-1,j))/4;
% end
% end


Q;

fprintf(1,'Press enter and type T for Temp. Distribution');pause;

 

[mmwave]

fprintf(1,'Press enter and type Q for Rate of Heat Flow');pause;

%T=Temperature Distribution Matrix
%Q=Rate of Heat Flow Matrix

 

[tacman]

disp(' ')
fprintf(1,'Press enter and type Q for Rate of Heat Flow');pause;disp(' ')
fprintf(1,'Press enter and type ALL for both T and Q');pause;disp(' ')

input('Enter your choice:','s');

if ans== 'T'
disp(' ');
disp('Temperature Distribution: ');
disp(T);
elseif ans== 'Q'
disp(' ');
disp('Rate of Heat Flow: ');
disp(Q);
elseif ans== 'ALL'
disp(' ');
disp('Temperature Distribution: ');
disp(T);
disp(' ');
disp('Rate of Heat Flow: ');
disp(Q);
else
disp('Invalid Choice!');
end



Hi there,

Thank you for sharing your code and the problem statement. From what I can see, you are on the right track in terms of applying the method of finite difference to solve for the steady state temperature distribution and heat flow rate in an I Beam. However, there are a few areas that could be improved upon to ensure the accuracy of your results.

Firstly, I recommend double checking your input values for the dimensions of the I Beam and the thermal conductivity. In your code, you have defined the thermal conductivity as 164 W/m deg C, but in the problem statement, it is listed as 164 W/m K. This may be causing errors in your calculations. Also, make sure that the dimensions of your I Beam are consistent with the problem statement. In your code, the width is listed as 0.024m, but in the problem statement it is listed as 2.4cm, which would be 0.024m. These small discrepancies can have a big impact on your results.

Next, I would suggest using a more accurate equation for the interior nodes. The equation you have used, T(m,n+1) + T(m,n-1) + T(m+1,n) +T(m-1,n) -4*T(m,n)=0, is accurate for a 2D square grid, but since the I Beam is not a perfect square, a more accurate equation would be T(m,n+1) + T(m,n-1) + T(m+1,n) +T(m-1,n) -4*T(m,n) = (dx^2 * dy^2 * Q(m,n)) / (2 * k * (dx^2 + dy^2)), where Q