Find intersection between two llines

Hi everyone,

I need to find the exact intersection between the vectors m_compressor and mass_flow_corrected.

clear
clc
close all

%% 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')

any (mass_flow_corrected == m_compressor)
intersection = find(mass_flow_corrected == m_compressor);
x_ER = E_R(intersection)
y_ER = mass_flow_corrected(intersection)

I am trying as I read on the forum, but it's not working, I wonder if there is any other way of dooing so.

Answer 1

[Star Strider]

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.

.

[Paul Rogers]

when I'll finish this endless prooject you'll be one I have to buy a beer to.

Thanks a lot, it's not the firrst time you give me help!

[Star Strider]

That would be nice!

As always, my pleasure!