MATLAB: Loop to create matrix with 2 inputs and 2 outputs

>2 variablesfor loopmatrixvector

Hello,

I recently had an assignment to create a code and got it right, but afterwards it was recommended to me that it would have been easier to do it a different way. The way the code below is set up was to be more like a calculator. You input the 2 values and it displays 2 outputs. This was very tedious to do and a loop would have been much quicker.

For my situation I have 2 inputs and 2 outputs that are a function of the 2 inputs. I want to have the inputs change through a loop and save the output variables to a matrix.

The value of TINT needs to be set to range from 1100 to 1800 at intervals of 100 while the OPR needs to be 4 to 30 increments of 2. I have a basic understanding of how to do this if it was only one variable changing. I am not sure how to do 2 variables.

%------------------------------------------------------

%Given Values                                         |

Ra = 0.287;   %kJ/kg/k                                |

ga = 1.4;     %Gamma Air                              |

gg = 1.3333;   %Gamma gas                              |

Nd = 0.93;    %Efficiency of Diffuser                 |

Nc = 0.87;    %Efficiency of Compressor               |

Nb = 0.98;    %Efficiency of Combustion Chamber       |

Nm = 0.99;    %Efficiency of Mechanical               |

Nt = 0.90;    %Efficiency of Turbines                 |

Nn = 0.95;    %Efficiency of Nozzle                   |

Ma = 0.84;    %Mach                                   |

Alt = 5;      %km                                     |

%OPR = 9;      %Overal Pressure Ratio                 |   

prompt = 'Overall Pressure Ratio? ';  %               |

OPR = input(prompt);                  %               |
PRClp = sqrt(OPR); %Pressure Ratio over LP Compressor |
PRChp = sqrt(OPR); %Pressure Ratio over HP Compressor |
deltap0 = 0.04;   % 4% of delievery a loss            |

%-----------------------------------------------------|



%From Handout 5 based on Altitude of 5km--------------|

Ta = 255.7;    %K                                     |

Pa = 54.05;    %kPa                                   |

rho = 0.736;  %kg/m^3                                 |

u = 1.63;     %N*s/m^2                                |

%-----------------------------------------------------|
Va = Ma * sqrt(ga*Ra*Ta*1000);   %in m/s

Cpoa = (ga*Ra)/(ga-1);          %in kJ/kg/K



T01 = Ta*(1 + ((ga-1)/2)*Ma^2);
P01 = Pa*(1 + Nd*((ga-1)/2)*(Ma^2))^(ga/(ga-1));
P02 = P01*PRClp;
T02 = T01 + ((T01*((PRClp)^((ga-1)/ga)-1)/Nc));
Wclp = Cpoa*(T01-T02);
P03 = P02*PRChp;
T03 = T02 + ((T02*((PRChp)^((ga-1)/ga)-1)/Nc));
Wchp = -Cpoa*(T03-T02);
Tr = T03;
   % Tp

    prompt = 'Temperature of Products? ';
    Tp = input(prompt);
    T04 = Tp;
      if (Tp <= 1600 && Tp>=400)
              COAp = 299180; COBp = 37.85; COCp = -4571.9; CO2Ap = 56835; 
              CO2Bp = 66.27; CO2Cp = -11634; H2OAp = 88923; H2OBp = 49.36; 
              H2OCp = -7940.8; N2Ap = 31317; N2Bp = 37.46; N2Cp = -4559.3;
              O2Ap = 43388; O2Bp = 42.27; O2Cp = -6635.4;
          %Tp > 1600K

          elseif (Tp > 1600);
              COAp = 309070; COBp = 39.29; COCp = -6201; CO2Ap = 93048; 
              CO2Bp = 68.58; CO2Cp = -16979; H2OAp = 154670; H2OBp = 60.43; 
              H2OCp = -19212; N2Ap = 44639; N2Bp = 39.32; N2Cp = -6753.4; 
              O2Ap = 127010; O2Bp = 46.25; O2Cp = -18798;
      end
