% Calculate the vapor Pressure of n-heptane using SRK-EOS(adsbygoogle = window.adsbygoogle || []).push({});

% 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

**Physics Forums - The Fusion of Science and Community**

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Seek Matlab Help in this problem

Loading...

Similar Threads - Seek Matlab Help | Date |
---|---|

How to program this in Matlab | Feb 28, 2018 |

Matlab Making a short test for self-adjointness | Feb 23, 2018 |

Matlab Plotting Coordinate Transformations in Matlab | Feb 12, 2018 |

Matlab MATLab: Not enough inputs for nlinfit | Nov 16, 2017 |

**Physics Forums - The Fusion of Science and Community**