MATLAB: Window Heat Transfer Problem

Hello I am currently working on a project for my heat transfer class where we need to calculate the heat transfer through a window for both summer and winter. Here is my current code

clear
AirTemperature = 100:50:500;
AirDynamicViscosityMatrix = [0.6924 1.0283 1.3289 1.488 1.893 2.075 2.286 2.484 2.671]; 
AirKinematicViscosityMatrix = [1.923 4.343 4.49 9.49 15.68 20.76 25.9 28.86 37.9];
AirDensityMatrix = [3.601 2.3675 1.7684 1.4128 1.774 0.998 0.8826 0.7833 0.7048];
KAirMatrix = [0.009246 0.013735 0.01809 0.02227 0.02624 0.03003 0.03365 0.03707 0.04038]; 
PrAirMatrix = [0.770 0.753 0.739 0.722 0.708 0.697 0.689 0.683 0.680];

WindowHeight = 1.8542;
WindowWidth = 0.762;
WindowArea= WindowHeight*WindowWidth;
WindowThickness = 0.00635;
Ti = 294.261111;
KWindow = 0.96;
OutsideWind = 6.7056;
Gravity = 9.81;
InsideAirTemperature = 294.26111;
i = 0;
NumberOfWindows = 11;
Rate = 0.09;

Season = input ('Please enter 1 for summer and 2 for winter ');
if (Season == 1)
    OutsideAirTemperature = 308.15;
else
    OutsideAirTemperature = 266.4833;
end

InsideTemperatureWindow = input ('Guess the temperature of the inside of the window (K) : ');
OutsideTemperatureWindow = input ('Guess the temperature of the outside of the window (K) : ');
TestTemperature1 =1;
TestTemperature2 = 1;

