Multiple For Loops to save data and then repeat again - Matrix Errors

Hi,

I am investigating the effect that certain parameters have on others, and in doing so I need to run the program at different variables, save the data then run again but I am getting lots of Matrix Errors. The code is very long but it is fairly clear of whats happening.

I would like to run the code at different CompOPR values to plot different Specific Thrust Values and Temperatures, save that data, then run again. This is to eventually plot the 4 sets of data on a graph (this part is fine but just getting the 4 sets of data and the errors I am running in to)

Any help would be massively appreciated :)

load('Thermo CW.mat')

P1 = 24990;%alt 40000ft
V1 =  238.32; %Mach 0.9
T1 = 220.79;


%Changing Parameters
CompStage = 9;
CompOPR = linspace(15,30,4)


for n = 1:length(CompOPR)
    
    
    CompSPR = CompOPR.^(1/CompStage);

    %Parameters
    CompEff = Values(5);
    CPair = [Values(4)];
    kAir = Values(28);
    CombdP = 0.02;
    TOPR = Values(24);
    TStages = 3;
    TSPR = TOPR^(1/TStages);
    kGas = [1.369863013698630 1.369863013698630 1.369863013698630 1.369863013698630];
    TEff = Values(9);
    CPgas = [1100 1100 1100 1100];


    
    T2s = T1 + ((V1^2)/(2*CPair));
    P2s = P1 * ((T2s/T1)^(kAir/(kAir-1)));
    T2 = T2s;
    P2 = P2s;

    CompT1 = T2;
    CompP1 = P2;
    CompT1s = T2;
    CompP2s = P2;

    CompT2s = CompT1*((CompSPR).^((kAir-1)/kAir));
    CompT2 = (((CompT2s - CompT1)/CompEff)+CompT1);
    ComP2s = CompP1 * ((CompT2s/CompT1).^(kAir/(kAir-1)));
    CompP2 = CompP1 * ((CompT2/CompT1).^(kAir/(kAir-1)));

    CompT3s = CompT2.*((CompSPR).^((kAir-1)/kAir));
    CompT3 = (((CompT3s - CompT2)/CompEff)+CompT2);
    ComP3s = CompP2 * ((CompT3s/CompT2).^(kAir/(kAir-1)));
    CompP3 = CompP2 * ((CompT3/CompT2).^(kAir/(kAir-1)));

    CompT4s = CompT3.*((CompSPR).^((kAir-1)/kAir));
    CompT4 = (((CompT4s - CompT3)/CompEff)+CompT3);
    CompP4s = CompP3 * ((CompT4s/CompT3).^(kAir/(kAir-1)));
    CompP4 = CompP3 * ((CompT4/CompT3).^(kAir/(kAir-1)));

    CompT5s = CompT4.*((CompSPR).^((kAir-1)/kAir));
    CompT5 = (((CompT5s - CompT4)/CompEff)+CompT4);
    CompP5s = CompP4 * ((CompT5s/CompT4).^(kAir/(kAir-1)));
    CompP5 = CompP4 * ((CompT5/CompT4).^(kAir/(kAir-1)));

    CompT6s = CompT5.*((CompSPR).^((kAir-1)/kAir));
    CompT6= (((CompT6s - CompT5)/CompEff)+CompT5);
    CompP6s = CompP5 * ((CompT6s/CompT5).^(kAir/(kAir-1)));
    CompP6 = CompP5 * ((CompT6/CompT5).^(kAir/(kAir-1)));

    CompT7s = CompT6.*((CompSPR).^((kAir-1)/kAir));
    CompT7= (((CompT7s - CompT6)/CompEff)+CompT6);
    CompP7s = CompP6 * ((CompT7s/CompT6).^(kAir/(kAir-1)));
    CompP7 = CompP6 * ((CompT7/CompT6).^(kAir/(kAir-1)));

    CompT8s = CompT7.*((CompSPR).^((kAir-1)/kAir));
    CompT8= (((CompT8s - CompT7)/CompEff)+CompT7);
    CompP8s = CompP7 * ((CompT8s/CompT7).^(kAir/(kAir-1)));
    CompP8 = CompP7 * ((CompT8/CompT7).^(kAir/(kAir-1)));

    CompT9s = CompT8.*((CompSPR).^((kAir-1)/kAir));
    CompT9= (((CompT9s - CompT8)/CompEff)+CompT8);
    CompP9s = CompP8 * ((CompT9s/CompT8).^(kAir/(kAir-1)));
    CompP9 = CompP8 * ((CompT9/CompT8).^(kAir/(kAir-1)));

    CompT10s = CompT9.*((CompSPR).^((kAir-1)/kAir));
    CompT10= (((CompT10s - CompT9)/CompEff)+CompT9);
    CompP10s = CompP9 * ((CompT10s/CompT9).^(kAir/(kAir-1)));
    CompP10 = CompP9 * ((CompT10/CompT9).^(kAir/(kAir-1)));

    T3s = T2s*((CompOPR).^((kAir-1)/kAir))
    T3 = CompT10;
    P3s = CompP10s;
    P3 = CompP10;

    CompWorks = CPair *(T3s-T2s);
    CompWork = CPair *(T3 - T2);

    T4 = linspace(1500,1900,20);


    for k = 1:length(T4)
        T4s = T4;
        P4s = P3;
        P4 = P4s * (1-CombdP);

        TT1 = T4;
        TT1s = T4s;
        TP1s = P4;
        TP1 = P4;

        TT2s = TT1 * ((TSPR^((kGas-1)/kGas)));
        TT2 = (((TT2s - TT1)*TEff)+TT1);
        TP2s = TP1 *((TT2s/TT1)^(kGas/(kGas-1)));
        TP2 = TP1 * ((TT2/TT1)^(kGas/(kGas-1)));

        TT3s = TT2 * ((TSPR^((kGas-1)/kGas)));
        TT3 = (((TT3s - TT2)*TEff)+TT2);
        TP3s = TP2 *((TT3s/TT2)^(kGas/(kGas-1)));
        TP3 = TP2 * ((TT3/TT2)^(kGas/(kGas-1)));

        TT4s = TT3 * ((TSPR^((kGas-1)/kGas)));
        TT4 = (((TT4s - TT3)*TEff)+TT3);
        TP4s = TP3 *((TT4s/TT3)^(kGas/(kGas-1)));
        TP4 = TP3 * ((TT4/TT3)^(kGas/(kGas-1)));

        AuxPwr = (0.08) * CompWork;

        TrueTurbOutT = AuxPwr + CompWork;

        TurbWorkMechLoss = TrueTurbOutT ./ (0.98);


        T5s = (-(TurbWorkMechLoss / CPgas)) + T4;
        T5 = T5s;
        P5 = P4 * ((T5/T4)^(kGas/(kGas - 1)));
        P5s = P5;

        P6s = P1;
        P6 = P1;
        T6s = T5s .*((P6s ./ P5s).^((kGas - 1)./ kGas));
        T6 = T5 .* ((P6 ./ P5).^((kGas - 1)./ kGas));

        V6 = (0.98*(((T5-T6)*(2*CPgas)).^(1/2)).^(2)).^(1/2);

        mDotair = Values(31);

        CPavg = (CPair + CPgas)/2;

        Qin = mDotair * CPavg *(T4-T3);

        LHV = Values(17)*(10^6);

        mDotfuel = Qin/LHV;

        Thrust = ((mDotair + mDotfuel).*V6)-(mDotair.*V1);

        mDotexhaust = mDotair + mDotfuel;

        SpecificThrust = Thrust/(mDotair);


    end


    plot(T4,mDotfuel)
    xlabel('Mass Flow Rate of Fuel, kg / s')
    ylabel('Specific Thrust, N s / kg')

