how to solve non-linear equations in a nozzle

Dear friends

I want to solve these non-linear equations for inlet of a de-laval nozzle(8 equations and 8 variables). All the equations are both non-linear and parametric. I tried to use fsolve to solve them but I got many errors. Thank you for guiding me in this regard.

inputs:

di = 0.04; % k

alpha = 12.67; % degree

Ti = 291.65; % K

Pi = 90; % bar

Vi = 43.5; % m/s

Ci = 443.8; % m/s

Mi = 0.098;

Zi = 1;

rhoi = 60.17; % kg/m^3

k = 1.32;

R = 8.3114; % kj/kmol.K

Ai = 0.0013; % m^2

equations:

F(1) = Pt/rhot - Zt*R*Tt;

F(2) = sqrt(k*R*Zt*Tt) - Vt;

F(3) = (k*R*Zt*Tt) - (Vi^2+((2*k)/(k-1))*(Zi*R*Ti-Zt*R*Tt));

F(4) = ((Zi*Ti)/(Zt*Tt)) - (rhoi/rhot)^(k-1);

F(5) = (rhoi/rhot)^(k-1) - (Pi/Pt)^((k-1)/k);

F(6) = Ai*(rhoi/rhot)*(Vi/Vt) - At;

F(7) = (Ai^0.5 - At^0.5)/((pi^0.5)*tand(alpha)) - xt;

output:

[ Pt , Tt , rhot , Zt , Vt , At , xt ]

Answer 1

[Torsten]

Maybe you can give better initial guesses in x0 for the solution than I can ...

    x0 = 10*rand(7,1);
options = optimset('MaxFunEvals',1000000,'MaxIter',1000000);
norm(fun(x0))

ans = 1.8557e+04

x = fsolve(@fun,x0,options)

Solver stopped prematurely.

fsolve stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 1.000000e+06.

x = 7×1

10.9152 12.7409 0.0048 21.5998 6.5697 8.0980 6.6604

norm(fun(x))

ans = 115.4423

function F = fun(x)
Pt = x(1);
Tt = x(2);
rhot = x(3);
Zt = x(4);
Vt = x(5);
At = x(6);
xt = x(7);
di = 0.04; % k
alpha = 12.67; % degree
Ti = 291.65; % K
Pi = 90; % bar
Vi = 43.5; % m/s
Ci = 443.8; % m/s
Mi = 0.098; 
Zi = 1;
rhoi = 60.17; % kg/m^3
k = 1.32;
R = 8.3114; % kj/kmol.K
Ai = 0.0013; % m^2
%equations:
F(1) = Pt/rhot - Zt*R*Tt;
F(2) = sqrt(k*R*Zt*Tt) - Vt;
F(3) = (k*R*Zt*Tt) - (Vi^2+((2*k)/(k-1))*(Zi*R*Ti-Zt*R*Tt));
F(4) = ((Zi*Ti)/(Zt*Tt)) - (rhoi/rhot)^(k-1);
F(5) = (rhoi/rhot)^(k-1) - (Pi/Pt)^((k-1)/k);
F(6) = Ai*(rhoi/rhot)*(Vi/Vt) - At;
F(7) = (Ai^0.5 - At^0.5)/((pi^0.5)*tand(alpha)) - xt;
end

[M.S.R]

Thanks for your reply.

You solved one of my most common problems on how to write the fsolve function !

But my answers won't converge to the desired answers unless I enter the initial guesses very precisely, which is not possible in solving a real problem!

Do you know a way around this?

The output is:

Pt = 49.42; % bar

Tt = 253.97; % k

rhot = 37.94; % kg/m^3

Zt = 1;

Vt = 414.2; % m/s

dt = 0.0163; % m

xt = 0.0527; % m

[Torsten]

Solve equation (3) for Zt*Tt.

Insert Zt*Tt in equation (2) and solve for Vt.

Insert Zt*Tt in equation (4) and solve for rhot.

Insert rhot in equation (5) and solve for Pt.

Insert rhot and Vt in equation (6) and solve for At.

Insert At in equation (7) and solve for xt.

Equation (1) is either right or wrong when you insert Pt, rhot and Zt*Tt from above - at least you cannot solve for Zt and Tt separately.

This shows that your system is either inconsistent or underdetermined.

syms ZtTt rhot Pt At xt Vt 
di = 0.04; % k
alpha = 12.67; % degree
Ti = 291.65; % K
Pi = 90; % bar
Vi = 43.5; % m/s
Ci = 443.8; % m/s
Mi = 0.098; 
Zi = 1;
rhoi = 60.17; % kg/m^3
k = 1.32;
R = 8.3114; % kj/kmol.K
Ai = 0.0013; % m^2
F(1) = Pt/rhot - ZtTt*R;
F(2) = sqrt(k*R*ZtTt) - Vt;
F(3) = (k*R*ZtTt) - (Vi^2+((2*k)/(k-1))*(Zi*R*Ti-ZtTt*R));
F(4) = ((Zi*Ti)/(ZtTt)) - (rhoi/rhot)^(k-1);
F(5) = (rhoi/rhot)^(k-1) - (Pi/Pt)^((k-1)/k);
F(6) = Ai*(rhoi/rhot)*(Vi/Vt) - At;
F(7) = (Ai^0.5 - At^0.5)/(sqrt(Pi)*tand(alpha)) - xt;
ZtTt1 = double(simplify(solve(F(3),ZtTt)))

ZtTt1 = 275.2123

Vt1 = double(simplify(solve(subs(F(2),ZtTt,ZtTt1),Vt)))

Vt1 = 54.9488

rhot1 = double(simplify(solve(subs(F(4),ZtTt,ZtTt1),rhot)))

rhot1 = 50.1936

Pt1 = double(simplify(solve(subs(F(5),rhot,rhot1),Pt)))

Pt1 = 70.8462

At1 = double(simplify(solve(subs(F(6),[rhot,Vt],[rhot1,Vt1]),At)))

At1 = 0.0012

xt1 = double(simplify(solve(subs(F(7),At,At1),xt)))

xt1 = 4.3680e-04

%Results should be 0, but isn't 
%This shows inconsistency of your system of equations
double(simplify(subs(F(1),[Pt rhot,ZtTt],[Pt1,rhot1,ZtTt1]))) 

ans = -2.2860e+03

[Alex Sha]

I think Torsten is right， there is no exact numerical solution， the approximate solution is:

pt: 596435.381381668
tt: -21.0901443320853
rhot: 260.748232462543
zt: -13.0493323055227
vt: 54.9487704256482
at: 0.000237481729284569
xt: 0.00968011270311718