2D heat transfer problem

Question

Tom Gayler 2016-4-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/277519-2d-heat-transfer-problem

评论： Walter Roberson 2020-5-9

Hello,

I'm doing a heat transfer problem as part of a course I am on. This is the basic assignment.

"This project concerns the calculation of temperatures inside a plate of size 0.8m× 0.4m as shown in figure 4. Temperature distributions T(x, y) satisfy Laplace’s equation in stationary equilibrium,

The plate is heated along parts of its top and left edges and cooled along parts of its lower and right edges. The faces and remainder of the edges are insulated. There is also a circular hole in the plate. The hole contains a highly conducting fluid, such that all points around the boundary of the hole have equal temperature."

I've written most of the code but cannot get it to display the answer. Here is the code; any help would be greatly appreciated.

   %set constants
  nx=320;
  ny=160;
  width = 0.8;
  length = 0.4;
  dx=width/nx;
  dy=length/ny;
phinew=zeros([ny nx]);
phiold=zeros([ny nx]);
%set up while loop
err=1;
tol=1e-6;
while err>tol
    %enforce boundary temperature
    phiold(2:40,160)=100;
    phiold(1,80:160)=100;
    phiold(280:320,1)=0;
    phiold(320,2:80)=0;
%calculate updated temp
for i=2:nx-1
    for j=2:ny-1
        phinew(i,j)=0.25*(phiold(i+1,j)+phiold(i-1,j)+phiold(i,j+1)+phiold(i,j-1));
    end
end
%top side
for i=41:319
    for j = 160
    phinew(i,j)=(1/3)*(phiold(i+1,j)+phiold(i-1,j)+phiold(i,j-1));
    end
end
%bottom side
for i=2:279
    for j=1
        phinew(i,j)=(1/3)*(phiold(i+1,j)+phiold(i-1,j)+phiold(i,j+1));
    end
end
%left side
for i=1
    for j=2:79
        phinew(i,j)=(1/3)*(phiold(i+1,j)+phiold(i,j+1)+phiold(i,j-1));
    end
end
%right side
for i=320
    for j=81:159
        phinew(i,j)= (1/3)*(phiold(i,j+1)+phiold(i-1,j)+phiold(i,j-1));
    end
end
%bottom left corner
for i=1
    for j=1
        phinew(i,j)=0.5*(phiold(i+1,j)+phiold(i,j+1));
    end
end
%top right corner
for i=320
    for j=160
        phinew(i,j)=0.5*(phiold(i-1,j)+phiold(i,j-1));
    end
end
err=99;
for i=1:nx
    for j=1:ny
        err=err+abs(phinew(i,j)-phiold(i,j));
    end
end
%swap old and new for next iteration after calculating error
for i=1:nx
    for j=1:ny
        phiold(i,j)=phinew(i,j);
end
end
end
surf(0:dx:(width-dx),0:dy:(length-dy),phinew);

2 个评论
显示无隐藏无

Walter Roberson 2016-4-7

You could improve the performance a lot by vectorizing your code.

Rena Berman 2017-1-20

(Answers Dev) Restored Question.

请先登录，再进行评论。

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

Answer 1

Walter Roberson 2016-4-21

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/277519-2d-heat-transfer-problem#answer_218879

You un-accepted my answer saying that it does not work for you, but you edited out the question, and do not say anything about what went wrong with the solution. We could make random guesses like, "Did you try using ZZ9 plural Z Alpha on line 42?", but it would be much easier to find a solution for you if you were to restore the question, post your current code, and describe the difficulties you are now encountering.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Walter Roberson 2016-4-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/277519-2d-heat-transfer-problem#answer_216820

Your code has an infinite loop. You have

err=1;
tol=1e-6;
while err>tol

that initial err is greater than that initial tol, so you will start the loop (as expected.)

Inside the loop, you never modify tol (which is not surprising.) Inside the loop, you modify err as follows:

