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?
Mostra commenti meno recenti
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
11 Commenti
Miriam
il 16 Nov 2018
Could you attach the function atmosi_100?
Ikenna Okoye
il 16 Nov 2018
Miriam
il 17 Nov 2018
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
il 17 Nov 2018
Walter Roberson
il 17 Nov 2018
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
il 17 Nov 2018
Star Strider
il 17 Nov 2018
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
il 17 Nov 2018
Walter Roberson
il 17 Nov 2018
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
il 17 Nov 2018
Walter Roberson
il 17 Nov 2018
Modificato: Walter Roberson
il 18 Nov 2018
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 .
Risposte (2)
Ikenna Okoye
il 18 Nov 2018
6 Commenti
Walter Roberson
il 18 Nov 2018
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
il 18 Nov 2018
Ikenna Okoye
il 18 Nov 2018
Walter Roberson
il 18 Nov 2018
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
il 19 Nov 2018
Modificato: Ikenna Okoye
il 19 Nov 2018
Walter Roberson
il 19 Nov 2018
"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.
Ikenna Okoye
il 19 Nov 2018
0 voti
Categorie
Scopri di più su Classification Ensembles in Centro assistenza e File Exchange
il 16 Nov 2018
il 19 Nov 2018
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
