Hey guys, I have a problem. Im writing a code that models engine thrust over an altitude range. It runs ok when i keep a variable constant, 'gamma', but when its varied, it dosent run. Can someone help?

clc
clear
z = 0:1:100;
[T,P,rho,del,sig,the]= atmosi_100(z);
 %% constants       
    r_o = 6.371e7;                                                          % earth radius, m
    g_o = 9.81;                                                             % gravity, m/s^2
    %T/W = F_req/(m_e*9.81);
    D = 100;                                                                % Drag, N, from aero
    m_i = 12000;                                                            % initial wet mass, kg
    P_c = 2.233e7;                                                          % ambient pressure, Pa
    g = 9.81;                                                               % gravitational acceleration, m/s^2
    R = 287;                                                                % gas constatnt, J/kg*K
    A_e = 1.1;                                                              % exit area, m^2
    % From App.B
    T_c = 3400;                                                             % temp, K
    gamma = 1.2:0.05:1.4;
    Mol = 13.5;                                                             % molecular mass, kg/kmol
    v_o = 10;                                                               % inicial velocity, m/s
    time = 140;                                                             % time to seperation, s
    
    %% Inputs
    F_req = 1334466;                                                        % required thrust, N (from Performance)

%% Engine characteristics/SSME design point    
% assuming combustion efficiency of 1.0
    ex = 77.5;
    P_c = 2.2339e+7;                                                        % chamber pressure, Pa
    A_t = (ex.^-1).*A_e;                                                       % throat area, m^2
    
     epsilon = sqrt(gamma./(((gamma+1)./2).^((gamma+1)./(gamma-1))));
     m_dot = P_c.*A_t.*g.*(epsilon./sqrt(R.*g.*T_c));

% assume f and I_sp
    f = 6;
    I_sp = 363;
% Newton-Raphson method for finding exit mach number    
    b = 2.*(gamma-1)./(gamma+1);
    c = ((gamma+1)./(gamma-1)).*((ex).^(b));
    d = 2./(gamma-1);
    
    Me_0 = 7.5;
    original = (c .* (Me_0 .^ b)) - (Me_0 .^ 2) - d;
    derivative = (b .* c .* (Me_0 .^ (b - 1))) - (2 .* Me_0);
    
     Me_1 = Me_0 - (original ./ derivative);
    
    while abs(original) > 10^-2
        Me_0 = Me_1;
        
        original = (c .* (Me_0 .^ b)) - (Me_0 .^ 2) - d;
        derivative = (b .* c .* (Me_0 .^ (b - 1))) - (2 .* Me_0);
        Me_1 = Me_0 - (original ./ derivative);
        
        fprintf('Iteration: M_e = %.20f, error = %.20f \n', Me_1, original);
    end
    
    mf_f = 1-exp(-(sqrt(g_o*r_o))./(I_sp*v_o*((1-D)/F_req)));                % fuel mass fraction
    mf = mf_f *m_i;                                                         % mass of fuel, kg,


P_e = P_c.*(1+((gamma-1)./2).*(Me_1.^2)).^(gamma./(1-gamma));               % exit static pressure, Pa
P_te = P_e./(1+((gamma-1)./(Me_1.^2))).^(-gamma./(gamma-1));                % exit total pressure, Pa

% using isentropic relations
    T_e = ((P_e./P_te).*T_c.^(gamma./(gamma-1))).^(1./(gamma./(gamma-1)));
    v_e = Me_0.*sqrt(gamma.*R.*T_e);
    F_ava = (m_dot.*v_e) + (P_e-P).*(A_t.*ex);
    S = (m_dot.*f)./F_ava;                                                   % specific fuel consumption, kg/s*N
    S_t = (S.*F_ava).*time;
    
%% engine sizing   
    D_e = 0.00357.*F_ava+14.48;                                              % engine diameter, m
    m_e = F_ava./(g*0.0006098.*F_ava+13.44);                                 % engine mass, kg
    L_e = 0.000030*F_ava+327.7;                                              % engine length, m
    Lstar = 1;                                                               % characteristic length, m
    A_c = A_t*(8*((2*sqrt(A_t/pi))^(-0.6)) + 1.25);                          % average combustion chamber cx area, m^2
    L_c = Lstar*(A_t/A_c);                                                   % combusiton chamber length, m            
    V_c = Lstar*A_t;                                                         % combustion chamber volume, m^3
    d_e  = sqrt((4.*ex.*A_t)/pi);                                            % exit diameter, m


