fzero instead of Newton-Raphson method in flash calculations!

3 views (last 30 days)
Dennis Qadri
Dennis Qadri on 23 Oct 2013
Answered: MEL on 7 May 2015
Hi
I am trying to do some flash calculations using Newtons method to define alpha.
But if i change the pressure i need to update the guess value to something very close (+/- 0.1), this is not allways possible! Is there a way to the following in another way like using fzero?
%Step two Newton-Raphson--------------------------------------------------%
alpha = 0.1; %Initial guess
epsilon2 = 1.0; %initial errors to start loops
while (epsilon2 >= 1.e-12); %Newton-Raphson while loop
h = 0; dh = 0; %Starting value for h and dh
for i=1:nc
h = z(i)*(K(i)-1)/(1+alpha*(K(i)-1))+h; %Iteration of alpha
dh = -(z(i)*(K(i)-1)^2/(1+alpha*(K(i)-1))^2+dh);%Iteration for next alpha
end
alphan = alpha - (h / dh); %Update alpha
epsilon2 = abs( (alphan-alpha)/alpha ); %Update error value
alpha = alphan; %Update alpha
end
  4 Comments
Show 2 older comments

Sign in to comment.

Sign in to answer this question.

Answers (2)

Dennis Qadri
Dennis Qadri on 23 Oct 2013
Edited: Dennis Qadri on 23 Oct 2013
This is going to be hard to understand!
%Start values*************************************************************%
%Operating temperature [Kelvin]
T = 363.15;
%Operating pressure [Atm]
P = 44.411547; %Lader til at udgaven på nettet bruge reduceret tryk!!
%Ideal gas universal constant [atm cc/mol K]
R = 82.0600;
%-------------------------------------------------------------------------%
%Feed mixture
z = [0.229 0.679 0.092];
z=z/sum(z)
%Critical pressures [Atm]
Pc = [4604*0.00986923267 4880*0.00986923267 4249*0.00986923267];
%Critical temperatures [Kelvin]
Tc = [190.6 305.3 369.8];
%Acentric factors
w = [0.0126 0.0978 0.1541];
%w = [0.013 0.099 0.150];
%-------------------------------------------------------------------------%
nc = length(z); %Lengt of z
Kij(1:nc,1:nc)=0; %Fitting parameter known as the coupling parameter
%Values for Peng-Robinson-------------------------------------------------%
for i=1:nc
Tr(i) = T/Tc(i); Pr(i) = P/Pc(i); %Reduced T and P
%Physical constants for Peng-Robinson EOS
a(i) = (0.457235*(R*Tc(i))^2)/Pc(i); %a from Peng-Robinson
b(i) = (0.077796*R*Tc(i))/Pc(i); %b from Peng-Robinson
ka(i) = 0.37464+1.54226*w(i)-0.26992*w(i)^2; %kappa from Peng-Robinson
al(i) = (1+ka(i)*(1.0 - sqrt(Tr(i)))^2); %alpha from Peng-Robinson
end
%*************************************************************************%
%Flash calculations*******************************************************%
%Values for Peng-Robinson-------------------------------------------------%
for i=1:nc
K(i) = (Pc(i)/P)*exp(5.37*(1+w(i))*(1-Tc(i)/T)) % k-values
end
%Step two Newton-Raphson--------------------------------------------------%
alpha = 0.1; iter = 0; %Initial guess
episilon1 = 1.0; epsilon2 = 1.0; %initial errors to start loops
%while (episilon1 >= 1.e-12); %K-value while loop
while (epsilon2 >= 1.e-12); %Newton-Raphson while loop
h = 0; dh = 0; %Starting value for h and dh, to start loop
for i=1:nc
h = z(i)*(K(i)-1)/(1+alpha*(K(i)-1))+h;%Iteration of alpha
dh = (z(i)*(K(i)-1)^2/(1+alpha*(K(i)-1))^2+dh);% next alpha
if (dh<0.0); dh = 1.0e-05;
end
dh=-dh;
alphan = alpha - (h / dh); %Update alpha
epsilon2 = abs( (alphan-alpha)/alpha ); %Update error value
alpha = alphan; %Update alpha
end
end
%Step three xi, yi--------------------------------------------------------%
for i=1:nc %For loop for xi and yi
x(i) = z(i)/(1.0 + alpha*(K(i)-1.0)); %Calculate xi
y(i) = K(i)*x(i); %calculate yi
if (x(i)<0.0); x(i) = 1.0e-05; end; %Stop xi from being <0
if (y(i)<0.0); y(i) = 1.0e-05; end; %Stop yi from being <0
end
%Step four Finding phiV and phiL------------------------------------------%
%Mixing and combining rules-----------------------------------------------%
AmixtureL = 0; BmixtureL = 0; %Initial guess for mixture constants
AmixtureV = 0; BmixtureV = 0; %Initial guess for mixture constants
BB(1) = 0;
for i=1:nc
A(i) = 0.45724*(a(i)*al(i))*(Pc/Tc);
B(i) = 0.07780*(Pc/Tc);
end
for i=1:nc %Loop for real values of mixture constants
for j=1:nc
Aij(i,j) = (1-Kij(i,j))*(A(i)*A(j))^0.5;
AmixtureL = x(i)*x(j)*Aij(i,j)+AmixtureL; %Mixture constant AmixtureL
AmixtureV = y(i)*y(j)*Aij(i,j)+AmixtureV; %Mixture constant AmixtureV
end
BmixtureL = x(i)*B(i)+BmixtureL; %Mixture constant BmixtureL
BmixtureV = y(i)*B(i)+BmixtureV; %Mixture constant BmixtureV
end
%Calculating the compressibillity [Z]-------------------------------------%
Avsum = zeros(1,nc) ; Alsum = zeros(1,nc); %Vector of zeros
for i=1:nc
for j=1:nc
Alsum(i) = Alsum(i)+x(j)*Aij(i,j);%Summation compressibility equation L
Avsum(i) = Avsum(i)+y(j)*Aij(i,j);%Summation compressibility equation V
end
end
Zl=roots([1 -1+BmixtureL AmixtureL-3*BmixtureL^2-2*BmixtureL -AmixtureL*BmixtureL+BmixtureL^2+BmixtureL^3]);
Zl=min(Zl(imag(Zl)==0)); % min of real values in vector
Zv=roots([1 -1+BmixtureV AmixtureV-3*BmixtureV^2-2*BmixtureV -AmixtureV*BmixtureV+BmixtureV^2+BmixtureV^3]);
Zv=max(Zv(imag(Zv)==0)); % max of real values in vector
%-------------------------------------------------------------------------%
phiL=exp((Zl-1).*B/BmixtureL-log(Zl-BmixtureL)-AmixtureL/(2*sqrt(2)*BmixtureL)*log((Zl+(1+sqrt(2))*BmixtureL)/(Zl+(1-sqrt(2))*BmixtureL)).*(2.*Alsum./AmixtureL-B./BmixtureL))
phiV=exp((Zv-1).*B/BmixtureV-log(Zv-BmixtureV)-AmixtureV/(2*sqrt(2)*BmixtureV)*log((Zv+(1+sqrt(2))*BmixtureV)/(Zv+(1-sqrt(2))*BmixtureV)).*(2.*Avsum./AmixtureV-B./BmixtureV))
MEL
MEL on 7 May 2015
Dear Dennis,
I used your code some time ago for a flash calculation. I think that you meant (in your comment) the following line:
while (episilon1 >= 1.e-12); %K-value while loop
instead of:
%while (episilon1 >= 1.e-12); %K-value while loop
because you have two "end"s after... However, I concluded that the use of two nested cycle is useless if you calculate the new K with a separate function that takes the current composition to calculate the new "K" inside the Newton-Raphson while. Tell me your opinion about this:
- you estimate the first K
- you calculate alpha
- you calculate a new K with the new composition (given by the new alpha)
and so iteratively until you achieve the convergent values of alpha and K.
Of course, this is valid if you calculate K by means of fugacities. If your K is indipendent from your current composition (for example taking K = vapur pressure_i /total pressure) my, but even more yours, argument is completely useless.
Thanks!

Categories

Find more on Newton-Raphson Method in Help Center and File Exchange

Tags

Products

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!