MC = 14.4;
MH = 24.9;
MO = 0;
Hrp = -8561991.6;
Ycc = MC + (MH/4) + (MO/2);
Ymin = Ycc - (MC/2); 
Mfuel = MC*12 + MH + MO*16;
Mair = 28.97;
%Find f

    %based on Tr < 1600K

    COAr = 299180; COBr = 37.85; COCr = -4571.9; CO2Ar = 56835;
    CO2Br = 66.27; CO2Cr = -11634; H2OAr = 88923; H2OBr = 49.36;
    H2OCr = -7940.8; N2Ar = 31317; N2Br = 37.46; N2Cr = -4559.3;
    O2Ar = 43388; O2Br = 42.27; O2Cr = -6635.4;
      %based on Excess Air

      N1 = 0;
      N2 = MC;
      N3 = (MH/2); 
      HbarCO = 0;
      HbarCO2 = (CO2Ap + CO2Bp*Tp + CO2Cp*log(Tp)) - (CO2Ar + CO2Br*Tr + CO2Cr*log(Tr));
      HbarH2O = (H2OAp + H2OBp*Tp + H2OCp*log(Tp)) - (H2OAr + H2OBr*Tr + H2OCr*log(Tr));
      HbarO2 = (O2Ap + O2Bp*Tp + O2Cp*log(Tp)) - (O2Ar + O2Br*Tr + O2Cr*log(Tr));
      HbarN2 = (N2Ap + N2Bp*Tp + N2Cp*log(Tp)) - (N2Ar + N2Br*Tr + N2Cr*log(Tr));
      syms Y
      eqn = Hrp + (N1 * 282800) + (N2 * HbarCO2) + (N3 * HbarH2O) + ((Y - Ycc)* HbarO2) + ((3.76 * Y)* HbarN2)==0;
      Y = solve(eqn,Y);
      %Finding f   

      fideal= ((1/(4.76*Y)))*((Mfuel)/(Mair));
      f = fideal/0.98;  %Getting f actual from ideal and efficiency

     P04 = (1- deltap0)*P03;
     Cp0g = (gg*Ra)/(gg-1);      %in kJ/kg/K
     T05 = T04 +((Wchp)/(Nm*(1+0.0168170)*Cp0g));  %K

     P05 = (P04)*(1-((1-(T05/T04))/Nt))^(gg/(gg-1));  %kPa

     T06 = T05 + (Wclp/(Nm*(1+f)*Cp0g));
     %double(T06);

     P06 = P05*(1-((1-(T06/T05))/Nt))^(gg/(gg-1));  %Page 26!!!!

    double(P06);
    %Return to the particular converging nozzle of the turbojet problem

    P0i = P06;
    T0i = T06;
    T0e = T06;
    PstarovP06 = (1-(1/Nn)*(1-(2/(gg+1))))^(gg/(gg-1));
    PaovP06 = Pa/P06;
    %Choke Test

    if(PstarovP06 >= PaovP06)
        choked = 1;
    elseif (PstarovP06 > PaovP06)
        choked =2;
    end
    %If Choked

    if(choked == 1)
        Me = 1;
        Pe = P06*(PstarovP06);
        Te = T06*(2/(gg+1));
        rhoe = Pe/(Ra*Te);   %kg/m^3



        Ve = sqrt(gg*Ra*1000*Te);
        double(Ve);   
    %If not choked

    elseif (choked == 2)
        Pe = Pa;  
        Te = T06*(1-(Nn(1-(Pe/P06)^((gg-1)/gg))));
        rhoe = Pe/(Ra*Te);   %kg/m^3
        Me = sqrt(((T06/Te)-1)*(2/(gg-1)));   %Me<1 if not choked

        Ve = Me*sqrt(gg*Ra*1000*Te);
    end
%Compute Fs (Specific Thrust)

Fs = ((1+f)*Ve)-Va +(Pe*1000-Pa*1000)*((1+f)/(rhoe*Ve));  %N/(kg/sec)

double(Fs)
%Compute tsfc

tsfc = f/Fs*3600;
double(tsfc)             %kgfuel/hour/N

Below is the way that I tried to edit the code to accomplish this and it has been running for about 5 minutes with no results.

clc clear format long g

