2D heat transfer problem

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);

[Walter Roberson]

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

[Rena Berman]

(Answers Dev) Restored Question.

Answer 1

[Walter Roberson]

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.

Answer 2

[Walter Roberson]

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.

Answer 3

[Michael Omodara]

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);

[Walter Roberson]

% To Fix: code doesn't work after here because the matrix is a 3x3 metrix

I do not see any 3 x 3 matrix?

% [C,h]= contourf(xx,yy);

xx and yy are 5 x 5 x (iters+1), so 5 x 5 x 721, and xx and yy contain only x and y coordinates and not temperature values. If you add in a parameter into the call so you are plotting xx, yy, and temperature, then you would be trying to contourf() a 3D matrix representing 4D data (3 independent variables + 1 result). You cannot contourf() that.

Potentially you could isosurface():

for it = 5:5:25
    isosurface(xx,yy,tt,Tt,it);     %you do not need hold on for this
end
set(findobj(gca,'type','patch'),'FaceAlpha',0.5)

By the way, the trick to plotting the temperature at the middle is

plot(squeze(T_plot))

[Michael Omodara]

Thank you for you suggestions. However, I don't understand completely what you are suggesting for me to do. May be I was not clear enough. I want a conttour plots of temperature in all the 5X5 points at a specified.

I used this :

for it = 5:5:25
    isosurface(xx,yy,tt,Tt,it);     %you do not need hold on for this
end
set(findobj(gca,'type','patch'),'FaceAlpha',0.5)

but it is not outputing any plot.

Also what does squeze stands for as the code is given the following error

Unrecognized function or variable 'squeze'.
Error in Bagged_grain_Modified (line 105)
plot(squeze(T_plot))

Thanks for your help

[Walter Roberson]

Sorry, should have been squeeze instead of squeze

I want a conttour plots of temperature in all the 5X5 points at a specified.

for it = 1 : size(Tt,3)
    contourf(xx(:,:,it), yy(:,:,it), Tt(:,:,it));
    title(sprintf('iter = %d', it));
    drawnow();
end

[Michael Omodara]

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]

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]

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]

[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).');