2D steady heat transfer output returning NaN

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

Answer 1

[Divyajyoti Nayak]

Hi @Jai Harnas,

I ran your code for heat transfer and did get the NaN values like you mentioned. On debugging I found that the first occurrence of NaN happens in the 13th iteration of the while loop. Looking further into the exact point where ‘T(j,i)’ becomes NaN, I found the cause for it.

Some of the values of ‘T’ are increasing exponentially and after the 12th iteration, they become so high that MATLAB treats them as infinity (inf). Any further algebraic operations done using both operands as inf gives NaN.

inf - inf

ans = NaN

Hope this helps.