MATLAB: Complex results instead of real for heat transfer calculation

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!

Answer 1

[Walter Roberson]

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

[Hemdeep]

Thanks for the help.

I used to display function to find out where the value was becoming negative.