%------------------------------------------------------
%Given Values                                         |
Ra = 0.287;   %kJ/kg/k                                |
ga = 1.4;     %Gamma Air                              |
gg = 1.3333;   %Gamma gas                              |
Nd = 0.93;    %Efficiency of Diffuser                 |
Nc = 0.87;    %Efficiency of Compressor               |
Nb = 0.98;    %Efficiency of Combustion Chamber       |
Nm = 0.99;    %Efficiency of Mechanical               |
Nt = 0.90;    %Efficiency of Turbines                 |
Nn = 0.95;    %Efficiency of Nozzle                   |
Ma = 0.84;    %Mach                                   |
Alt = 5;      %km                                     |
%OPR = 9;      %Overal Pressure Ratio                 |   
%prompt = 'Overall Pressure Ratio? ';  %               |
%OPR = input(prompt);                  %               |
%PRClp = sqrt(OPR); %Pressure Ratio over LP Compressor |
%PRChp = sqrt(OPR); %Pressure Ratio over HP Compressor |
deltap0 = 0.04;   % 4% of delievery a loss            |
%-----------------------------------------------------|
%From Handout 5 based on Altitude of 5km--------------|
Ta = 255.7;    %K                                     |
Pa = 54.05;    %kPa                                   |
rho = 0.736;  %kg/m^3                                 |
u = 1.63;     %N*s/m^2                                |
%-----------------------------------------------------|
     % Tp
      %prompt = 'Temperature of Products? ';
      %Tp = input(prompt);
      %T04 = Tp;
     for TINT = 1000:1700
         TINT = TINT + 100;
         for OPR = 2:28
             OPR = OPR + 2;
             PRClp = sqrt(OPR); %Pressure Ratio over LP Compressor 
             PRChp = sqrt(OPR); %Pressure Ratio over HP Compressor
              Va = Ma * sqrt(ga*Ra*Ta*1000);   %in m/s
              Cpoa = (ga*Ra)/(ga-1);          %in kJ/kg/K
              T01 = Ta*(1 + ((ga-1)/2)*Ma^2);
              P01 = Pa*(1 + Nd*((ga-1)/2)*(Ma^2))^(ga/(ga-1));
              P02 = P01*PRClp;
              T02 = T01 + ((T01*((PRClp)^((ga-1)/ga)-1)/Nc));
              Wclp = Cpoa*(T01-T02);
              P03 = P02*PRChp;
              T03 = T02 + ((T02*((PRChp)^((ga-1)/ga)-1)/Nc));
              Wchp = -Cpoa*(T03-T02);
              Tr = T03;
              Tp = TINT;
              T04 = Tp;
      if (Tp <= 1600 && Tp>=400)
              COAp = 299180; COBp = 37.85; COCp = -4571.9; CO2Ap = 56835; 
              CO2Bp = 66.27; CO2Cp = -11634; H2OAp = 88923; H2OBp = 49.36; 
              H2OCp = -7940.8; N2Ap = 31317; N2Bp = 37.46; N2Cp = -4559.3;
              O2Ap = 43388; O2Bp = 42.27; O2Cp = -6635.4;
          %Tp > 1600K
          elseif (Tp > 1600);
              COAp = 309070; COBp = 39.29; COCp = -6201; CO2Ap = 93048; 
              CO2Bp = 68.58; CO2Cp = -16979; H2OAp = 154670; H2OBp = 60.43; 
              H2OCp = -19212; N2Ap = 44639; N2Bp = 39.32; N2Cp = -6753.4; 
              O2Ap = 127010; O2Bp = 46.25; O2Cp = -18798;
      end
