Add these lines to find the intersection:
x_ER = interp1(mass_flow_corrected - m_compressor, 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) * p01;
so the full code is now:
%% 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)
%% Turbine's parameters
A_eff = .0089; %[m^2]
A_0 = A_eff/.6;
dx = 0.001;
E_R = [1:dx:4]; %Expansion ratio
T_3 = 1173.15; %Discharge temperature [K]
%% Compressor's parameters
U = 164.17; %[m/s]
Ac = 0.0026; %[m^2]
m_c_adim = 0.273675009794891; %adimensional mass flow from Greitzer steady-state
m_c = m_c_adim*[ro1*Ac*100*(T01^0.5)]; %actual compressor's mass flow
p2 = 2.591609008885050*p01;
P_R = (p2)/(0.5*ro1*(U^2));
%%
m_compressor = m_c*ones(1,length(E_R));
mass_flow = 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;
%% Plots
figure()
plot(E_R,mass_flow_corrected,...
E_R,m_compressor)
grid on
grid minor
xlabel('ER')
legend('turbine','compressor')
ylabel('corrected mass flow')
title('ER vs mass flow')
x_ER = interp1(mass_flow_corrected - m_compressor, 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) * p01;
hold on
plot(x_ER, y_ER, 'pg')
hold off
text(x_ER, y_ER, sprintf('\\uparrow\nIntersection:\nx\\_ER = %.3f\ny\\_ER = %.3f', x_ER, y_ER), 'HorizontalAlignment','left', 'VerticalAlignment','top')
and the intersection is plotted and labeled as well (if necessary).
The idea is to interpolate the difference and return the associated value of ‘x_ER’, then plug that into the ‘mass_flow_corrected’ expression to get the associated value for ‘y_ER’. It works here because the ‘mass_flow_corrected’ values are monotonically increasing.
.
