# Efficient implementation of cubic equation of state

...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?

• Note that for cubic polynomials with all roots real, Root[]'s numbering system is somewhat clear: Root[poly, 1] is the smallest root, and Root[poly, 3] is the biggest one. – J. M.'s discontentment Aug 22 '12 at 15:50
• Okay, so since this mixture input has very low boiling points, we're not seeing any liquid phase, therefore, there is only one real root (the vapor root). What I want FindRoot to do is start searching at a value of 1. – kale Aug 22 '12 at 15:53
• In the case of you seeking the "vapor root", then yes, the ideal gas approximation as the starting point is a logical choice... – J. M.'s discontentment Aug 22 '12 at 15:54

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*)
ω = {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*ω - 0.26992*ω^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}];

• Damn, this is actually some of the stuff I'm currently working on... :o I use NSolve[] myself instead of FindRoot[], as the equation is a polynomial, after all. EOSs with transcendental parts, like Dieterici, are where you'll need the services of FindRoot[]... – J. M.'s discontentment Aug 22 '12 at 15:48
• I do believe the code would be more efficient if you leverage dot products (via Dot[]) to compute your sums... – J. M.'s discontentment Aug 22 '12 at 15:58
• I chose FindRoot because I wanted it to only return one root near x_0. However if NSolve can be implemented in a quicker manner, I'm all game. Feel free to post an answer showing me shortcuts. I really want to see how inefficient my code is. – kale Aug 22 '12 at 16:03