%% Plots/outputs

    fprintf('m_dot = %f [kg/s] \n',m_dot)
    fprintf('v_e = %f [m/s] \n',v_e)
    fprintf('V_c = %f [m^3] \n,',V_c)
    fprintf('A_c = %f [m^2] \n', A_c)
    fprintf('m_e = %f [kg] \n',m_e)
    fprintf('L_e = %f [m] \n', L_e)
    fprintf('D_e = %f [m]',D_e)
    fprintf('L_c = %f [m] \n', L_c)
    fprintf('S_t = %f [kg]\n', S_t)
%   fprintf('m_f = %f [kg]\n', m_f)
 
    figure(1)
    plot(z,F_ava)
    title('Avalible thrust vs. Altitude')
    xlabel('Altitude (km)')
    ylabel('F_ava (N)')
    
    figure(2)
    plot(z,S)
    title('Specific fuel consimption vs. Altitiude')
    xlabel('Altitude (km)')
    ylabel('S (kg/s*N)')
    
 %% Trade Studies
% plot(F_ava,gamma)  % vary gamma from 1.2-1.4

[Miriam]

Could you attach the function atmosi_100?

[Ikenna Okoye]

function[T,P,rho,del,sig,the]= atmosi_100(z)
global r R_o
%% Constants
% Sea level conditions
P_o   = 101325; % (Pa)
rho_o = 1.225; % (Kg/m^3)
T_o   = 288.15; % (K)
mu_o  = 1.789*10^-5; % (Kg/m/s)
nu_o  = 1.46*10^-5; % (m^2/s)
tc_o  = 0.02596; % (W/m/K)
R_o   = 287.053; % (J/Kg/K)
g_o   = 9.80665; % (m/s^2)
Cp    = 1005; % (J/Kg/K)
r_o   = 6356.766; % (km) 
r     = 1.4;
%% Process
for i = 1:length(z)    
  
h(i) = (r_o*z(i))/(r_o+z(i));

if h(i) == 0
    P(i) = P_o;
    T(i) = T_o;
  rho(i) = rho_o;
  
elseif 0<h(i) && h(i)<=11
    P(i) = P_o*(T_o/(T_o-6.5*h(i)))^(34.1632/-6.5);
    T(i) = T_o-6.5*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = z(i)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 11<h(i) && h(i)<=20
    T(i) = 216.65;
    P(i) = 22632.06*exp(-34.1632*(h(i)-11)/T(i));
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-11)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 20<h(i) && h(i)<=32
    P(i) = 5474.889*(216.65/(216.65+(h(i)-20)))^34.1632;
    T(i) = 196.65+h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-20)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
 
elseif 32<h(i) && h(i)<=47
    P(i) = 868.0187*(228.65/(228.65+2*r*(h(i)-32)))^(34.1632/(2*r));
    T(i) = 139.05+2*r*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-32)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 47<h(i) && h(i)<=51
    T(i) = 270.65;
    P(i) = 110.9063*exp(-34.1632*(h(i)-47)/T(i));
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-47)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 51<h(i) && h(i)<=71
    P(i) = 66.93884*(270.65/(270.65-2*r*(h(i)-51)))^(34.1632/(-2*r));
    T(i) = 413.45-2*r*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-51)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 71<h(i) && h(i)<=85
    P(i) = 3.956420*(214.65/(214.65-2*(h(i)-71)))^(34.1632/-2);
    T(i) = 356.65-2*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-71)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 85<z(i) && z(i)<=91
    T(i) = 186.8673;
    P(i) = exp(2.159582*10^(-6)*z(i)^3 + -4.836957*10^(-4)*z(i)^2 + -.14251928*z(i)  +  13.47530);
  rho(i) = exp(-3.322622*10^(-6)*z(i)^3 + 9.111460*10^(-4)*z(i)^2  +  -.2609971*z(i)  +  5.944694);
  eps(i) = (z(i)-85)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
else 
    T(i) = 263.1905-76.3232*sqrt(1-((z(i)-91)/-19.9429)^2);
    P(i) = exp(3.304895*10^-5*z(i)^3  +  -.009062730*z(i)^2  +  .6516698*z(i)  +  -11.03037);
  rho(i) = exp(2.873405*10^-5*z(i)^3  +  -.008492037*z(i)^2  +  .6541179*z(i)  +  -23.6201);
  eps(i) = (z(i)-85)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% run from this code below