MC = 14.4;
MH = 24.9;
MO = 0;
Hrp = -8561991.6;
Ycc = MC + (MH/4) + (MO/2);
Ymin = Ycc - (MC/2); 
Mfuel = MC*12 + MH + MO*16;
Mair = 28.97;
%Find f
    %based on Tr < 1600K
    COAr = 299180; COBr = 37.85; COCr = -4571.9; CO2Ar = 56835;
    CO2Br = 66.27; CO2Cr = -11634; H2OAr = 88923; H2OBr = 49.36;
    H2OCr = -7940.8; N2Ar = 31317; N2Br = 37.46; N2Cr = -4559.3;
    O2Ar = 43388; O2Br = 42.27; O2Cr = -6635.4;
      %based on Excess Air
      N1 = 0;
      N2 = MC;
      N3 = (MH/2); 
      HbarCO = 0;
      HbarCO2 = (CO2Ap + CO2Bp*Tp + CO2Cp*log(Tp)) - (CO2Ar + CO2Br*Tr + CO2Cr*log(Tr));
      HbarH2O = (H2OAp + H2OBp*Tp + H2OCp*log(Tp)) - (H2OAr + H2OBr*Tr + H2OCr*log(Tr));
      HbarO2 = (O2Ap + O2Bp*Tp + O2Cp*log(Tp)) - (O2Ar + O2Br*Tr + O2Cr*log(Tr));
      HbarN2 = (N2Ap + N2Bp*Tp + N2Cp*log(Tp)) - (N2Ar + N2Br*Tr + N2Cr*log(Tr));
      syms Y
      eqn = Hrp + (N1 * 282800) + (N2 * HbarCO2) + (N3 * HbarH2O) + ((Y - Ycc)* HbarO2) + ((3.76 * Y)* HbarN2)==0;
      Y = solve(eqn,Y);
      %Finding f   
      fideal= ((1/(4.76*Y)))*((Mfuel)/(Mair));
      f = fideal/0.98;  %Getting f actual from ideal and efficiency
     P04 = (1- deltap0)*P03;
     Cp0g = (gg*Ra)/(gg-1);      %in kJ/kg/K
     T05 = T04 +((Wchp)/(Nm*(1+0.0168170)*Cp0g));  %K
     P05 = (P04)*(1-((1-(T05/T04))/Nt))^(gg/(gg-1));  %kPa
     T06 = T05 + (Wclp/(Nm*(1+f)*Cp0g));
     %double(T06);
     P06 = P05*(1-((1-(T06/T05))/Nt))^(gg/(gg-1));  %Page 26!!!!
    double(P06);
    %Return to the particular converging nozzle of the turbojet problem
    P0i = P06;
    T0i = T06;
    T0e = T06;
    PstarovP06 = (1-(1/Nn)*(1-(2/(gg+1))))^(gg/(gg-1));
    PaovP06 = Pa/P06;
    %Choke Test
    if(PstarovP06 >= PaovP06)
        choked = 1;
    elseif (PstarovP06 > PaovP06)
        choked =2;
    end
    %If Choked
    if(choked == 1)
        Me = 1;
        Pe = P06*(PstarovP06);
        Te = T06*(2/(gg+1));
        rhoe = Pe/(Ra*Te);   %kg/m^3
        Ve = sqrt(gg*Ra*1000*Te);
        double(Ve);   
    %If not choked
    elseif (choked == 2)
        Pe = Pa;  
        Te = T06*(1-(Nn(1-(Pe/P06)^((gg-1)/gg))));
        rhoe = Pe/(Ra*Te);   %kg/m^3
        Me = sqrt(((T06/Te)-1)*(2/(gg-1)));   %Me<1 if not choked
        Ve = Me*sqrt(gg*Ra*1000*Te);
    end
%Compute Fs (Specific Thrust)
Fs = ((1+f)*Ve)-Va +(Pe*1000-Pa*1000)*((1+f)/(rhoe*Ve));  %N/(kg/sec)
double(Fs);
%Compute tsfc
tsfc = f/Fs*3600;
double(tsfc);            %kgfuel/hour/N
Fs(TINT) = Fs;
       end
   end
   Fs

Thank you for any help you may be able to provide,

Colby

Best Answer

Colby - you almost have it. Your code for the two for loops is as follows

     for TINT = 1000:1700
         TINT = TINT + 100;
         for OPR = 2:28
             OPR = OPR + 2;

Since TINT and OPR are your indexing variables, you don't need to increment either through code. This will happen automatically so long as you set your step value correctly in the for loop. The above can be rewritten as

     for TINT = 1100:100:1800
         for OPR = 4:2:30
             % etc.
         end
     end

So on the first iteration of the outer for loop, TINT will be assigned a value on 1100, and subsequent iterations of this loop will assign values of 1200, 1300, 1400,...,1800 to TINT.

The inner for loop will initialize *OPR to 4 with subsequent iterations assigning it values of 6,8,10,...,30.

Related Solutions

