Hi, I have a problem with the function solve (I tried vpasolve and numeric::solve MuPad toolbox). I have a system with one unknown and one equation. It's the equation of isentropic flow in Laval nozzle. I want to determine the outlet pressure. but using solve I have an empty system, and the response is complex using numeric::solve. You can found my code, perhaps I made a mistake, or can't Matlab do it?
Under a piece of the code with the solve function, and following my entire function.
Thank you.
Solange
The entrance values
test = Leakage_LabyrinthSeal(15*10^-6, 100+273, 10*10^5, 6*10^5,81.3 ,4.8*10^-3, 0*10^-3, 'R134a', 20*10^-3)
Solve function
function [P1_sol] = Leakage_LabyrinthSeal(Gap, Temp, Pa, Pb,R,Radius_first_tooth, pitch, fluide, leakage_initialization)
S2 = 2*pi*(Radius_first_tooth+pitch)*Gap iteration = 1;
if Pa > Pb
Pin = Pa;
Pout = Pb;
S3 = 2*pi*Radius_first_tooth*Gap;
S1 = 2*pi*(Radius_first_tooth+2*pitch)*Gap;
else
Pin = Pb;
Pout = Pa;
S1 = 2*pi*Radius_first_tooth*Gap;
S3 = 2*pi*(Radius_first_tooth+2*pitch)*Gap;
end
Gamma= refpropm('K', 'T', Temp,'P', Pin*10^-3,fluide)
mdot_1(1) = leakage_initialization;
for i=1:iteration
syms P1
P1_sol = solve( 0==S1^2*Pin^2*(2*Gamma/((Gamma-1)*R*Temp))*(P1/Pin)^(2/Gamma)*(1-(P1/Pin)^((Gamma-1)/Gamma))-mdot_1(i)^2, P1)
P1(i) = P1_sol;
end
Entire code
%leakage massflow on a labyrith seal in axial orientation %Sharp knif model, Tuyère de Laval %Take in consideration the direction on the flow, P2>P1 --> negativ value % Unit [SI] % Radius = position of the seal % Gap = width between rotor and knif % Gamma = Cp/Cv, and r_cst = R/M, fluid properties % Temp = temperature of the fluid around the seal % Pa, Pb = Pressure on both side of the seal, Pa highest radius
% pour trois dents
function [MassFlow, Erreur] = Leakage_LabyrinthSeal(Gap, Temp, Pa, Pb,R,Radius_first_tooth, pitch, fluide, leakage_initialization)
S2 = 2*pi*(Radius_first_tooth+pitch)*Gap
iteration = 1;
if Pa > Pb
Pin = Pa;
Pout = Pb;
S3 = 2*pi*Radius_first_tooth*Gap;
S1 = 2*pi*(Radius_first_tooth+2*pitch)*Gap;
else
Pin = Pb;
Pout = Pa;
S1 = 2*pi*Radius_first_tooth*Gap;
S3 = 2*pi*(Radius_first_tooth+2*pitch)*Gap;
end
Gamma= refpropm('K', 'T', Temp,'P', Pin*10^-3,fluide)
mdot_1 = zeros(1,iteration);
mdot_2 = zeros(1,iteration);
mdot_3 = zeros(1,iteration);
P_1 = zeros(1,iteration);
P_2 = zeros(1,iteration);
E = zeros(1,iteration);
mdot_1(1) = leakage_initialization;
for i=1:iteration
syms P1
P1_sol = solve( 0==S1^2*Pin^2*(2*Gamma/((Gamma-1)*R*Temp))*(P1/Pin)^(2/Gamma)*(1-(P1/Pin)^((Gamma-1)/Gamma))-mdot_1(i)^2, P1)
P_1(i) = P1_sol;
mdot_2(i) = mdot_1(i);
mdot2=mdot_2(i);
P2_sol = vpasolve( mdot2^2 == S2^2* P1_sol^2*(2*Gamma/((Gamma-1)*R*Temp))*(P2/P1_sol)^(2/Gamma)*(1-(P2/P1_sol)^((Gamma-1)/Gamma)),P2);
P_2(i) = P2_sol;
PL = P_2(i).*(2/(Gamma+1)).^(Gamma/(Gamma-1));
if Pout > PL
Pout = Pout;
else
Pout = PL;
end
mdot_3(i) = S3.* P_2(i).*sqrt((2*Gamma/((Gamma-1)*R.*Temp)).*(Pout./P_2(i)).^(2/Gamma).*(1-(Pout./P_2(i)).^((Gamma-1)/Gamma)));
E(i) = (mdot_3(i)-mdot_1(i))./mdot_3(i);
if(abs(E(i))<0.001)
MassFlow=mdot_3(i);
break;
end
mdot_1(i+1) = mdot_1(i) + (mdot_3(i)-mdot_1(i))./2;
end
Erreur = E(iteration);
MassFlow = mdot_3(iteration);