I'm trying to use the Peng-Robinson-Stryjek-Vera 2 (PRSV2) Equation of State to calculate the density of a compound, in this example: methanol. At ambient temperature and pressure, shouldn't I be able to solve the cubic equation of state for the molar volume and calculate a density corresponding reasonably well to the density of the liquid given in the literature at those conditions of temperature and pressure? The values are considerably off, at least farther off than I would expect from such a well tested equation of state. For instance, for methanol, I obtain a density of 0.672 g/mL rather than the literature value of 0.791 g/mL. For acetonitrile, it's even worse: 0.480 g/mL rather than 0.786 g/mL. Does anyone have experience with such cubic equations of state? Are they simply not that accurate? Am I doing something wrong? I've checked my implementation numerous times for typos, but cannot find any.
This implementation in MATLAB comes from the following two papers:
Stryjek, R. and J. H. Vera, (1986). "PRSV: An Improved Peng-Robinson Equation of State for Pure Compounds and Mixtures", Can. J. Chem. Eng., 64, 323-333.
Stryjek, R. and J. H. Vera, (1986). "PRSV2: A Cubic Equation of State for Accurate Vapor-Liquid Equilibria Calculations", Can. J. Chem. Eng., 64, 820-826.
T = 298; % Temperature, in K
P = 1; % Pressure, in atm
% Values for methanol
M = 32.04; % Molar mass, in g mol^-1
T_c = 512.58; % Critical temperature, in K
P_c = 8095.79; % Critical pressure, in kPa
omega = 0.56533; % Acentric factor, from PRSV2
kappa1 = -0.16816; % From PRSV
kappa2 = -1.3400; % From PRSV2
kappa3 = 0.588; % From PRSV2
R = 8.3144598; % Ideal gas constant, in Pa m^3 mol^-1 K^-1
R = R*1000; % Ideal gas constant, in kPa cm^3 mol^-1 K^-1
T_R = T/T_c; % Reduced Temperature
kappa0 = 0.378893 + 1.4897153*omega - 0.17131848*omega^2 + ...
0.0196544*omega^3;
kappa = kappa0 + (kappa1 + kappa2*(kappa3 - T_R)*(1 - (T_R)^0.5))*...
(1 + (T_R)^0.5)*(0.7 - T_R);
alpha = (1 + kappa*(1 - (T_R)^0.5))^2;
a = (0.457235 * R^2 * (T_c)^2 / P_c)*alpha;
b = 0.077796*R*(T_c)/(P_c);
P = P*101.325; % Pressure, in kPa
% Using the PRSV EOS cubic expansion for molar volume
C3 = 1;
C2 = b - R*T/P;
C1 = a/P - 2*R*T*b/P - 3*b^2;
C0 = (b^2 + R*T*b/P - a/P)*b;
cubic = [C3 C2 C1 C0];
v = roots(cubic); % Molar volumes, in cm^3 mol^-1
rho1 = M/v(1) % Density, in g cm^-3
rho2 = M/v(2) % Density, in g cm^-3
rho3 = M/v(3) % Density, in g cm^-3