err=99;
for i=1:nx
    for j=1:ny
        err=err+abs(phinew(i,j)-phiold(i,j));
    end
end

abs() anything is never negative, so your loop can only leave err the same (if the difference is 0) or increase it. Your final err cannot be any smaller than the "err=99" that you assigned. And 99 > tol so you are never going to leave the loop.

You should be using

err=0;

You should also vectorize your code. For example, your lines

for i=320
    for j=81:159
        phinew(i,j)= (1/3)*(phiold(i,j+1)+phiold(i-1,j)+phiold(i,j-1));
    end
end

can be rewritten as

i=320
j=81:159
phinew(i,j) = (1/3)*(phiold(i,j+1)+phiold(i-1,j)+phiold(i,j-1));

which will calculate for all of those j values at the same time.

Likewise, your code

for i=1:nx
    for j=1:ny
        phiold(i,j)=phinew(i,j);
end
end

is the same as

phiold(1:nx, 1:ny) = phinew(1:nx, 1:ny);

with no loop at all.

The natural impulse would be to further optimize that as

phiold = phinew;

but that is not equivalent. Look at the way you initialized:

phinew=zeros([ny nx]);
phiold=zeros([ny nx]);

so your arrays are ny by nx, but the work of the loop is to copy an nx by ny array because you hae "for i=1:nx" and you use nx as your first index, when it would be expected to be your second index for consistency with initializing to ny by nx. Chances are that there are other places in the code where you have made a similar mistake. The first index does correspond to y, so it is not clear why you would be using 1:nx for the first index.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 3

Michael Omodara 2020-3-27

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/277519-2d-heat-transfer-problem#answer_422392

I am trying to solve a 2D heat ttransfer in a box of grain. I have the code working as much as I would think but I could not plot the contour because of out put temperature has three coluumns. Can any body help me with this?

Also I want to extract the temperature values at the center of the grid and plot the values against time but I have no idea how to do that.

This is my code

%% Input Grain properties
MC= 15;                                         
k = 0.1409+0.00112*MC;                          
Cp   = 1000*(1.465+0.0356*MC);                  
rho  = 721;                                     
alpha = k/(Cp*rho);                             
%% Estimate Fourier number and grid space
Lx= 0.686;                                      
Ly= 0.267;                                      
% Lx= 1.371;                                     
% Ly= 2.136;                                    
% Lx=5.488;
% Ly=3.204;
% Lx= 2.744;                                     
% Ly= 2.136;                                    
Nx=5;                                            
Ny=5;                                            
dx=Lx/Nx;                                      
dy=Ly/Ny;                                        
x = linspace(0,Lx,Nx);
y = linspace(0,Ly,Ny);
%% Simulation parameter (difff simulation periods)
t_end=2592000;                                                           
dt=3600;                                       
iter=t_end/dt; 
                        
k1=alpha*dt/(dx^2);                             % Fourier number in x direction
k2=alpha*dt/(dy^2);                             % Fourier number in y direction
error = 6e9;
tol = 1e-4;
M=(dx^2+dy^2)/(2*alpha*dt);
%% Check stability of solution
if M>=4                                         % test the stability condition 
else 
fprintf('Error, the stability condition is not met\n Please return to "Inputs Section" and choose a "dt" less than %f \n',dx^2+dy^2/(8*alpha)) 
return 
end 
%% Initial and Boudary conditions 
                   
