2D steady heat transfer output returning NaN

4 visualizaciones (últimos 30 días)
Jai Harnas
Jai Harnas el 20 de Sept. de 2021
Hi, i am currently writing code for a 2D steady heat transfer problem. The problem is a metal plate with a heated coil in the middle and a wet timber block being dried on the end of the plate. I have set up the matrices and all the equations for the gemoetry with dx=dy. However, when i run the counting k variable shows that the code has been run 111 times but each cell apart from the coil is returning NaN and im not sure why. Im aware the if statements is a messy way to write the code, but its how ive understood it and i have done for loops to separate the gemoetry into nodes around the coil, nodes in the plate passed the coil and nodes in the timber block. Im essentially after if someone can tell me why all my points in the matrix is returning NaN if the code is being run 111 times.
I have also tried to include it all in one for loop, keeping the if statements, and the same problem occurs. Thank you, here is my code:
clear all;
clc;
nx=200;
ny=91;
w=0.2;
H=0.09;
dx=0.001;
dy=0.001;
Tinf=20;
h=7.7;
kt=0.6;
kp=14.4;
Tcoil=300;
x=linspace(0,dx,w);
y=linspace(0,dy,H);
T=zeros(ny,nx);
%set coil temperature
T(87,1:125)=Tcoil;
T(85,1:125)=Tcoil;
T(86,1:125)=Tcoil;
tol=1e-6;
error=1;
k=0;
while error > tol
k=k+1;
Told=T;
%following for loop is for the nodes surrounding the coil and
%insulation
for i=1:125
for j=81:ny
if j==81 && i==1 %plate top exterior node next to insulation
T(j,i)=h*dx*(Tinf-T(j,i))+kp*(T(j+1,i)-T(j,i))+kp*0.5*(T(j,i+1)-T(j,i));
elseif i==1&&j>81&&j<85 || i==1&&j>87&&j<ny %interior nodes next to insulation
T(j,i)=kp*(T(j+1,i)-T(j,i))+kp*(T(j-1,i)-T(j,i))+kp*(T(j,i+1)-T(j,i));
elseif i==1&& j==ny %lower exterior node next to insulation
T(j,i)=h*dx*(Tinf-T(j,i))+kp*(T(j-1,i)-T(j,i))+kp*0.5*(T(j,i+1)-T(j,i));
elseif i>1 && j==81
%top exterior nodes without insulation
T(j,i)=h*dx*(Tinf-T(j,i))+kp*(T(j+1,i)-T(j,i))+kp*0.5*(T(j,i+1)-T(j,i))+kp*0.5*(T(j,i-1)-T(j,i));
elseif i>1 && j==ny %lower exterior nodes without insulation
T(j,i)=h*dx*(Tinf-T(j,i))+kp*(T(j-1,i)-T(j,i))+kp*0.5*(T(j,i+1)-T(j,i))+kp*0.5*(T(j,i-1)-T(j,i));
elseif i>1 && j>81&&j<85 || i>1&&j>87&&j<ny %interior nodes without insulation
T(j,i)=kp*(T(j+1,i)-T(j,i))+kp*(T(j-1,i)-T(j,i))+kp*(T(j,i+1)-T(j,i))+kp*(T(j,i-1)-T(j,i));
end
end
end
%following for loop is for the plate nodes passed the coil (including
%inbetween nodes)
for i=126:nx
for j=81:ny
if j==81 && i<150 %upper exterior nodes passed coil
T(j,i)=h*dx*(Tinf-T(j,i))+kp*(T(j+1,i)-T(j,i))+kp*0.5*(T(j,i+1)-T(j,i))+kp*0.5*(T(j,i-1)-T(j,i));
elseif j>81 && j<ny && i>125 && i<nx %interior nodes passed coil
T(j,i)=kp*(T(j+1,i)-T(j,i))+kp*(T(j-1,i)-T(j,i))+kp*(T(j,i+1)-T(j,i))+kp*(T(j,i-1)-T(j,i));
elseif j==ny && i<nx %lower exterior nodes passed coil
T(j,i)=h*dx*(Tinf-T(j,i))+kp*(T(j-1,i)-T(j,i))+kp*0.5*(T(j,i+1)-T(j,i))+kp*0.5*(T(j,i-1)-T(j,i));
elseif j>81 && j<ny && i==nx %side exterior nodes
T(j,i)=h*dy*(Tinf-T(j,i))+kp*0.5*(T(j-1,i)-T(j,i))+kp*(T(j,i-1)-T(j,i))+kp*0.5*(T(j+1,i)-T(j,i));
elseif j==ny && i==nx %bottom right corner node
T(j,i)=h*(dy/2)*(Tinf-T(j,i))+h*(dx/2)*(Tinf-T(j,i))+kp*0.5*(T(j,i-1)-T(j,i))+kp*0.5*(T(j-1,i)-T(j,i));
elseif j==81 && i==150 %inbetween elbow node
T(j,i)=kp*0.5*(T(j,i+1)-T(j,i))+kp*(T(j-1,i)-T(j,i))+kp*0.5*(T(j,i-1)-T(j,i))+h*(dy/2)*(Tinf-T(j,i))+h*(dx/2)*(Tinf-T(j,i))+kt*0.5*(T(j+1,i)-T(j,i))+kt*0.5*(T(j,i+1)-T(j,i));
elseif j==81 && i>150 && i<nx %inbetween interior node
T(j,i)=kp*0.5*(T(j,i+1)-T(j,i))+kp*(T(j-1,i)-T(j,i))+kp*0.5*(T(j,i-1)-T(j,i))+kt*0.5*(T(j,i-1)-T(j,i))+kt*(T(j-1,i)-T(j,i))+kt*0.5*(T(j,i+1)-T(j,i));
elseif j==81 && i==nx %inbetween exterior node
T(j,i)=h*dy*(Tinf-T(j,i))+kp*0.5*(T(j+1,i)-T(j,i))+kp*0.5*(T(j,i-1)-T(j,i))+kt*0.5*(T(j,i-1)-T(j,i))+kt*0.5*(T(j-1,i)-T(j,i));
end
end
end
%following for loop is for the timber nodes
for i=150:nx
for j=1:81
if i==150 && j>1 %left exterior node (timber)
T(j,i)=kt*(T(j,i+1)-T(j,i))+kt*0.5*(T(j-1,i)-T(j,i))+h*dy*(Tinf-T(j,i))+kt*0.5*(T(j+1,i)-T(j,i));
elseif i==150 && j==1 %left corner node (timber)
T(j,i)=kt*0.5*(T(j,i+1)-T(j,i))+kt*0.5*(T(j+1,i)-T(j,i))+h*(dy/2)*(Tinf-T(j,i))+h*(dx/2)*(Tinf-T(j,i));
elseif i==nx && j==1 %right corner node (timber)
T(j,i)=h*(dy/2)*(Tinf-T(j,i))+kt*0.5*(T(j+1,i)-T(j,i))+kt*0.5*(T(j,i-1)-T(j,i))+h*(dx/2)*(Tinf-T(j,i));
elseif i>150 && i<nx && j==1 %upper exterior node (timber)
T(j,i)=h*dx*(Tinf-T(j,i))+kt*(T(j+1,i)-T(j,i))+kt*0.5*(T(j,i+1)-T(j,i))+kt*0.5*(T(j,i-1)-T(j,i));
elseif i==nx && j>1 %right exterior node (timber)
T(j,i)=h*dy*(Tinf-T(j,i))+kt*0.5*(T(j-1,i)-T(j,i))+kt*(T(j,i-1)-T(j,i))+kt*0.5*(T(j+1,i)-T(j,i));
elseif i>150 && i<nx && j>1 %interior node (timber)
T(j,i)=kt*(T(j+1,i)-T(j,i))+kt*(T(j-1,i)-T(j,i))+kt*(T(j,i+1)-T(j,i))+kt*(T(j,i-1)-T(j,i));
end
end
end
error=max(max(abs(Told-T)));
end

Respuestas (0)

Categorías

Más información sobre Unit Conversions en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by