Problem with time in temperature distribution

4 次查看（过去 30 天）
Green Sal
Green Sal 2011-11-29
Hello guys, i am trying to find the temperature distribution in a duct that is insulated from both sides, constant temperature on one side, and const. heat flux from the top. in the inner duct steam flows at constant temperature. The following is the code that i used but i know something is not working. I am trying to validate my code by changing the values of the heat flux but the temperature distribution seems to stay the same. Also, if i run the code at 10 seconds and 100 seconds, the results are the same. If anybody understand such problems please help. Thank you.
clear
clc
T = zeros(20,20);
% Initial condition
% Set all material to 20C ambient temp
for y = 1:20;
for x = 1:20;
T(y,x) = 20;
end
end
h=[1000000;0;0;0];
dt = .1; %(Time increment, Seconds)
p = 9000; %(Density of Copper, kg/m^3)
Cp = 380; %(Specific heat of Copper, J/kg*K)
qT = 500; %(Solar heat flux, W/m^2)
TL = 50; %(Temperature of the left side of the block, Celsius)
Tsteam = 400; % (Temperature of steam inside the duct, Celsius)
l = 0.005; % 0.005 Meters
K = 400; % (Thermal conductivity of Coppy,W/m*C)
Alpha = K/(p*Cp); % (Thermal diffusivity)
Fo = (Alpha*dt)/(l^2); % Seconds
Bi=(h*l)/K;
b(1)=det(T);
fprintf('Enter Time (s):\n')
t = input(' ');
A=t/dt; %time/step increment size
%instant steam at center
for y = 6:15;
for x = 6:15;
T(y,x) = Tsteam;
end;
end;
%for whole matrix
for n = 1:A
% BOUNDRY CONDITIONS
% Node 1 - Top left node
T(1,1) = T(1,1)*(1-4*Fo)+2*Fo*(T(1,2) + T(2,1) + (2/K)*TL+ qT*l/K);
% Node 21 - Top right node
T(1,20) = T(1,20)*(1-4*Fo)+2*Fo*(T(1,19) + T(2,20)+ qT*l/K);
% Left nodes - Left side BC
for x = 1;
for y = 2:19;
T(y,x) = TL;
end
end
% Right nodes - Right surface BC (Insulated section (no qT))
for x = 20;
for y = 3:19;
T(y,x) = T(y,x)*(1-4*Fo) + Fo*(T(y-1,x) + T(y+1,x) + 2*T(y,x-1));
end;
end;
% Bottom left Node - insulated with left BC
T(20,1) = T(20,1)*(1-4*Fo) + 2*Fo*(T(19,1) + T(20,2) + (2*TL)/K);
% Bottom right Node - insulated with right BC
T(20,20) = T(20,20)*(1-4*Fo) + 2*Fo*(T(19,20) + T(20,19));
% Top nodes calcs - Top surface (Heat Flux addition)
for y = 1;
for x = 2:19;
T(y,x) = T(y,x)*(1-3*Fo) + Fo*((qT*l/K) + T(y,x-1) + T(y,x+1) + T(y+1,x));
end;
end;
% Calculating the in between nodal temperatures at each delta(time)
for y = 2:19; %inside left slab and right slab
for x = [2:5,16:19]
T(y,x) = T(y,x)*(1-4*Fo) + Fo*(T(y,x-1) + T(y,x+1) + T(y-1,x) + T(y+1,x));
end;
end;
%central nodes
for x = 6:19;%center slab portions above and below duct
for y = [2:5,16:19];
T(y,x) = T(y,x)*(1-4*Fo) + Fo*(T(y,x-1) + T(y,x+1) + T(y-1,x) + T(y+1,x));
end;
end;
% Bottom Nodes - Insulated section (no qT)
for y = 20;
for x = 2:19;
T(y,x) = T(y,x)*(1-4*Fo) + Fo*(T(y,x-1) + T(y,x+1) + 2*T(y-1,x));
end;
end;
%this part below kills the initial for loop if the duct has reached steady
%state, and displays the time in seconds
b(n+1)=det(T);
if b(n)== b(n+1);
s=(n-1)*.1;
fprintf('The Duct reached Steady State at: %4.1f seconds \n',s);
break
end
end;
Tfinal = flipud(T);
% Graphing the temperature isotherm
Spacing = (20:5:400);
[C,h]=contourf(Tfinal,Spacing);
clabel(C,h)
[val,ind] = max(Tfinal(:));
[max_x,max_y] = ind2sub(size(Tfinal),ind);
  1 个评论
Walter Roberson
Walter Roberson 2011-11-29
Image Analyst asked approximately: Have you stepped through the code with the debugger?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Grzegorz Knor
Grzegorz Knor 2011-11-29
Probably the temperature reaches a steady state before this time. Put these lines:
Tfinal = flipud(T);
Spacing = (20:5:400);
[C,h]=contourf(Tfinal,Spacing);
clabel(C,h)
drawnow
into loop and you will see changes in temperature.
  3 个评论
显示 1更早的评论
Grzegorz Knor
Grzegorz Knor 2011-11-30
delete clabel function, and increase step in Spacing:
Spacing = (20:10:400);
[C,h]=contourf(Tfinal,Spacing);
title(num2str((n-1)*.1))
drawnow
Green Sal
Green Sal 2011-11-30
Grzegorz thank you for your support. i have one last question for you. So now i know the time at which the duct reach steady state, now i would like to plot the isotherms at only t=0.25*t(ss), t=0.50*t(ss) and t(ss) is the time at which it reaches steady state.
Thank you!

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Aerodynamics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by