prompt = 'For SI select 1 for imperial select 2?';
unit = input(prompt);
if unit==1
prompt = 'What is your height value in kilometers?';
z = input(prompt);
else unit==2
prompt = 'What is your height value in feet';
zi = input(prompt);
z = convlength([zi],'ft','km');
end  
[T,P,rho,zIm,TIm,PIm,rhoIm,a,del,sig,the,aIm]= atmos0_100(z);
fprintf('For altitude of   z = %.3f [km]     z = %.3f [ft] \n',z,zIm);
fprintf('Speed of sound is a = %.3f [m/s]    a = %.3f [ft/s] \n',a,aIm);
fprintf('        SI                    Imperial                          Ratios   \n');
fprintf('------------------------------------------------------------------------- \n');
fprintf('P =   %.3f [Pa]       P = %.3f [psf]                del = %.5f  \n',P,PIm,del);
fprintf('T =   %.3f [K]          T = %.3f [R]                   the = %.5f  \n',T,TIm,the);
fprintf('rho = %f [kg/m^3]    rho = %f [slug/ft^3]           sig = %f  \n',rho,rhoIm,sig);

[Miriam]

Would you want to generate the plots and print outputs for each value of gamma? In that case you could use a for loop.

[Ikenna Okoye]

No, my plan was to make a carpet plot of thrust vursus altitude, by varying gamma and expantion ratio,e.

[Walter Roberson]

long code and you did not give us any error message to work with and you did not describe what we should pay attention to in order to know if it is working or not . Are you expecting people to be able to just glance at the code and read your mind and tell you what change needs to be made to have it do whatever is in your imagination ?

[Ikenna Okoye]

There was only one error that came up when run in MATLAB with the above code and the above comment detial what I wanted to do. Original question is first code file, the comments have an atmosphere code that finds the perameters to 100km altitude, one is the function file the other below is where its run from.

Here's the error that arises:

Matrix dimensions must agree.

Error in Pro_rev3 (line 68)

F_ava = (m_dot.*v_e) + (P_e-P).*(A_t.*ex);

[Star Strider]

You need to post your ‘atmosi_100’ function.

I do not see it in what you posted, and your code will noit work without it.

[Ikenna Okoye]

It's in the second comment above.

[Walter Roberson]

We need function atmos0_100 as well. It is called at line 120 of atmosi_100()

It would also help to know what inputs to use to reproduce this.

[Ikenna Okoye]

What inputs do you mean? As in equations, or do you mean the values already hard coded in? The code runs ok when 'gamma' is kept constant, everything works ok then, so im not sure why when run by itself it dosent work for me.

[Walter Roberson]

the code prompts 1 for si and 2 for imperial, and prompts for aa numeric height value . Then it tries to call atmos0_100 which we do not have the code for .

Answer 1

[Ikenna Okoye]

Sorry, the script that i first posted should already be atmos0_100, the 'i' should be a '0':

function[T,P,rho,del,sig,the]= atmos0_100(z)
global r R_o
%% Constants
% Sea level conditions
P_o   = 101325; % (Pa)
rho_o = 1.225; % (Kg/m^3)
T_o   = 288.15; % (K)
mu_o  = 1.789*10^-5; % (Kg/m/s)
nu_o  = 1.46*10^-5; % (m^2/s)
tc_o  = 0.02596; % (W/m/K)
R_o   = 287.053; % (J/Kg/K)
g_o   = 9.80665; % (m/s^2)
Cp    = 1005; % (J/Kg/K)
r_o   = 6356.766; % (km) 
r     = 1.4;
%% Process
for i = 1:length(z)    
  
h(i) = (r_o*z(i))/(r_o+z(i));
if h(i) == 0
    P(i) = P_o;
    T(i) = T_o;
  rho(i) = rho_o;
  
elseif 0<h(i) && h(i)<=11
    P(i) = P_o*(T_o/(T_o-6.5*h(i)))^(34.1632/-6.5);
    T(i) = T_o-6.5*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = z(i)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 11<h(i) && h(i)<=20
    T(i) = 216.65;
    P(i) = 22632.06*exp(-34.1632*(h(i)-11)/T(i));
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-11)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 20<h(i) && h(i)<=32
    P(i) = 5474.889*(216.65/(216.65+(h(i)-20)))^34.1632;
    T(i) = 196.65+h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-20)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
 
elseif 32<h(i) && h(i)<=47
    P(i) = 868.0187*(228.65/(228.65+2*r*(h(i)-32)))^(34.1632/(2*r));
    T(i) = 139.05+2*r*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-32)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 47<h(i) && h(i)<=51
    T(i) = 270.65;
    P(i) = 110.9063*exp(-34.1632*(h(i)-47)/T(i));
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-47)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 51<h(i) && h(i)<=71
    P(i) = 66.93884*(270.65/(270.65-2*r*(h(i)-51)))^(34.1632/(-2*r));
    T(i) = 413.45-2*r*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-51)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 71<h(i) && h(i)<=85
    P(i) = 3.956420*(214.65/(214.65-2*(h(i)-71)))^(34.1632/-2);
    T(i) = 356.65-2*h(i);
  rho(i) = P(i)/(R_o*T(i));
  eps(i) = (z(i)-71)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