T_Top = 25;                                      
T_Left = 25;                                    
% T_Left = 5;                                    
T_Bottom = 25;                                   
T_Right = 25;                                                        
%% Face grid points
T=Inf(iter+1,Nx,Ny);                              
T(1,:,:)=5;                                       
T(:,1,2:end-1)=T_Top;                              
T(:,2:end-1,1)= T_Left;                                                  
T(:,Nx,2:end-1)=T_Bottom;                         
T(:,2:end-1,Ny)=T_Right;                          
%% Corner grid points
T(:,1,1) = 0.5*(T_Top + T_Left);                  
T(:,Nx,1) = 0.5*(T_Bottom + T_Left);              
T(:,1,Ny) = 0.5*(T_Top + T_Right);                
T(:,Nx,Ny) = 0.5*(T_Bottom + T_Right);            
%%Calculate Temperature Distribution using explicit method
[xx,yy,tt]= meshgrid (x,y,0:iter);
tic;
for t = 2:iter+1
 for i= 2:Nx-1
     for j = 2:Ny-1
         term1= T(t-1,i,j);
         term2= k1*(T(t-1,i-1,j)-2*T(t-1,i,j)+T(t-1,i+1,j));
         term3= k2*(T(t-1,i,j-1)-2*T(t-1,i,j)+T(t-1,i,j+1));
         T(t,i,j)= term1+term2+term3;
     end
 end
end
% time_counter=toc;
Tt=permute(T,[3,2,1]);
T_plot=Tt(3,3,:);
% 
% T_plot= T.';
% Tp=transpose(T);
% newarr = permute(sysarr,[3 2 1]);
%% Plot results
% To Fix: code doesn't work after here because the matrix is a 3x3 metrix
figure(1)
% [C,h]= contourf(xx,yy);
% colorbar;
colormap(jet);
% clabel(C,h);
% title(sprintf('2D Transient state heat conduction n No of iterations: %d(Time: %0.2fs) n Computation Time: %0.4f',counter-1,z*dt,time_counter));
% xlim([0 Lx]);
% ylim([0 Ly]);
% xlabel('Width of Stack (m)','FontWeight','bold','FontSize',16);
% ylabel('Height of Stack (m)','FontWeight','bold','FontSize',16);
% zlabel('Temperature ({\circC})','FontWeight','bold','FontSize',16);
% hold on
figure (2)
surf(tt,Tt);
xlim([0 Lx]);
ylim([0 Ly]);
zlim([5 25]);
xlabel('Width of Stack (m)','FontWeight','bold','FontSize',16);
ylabel('Height of Stack (m)','FontWeight','bold','FontSize',16);
zlabel('Temperature ({\circC})','FontWeight','bold','FontSize',16);

7 个评论
显示 5更早的评论隐藏 5更早的评论

Michael Omodara 2020-3-31

编辑：Walter Roberson 2020-4-1

Thank you for your help.

I modified the code to get the contour plots. And I want to make a movie of the hoourly plots.

However, it only plots the different countours without showing the movie

time_counter=toc;                                   
Tt=permute(T,[3,2,1]);                                                             
tall=0:iter;                                       
n = length(T_plot);                                
newValues = zeros(1,n);                             
for i =1:n
    
    newValues(i) = T_plot(:,:,i);                   
end
%% Plot results
figure(1)
[C,h]= contourf(xx(:,:,1),yy(:,:,1),Tt(:,:,1));     
colorbar;                                           
title_var_init = 1;                                 
colormap(jet);
clabel(C,h);
title(sprintf('2D Transient state heat conduction for %d days: time stop of 1 hour interval n Computation Time: %0.4f',title_var_init,time_counter));
xlim([0 Lx]);
ylim([0 Ly]);
xlabel('Width of Stack (m)','FontWeight','bold','FontSize',16);
ylabel('Height of Stack (m)','FontWeight','bold','FontSize',16);
zlabel('Temperature ({\circC})','FontWeight','bold','FontSize',16);
g = 24;                                            
g_end = iter+1;                                     
ff = 2;                                             
title_var = 2;
while g < g_end
figure(ff)                                          
[C,h]= contourf(xx(:,:,g),yy(:,:,g),Tt(:,:,g));     
g = g + 24;                                         
ff = ff + 1;                                        
colorbar;
colormap(jet);
clabel(C,h);
title(sprintf('2D Transient state heat conduction for %d days:  at 1 hour interval Computation time: %0.4f',title_var,time_counter));
title_var = title_var + 1;                          
xlim([0 Lx]);
ylim([0 Ly]);
xlabel('Width of Stack (m)','FontWeight','bold','FontSize',16);
ylabel('Height of Stack (m)','FontWeight','bold','FontSize',16);
zlabel('Temperature ({\circC})','FontWeight','bold','FontSize',16);
end
%% Animation
axis tight
set(gca,'nextplot','replacechildren');
for i = 1:n
    F(i) = getframe;
