View Single Post
 P: 5 % Calculate the vapor Pressure of n-heptane using SRK-EOS % P = RT / V - b - alpha*a / V(V-b) R = 8.134e-3; % MPa m3 / mole.K n = 0; n = n+1; for i=1:5 Pn = 1.6714; %(Assume Pressure)> this needs updating Tr = 0.9; Tc = 540.13; % K Pc = 2.74; % Mpa; T = Tr*Tc; w = 0.350; m = 0.480+1.574*w -0.176*w^2; al=(1+m*(1-Tr^0.5))^2; a = 0.42748*(((R*Tc)^2) /Pc); b = 0.08664*(R*Tc / Pc); % Z3 - Z2 + (A-B-B^2)Z-AB=0 A = (a*Pn*al)/(R*T)^2; B = (b*Pn)/(R*T); Zv=max(roots([1 -1 A-B-B^2 -A*B])); Zl=min(roots([1 -1 A-B-B^2 -A*B])); disp(Zv) disp(Zl) phiv(i)=exp((Zv -1)-log(Zv -B)-(A/B) *(log(Zv+B/Zv ))); phil(i)=exp((Zl -1)-log(Zl -B)-(A/B) *(log(Zl+B/Zl ))); Fv=Pn*phiv; Fl=Pn*phil; disp(phiv) disp(phil); disp(Fv); disp(Fl); Error=Fl-Fv; if (Error < 10^-4); Pnew=Pn; else P(n+1) =Pn * (phil/phiv); Pn = P(n+1); end end Basically in this code I am trying to assume a pressure 'Pn' to calculate the results which satisfy this relation (Error<10^-4) if this satisfy's then our assume pressure is right else if not then we use the relation P(n+1) = Pn*(Phil/Phiv) to caluclate another pressure and update the above pressure (previously assume) with this one. I have error with this and need advice please help