Info
This question is closed. Reopen it to edit or answer.
How to solve non-linear equations with several unknowns and several equations?
1 view (last 30 days)
Show older comments
I have a problem with xz. After I run this code I get the same value of "xz" for 500 times.
"xz" should not get the same value. It should vary with P, T_sat and q.
I think the problem is "eps" in the second round cannot get beyond "eps>0.0001" then the rest of "xz" has the same values.
T0=-90;
%Tz=1*ones(500,1);
A0 = 5*10^-6;
z= -500:0;
gamma_0 = 0.4 ;
D =1200;
q0= 150*10^-3;
a= 0.001;
m_fluid = 150;
m_steam = 30;
d = 0.4318;
rho_l= 1000;
rho_g= 0.804;
G= (4*(m_fluid+m_steam))/(pi()*d^2);
u = (G/rho_l);
g=9.81;
alpha0=0.1;
xz(1) = rho_g*alpha0/((rho_l*(1-alpha0))+(rho_g*alpha0));
for i=2:501
eps=1;
Tz(i) = -(1/a)*((1+a*T0)*exp(((a*A0*D.^2/gamma_0)*(1-exp(-z(i)/D))+((a*q0*z(i))/(gamma_0+gamma_0*a*T0)))-((a*A0*D*z(i))/gamma_0))-1);
sigma(i)=((-7.5*(10^(-12)))*Tz(i).^4)+(5.8*(10^(-9))*(Tz(i).^3))-1.8*(10^(-6))*(Tz(i).^2)-(1.9*(10^(-5))*Tz(i))+0.07;
h_lg(i) = ((-5.42*(10^(-5)))*Tz(i).^4)+((0.013)*(Tz(i).^3))-4.5*(Tz(i).^2)-(1.9*(10^(3))*Tz(i))+2.5*(10^(6));
Cp(i)= (0.0001*Tz(i).^3)-(0.027*Tz(i).^2)+ 3.9*Tz(i)+(4*10.^3);
k(i)=((-1.3853*(10^(-10)))*Tz(i).^4)+(8.7543*(10^(-8)))*(Tz(i).^3)-(2.43*(10^(-5))*(Tz(i).^2)-(0.0030633)*Tz(i))+0.54455;
Tz_k(i)=Tz(i)+273.15;
mu_l(i)= 2.414*10.^-5*10.^(247.8/(Tz_k(i)-140));
mu_g(i) = (18.27*10.^-6)*((291.15+120)/(Tz_k(i)+120))*((Tz_k(i)/291.15).^(3/2));
Pr(i)= ((Cp(i)*mu_l(i))/k(i));
Tz_new=Tz(i);
sigma_new=sigma(i);
h_lg_new=h_lg(i);
Pr_new=Pr(i);
k_new=k(i);
mu_l_new=mu_l(i);
mu_g_new=mu_g(i);
xz_new=xz(i-1);
q_new = 150*10^-3;%Initial condition
h_new = 1*10^5; %initial condition it's here$ the inital condition is wrong
while eps>0.0001
xz_ref= xz_new;
sigma_ref=sigma_new;
h_lg_ref= h_lg_new;
Tz_ref=Tz_new;
Pr_ref=Pr_new;
mu_l_ref=mu_l_new;
mu_g_ref=mu_g_new;
k_ref=k_new;
q_ref=q_new;
h_ref=h_new;% it copy here
%%Pressure%%
%%Firction Pressure%%
rho_x = ((xz_ref/rho_g)+((1-xz_ref)/rho_l)).^-1;
WE = ((rho_l^2)*(u^2)*d)/(g*sigma_ref*rho_x);
Fr=((rho_l^2)*(u^2))/(g*(rho_x.^2));
H=((rho_l/rho_g).^0.91)*((mu_g_ref/mu_l_ref).^0.19)*((1-(rho_g/rho_l)).^0.7);
F=(xz_ref.^0.78)*((1-(xz_ref.^2)).^0.24);
Re_l =((G*d)/mu_l_ref);
Re_g =((G*d)/mu_g_ref);
fl =0.079/(Re_l.^0.25);
fg = 0.079/(Re_g.^0.25);
E=(1-(xz_ref.^2))+((xz_ref.^2)*((rho_l*fg)/(rho_g*fl)));
frict2 = E+((3.24*F*H)/((Fr.^0.045)*(WE.^0.035)));
PL = 4*fl*(-z(i)/d)*(G.^2)*(1/(2*rho_l));
P_frict = PL*frict2;
alpha = (1+0.28*((1+xz_ref)/xz_ref)^0.64*(rho_g/rho_l)*(mu_l_ref/mu_g_ref).^0.07)^(-1);
rho_h = (rho_l*(1-alpha))+(rho_g*alpha);
P_static = rho_h*g;
P_mom = (G^2)*((((1-xz_ref).^2)/(rho_l*(1-alpha)))+ ((xz_ref).^2)/(rho_g*alpha))-(G^2)*((((1-xz(i-1)).^2)/(rho_l*(1-alpha)))+ ((xz(i-1).^2)/(rho_g*alpha)));
P= P_frict+P_static-P_mom;
Tsat = -1.8*10.^(-24)*P.^4+2*10.^(-17)*P.^3-8.1*10.^(-11)*P.^2+0.00016*P+88;
h_l = 0.023*(k_ref/d)*((G*(1-xz_ref)*d)/(mu_l_ref).^(0.8))*(Pr_ref.^(1/3));
q_new = h_ref*(Tz_ref-Tsat);
h_new = h_l*((1+(3000*((q_ref/G*h_lg_ref).^0.86)+(((xz_ref/(1-xz_ref)).^(3/4))*((rho_l/rho_g).^0.41)))));
xz_new= xz_ref+((4*q_ref)/(d*G*h_lg_ref))*(-z(i-1)+z(i));
eps = abs(xz_ref - xz_new);
end
xz(i) = xz_new;
Tz(i) = Tz_new;
sigma(i)=sigma_new;
%Tz_k=Tz_k_new;
mu_l(i)= mu_l_new;
mu_g(i) = mu_g_new;
h_lg(i)= h_lg_new;
Pr(i)=Pr_new;
q=q_new;
h=h_new;
end
Thank you
2 Comments
John D'Errico
on 9 Mar 2015
It is a really bad idea to name a variable eps. This is a useful function in MATLAB.
For this code though, only you know what it is doing. Learn to use the debugger.
Answers (0)
See Also
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!