MATLAB: How to run script iteratively

 clc; clear; close all;
          format long 
          Ma=.84
          Ra=.287
          Rg=.287
          siga=1.4
          sigg=1.3333
          Pa=54.05
          Ta=225.7
          OPR=2:3:200
          tsfc=cell(1,67) ;      %PRE-ALLOCATION
          for i = 1:numel(OPR)
          Va=Ma*sqrt(siga*Ra*Ta)
          Cpoa= (siga*Ra)./(siga - 1)
          TINT = 1100:100:1800
          PRC = sqrt(OPR(i))
          % Diffuser
          efd=.93
          To1 = Ta*(1+((siga-1)/2)*(Ma^2))
          Toa=To1
          Po1 = Pa*(1 + efd*((siga - 1)/2)*(Ma^2))^(siga/siga-1)
          Poa=Pa*(Toa/Ta)^(siga/siga-1)
          %Compressors  
          efc=.87
          %Low Pressure
          Po2=Po1.*(PRC)
          To2 = To1 + ((To1*((PRC.^((siga-1)/siga)-1)))/efc)
          WC_LP = Cpoa*(To1-To2)
          %High Pressure
          To3 = To2 + ((To2.*((PRC.^((siga-1)/siga)-1)))/efc)
          Po3=Po2.*(PRC)
          WC_HP = Cpoa*(To2-To3)
          %Combustion
          efb =.98
          deltPob = .4 
          Po4= (1-deltPob)*Po3
          To4=TINT
          %Combustion Calculations for f
          MC=14.4
          MH=24.9
          MO=0
          Ycc=MC + (MH/4)- (MO/2)
          Tp=To4
          Tr=To3
          To= 298
          HRP_to = -8561991.6
          if Tp <= 1600
          hco2p = 56835 + 66.27*Tp + -11634.0*log(Tp)
          hh2op = 88923 + 49.36*Tp + -7940.8*log(Tp)
          hn2p= 31317 + 37.46*Tp + -4559.3*log(Tp)
          ho2p = 43388 + 42.27*Tp + -6635.4*log(Tp)
          else
          hco2p = 93048 + 68.58*Tp + -16979.0*log(Tp)
          hh2op = 154670 + 60.43*Tp + -19212.0*log(Tp)
          hn2p = 44639 + 39.32*Tp + -6753.4*log(Tp)
          ho2p =127010 + 46.25*Tp + -18798*log(Tp)   
          end
          if Tr <= 1600
          hco2r = 56835 + 66.27*Tr + -11634.0*log(Tr)
          hh2or = 88923 + 49.36*Tr + -7940.8*log(Tr)
          hn2r= 31317 + 37.46*Tr + -4559.3*log(Tr)
          ho2r = 43388 + 42.27*Tr + -6635.4*log(Tr)
          else
          hco2r = 93048 + 68.58*Tr + -16979.0*log(Tr)
          hh2or = 154670 + 60.43*Tr + -19212.0*log(Tr)
          hn2r = 44639 + 39.32*Tr + -6753.4*log(Tr)
          ho2r =127010 + 46.25*Tr + -18798*log(Tr)
          end
          hco2 = hco2p-hco2r
          hh2o = hh2op-hh2or
          hn2 = hn2p-hn2r
          ho2 =ho2p-ho2r
          h2oa=30.5
          h2ob=.0103
          co2a=28.8
          co2b=.028
          o2a=27
          o2b=.0079
          fuela= 38.1
          fuelb=.656
          Delt_HRP = (MH/2)*(h2oa*Tr +.5*h2ob*Tr.^2) - (h2oa*To +.5*h2ob*To.^2)+MC*(co2a*Tr +.5*co2b*Tr.^2) - (co2a*To +.5*co2b*To.^2)- Ycc*(o2a*Tr +.5*o2b*Tr.^2) - (o2a*To +.5*o2b*To.^2)-(fuela*Tr +.5*fuelb*Tr.^2) - (fuela*To +.5*fuelb*To.^2)
          HRP =HRP_to + Delt_HRP
          y = (-HRP - MC*hco2 - (MH/2).*hh2o + Ycc.*ho2)./(3.76.*hn2+ho2)
          fideal=197.7./(4.76*y*28.97)
          f=fideal*efb
          %Turbines | HIgh and Low have equal efficiencies 
          eft=.90
          efm=.99
          Cpog = (sigg*Rg)/(sigg-1)
          % High Pressure
          To5= To4 + (WC_HP./(efm*(1+f)*Cpog))
          Po5 = Po4*(1-((1-(To5/To4))/eft))^(sigg/(sigg-1))
          % Low Pressure
          To6= To5 + (WC_LP./(efm*(1+f)*Cpog))
          Po6 = Po5*(1-((1-(To6/To5))/eft))^(sigg/(sigg-1))
          %Nozzle
          %Choke Test
          efn = .95
          P_Po6 = (1-((1/efn)*(1-(2/(sigg+1)))))^(sigg/(sigg-1))
          R=Pa./Po6
          Me=1
          Pe=Po6*(P_Po6)
          Te=To6*(2/(sigg+1))
          roe= Pe./(Rg*Te)
          Ve= Me*sqrt(sigg*Rg*1000*Te)
          % Fs and T.S.F.C Calculations
          Fs = ((1+f).*Ve-Va) + (Pe-Pa).*((1+f)./(roe.*Ve))
          tsfc{i}= (f./Fs).*3600;
          end
          celldisp(tsfc)        %each subcell represents each value of OPR
