MATLAB: MATLAB: Complex results instead of real for heat transfer calculation

complex resultsengineeringheat transferwindow pane

Hi I am writing a program to calculate the heat transfer across a window pane. Here is the 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.78;
OutsideWind = 6.7056;
Gravity = 9.81;
InsideAirTemperature = 294.26111;
i = 1;
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) : ');
NewOutsideTemperatureWindow =1;
NewInsideTemperatureWindow = 1;
while abs(NewOutsideTemperatureWindow-OutsideTemperatureWindow) > 0.001 && abs(NewInsideTemperatureWindow-InsideTemperatureWindow) > 0.001
  i = i + 1;   
  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);
  InsideGrashofPr = Gravity*InsideAirBeta*(InsideTemperatureWindow-InsideAirTemperature)*(WindowHeight^3)*InsidePrAir/(InsideAirKinematicViscosity^2);
  if (OutsideGrashofPr < 10000)
      c = 0.59;
      m = 0.25;
  elseif (10000 <= OutsideGrashofPr < 1000000000)
      c = 0.021;
      m =0.4;
  else
      disp('GrPr value is outside of range');
      Stop
  end
  OutsideNusselt = c*(OutsideGrashofPr^m);
  InsideNusselt = c*(InsideGrashofPr^m);
  ho = OutsideNusselt*OutsideKAir/WindowHeight;
  hi = InsideNusselt*InsideKAir/WindowHeight;
  Sigma = 0;
  OutsideEmissivity = 0;
  InsideEmissivity = 10000000;
  R1 = (1/((ho*WindowArea)+(OutsideEmissivity*Sigma*(OutsideAirTemperature+OutsideTemperatureWindow)*((OutsideAirTemperature^2)+(OutsideTemperatureWindow^2))*WindowArea)));
  R2 = WindowThickness/(KWindow*WindowArea);
  R3 = (1/((hi*WindowArea)+(OutsideEmissivity*Sigma*(InsideAirTemperature+InsideTemperatureWindow)*((InsideAirTemperature^2)-(InsideTemperatureWindow^2))*WindowArea)));
  q = (OutsideAirTemperature-InsideAirTemperature)/(R1+R2+R3);
  NewOutsideTemperatureWindow = (OutsideAirTemperature-(q*R1));
  NewInsideTemperatureWindow = (q*R3)+InsideAirTemperature;
  OutsideTemperatureWindow = NewOutsideTemperatureWindow;
  InsideTemperatureWindow = NewInsideTemperatureWindow;
  NewInsideTemperatureWindow = 1;
  NewOutsideTemperatureWindow = 1;
end
q1 = (OutsideAirTemperature-OutsideTemperatureWindow)/R1;
q2 = (OutsideTemperatureWindow-InsideTemperatureWindow)/R2;
q3 = (InsideTemperatureWindow-InsideAirTemperature)/R3;
Rv = (OutsideAirTemperature-InsideAirTemperature)/(q/WindowArea);
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 value of q4 is ' num2str(q) ' W'])
Unfortunately, I am getting a complex number during the second iteration as shown below
Please enter 1 for summer and 2 for winter 2
Guess the temperature of the inside of the window (K) : 300
Guess the temperature of the outside of the window (K) : 200
The value of Tfo is 233.2416 K
The value of Tfi is 297.1306 K
The value of Tfo is 270.3712-2.388748i K
The value of Tfi is 284.423-2.349841i K
??? Error using ==> interp1 at 166
The interpolation points XI should be real.
Error in ==> SinglePane at 39
    InsideAirKinematicViscosity = (interp1(AirTemperature,
    AirKinematicViscosityMatrix, Tfi))*(10^(-6));
Any help would be great. Thanks!

Best Answer

The first place I would look would be at your statements,
OutsideNusselt = c*(OutsideGrashofPr^m);
InsideNusselt = c*(InsideGrashofPr^m);
Your m is floating point (0.4) so this will be calculated using logs. If the value being taken to a power was negative, the result will be complex.
Remember that x^(0.4) is not defined to be the same as (x^2)^(1/5) where the latter would come out positive for negative x.
The usual heat transfer equations are order 16 polynomials, and commonly when they are worked through, 12 of the roots are imaginary and four are real. The real solutions then get brought in to practicality by noticing that three of the four real solutions involve temperatures below Absolute 0. See this previous discussion
Related Question