end



%% These are the errors I am getting
Matrix dimensions must agree.

Error in Untitled2 (line 137)
        T6s = T5s .*((P6s ./ P5s).^((kGas - 1)./ kGas));
 

Answer 1

[SandeepKumar R]

T6s is the expression that has this error. Split it up into sub expressions and evaluate. When the matrix dimension error pops up for that sub expression, inspect the size of the matrices ( size(A) ). This will show where incorrect matrix operation is taking place.

[Matthew Hole]

Hi SandeepKumar, this helped solve that error thank you, and I'm applying it to the similar ones that come up. In your opinion is what I'm attempt possible, or is there a way to streamline the process?

Thanks

Answer 2

[SandeepKumar R]

In my experirence, careful bookeeping of the matrix dimensions in your mind while multiplying them is the only fix. If at all you're not sure of the dimension you can look at the workspace.

[Matthew Hole]

I've managed to fix the errors and the program runs, but its not functioning as intended. Code is shown below.

I want it to run for each value of the first for loop, which runs for all the values of the second, saves the values, then runs again until the length of CompOPR is complete, ie the 4 times. Sorry if this doesnt make sense

load('Thermo CW.mat')

P1 = 24990;%alt 40000ft
V1 =  238.32; %Mach 0.9
T1 = 220.79;


