2D heat transfer problem
14 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(3 个)
Walter Roberson
2016-4-21
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 个评论
Walter Roberson
2016-4-7
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 个评论
Michael Omodara
2020-3-27
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 个评论
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).');
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Geometry and Mesh 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!