Empty system, 1 equation 1 unknow with solve

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;

    %Laval Pressur

    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);

[Solange Cevat]

I solve my problem with a for iteration. I initialize my pressure, calculate the mass flow. I change the pressure until I become the right mass flow.

Code

    P1 = Pin-(Pin-Pout_in)*1/10;
    for p=1:100
        mdot1_p = S1*Pin*sqrt((2*Gamma/((Gamma-1)*R*Temp))*((P1/Pin)^(2/Gamma)-(P1/Pin)^((Gamma+1)/Gamma)));
        Erreur_mdot1_p = (mdot1_p-mdot_1)/mdot_1;
        if abs(Erreur_mdot1_p) < 0.005
            break
        end
        P1 = P1+ sign(Erreur_mdot1_p)* ((Pin-Pout_in)*1/1000);
    end

Answer 1

[Roger Stafford]

The function 'solve' is intended to find symbolic solutions to equations. There is no guarantee that it will succeed. If it fails, you should use 'fzero' for one unknown and one equation to get a numeric solution. All other quantities should be known. For more than one unknown and one equation, use 'fsolve'.

[Solange Cevat]

I tried fzero and fsolve but with both functions it didn't work. for fzero it became a complex number and with fsolve, it said that the initial point is a solution. I tried with both functions to change the initial value. following you can see how I implement my functions:

Code

    P1_sol = fsolve ( @(P1) S1^2*Pin^2*(2*Gamma/((Gamma-1)*R*Temp))*(P1/Pin)^(2/Gamma)*(1-(P1/Pin)^((Gamma-1)/Gamma))-mdot_1(i)^2,20*10^-3);

    P1_sol = fzero ( @(P1) S1^2*Pin^2*(2*Gamma/((Gamma-1)*R*Temp))*(P1/Pin)^(2/Gamma)*(1-(P1/Pin)^((Gamma-1)/Gamma))-mdot_1(i)^2,20*10^-3);

Answer 2

[Michael Haderlein]

- What is the value of S1, Gamma, R, Temp, Pin and mdot_1? It will strongly depend on these values if there is a solution.

- If you can estimate the order of magnitude of the solution, you can plot the function with

fun=@(P1) ...
  S1^2*Pin^2*(2*Gamma/((Gamma-1)*R*Temp))...
  *(P1/Pin).^(2/Gamma)...
  .*(1-(P1/Pin).^((Gamma-1)/Gamma))-...
  mdot_1^2;
ezplot(fun,[1e5,1e6])

- or get the values of the function with

y=feval(fun,1e5:1e6);

I just tried it with fsolve and S1=1; Gamma=2; R=8.314; Temp=300; Pin=2e5; mdot_1=1;

and got quickly the solution P1_sol=0.0031. However, I guess P1_sol is some pressure and should be in the range of 1 bar resp. 1e5 Pa. So you're not interested in the solution at 0.0031 Pa but in a solution in the range of 1e5. If you initialize with, for instance, 1e5, you'll get the result 2e5 which maches beautifully with the outcome of ezplot. So, if you have multiple roots, the solution will depend on the initial guess.

- One remark on your equation: Given that mdot_1 is small, there are the two roots P1=0 and P1=Pin. Depending on which root you want to find, this estimation should well serve as initial guess.

[Solange Cevat]

My values are much smaller, I tried ezplot but I have just negative value. So I tried the same equation to found a solution with possible value P1 = 8e5. And I found MassFlow = 0.0014. So there is a solution. The value I used are given just under. And the entire code in the question above. But I don't undertand why it doesn't work, perhaps I define bad my function.

Thank you very much.

Value

Leakage_LabyrinthSeal(15*10^-6, 100+273, 10*10^5, 6*10^5,81.3 ,4.8*10^-3, 0*10^-3, 'R134a', 10*10^-3) Leakage_LabyrinthSeal(Gap, Temp, Pa, Pb,R,Radius_first_tooth, pitch, fluide, leakage_initialization)

Code

fun=@(P1)... S1^2*Pin^2*(2*Gamma/((Gamma-1)*R*Temp))*(P1/Pin).^(2/Gamma).*(1-(P1/Pin).^((Gamma-1)/Gamma))- mdot_1^2;

ezplot(fun,[1e5,1e6])

P1_test = 8*10^5;

mdot2 = S1^2* Pin^2*(2*Gamma/((Gamma-1)*R.*Temp))*(P1_test/Pin)^(2/Gamma)*(1-(P1_test/Pin)^((Gamma-1)/Gamma));

MassFlow = sqrt(mdot2);

[Michael Haderlein]

Still Gamma is not known.

[Solange Cevat]

I could write it 1.1462 and is calculate in the code given in the question. Sorry.

[Michael Haderlein]

Well, I don't get any positive value. There's a maximum at 5.75e5, but its value is -9.73e-5. What makes you think that there's a solution at 8e5? You should check if there's a typo in your equations, maybe a bracket closed at the wrong position.