end
movie(F,8) % Play the movie twenty times

Can anyone show me what I am not doing right with the animation. Thank you

Walter Roberson 2020-4-1

time_counter=toc;                                   
Tt=permute(T,[3,2,1]);                                                             
tall=0:iter;                                       
n = length(T_plot);                                
newValues = zeros(1,n);                             
for i =1:n
    
    newValues(i) = T_plot(:,:,i);  %what is the purpose of this variable??                
    %% Plot results
    figure(1)
    [C,h]= contourf(xx(:,:,1),yy(:,:,1),Tt(:,:,1));     
    colorbar;                                           
    title_var_init = 1;                                 
    colormap(jet);
    clabel(C,h);
    title(sprintf('2D Transient state heat conduction for %d days: time stop of 1 hour interval n Computation Time: %0.4f',title_var_init,time_counter));
    xlim([0 Lx]);
    ylim([0 Ly]);
    xlabel('Width of Stack (m)','FontWeight','bold','FontSize',16);
    ylabel('Height of Stack (m)','FontWeight','bold','FontSize',16);
    zlabel('Temperature ({\circC})','FontWeight','bold','FontSize',16);
    g = 24;                                            
    g_end = iter+1;                                     
    ff = 2;                                             
    title_var = 2;
    while g < g_end
    figure(ff)                                          
    [C,h]= contourf(xx(:,:,g),yy(:,:,g),Tt(:,:,g));     
    g = g + 24;                                         
    ff = ff + 1;                                        
    colorbar;
    colormap(jet);
    clabel(C,h);
    title(sprintf('2D Transient state heat conduction for %d days:  at 1 hour interval Computation time: %0.4f',title_var,time_counter));
    title_var = title_var + 1;                          
    xlim([0 Lx]);
    ylim([0 Ly]);
    xlabel('Width of Stack (m)','FontWeight','bold','FontSize',16);
    ylabel('Height of Stack (m)','FontWeight','bold','FontSize',16);
    zlabel('Temperature ({\circC})','FontWeight','bold','FontSize',16);
    axis tight
    sleep(1);     %give time to paint screen
    F(i) = getframe;
end
%% Animation
movie(F,20) % Play the movie twenty times

Michael Omodara 2020-4-1

编辑：Walter Roberson 2020-5-9

Hi, I observed that the plot from

[C,h]= contourf(xx(:,:,g),yy(:,:,g),Tt(:,:,g));

interchanged the bottom side with the left side and the topside with the right.

However, when I used

contourf(x,y,T(:,:,i))

instead, the plot come out right but without contour labels.

What can I do to resolve this problem. i have spent weeks on the issue without success.

Thank you that works perfectly. I used pauce (1) in place of sleep (1) when it gave error that function sleep is not recognized.

Thankk you so much

Walter Roberson 2020-5-9

[C,h]= contourf(xx(:,:,g),yy(:,:,g),Tt(:,:,g));

Try

[C,h]= contourf(yy(:,:,g),xx(:,:,g),Tt(:,:,g).');

and earlier the corresponding

[C,h]= contourf(yy(:,:,1),xx(:,:,1),Tt(:,:,1).');

请先登录，再进行评论。

2D heat transfer problem

2 个评论
显示无隐藏无

回答（3 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论

另请参阅

类别

标签

Community Treasure Hunt

2D heat transfer problem

2 个评论 显示 无隐藏 无

回答（3 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

7 个评论 显示 5更早的评论隐藏 5更早的评论

另请参阅

类别

标签

Community Treasure Hunt

WeChat

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论