tsfc = vertcat(tsfc{:})    %double array

MATLAB: Differential Equation Problem using ode45

[W,Fa] = PBR_Isothermal;
% plot results
subplot(3,2,1)
plot(W,Fa(:,1))
subplot(3,2,2)
plot(W,Fa(:,2))
subplot(3,2,3)
plot(W,Fa(:,3))
subplot(3,2,4)
plot(W,Fa(:,4))
subplot(3,2,5)
plot(W,Fa(:,5))
subplot(3,2,6)
plot(W,Fa(:,6))
function [W,Fa] = PBR_Isothermal
  clc; clear;
  %A = Dammitol
  %B = Oxygen
  %C = Valualdehyde
  %D = Carbon Dioxide
  %E = Water
  %I = Nitrogen
  
  %Feed based on inlet as 100 kmol/h
  
  Fa0 = 0.1*100; 
  Fb0 = 0.07*100;
  Fc0 = 0;
  Fd0 = 0;
  Fe0 = 0.02*100;
  Fi0 = 0.81*100;
 
  
  [W, Fa]  = ode45(@FunA,[0 1000],[Fa0 Fb0 Fc0 Fd0 Fe0 Fi0]);
%   [Fb W] = ode45(FunA,[0 1000],Fb0);
%   [Fc W] = ode45(FunA,[0 1000],Fc0);
%   [Fd W] = ode45(FunA,[0 1000],Fd0);
%   [Fe W] = ode45(FunA,[0 1000],Fe0);
end
function A = FunA(W,F)
%Total Pressure
Ptot = 114.7;
%Defining Constant
FT   = F(1)+F(2)+F(3)+F(4)+F(5)+F(6);
RT   = (1/353-1/373)/1.987;
%Partial Pressure of Each
PDT  = (F(1)/FT)*Ptot;
PO2  = (F(2)/FT)*Ptot;
PVA  = (F(3)/FT)*Ptot;
PCO  = (F(4)/FT)*Ptot;
PWA  = (F(5)/FT)*Ptot;
PN2  = (F(6)/FT)*Ptot;
%Constants
k1 = 1.771*10^-3; 
k2 = 23295; 
k3 = 0.5; 
k4 = 1.0; 
k5 = 0.8184; 
k6 = 0.5; 
k7 = 0.2314;
k8 = 1.0; 
k9 = 1.25; 
k10 = 2.0; 
k11 = 2.795*10^-4; 
k12 = 33000; 
k13 = 0.5; 
k14 = 2.0;
%Rate of Reaction
r1 = (k1*exp(-k2*RT)*PO2^k3*PDT^k4)/(1+k5*PO2^k6+k7*PDT^k8+k9*PVA^k10); 
r2 = (k11*exp(-k12*RT)*PO2^k13*PVA^k14)/(1+k5*PO2^k6+k7*PDT^k8+k9*PVA^k10);
%Molar Flow Rate of Each Component
A(1) = -r1;
A(2) = -1/2*r1;
A(3) = r1;
A(4) = 2*r2;
A(5) = r1+2*r2;
Fi0 = 0.81*100;
A(6) = Fi0;
A=A';
end

Best Answer

Related Solutions

MATLAB: How to run script iteratively

MATLAB: Differential Equation Problem using ode45

Related Question