PNG 520
Phase Relations in Reservoir Engineering

Solution Techniques for Cubic Expressions and Root Finding

PrintPrint

Hopefully, we have convinced you that the use of cubic equations of state can represent a very meaningful and advantageous way of modeling the PVT behavior of petroleum fluids. What we need now are the tools that will allow us to get the information that we want out of them. Even though cubic equations of state are explicit in pressure, pressure is not the common unknown to be calculated in the typical problem. In the most common problem, pressure and temperature are known and we want either molar volume (or its reciprocal, molar density) or compressibility factor (the most likely case). Therefore, we are faced very often with the need to solve for the roots of a cubic expression. Here we present a number of approaches that may be followed.

Analytical scheme

Given the cubic polynomial with real coefficients: x3 + ax2 + bx + c = 0, the first step is to calculate the parameters:

Q= a 2 3b 9  and R= 2 a 3 9ab+27c 54 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.5)

Now let M = R2 – Q3 be the discriminant. We then consider the following cases:

  1. If M < 0 (R2 < Q3), the polynomial has three real roots. For this case, compute θ=arccos( R Q 3 ) This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. and calculate the three distinct real roots as:

    x 1 =( 2 Q cos θ 3 ) a 3 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.6a)

    x 2 =( 2 Q cos θ+2π 3 ) a 3 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.6b)

    x 3 =( 2 Q cos θ2π 3 ) a 3 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.6c)

    Note that x1, x2, x3 are not given in any special order, and that θ This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. has to be calculated in radians.

  2. If M > 0 (R2 > Q3), the polynomial has only one real root. Compute:

    S= R+ M 3 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.7a)

    T= R M 3 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.7b)

    and calculate the real root as follows:

    x 1 =S+T a 3 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.7c)

Two complex roots (complex conjugates) may be found as well. However, they are of no interest for our purposes, and thus no formulas are provided. Such formulas can be found in the following suggested readings:

  • W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran, 2nd Edition, Cambridge Univ. Press, (1992), p. 179.
  • Spiegel, M., Liu, J., Mathematical handbook of formulas and tables, 2nd Edition, Schaum’s Outline Series, McGraw Hill, p.10.

Sometimes, the equations for S and T listed above cause problems while programming. This usually happens whenever the computer/calculator performs the cubic root of a negative quantity. If you want to avoid such a situation, you may compute S’ and T’ instead:

S'=sign( R ) abs( R )+ M 3 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. T'=Q/S' This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (making T’=0 when S’=0)

where:

abs(R) = Absolute value of R

and

sign(R) = (+1) or (– 1) if R is positive or negative respectively.

It may be defined as:

sign(R) = R/Abs(R)

and then the real root is:

x 1 =S'+T' a 3 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers.

Keep in mind the following useful relationships among the roots of any cubic expression:

x 1 + x 2 + x 3 =a This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.8a)

x 1 x 2 + x 2 x 3 + x 3 x 1 =+b This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.8b)

x 1 x 2 x 3 =c This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.8c)

Self-Check

Test your understanding of cubic root calculations by analytical means by solving the following examples.

With general cubic expressions,

x 3 6 x 2 +11x6=0 x 1 =1 x 2 =2 x 3 =3 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers.
x 3 +7 x 2 +49x+343=0 x 1 =7 x 2 =0+7i x 3 =07i This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers.
x 3 +2 x 2 +3x+4=0 x 1 =1.65063 x 2 =0.174685+1.54687i x 3 =0.1746851.54687i This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers.

with EOS in terms of volume (v), for a pure component,

v 3 7.8693 v 2 +13.3771v6.5354=0 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers.
=> One real root (one phase),
v 1 =5.7357 v 3 15.6368 v 2 +30.315v14.8104=0 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers.
=> Two possible phases (three real roots),
v1 = 0.807582 (liquid phase)
v2 = 1.36174  (rejected)
v3 = 13.4675(gas phase)

with EOS in terms of the compressibility factor (z), for a pure component,

z 3 1.0595 z 2 +0.2215z0.01317=0 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers.
=> One phase,
z=0.8045 z 3 z 2 +0.089z0.0013=0 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. => Two possible phases,
z liquid =0.0183012 z x =0.0786609 (useless) z vapor =0.903038 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers.

Numerical Scheme

The Newton-Raphson method provides a useful scheme for solving for a non-explicit variable from any form of equation (not only cubic ones). Newton-Raphson is an iterative procedure with a fast convergence, although it is not always capable of providing an answer — because a first guess close enough to the actual answer must be provided.

In solving for “x” in any equation of the type f(x)=0, the method provides a new estimate (new guess) closer to the actual answer, based on the previous estimate (old guess). This is written as follows:

x new = x old f( x old ) f'( x old ) This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.9)

Considering a cubic equation f( x )= x 3 +a x 2 +bx+c=0 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. , the previous equation takes the form:

x new = x old x old 3 +a x old 2 +b x old +c 3 x old 2 +2a x old +b This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.10)

The iterations continue until no significant improvement for “xnew” is achieved, i.e., | xnew – xold | < tolerance. An educated guess must be provided as the starting value for the iterations. If you are solving a cubic equation in Z (compressibility factor), it is usually recommend to take Z = bP/RT as the starting guess for the compressibility of the liquid phase and Z = 1 for the vapor root.

Semi-analytical Scheme

If you used the previous numerical approach to calculate one of the roots of the cubic expression, the semi-analytical scheme can give you the other two real roots (when they exist). By using the relationships given before, with the value ‘x1’ as the root already known, the other two roots are calculated by solving the system of equations:

x 2 + x 3 =a x 1 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.11a)

x 2 x 3 =c/ x 1 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.11b)

which leads to a quadratic expression.

This procedure can be reduced to the following steps:

  1. Let x 3 +a x 2 +bx+c=0 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. be the original cubic polynomial and “E” the root which is already known (x1 = E). Then, we may factorize such a cubic expression as:

    ( xE )( x 2 +Fx+G )=0 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.12a)

    where:

    F=a+E This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.12b)

    G=c/E This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.12c)

  2. Solve for x2, x3 by using the quadratic expression formulae,

    x 1 =E This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.13a)

    x 2 = F+ F 2 4G 2 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.13b)

    x 3 = F F 2 4G 2 This equation is not rendering properly due to an incompatible browser. See Technical Requirements in the Orientation for a list of compatible browsers. (11.13c)