What is this error in my iterative program?

Im to solve for the temperature of the plate Tw, but to do so I need to know coefficient of convection h which is dependent on Tw. This is done via an iterative solution, only,

I'm not sure what the errors are that I am seeing here.

In addition to troubleshooting the error here, I am open to any advice of how the program can be improved (better method, clarity, formatting, etc.).

Thanks yall, I will respond to any questions as soon as possible.

flatplate_iteration.m

clear;
clc;
airProp=readtable('airProp_table.txt');
Tw= 30;
Ti=25;
ui=3;
L=.25;
A=L^2;
h=800/(Tw-Ti);
Nu=0;
k=0;
h_calc=Nu*k/L;
S=10;

while S>.001
i=1;
T1=airProp.Temperature(i);
Tf=(Tw-Ti)/2;
difference=Tf-T1;

while difference>=5
    i=i+1;
    T1=airProp.Temperature(i);
    difference=Tf-T1;
end
if difference==0
    d=airProp.Density(i);
    mu=airProp.DynamicViscosity(i);
    Pr=airProp.PrandtlNumber(i);
    k=airProp.ThermalConductivity(i);
else
    clear T2 d1 d2 mu1 mu2 Pr1 Pr2 k1 k2;
    T2=airProp.Temperature(i+1);
    d1=airProp.Density(i);
    d2=airProp.Density(i+1);
    d=d1+(Tf-T1)*(d2-d1)/(T2-T1);
    mu1=airProp.DynamicViscosity(i);
    mu2=airProp.DynamicViscosity(i+1);
    mu=mu1+(Tf-T1)*(mu2-mu1)/(T2-T1);
    Pr1=airProp.PrandtlNumber(i);
    Pr2=airProp.PrandtlNumber(i+1);
    Pr=Pr1+(Tf-T1)*(Pr2-Pr1)/(T2-T1);
    k1=airProp.ThermalConductivity(i);
    k2=airProp.ThermalConductivity(i+1);
    k=k1+(Tf-T1)*(k2-k1)/(T2-T1);
end
Re=d*ui*L/mu;
Nu=.664*Re^(1/2)*Pr^(1/3);
h_calc=Nu*k/L;
h=800/(Tw-Ti);
T_calc=Tw;
Tw=Tw+.1;
S=abs(h_calc-h);
end

[Tommy]

The error I get occurs when i equals 30. This line:

T2=airProp.Temperature(i+1);

fails because your table airProp only has 30 rows and you are trying to access row #31. This happens after many iterations. I see that you are updating Tw on every iteration of the (outer) loop. How can you be sure that this code:

i=1;
T1=airProp.Temperature(i);
Tf=(Tw-Ti)/2;
difference=Tf-T1;
while difference>=5
    i=i+1;
    T1=airProp.Temperature(i);
    difference=Tf-T1;
end

always leaves i with a value equal to 29 or less? (At the moment, it doesn't.)

[Rik]

You could also replace this entire loop by an array comparison.

Answer 1

[Rik]

Your code has some minor issues. I fixed some of them below. I also added a block of code where you need to add either an error or a mitigation. I also changed you copy-paste code to a function. Whenever you find yourself repeating code you should consider making it a local function. Make sure to give it a descriptive name and write documentation.

The most important thing you need to do is adding comments about what is happening. Then you can fill in the code that makes sure i has a valid value.

clear %you are setting all variables, so you should only need this during debugging
clc %you are not writing any text to the command window, so you only need this to clear the previous error during debugging
airProp=readtable('airProp_table.txt');
Tw= 30;
Ti=25;
ui=3;
L=.25;
A=L^2;
h=800/(Tw-Ti);
Nu=0;
k=0;
h_calc=Nu*k/L;
S=10;

while S>.001
    Tf=(Tw-Ti)/2;
    
    differences=Tf-airProp.Temperature;
    i=find(differences<5,1);
    difference=differences(i);
    if isempty(i)
        %(valid difference not found)
        %throw error or set some value for i that will work
    elseif i==numel(airProp.Temperature)
        %(later code will only work if difference==0)
        %throw error or set some value for i that will work
    end
    
    if difference==0
        d=airProp.Density(i);
        mu=airProp.DynamicViscosity(i);
        Pr=airProp.PrandtlNumber(i);
        k=airProp.ThermalConductivity(i);
    else
        T1=airProp.Temperature(i);
        T2=airProp.Temperature(i+1);
        
        d=FindMidPoint(airProp.Density,i,T1,T2,Tf);
        mu=FindMidPoint(airProp.DynamicViscosity,i,T1,T2,Tf);
        Pr=FindMidPoint(airProp.PrandtlNumber,i,T1,T2,Tf);
        k=FindMidPoint(airProp.ThermalConductivity,i,T1,T2,Tf);
    end
    Re=d*ui*L/mu;
    Nu=.664*Re^(1/2)*Pr^(1/3);
    h_calc=Nu*k/L;
    h=800/(Tw-Ti);
    T_calc=Tw;
    Tw=Tw+.1;
    S=abs(h_calc-h);
end
function out=FindMidPoint(tablevar,index,T1,T2,Tf)
%one-line description goes here
%longer description (if needed) goes here
%syntax description goes here (what does each variable mean)
val1=tablevar(index);
val2=tablevar(index+1);
out=val1+(Tf-T1)*(val2-val1)/(T2-T1);
end