Simple For cycle (a small problem with the indexes)

Question

Paul Rogers 2021-1-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/709328-simple-for-cycle-a-small-problem-with-the-indexes

评论： Mischa Kim 2021-1-5

采纳的回答： Mischa Kim

matlab_10.mat

在 MATLAB Online 中打开

I have a small problem with the last loop in the code:

clear
clc
close all
load matlab_10.mat
init= 24500;
%% Gas properties
k = 1.4;          %Gas constant ratio, cp/cv()
R = 287.05;       %Ideal gas constant [J/Kg*K]
p01 = 1e5;        %Compressor inlet pressure(Pa)
T01 = 303.35;     %Compressor inlet temperature(K)
cp = 1005;        %Specific heat capacity for constant pressure, property of gas(J/(kg*K))
ro1 = 1.15;       %Gas density(kg/m^3)
AFR = 14.71;      %Air Fuel Ratio
%% Turbine's parameters
A_eff = .03;       %[m^2]
A_0 = A_eff/.6;
E_R = [1:0.001:4]; %Expansion ratio
T_3 = 1173.15;     %Discharge temperature [K]
eta_t = 0.6;
%% Compressor's parameters
eta_c = 0.95;
I = 0.001;                                  %Moment of inertia for the compressor spool(kg*m^2)
mass_flow_compressor = 0.273675009794891;   %adimensional mass flow from Greitzer steady-state [kg/s]
% m_c = mass_flow_compressor*(ro1*Ac*U);    %actual compressor's mass flow [kg/s]  
mass_flow_compressor_corrected = (mass_flow_compressor*(T01^0.5))./p01; %[m*s*(k)^0.5]
P_R = 2.591609007038750; %pressure ratio from Greeitzer
p2 = P_R*p01;  %[Pa] 
%% 
mass_flow_compressor_corrected = (mass_flow_compressor_corrected*ones(1,length(E_R))).*(1+(1./AFR));
mass_flow_turbine_corrected = A_eff.*((k./R).^0.5).*((1./E_R).^(1/k)).*(((2./k-1).*(1-(1./E_R).^((k-1)./k))).*0.5);
% mass_flow_corrected = (mass_flow).*p01;
x_ER = interp1(mass_flow_turbine_corrected - mass_flow_compressor_corrected, E_R, 0);
y_ER = A_eff.*((k./R).^0.5).*((1./x_ER).^(1/k)).*(((2./k-1).*(1-(1./x_ER).^((k-1)./k))).*0.5);
figure()
plot(E_R,mass_flow_turbine_corrected,...
    E_R,mass_flow_compressor_corrected,'linewidth',3)
grid on
grid minor
xlabel('ER')
legend('turbine','compressor')
ylabel('m*s*K^{0.5}')
title('ER vs mass flow')
hold on
plot(x_ER, y_ER, 'pg','linewidth',3)
hold off
text(x_ER, y_ER, sprintf('\\uparrow\nIntersection:\nx\\_ER = %.3f\ny\\_ER = %.3f', x_ER, y_ER),...
    'HorizontalAlignment','left', 'VerticalAlignment','top')
%% Power Steady Case
Pt = ((y_ER.*p01)./(T01^0.5)).*eta_t.*(k./(k-1)).*R.*T_3.*(1-(1./(x_ER.^((k-1)./k))));
Pc = ((y_ER.*p01)./(T01^0.5)).*(1./eta_c).*(k./(k-1)).*R.*T01.*((p2./p01).^((k-1)./k)-1);
%% Pulsating
mass_flow_compressor_vet = phi_10(init:end);
mass_flow_compressor_vet_corrected = (mass_flow_compressor_vet*(T01^0.5))./p01; %[m*s*(k)^0.5]
mass_flow_compressor_vet_corrected = mass_flow_compressor_vet_corrected*ones(1,length(E_R)).*(1+(1./AFR));
for i=1:length(phi_10(init:end))
%        
   x_ER_vet(i,:) = interp1(mass_flow_turbine_corrected - mass_flow_compressor_vet_corrected(i,:), E_R, 0);
   y_ER_vet(i,:) = A_eff.*((k./R).^0.5).*((1./x_ER_vet(i,:)).^(1/k)).*(((2./k-1).*(1-(1./x_ER_vet(i,:)).^((k-1)./k))).*0.5);
% 
%    Pt_vet = ((y_ER_vet.*p01)./(T01^0.5)).*eta_t.*(k./(k-1)).*R.*T_3.*(1-(1./(x_ER_vet.^((k-1)./k))));
%    Pc_vet = ((y_ER_vet.*p01)./(T01^0.5)).*(1./eta_c).*(k./(k-1)).*R.*T01.*((p2./p01).^((k-1)./k)-1);
%    omega_vet = (2.*((Pt_vet-Pc_vet)./(I))+55000.^2).^0.5;
end
figure()
plot(t_10(init:end),omega_vet,'linewidth',3)
grid on
grid minor
xlabel('t [s]')
legend('\omega')
ylabel('rpm')
hold on
title('\omega')

my problem is I have the first omega, let's call it omega(1), and by that iteration I need to find the others.

The algorithm I wrote on paper foor the first two steps is in the following figure: (anyway I've already written the code in the comment inside the for-loop, i just miss the indexes)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mischa Kim 2021-1-5

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/709328-simple-for-cycle-a-small-problem-with-the-indexes#answer_591463

编辑：Mischa Kim 2021-1-5

在 MATLAB Online 中打开

There you go:

omega_vet(1) = 55000;
for i=1:length(phi_10(init:end))
       
   x_ER_vet(i) = interp1(mass_flow_turbine_corrected - mass_flow_compressor_vet_corrected(i), E_R, 0);
   y_ER_vet(i) = A_eff.*((k./R).^0.5).*((1./x_ER_vet(i)).^(1/k)).*(((2./k-1).*(1-(1./x_ER_vet(i)).^((k-1)./k))).*0.5);
   Pt_vet(i)   = ((y_ER_vet(i).*p01)./(T01^0.5)).*eta_t.*(k./(k-1)).*R.*T_3.*(1-(1./(x_ER_vet(i).^((k-1)./k))));
   Pc_vet(i)   = ((y_ER_vet(i).*p01)./(T01^0.5)).*(1./eta_c).*(k./(k-1)).*R.*T01.*((p2./p01).^((k-1)./k)-1);
   
   omega_vet(i+1) = (2.*((Pt_vet(i)-Pc_vet(i))./(I)) + omega_vet(i).^2).^0.5;
end
figure()
plot(t_10(init:end),omega_vet(1:end-1),'linewidth',3)

2 个评论
显示无隐藏无

Paul Rogers 2021-1-5

thanks a lot, but nope, since the first value of the omega_vet is 55000.

Mischa Kim 2021-1-5

Right you are. For some reason I missed the info on your paper. See the updated code above. Note that in the final plot command you need to account for the fact that omega is always one index ahead. So you need to adjust the length of the time or the omega vector.

请先登录，再进行评论。