elseif 85<z(i) && z(i)<=91
    T(i) = 186.8673;
    P(i) = exp(2.159582*10^(-6)*z(i)^3 + -4.836957*10^(-4)*z(i)^2 + -.14251928*z(i)  +  13.47530);
  rho(i) = exp(-3.322622*10^(-6)*z(i)^3 + 9.111460*10^(-4)*z(i)^2  +  -.2609971*z(i)  +  5.944694);
  eps(i) = (z(i)-85)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
else 
    T(i) = 263.1905-76.3232*sqrt(1-((z(i)-91)/-19.9429)^2);
    P(i) = exp(3.304895*10^-5*z(i)^3  +  -.009062730*z(i)^2  +  .6516698*z(i)  +  -11.03037);
  rho(i) = exp(2.873405*10^-5*z(i)^3  +  -.008492037*z(i)^2  +  .6541179*z(i)  +  -23.6201);
  eps(i) = (z(i)-85)*10^3/(R_o*T(i));
  del(i) = P(i)/P_o;
  sig(i) = rho(i)/rho_o;
  the(i) = T(i)/T_o;
  
end
end

[Walter Roberson]

Inside atmosi_100 you have the call

[T,P,rho,zIm,TIm,PIm,rhoIm,a,del,sig,the,aIm]= atmos0_100(z);

Therefore your function atmos0_100 must be able to return at least 12 outputs.

The code for atmos0_100 that you have posted here starts with

function[T,P,rho,del,sig,the]= atmos0_100(z)

so it returns 6 outputs. This conflicts with the call to the function. Either your atmosi_100 is making the wrong call or else you have posted the wrong atmos0_100

[Ikenna Okoye]

I have to see what the deal is, the atmospere is a team members code that im imbeding into my propulsion code. I'll get back to you.

[Ikenna Okoye]

Revised code, works with atmos0_100 above:

prompt = 'For SI select 1 for imperial select 2?';
unit = input(prompt);
if unit==1
prompt = 'What is your height value in kilometers?';
z = input(prompt);
else unit==2
prompt = 'What is your height value in feet';
zi = input(prompt);
z = convlength([zi],'ft','km');
end  
[T,P,rho,del,sig,the]= atmos0_100(z)
fprintf('For altitude of   z = %.3f [km]     ',z);

fprintf('        SI                    Imperial                          Ratios   \n');
fprintf('------------------------------------------------------------------------- \n');
fprintf('P =   %.3f [Pa]                      del = %.5f  \n',P,del);
fprintf('T =   %.3f [K]                         the = %.5f  \n',T,the);
fprintf('rho = %f [kg/m^3]             sig = %f  \n',rho,sig);

[Walter Roberson]

Where is that code to be placed? In the version you posted before it goes immediately after the end of the loop

for i = 1:length(z)

in atmosi_100. The code

[T,P,rho,del,sig,the]= atmos0_100(z)

assigns to P, overwriting the P that was constructed in the for loop, leading to trouble.

[Ikenna Okoye]

The code i posted, 2 comments up, would be in its own file. Total there are three code files: the one just mentioned, the atmosphere program it calls, atmosi0_100 and the propulsion code that I first posted in my question.

[Walter Roberson]

"Total there are three code files: the one just mentioned, the atmosphere program it calls, atmosi0_100 and the propulsion code that I first posted in my question. "

"The one just mentioned" appears to be the script starting from "prompt =".

You seem to be saying that the atmosphere program it calls is the routine named atmosi0_100, which would be consistent for that script.

The propulsion code that you first posted appears to refer to the block of code that starts with

clc
clear
z = 0:1:100;
[T,P,rho,del,sig,the]= atmosi_100(z);

which calls atmosi_100 which is a different function, a fourth function.

Considering that the code posted first starts with "clear" then it would not appear to make sense to call the script first to establish any variables, as any variables would be erase by the clear. This suggests that the script that does the prompt is to be called at some point after the propulsion code, but the relationship between them is not obvious.

Answer 2

[Ikenna Okoye]

Right the propulsion code calls the funciton file, that starts with prompt, that then calls atmos0_100. Does that clarify?