Info

This question is closed. Reopen it to edit or answer.

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

1 view (last 30 days)
Matthew Hole
Matthew Hole on 6 Mar 2019
Closed: MATLAB Answer Bot on 20 Aug 2021
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));

Answers (2)

SandeepKumar R
SandeepKumar R on 6 Mar 2019
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.
  1 Comment
Matthew Hole
Matthew Hole on 6 Mar 2019
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

SandeepKumar R
SandeepKumar R on 6 Mar 2019
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.
  1 Comment
Matthew Hole
Matthew Hole on 6 Mar 2019
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')

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by