while abs(TestTemperature1-OutsideTemperatureWindow) > 0.001 || abs(TestTemperature2-InsideTemperatureWindow) > 0.001
    i = i + 1;
    disp (['Loop number ' num2str(i)])
    Tfo = (OutsideAirTemperature+OutsideTemperatureWindow)/2;
    Tfi = (InsideAirTemperature+InsideTemperatureWindow)/2;
    disp (['The value of Tfo is ' num2str(Tfo) ' K'])
    disp (['The value of Tfi is ' num2str(Tfi) ' K'])

      InsideAirKinematicViscosity = (interp1(AirTemperature, AirKinematicViscosityMatrix, Tfi))*(10^(-6)); 
      InsideAirDynamicViscosity = (interp1(AirTemperature, AirDynamicViscosityMatrix, Tfi))*(10^(-5)); 
      InsideAirDensity = interp1(AirTemperature, AirDensityMatrix, Tfi);
      InsideKAir = interp1(AirTemperature, KAirMatrix, Tfi);
      InsidePrAir = interp1(AirTemperature, PrAirMatrix, Tfi);
      OutsideAirKinematicViscosity = (interp1(AirTemperature, AirKinematicViscosityMatrix, Tfo))*(10^(-6)); 
      OutsideAirDynamicViscosity = (interp1(AirTemperature, AirDynamicViscosityMatrix, Tfo))*(10^(-5)); 
      OutsideAirDensity = interp1(AirTemperature, AirDensityMatrix, Tfo);
      OutsideKAir = interp1(AirTemperature, KAirMatrix, Tfo);
      OutsidePrAir = interp1(AirTemperature, PrAirMatrix, Tfo);

      OutsideReynold = OutsideWind*WindowHeight/OutsideAirKinematicViscosity;
      OutsideAirBeta = 1/Tfo;
      InsideAirBeta = 1/Tfi;

      OutsideGrashofPr = Gravity*OutsideAirBeta*(OutsideTemperatureWindow-OutsideAirTemperature)*(WindowHeight^3)*OutsidePrAir/(OutsideAirKinematicViscosity^2);
      disp (['The value of OutsideGrashofPr is ' num2str(OutsideGrashofPr)])
      InsideGrashofPr = Gravity*InsideAirBeta*(InsideTemperatureWindow-InsideAirTemperature)*(WindowHeight^3)*InsidePrAir/(InsideAirKinematicViscosity^2);
      disp (['The value of InsideGrashofPr is ' num2str(InsideGrashofPr)])
      if (OutsideGrashofPr < 10000)
          cOutside = 0.59;
          mOutside = 0.25;
      elseif (10000 <= OutsideGrashofPr < 1000000000)
          cOutside = 0.021;
          mOutside =0.4;
      else
          disp('GrPr value is outside of range');
          Stop
      end
      if (InsideGrashofPr < 10000)
          cInside = 0.59;
          mInside = 0.25;
      elseif (10000 <= OutsideGrashofPr < 1000000000)
          cInside = 0.021;
          mInside =0.4;
      else
          disp('GrPr value is outside of range');
          Stop
      end
      OutsideNusselt = cOutside*(OutsideGrashofPr^mOutside);
      disp (['The value of OutsideNusselt is ' num2str(OutsideNusselt)])
      InsideNusselt = cInside*(InsideGrashofPr^mInside);
      disp (['The value of InsideNusselt is ' num2str(InsideNusselt)])
      ho = OutsideNusselt*OutsideKAir/WindowHeight;
      hi = InsideNusselt*InsideKAir/WindowHeight;

      Sigma = 5.669*(10^-8);
      OutsideEmissivity = 0.96;
      InsideEmissivity = 0.96;

      R1 = (1/((ho*WindowArea)+(OutsideEmissivity*Sigma*(OutsideAirTemperature+OutsideTemperatureWindow)*((OutsideAirTemperature^2)+(OutsideTemperatureWindow^2))*WindowArea)));
      R2 = WindowThickness/(KWindow*WindowArea);
      R3 = (1/((hi*WindowArea)+(InsideEmissivity*Sigma*(InsideAirTemperature+InsideTemperatureWindow)*((InsideAirTemperature^2)+(InsideTemperatureWindow^2))*WindowArea)));
      q = (OutsideAirTemperature-InsideAirTemperature)/(R1+R2+R3);
      NewOutsideTemperatureWindow = (q*R1)+OutsideAirTemperature;
      NewInsideTemperatureWindow = (q*R3)+InsideAirTemperature;
      TestTemperature1 = OutsideTemperatureWindow;
      TestTemperature2 = InsideTemperatureWindow;
      OutsideTemperatureWindow = NewOutsideTemperatureWindow;
      InsideTemperatureWindow = NewInsideTemperatureWindow;
      q1 = (OutsideAirTemperature-OutsideTemperatureWindow)/R1;
      q2 = (OutsideTemperatureWindow-InsideTemperatureWindow)/R2;
      q3 = (InsideTemperatureWindow-InsideAirTemperature)/R3;
      q = (OutsideAirTemperature-InsideAirTemperature)/(R1+R2+R3);
      disp (['The value of q1 is ' num2str(q1) ' W'])
      disp (['The value of q2 is ' num2str(q2) ' W'])
      disp (['The value of q3 is ' num2str(q3) ' W'])
      disp (['The overall value of q is ' num2str(q) ' W'])
  end
  Rv = (OutsideAirTemperature-InsideAirTemperature)/(q/WindowArea);
  Cost = q*0.09*30*24*3600*NumberOfWindows/1000;
  disp (['The value of q1 is ' num2str(q1) ' W'])
  disp (['The value of q2 is ' num2str(q2) ' W'])
  disp (['The value of q3 is ' num2str(q3) ' W'])
  disp (['The overall value of q is ' num2str(q) ' W'])
  disp (['The money lost this month is $' num2str(Cost)])
  disp (['The value of Rv is ' num2str(Rv)])

Although I am getting it to run for summer, I should be getting all the values of q to be the same. I spent hours but cannot find out why the values are not the same. Also, when working with winter, it fails to complete the run. It says ??? Error using ==> interp1 at 166 The interpolation points XI should be real.

Error in ==> SinglePane at 42 InsideAirKinematicViscosity = (interp1(AirTemperature, AirKinematicViscosityMatrix, Tfi))*(10^(-6));

Any help would be appreciated. Thanks

Answer 1

[Walter Roberson]

if 10000 <= OutsideGrashofPr < 1000000000

is interpreted as

t = 10000 <= OutsideGrashofPr;  %a logical vector
if t < 1000000000   %always true since false = 0 and true = 1