Trial and error help is not giving me solution to that (still waiting ....

Question

Nasir Qazi 2012-2-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/30397-trial-and-error-help-is-not-giving-me-solution-to-that-still-waiting

% Calculate Bubble Point and Dew Point SRK-EOS
% P = RT / V - b  -  alpha*a / V(V-b)
%----------------------------------------------
% liquid phase composition
y = [0.86 0.05 0.05 0.02 0.01 0.0095 0.0005]';
% vapor phase composition
x = [0.45 0.05 0.05 0.03 0.01 0.01 0.4]';
% critical pressure
Pc =[4.5947 4.8711 4.2472 3.6397 3.3688 3.1233 1.965]';
% critical temperature
Tc =[190.7389 305.5111 370.0333 425.3444 469.8889 512.7778 644.4444]';
% Accentric Factor
w = xlsread('myfile.xlsx','Sheet1','E2:E8');
%---------------------------------------------
L = 1;
R = 8.314e-3; % MPa m3 / mole.K
T = 200; % K
Tr = T./Tc; % K
P(L) = 5.0; %MPa
m = 0.480+1.574.*w -0.176.*w.^2;
al=1+ m.*(1-sqrt(Tr)).^2;
a = 0.42748.*(((R.*Tc).^2)./Pc);
b = 0.08664.*(R.*Tc./Pc);
q=0;
while q==0
      % Z3 - Z2 + (A-B-B^2)Z-AB=0
      % for gas phase
      SA = 0;
      RR = zeros(7,1);
      for i = 1:7
          for j = 1:7
              SA = SA + y(i).*y(j).*sqrt(a(i).*a(j).*al(i).*al(j));
              RR(i) = RR(i) + y(j).*sqrt(a(i).*a(j).*al(i).*al(j));
          end
      end
      SB = sum(y.*b);
      % For Liquid Phase
      SAA = 0;
      MM = zeros(7,1);
      for i = 1:7
          for j = 1:7
              SAA = SAA + x(i).*x(j)*sqrt(a(i).*a(j)*al(i).*al(j));
              MM(i) = MM(i) +x(j).*sqrt(a(i).*a(j).*al(i).*al(j));
          end
      end
      SBB = sum(x.*b);
      %for Gas Phase
      A = SA.*P(L)./(R*T).^2;
      B = SB.*P(L)./(R*T);
      ZV=roots([1 -1 A-B-B^2 -A*B]);
      Zv=max(ZV);
      Zvc=isreal(Zv);
      if Zvc == false
          Zv=abs(Zv);
          disp('error 1')
      end
      %for liquid Phase
      AA = (SAA.*P(L))./(R*T).^2;
      BB = (SBB.*P(L)./(R*T));
      Z=roots([1 -1  AA-BB-BB^2 -AA*BB]);
      Zl=min(Z);
      Zvl=isreal(Zl);
      if Zvl == false
          Zl=abs(Zl);
          disp('error 2')
      end
      phiv =exp((b.*(Zv-1)/SB) - log(Zv-B)-(A/B).*((2.*RR/SA)-
(b/SB)).*log(1 + B/Zv));
      phil =exp((b.*(Zl-1)./SBB) - log(Zl-BB)-(AA./BB).*((2.*MM/SAA)-
(b./SBB)).*log(1 + BB./Zl));
      K = phil./phiv;
      RA1= sum(y./K)-1;
      if abs(RA1) < 1e-5 && abs((y./K) - x) < 1e-5
          break
      else
          if L==1
              L = L+1;
              P(L) = 1.01.*P(L-1);
              RA1;
          else
              if L==1
                  P(L)=P(L-1)-(RA1./(RA1-RA2))*(P(L-1)-P(L-2));
                  RA2 = RA1;
              end
              %adjustment
              x2 = (1-0.5)*(y./K) + 0.5*x;
              x = x2;
          end
      end
end
disp(P(L));

% I am trying to using trial and error to get the results of P(L) , but my solution is not converging and keep going on. I need to know what can I do to resolve the issue .

Help Please