Efficient implementation of cubic equation of state

Question

...In this case the Peng Robinson Equation of State. Equations of state are empirical equations with parameters derived from experimental data. They are used to predict pure component and mixture properties such as compressibility, fugacity, and mixture equilibriums. See Wikipedia. Although I'm dealing with a Peng-Robinson EOS here, the algorithm can be applied to any other cubic EOS such as Van der Waals or Redlich–Kwong.

The equation is defined as such:

$p=\frac{R T}{V_m-b}-\frac{a}{-b^2+2 b V_m+V_m^2}$

where

$a=\frac{0.457235 R^2 T_c^2}{p_c} \alpha(T)$

$b = 0.07780 \frac{R T_c}{p_c}$

$\alpha(T) = [1+\kappa(1-(\frac{T}{T_c})^\frac{1}{2})]^2$

$\kappa = 0.37464+1.54226\omega-0.26992\omega^2$.

Translated into polynomial form:

$A = \frac{a p}{R^2 T^2}$ and $B = \frac{b p}{R T}$ which leads to the cubic form:

$Z^3-(1-B)Z^2+(A-2B-3B^2)Z-(A B - B^2 -B^3)=0$

When dealing with multi-component mixtures, Van der Waals mixing rules can be implemented as such:

$amix =\sum_{i=1}^N \sum_{j=1}^N x_i x_j Aij $

$bmix = \sum_{i=1}^N x_i b_i$

where

$Aij = Aji = (a_i a_j)^\frac{1}{2} (1-kij)$

The experimental parameters provided from literature are $T_c$ (critical temperatures), $P_c$ (critical pressures), $R$ (universal gas constant), $\omega$ (acentric factor), and $kij$ (binary interaction parameters). The component molar compositions are $x_i$ and is $N$ in length.

When solving the cubic equation there will, of course, be three roots. The max real root is the vapor root, and if a liquid phase is present, the min real root is the liquid root. The middle root is a fake root. If there is no liquid phase present, the vapor root will be the only real root.

What is the most efficient (both in terms of memory and speed) way to implement this algorithm to calculate the vapor z-factor?

Answer 1

Although this may be a little esoteric of a question, I can't begin to tell you how many times I've searched the web for a good implementation of an EOS. Engineers (especially petroleum and chemical engineers) use them often.

Here is my implementation for a three component system: (1-methane, 2-ethane, 3-carbon dioxide):

Define experimental parameters in initialization cell:

xi = {0.966294, 0.015073, 0.018633}; (*mol percent methane, ethane, carbon dioxide*)
Tc = {190.4, 305.4,304.1}; (*critical temps in K*)
Pc = {46, 48.8, 73.8}; (*critical pressures in bar*)
\[Omega] = {0.011, 0.099, 0.239}; (*acentric factors*)
kij = {{0, -0.003, 0.09},{-0.003, 0, 0.13},{0.09, 0.13, 0}}; (*binary interaction parameters*)
rgas = 0.00008314; (*universal gas constant*)
comps = Length[xi]; (*number of components*)

Set up polynomial coefficients:

kappa = 0.37464 + 1.54226*\[Omega] - 0.26992*\[Omega]^2;    
PREOSpoly[T_, P_, x_] :=(
  alpha = Table[(1 + kappa[[i]]*(1 - (T/Tc[[i]])^.5))^2, {i, 1, comps}];
  aterm = Table[0.45724*rgas^2*Tc[[i]]^2/Pc[[i]]*alpha[[i]], {i, 1, comps}];
  amix = Sum[xi[[i]]*xi[[j]]*(aterm[[i]]*aterm[[j]])^.5*(1 - kij[[i, j]]), {i, 1, comps}, {j, 1, comps}];
  bmix = Sum[xi[[i]]*(0.0778*rgas*Tc[[i]]/Pc[[i]]), {i, 1, comps}];
  AA = amix*P/(rgas^2*T^2);
  BB = bmix*P/(rgas*T);
  x^3 - (1 - BB)*x^2 + (AA - 2*BB - 3*BB^2)*x - (AA*BB - BB^2 - BB^3)
);

Find vapor root of the cubic polynomial:

zfactor[T_, P_] := root /. FindRoot[PREOSpoly[T, P, root], {root, 1}];