%Changing Parameters
CompStage = 9;
CompOPR = linspace(15,30,4);


for n = 1:length(CompOPR)
    
    
    CompSPR = CompOPR.^(1/CompStage);

    %Parameters
    CompEff = Values(5);
    CPair = [1005 1005 1005 1005];
    kAir = Values(28);
    CombdP = 0.02;
    TOPR = Values(24);
    TStages = 3;
    TSPR = TOPR^(1/TStages);
    kGas = [1.369863013698630];
    TEff = Values(9);
    CPgas = [1100 1100 1100 1100];


    
    T2s = T1 + ((V1^2)./(2.*CPair));
    P2s = P1 * ((T2s/T1).^(kAir/(kAir-1)));
    T2 = T2s;
    P2 = P2s;

    CompT1 = T2;
    CompP1 = P2;
    CompT1s = T2;
    CompP2s = P2;

    CompT2s = CompT1.*((CompSPR).^((kAir-1)/kAir));
    CompT2 = (((CompT2s - CompT1)/CompEff)+CompT1);
    ComP2s = CompP1 * ((CompT2s/CompT1).^(kAir/(kAir-1)));
    CompP2 = CompP1 * ((CompT2/CompT1).^(kAir/(kAir-1)));

    CompT3s = CompT2.*((CompSPR).^((kAir-1)/kAir));
    CompT3 = (((CompT3s - CompT2)/CompEff)+CompT2);
    ComP3s = CompP2 * ((CompT3s/CompT2).^(kAir/(kAir-1)));
    CompP3 = CompP2 * ((CompT3/CompT2).^(kAir/(kAir-1)));

    CompT4s = CompT3.*((CompSPR).^((kAir-1)/kAir));
    CompT4 = (((CompT4s - CompT3)/CompEff)+CompT3);
    CompP4s = CompP3 * ((CompT4s/CompT3).^(kAir/(kAir-1)));
    CompP4 = CompP3 * ((CompT4/CompT3).^(kAir/(kAir-1)));

    CompT5s = CompT4.*((CompSPR).^((kAir-1)/kAir));
    CompT5 = (((CompT5s - CompT4)/CompEff)+CompT4);
    CompP5s = CompP4 * ((CompT5s/CompT4).^(kAir/(kAir-1)));
    CompP5 = CompP4 * ((CompT5/CompT4).^(kAir/(kAir-1)));

    CompT6s = CompT5.*((CompSPR).^((kAir-1)/kAir));
    CompT6= (((CompT6s - CompT5)/CompEff)+CompT5);
    CompP6s = CompP5 * ((CompT6s/CompT5).^(kAir/(kAir-1)));
    CompP6 = CompP5 * ((CompT6/CompT5).^(kAir/(kAir-1)));

    CompT7s = CompT6.*((CompSPR).^((kAir-1)/kAir));
    CompT7= (((CompT7s - CompT6)/CompEff)+CompT6);
    CompP7s = CompP6 * ((CompT7s/CompT6).^(kAir/(kAir-1)));
    CompP7 = CompP6 * ((CompT7/CompT6).^(kAir/(kAir-1)));

    CompT8s = CompT7.*((CompSPR).^((kAir-1)/kAir));
    CompT8= (((CompT8s - CompT7)/CompEff)+CompT7);
    CompP8s = CompP7 * ((CompT8s/CompT7).^(kAir/(kAir-1)));
    CompP8 = CompP7 * ((CompT8/CompT7).^(kAir/(kAir-1)));

    CompT9s = CompT8.*((CompSPR).^((kAir-1)/kAir));
    CompT9= (((CompT9s - CompT8)/CompEff)+CompT8);
    CompP9s = CompP8 * ((CompT9s/CompT8).^(kAir/(kAir-1)));
    CompP9 = CompP8 * ((CompT9/CompT8).^(kAir/(kAir-1)));

    CompT10s = CompT9.*((CompSPR).^((kAir-1)/kAir));
    CompT10= (((CompT10s - CompT9)/CompEff)+CompT9);
    CompP10s = CompP9 * ((CompT10s/CompT9).^(kAir/(kAir-1)));
    CompP10 = CompP9 * ((CompT10/CompT9).^(kAir/(kAir-1)));

    T3s = T2s.*((CompOPR).^((kAir-1)/kAir))
    T3 = CompT10;
    P3s = CompP10s;
    P3 = CompP10;

    CompWorks = CPair .*(T3s-T2s);
    CompWork = CPair .*(T3 - T2);

    T4 = linspace(1500,1900,20);


    for k = 1:length(T4)
        T4s = T4;
        P4s = P3;
        P4 = P4s * (1-CombdP);

        TT1 = T4;
        TT1s = T4s;
        TP1s = P4;
        TP1 = P4;

        TT2s = TT1 * ((TSPR^((kGas-1)/kGas)));
        TT2 = (((TT2s - TT1)*TEff)+TT1);
        TP2s = TP1 *((TT2s/TT1)^(kGas/(kGas-1)));
        TP2 = TP1 * ((TT2/TT1)^(kGas/(kGas-1)));

        TT3s = TT2 * ((TSPR^((kGas-1)/kGas)));
        TT3 = (((TT3s - TT2)*TEff)+TT2);
        TP3s = TP2 *((TT3s/TT2)^(kGas/(kGas-1)));
        TP3 = TP2 * ((TT3/TT2)^(kGas/(kGas-1)));

        TT4s = TT3 * ((TSPR^((kGas-1)/kGas)));
        TT4 = (((TT4s - TT3)*TEff)+TT3);
        TP4s = TP3 *((TT4s/TT3)^(kGas/(kGas-1)));
        TP4 = TP3 * ((TT4/TT3)^(kGas/(kGas-1)));

        AuxPwr = (0.08) * CompWork;

        TrueTurbOutT = AuxPwr + CompWork;

        TurbWorkMechLoss = TrueTurbOutT ./ (0.98);
        
        T51 = reshape(T5,1,20)


        T5s = (-(TurbWorkMechLoss / CPgas)) + T4;
       
        P5 = P4 .* ((T51/T4).^(kGas./(kGas - 1)));
        P5s = P5;

        P6s = P1;
        P6 = P1;
        
        T6s1 = size(((P6s ./ P5s) .^((kGas - 1) ./ kGas)))
        
        T5s1 = reshape(T5s,20,1)
        T5 = T5s1
        size(T5s1)
        
        T6s = (T5s1 .* ((P6s ./ P5s) .^((kGas - 1) ./ kGas)));
        T6 = T5 .* ((P6 ./ P5).^((kGas - 1)./ kGas));

        V6 = (0.98.*(((T5-T6).*(2.*CPgas)).^(1/2)).^(2)).^(1/2);

        mDotair = Values(31);

        CPavg = (CPair + CPgas)/2;
        
        T31 = reshape(T3,4,1)
        
        Qin = mDotair * CPavg *(T4-T31);

        LHV = Values(17)*(10^6);

        mDotfuel = Qin/LHV;
        
        
        mOut = (mDotair + mDotfuel)
        size(mOut)
        size(V6)
        momOut = (mOut)*V6
        momIn = (mDotair.*V1)
        
       
        
        
        size(mDotair)
        
        size(V1)
       

        Thrust = momOut-momIn;

        mDotexhaust = mDotair + mDotfuel;

        SpecificThrust = Thrust/(mDotair);


    end



end


    plot(T4,mDotfuel)
    xlabel('Mass Flow Rate of Fuel, kg / s')
    ylabel('Specific Thrust, N s / kg')