Introduction

Stability testing analysis has been carried out for pure systems as well as for mixtures by global minimization of the Gibbs energy. The tangent plane distance (TPD) function minimization method [8, 12] was among the earliest approaches and since then many researchers have improved on various solution methods [16]. According to these works, the stability status—whether stable or unstable—of a system is tested at specified temperature and pressure conditions The procedure involves the formulation of ‘mesh’ equations whose solution is a very complex task and often gives false solutions. Research efforts to date have gone into simplifying these solutions and eliminating the false ones. Stability testing analysis, as advocated by Michelsen [12, 13], has become accepted as a necessary preliminary step for phase equilibrium calculations. This is because the number of phases involved in a phase equilibrium calculation is not known beforehand. This number is, however, crucial for efficient calculations and without it, huge computation costs are incurred due to wrong initial guesses and false solutions. Once it is established through stability testing analysis that a phase is stable, the initial estimate required for efficient equilibrium calculations is automatically provided by the stability test results. Thus, stability analysis has become an area of concern not only in equilibrium calculations, but also in phase behavior modeling in general. This generalization will be quite complicated if it is to be applied to the management of phase behavior and phase split calculations in liquefied natural gas (LNG) systems, where volume requirements for storage, handling, and transportation must be minimized through accurate prediction of phase stability and phase split conditions.

Almost all of the stability testing methods are based on the minimization of the Gibbs free energy of the system at a given temperature and pressure condition. This choice of set of thermodynamic properties is made for convenience and ease of development of equations whose solution has proved to be quite complicated [10, 14].

In this work, the minimization of the Helmholtz free energy (A) at constant temperature (T) and volume (V) is chosen as an alternative set of thermodynamic properties. Though mention has been made of this set as a possible and preferable alternative in the open literature [9], very few researchers [15] have employed this approach to determine the stability limit of pure and binary systems. Danesh [9] and Michelsen [14] agree that the Helmholtz energy formulation would yield simpler equations that will enhance easier solutions. We note here that this is true for stability analysis testing for pure compounds and binary systems but would likely, as indicated by Modell and Reid [15], become more complex with mixtures. This complexity, however, may be reduced if a major modification of the attractive force parameter of the EOS model is carried out [7, 9].

Stability testing and stability limit determination

Stability testing solution methods, as the name implies, are used as point tests to determine whether a system at a specified condition is stable. Such methods are limited in the sense that they cannot be used to determine the point at which an unstable system would become stable and vice versa, neither can they serve to determine the degree (extent) of stability or instability of a system. The result may be different if a change occurs in the temperature and/or pressure conditions of the said system. This limits any stability test to a point test and must be repeated if any change in T or P takes place. The degree of this change is another issue yet to be addressed but suffice to state qualitatively here that such changes are not likely to be negligible. The Gibbs free energy minimization carried out at constant T and P is therefore a very convenient choice for stability testing but may not be the best choice for wide-range stability limit determination.

A stability limit determination model, on the other hand, ‘searches’ over a wide pressure range to determine the limit of intrinsic thermodynamic phase instability and tells at what point (of pressure) the liquid or vapor system in question will become stable. This approach is based on the Helmholtz energy stability criterion which is given in Eq. (1):

$$A_{\text{VV}} = \left( {\frac{\partial P}{\partial V}} \right)_{T,N} > 0.$$
(1)

The advantages of stability limit determination over stability testing are many and include the following:

  1. a.

    a single solution determines the phase state and its stability status over a range of pressure values;

  2. b.

    we do not merely discover whether a system is stable but we are able to know when a system will become stable (if it will) over the pressure range considered;

  3. c.

    computation efforts are conserved because we are told if the system will never become stable at the given conditions;

  4. d.

    in stability testing where T and P are held constant, any change in either T or P will require a fresh test but in stability limit determination, this set of two variables is reduced to one, i.e., temperature alone may vary;

  5. e.

    a single result gives the required value for both liquid and vapor phases.

Stability limit profile

It is very needful to have, in open literature, safe regions of meaningful equilibrium calculations for various compounds (both for pure systems and for complex mixtures) that are of practical importance, especially with respect to reservoir fluids. This will save a lot of computation efforts presently being lost due to false solutions, trivial solutions, local minima, and poor initial guesses. A novel work on stability limit trends for pure systems was reported for cyclohexane [6]. It is here extended to LNG systems. Firstly, stability limit profiles are determined using the new model for constituents of an LNG and then with improved robustness, the model is employed to determine the stability limit profile of LNG mixture over a wide range of pressure and temperature.

Methodology

We begin with the determination of stability limits of LNG constituents using a new model based on the stability criterion of the Helmholtz energy using the Peng–Robinson Equation of State (PR-EOS) [3]. This EOS model is then extended to cover a wide range of pressures and then generate a pressure–temperature stability profile for each of the pure compounds contained in LNG. Next, we apply the model to an LNG mixture using mixing rules for the parameters to generate its stability limit profile. Finally, we generalize the model for applicability to any cubic EOS and demonstrate its robustness by an algorithm that calculates the spinodal and binodal curves for all constituents of LNG. This is done for both liquid and vapor phases using two cubic EOS namely; the Peng–Robinson (PR-EOS) and the Soave–Redlich–Kwong (SRK-EOS). A comparison is done with the recent work of Aursand et al. [1] to both validate our method and demonstrate its robustness.

In most LNG systems, nitrogen is negligible both in percentage concentration and in its impact on the stability limit profile and phase behavior of the mixture. Where nitrogen content is significant, above 1.0 mol%, it cannot be neglected in phase split calculations for LNG systems and must therefore be included in the combined model for the LNG mixture [4, 11].

Mathematical formulation

Model for pure systems

The mathematical formulation begins with the criterion for intrinsic stability in the Helmholtz energy representation as stated in Eq. (1). The PR-EOS as given in Eq. (2) is employed for both the vapor and the liquid phases of the systems considered in this work:

$$P = \frac{RT}{v - b} - \frac{{a(\omega ,T_{\text{r}} )}}{v(v + b) + b(v - b)}.$$
(2)

At the limit of intrinsic stability, Avv = 0, i.e.,

$$\left( {\frac{\partial P}{\partial V}} \right)_{T,N} = 0.$$
(3)

For a pure material, we substitute Eq. (2) into Eq. (3) and differentiate, to obtain

$$PV^{4} + 4PBV^{3} + (2PB^{2} - A)V^{2} + (2AB - 4PB^{3} )V + PB^{4} + AB^{2} = 0,$$
(4)

where

$$\begin{aligned} & A = a \times a = a(\omega ,T_{\text{r}} ) \\ & B = b. \\ \end{aligned}$$

Equation (4) can be re-written as follows:

$$PV^{4} + EV^{3} + FV^{2} + GV + H = 0,$$
(5)

where

$$\begin{aligned} & E = 4PB \\ & F = 2PB^{2} \, {-} \, A \\ & G = 2AB \, {-} \, 4PB^{3} \\ & H = PB^{4} + AB^{2} . \\ \end{aligned}$$

The solution of Eq. (5) yields four roots. The first root is negative; the second is smaller than b. The third and fourth roots represent the volume for the liquid phase and vapor phase, respectively. The root with a negative value has no physical implication and so is discarded. The root whose value is less than ‘b’ is also discarded since ‘b’ represents the smallest molar volume of the system at infinite pressures [12].

The results from Eq. (5) give the limit-of-stability curve (on a P–V plane) for a pure system. Its point of stability at a selected temperature is the P value at which the T-isotherm intersects the limit-of-stability curve [15].

Hence, a T-isotherm can be obtained from the EOS (Eq. 2). When the Peng–Robinson EOS is rearranged, it yields Eq. (6) given as

$$PV^{3} + JV^{2} + KV + L = 0,$$
(6)

where

$$\begin{aligned} & RT = C \\ & J = PB -C \\ & K = A - 3PB^{2} - 2CB \\ & L = CB^{2} + PB^{3} - AB. \\ \end{aligned}$$

Equation (6) yields three roots when solved; the smallest (which is greater than ‘b’) represents the liquid molar volume, the largest is the vapor molar volume while the intermediate is irrelevant since it is taken to represent a value in the metastable region [15].

Equations (5) and (6) can be solved simultaneously to eliminate one of the set of variables (P, V, T) to determine the limit of stability. Alternatively, this can be accomplished by the intersection of the stability curve and the isotherm when displayed on the same graph for the determination of the limit of stability for the pure material at a given temperature. If such a limit exists, then the intersection of these two curves indicates the limit of stability (in pressure) as illustrated by Fig. 1. Below the pressure at the intersection of the two curves, the liquid phase is stable, while above the pressure of intersection, the liquid phase is unstable; if, however, the two curves do not meet within the region of interest, then there is no limit of stability for that substance at the temperature of interest and the system is then said to be in a metastable or unstable state (an example is shown in Fig. 2). In this work, Eqs. (5) and (6) were used to determine the stability limits for LNG constituents namely; methane, ethane, propane, and butane as given in Table 1.

Fig. 1
figure 1

Limit of stability for benzene (liquid) at 522 K (stable at 3.6 bar)

Full size image
Fig. 2
figure 2

Limit of stability for cyclohexane (liquid) at 420 K (undetermined)

Full size image
Table 1 Properties of LNG constituents
Full size table

Extension to mixtures

With the Peng–Robinson EOS given in Eq. (7) for mixtures and specifically for 1 mol of the mixture,

$$P = \frac{RT}{{V_{\text{m}} - b}} - \frac{a}{{V_{\text{m}}^{2} + 2bV_{\text{m}} - b^{2} }},$$
(7)

where \(V_{\text{m}}\) is the volume of the mixture, T is the temperature of the system, and P is the total pressure of the system. The application to mixtures was done using the following mixing rules:

$$a = \mathop \sum \limits_{ij} y_{i} y_{j} a_{ij} ,$$
$$a_{ij} = (a_{i} a_{j} )^{1/2} ,$$
$$a_{i} = a\left( {T_{{{\text{c}}_{i} }} } \right)a_{i} (T),$$
$$a_{i} (T) = [1 + \kappa_{i} (1 - T_{{{\text{r}}_{i} }}^{1/2} )]^{2} ,$$
$$\kappa_{i} = 0.37464 + 1.54226\omega_{i} - 0.26992\omega_{i}^{2} ,$$
$$a(T_{{{\text{c}}_{i} }} ) = 0.45724\frac{{(RT_{{{\text{c}}_{i} }} )^{2} }}{{P_{{{\text{c}}_{i} }} }},$$
$$b = \mathop \sum \limits_{i} y_{i} b_{i},$$

where \(b_{i}\) is the parameter ‘b’ for pure species i and \(y_{i}\) is the mole fraction of pure species i in the system.

These mixing rules were employed along with Eqs. (5) and (6) to determine stability limits at various T values for the LNG mixture.

Model for spinodal and binodal curves for pure systems

For generalized cubic EOS of the form:

$$P = \frac{RT}{v - b} - \frac{a}{v(v + b) + b(v - b)},$$
(8)

where P is pressure, v is molar volume, and T is temperature.

The parameters for PR-EOS are:

$$a = 0.45724 \propto R^{2} T_{\text{c}}^{2} /P_{\text{c}} ,$$
$$b = 0.0778RT_{\text{c}} /P_{\text{c}} ,$$
$$\propto \, = [1 + m(1 - T_{\text{r}}^{0.5} )]^{2} ,$$
$$m = 0.37464 + 1.5422\omega - 0.26992\omega^{2} ,$$
$$u = 2b,$$
$$w = b.$$

For SRK-EOS, the parameters are given by:

$$a = 0.42747 \propto R^{2} T_{\text{c}}^{2} /P_{\text{c}} ,$$
$$b = 0.0778RT_{\text{c}} /P_{\text{c}} ,$$
$$\propto \, = \, [1 + m(1 - T_{\text{r}}^{0.5} )]^{2} ,$$
$$m = 0.48508 + 1.55171\omega - 0.15613\omega^{2} .$$

Applying the Helmholtz energy criterion for stability given as:

$$\partial^{2} A = 0.$$
(9)

and from basic classical thermodynamics, the derivative of the Helmholtz free energy is:

$${\text{d}}A = - \,S{\text{d}}T + P{\text{d}}V,$$
$$({\text{d}}A)_{T} = P{\text{d}}V,$$
$$\left( {\frac{\partial A}{\partial V}} \right)_{T} = P,$$
$$\left( {\frac{{\partial^{2} A}}{{\partial V^{2} }}} \right)_{T} = \left( {\frac{\partial P}{\partial V}} \right)_{T} ,$$
$$\partial^{2} A = \left( {\frac{\partial P}{\partial V}} \right)_{T} \partial^{2} V.$$
(10)

Equating Eqs. (9) and (10), we obtain

$$\begin{aligned} \left( {\frac{\partial P}{\partial V}} \right)_{T} \partial^{2} V = 0 \hfill \\ \left( {\frac{\partial P}{\partial V}} \right)_{T} = 0. \hfill \\ \end{aligned}$$
(11)

Differentiating Eq. (8), we obtain

$$\left( {\frac{\partial P}{\partial V}} \right)_{T} = \frac{RT}{{(v - b)^{2} }} + \frac{ \propto (2v + u)}{{(v^{2} + uv - w^{2} )^{2} }}.$$
(12)

Equating Eqs. (11) and (12) we have

$$\begin{aligned} & 0 = \frac{RT}{{(v - b)^{2} }} + \frac{ \propto (2v + u)}{{(v^{2} + uv - w^{2} )^{2} }} \\ & 0 = \frac{{RT(v^{2} + uv - w^{2} )^{2} + \propto (2v + u)(v - b)^{2} }}{{(v - b)^{2} (v^{2} + uv - w^{2} )^{2} }}. \\ \end{aligned}$$

So that

$$RT(v^{2} + uv - w^{2} )^{2} + a(2v + u)(v - b)^{2} = 0.$$
(13)

Rearranging Eq. (7), we obtain

$$\begin{aligned} & v^{4} - \left( {\frac{2a}{RT} - 2v} \right)v^{3} + \left( {u^{2} + \frac{au}{RT} - 2w^{2} + \frac{4ab}{RT}} \right)v^{2} \\ & \quad - \left( {\frac{{2ab^{2} }}{RT} - \frac{2aub}{RT} + 2uw^{2} } \right)v + \left( {w^{4} - \frac{{aub^{2} }}{RT}} \right) = 0. \\ \end{aligned}$$
(14)

Equation (14) was solved for a pure system at a selected temperature to obtain four roots (volume). The values that are negative or less than ‘b’ were discarded, of the remaining two, the smaller value is for the liquid phase while the larger value is for the vapor phase. Substituting these two values into Eq. (8) yields the corresponding value of pressure (P).

Algorithm

Liquefied natural gas (LNG) is a natural gas that has been cooled to its two-phase region where it condenses to liquid; this occurs at a temperature of approximately 112 K (− 256 °F, − 161 °C) at atmospheric pressure. Methane is the major constituent of the natural gas mixture with small amounts of ethane, propane, and other hydrocarbons. The inorganic compound impurities in trace amounts are nitrogen, carbon dioxide, and hydrogen. The inorganic compounds are not desirable because they are not combustible and cause corrosion and other problems in gas production and processing systems. The algorithm used to determine the spinodal (curve) for methane, ethane, propane, butane, and nitrogen is shown using the flowchart in Fig. 3.

Fig. 3
figure 3

Algorithm for determination of the spinodal points of a pure component

Full size image

Results and discussion

We present the results here in three parts:

Firstly, the stability limit P–T profiles for constituents of LNG is given by the solution of Eqs. (5) and (6) as shown in Fig. 4 for the liquid phase, and in Fig. 5 for the vapor phase. Theoretically, the stability limit curve intersects the T-isotherm at the point of stability of the system. An illustration of this is given, as earlier mentioned, in Figs. 1 and 2. By such intersections (where they exist), the curves in Figs. 4 and 5 were all generated. This was demonstrated mathematically for pure and binary systems by Modell and Reid [15] and validated using the experimental works of Skripov and Ermakov [17]. The model equations used here were developed and first applied by Babalola and Susu [5] to pure n-hexane giving a near-perfect match to the results obtained mathematically by Modell and Reid [15], and that obtained experimentally by Skripov and Ermakov [17]. Here, this EOS model is for the first time, applied to LNG constituents and extended to LNG mixture. Beyond determining, by point tests, whether a phase of a system is stable at the prevailing P and T, this work presents the stability limit trends (profiles) of both the liquid and vapor phases of LNG systems over the temperature and pressure ranges encountered in the industry. Such a tool has been lacking in the thermodynamic modeling framework (TMF) for LNG systems according to Firoozabadi and Pan [10]. The TMF of complex LNG systems has various components which are being addressed by researchers to improve thermodynamic properties and phase split predictions [4]. This, therefore, is a significant contribution to facilitate improved design and operation of separation processes in the oil and gas industry especially at low-temperature zones. Hence, with the data information in Figs. 4 and 5, it is easy to identify the T and P conditions under which any of the four constituents of LNG would be stable. This knowledge will be helpful in making crucial decisions not only in design and operation of LNG facilities, but also in its handling and transportation. Obviously, liquid methane requires the highest pressure to be stable and the stability pressure range is narrow (42–52 Bars) with a fairly wider temperature range of 113–173 K. For the vapor phase, however, the pressure range of stability is much wider (10–50 Bars) for about the same temperature range.

Fig. 4
figure 4

Stability limit profiles for LNG constituents (liquid phase)

Full size image
Fig. 5
figure 5

Stability limit profiles for LNG constituents (vapor phase)

Full size image

Next, in Figs. 6 and 7 are shown the stability limit profiles for the LNG mixture whose composition is given in Table 1 for liquid phase and vapor phase, respectively. These profiles were generated for PR-EOS and SRK-EOS. The stability pressure range for the liquid phase is much shorter than that for the vapor phase. The range predicted by PR-EOS is somewhat wider than that predicted by the SRK-EOS for the liquid phase, while the reverse is the case for the vapor phase mixture.

Fig. 6
figure 6

Stability limit profiles for LNG mixture (liquid) using PR-EOS and SRK-EOS

Full size image
Fig. 7
figure 7

Stability limit profiles for LNG mixture (vapor) using PR-EOS and SRK-EOS

Full size image

Lastly, the spinodal curve is here for the first time, as far as we are aware, generated for LNG constituents in Figs. 8, 9, 10, 11, 12, 13, 14, 15, 16 and 17. Most results in literature determine the spinodal point; this is the limit of stability for a system at a specified T and P, for which the corresponding volume or density can be calculated. In contrast, our model here presented gives the spinodal curve and so covers a wide range of T while yielding the corresponding range of pressures. The binodal curve is the equilibrium curve for the liquid and vapor phases; while the spinodal represents the limit of stability with suppression of nucleation sites. At the limit of stability (spinodal), therefore, the system must spontaneously split into two phases. The region between the spinodal and the binodal curves represent the metastable region which tappers towards the critical point beyond which both phases are no more distinguishable. Our model accurately predicts this critical point P and T for each of the constituents of LNG and is hence validated to some extent.

Fig. 8
figure 8

Liquid/vapor spinodal and binodal curves for methane using PR-EOS

Full size image
Fig. 9
figure 9

Liquid/vapor spinodal and binodal curves for methane using SRK-EOS

Full size image
Fig. 10
figure 10

Liquid/vapor spinodal and binodal curves for ethane using PR-EOS

Full size image
Fig. 11
figure 11

Liquid/vapor spinodal and binodal curves for ethane using SRK-EOS

Full size image
Fig. 12
figure 12

Liquid/vapor spinodal and binodal curves for propane using PR-EOS

Full size image
Fig. 13
figure 13

Liquid/vapor spinodal and binodal curves for propane using SRK-EOS

Full size image
Fig. 14
figure 14

Liquid/vapor spinodal and binodal curves for butane using PR-EOS

Full size image
Fig. 15
figure 15

Liquid/vapor spinodal and binodal curves for butane using SRK-EOS

Full size image
Fig. 16
figure 16

Liquid/vapor spinodal and binodal curves for nitrogen using PR-EOS

Full size image
Fig. 17
figure 17

Liquid/vapor spinodal and binodal curves for nitrogen using SRK-EOS

Full size image

The spinodals were determined for liquid and vapor phases of methane using both PR-EOS (Fig. 8) and SRK-EOS (Fig. 9). Interestingly, this robust model also includes the binodal curve over the range of T and P so determined. An additional validation from the spinodal for methane is in comparison with the recent work of Aursand et al. [1]. In their work, they used the PR-EOS to determine the spinodal point by drawing an isotherm for methane liquid at 175 K and obtained a P value of 1.13 MPa at a density of 236.1 kg/m3. Our PR-EOS prediction of spinodal at the same T and density is 1.1626 MPa, as can be read off from Fig. 8, is almost in perfect agreement with theirs. Also, their spinodal point for vapor phase methane at T = 175 K and density of 87.5 kg/m3 was 3.243 MPa while the spinodal pressure predicted here at the same T and density is 3.2578 MPa. Furthermore, the robustness of our model is demonstrated in the fact that it predicts the spinodal T and P over a wide range. A robust model by Aursand et al. [1] was also able to generate the spinodal of methane over a range of T and P to which our results closely compares as follows: at 185 K, their spinodal for liquid methane was 3.74 MPa while ours is 3.58 MPa, also, for methane vapor at 180 K they had 3.68 MPa while we got 3.61 MPa. The spinodal and binodal curves are shown for methane using SRK-EOS in Fig. 9 for which the spinodal P = 1.2932 MPa for the liquid phase and 3.2627 MPa for the vapor phase at 175 K. Interestingly, both EOS models predict most closely at the critical points.

Figures 10 and 11 show the spinodal and binodal curves for ethane while those of propane are shown in Figs. 12 and 13. Figures 14 and 15 are for butane while those of nitrogen are shown in Figs. 16 and 17. For nitrogen, Aursand et al. [1] predicted the spinodal for liquid nitrogen at 120 K to be 1.96 MPa while our model gave a value of 2.04 MPa, and for nitrogen vapor at 120 K, they predicted 2.64 MPa while we obtained 2.73 MPa. Most uniquely for  the LNG mixture are the spinodals for the first time here given in Figs. 18 and 19. Here it is shown that for the LNG mixture, PR-EOS is more reliable than SRK-EOS for the prediction of spinodals. This is evident from the negative pressure values (suction) predicted by the SRK model.

Fig. 18
figure 18

Spinodal of liquid LNG mixture

Full size image
Fig. 19
figure 19

Spinodal of vapor LNG mixture

Full size image

The region between the spinodals and the binodal curve leaves much to be addressed by investigators; according to Aursand et al. [1], much of the metastable region of the phase diagram is experimentally unavailable at present. Few researchers, however, like Avedisian [2], have published some experimental works in this area. Therefore, reliable models such as this are the needed tools for skilled decision-making on regions of operation where phase stability is guaranteed. For the LNG system treated here, only two phases (liquid and vapor) were determined since the concentration of nitrogen and other compounds that could form a third phase at extreme conditions were negligible.

Conclusion

A model for stability of pure systems, by intersection point of stability limit curve and T-isotherm, is employed to determine the stability limit profiles of LNG constituents through generation of multiple intersection points.

Subsequent to this, an extension was made to multicomponent mixtures and was successfully applied to determine the stability limit profiles of an LNG system.

A new and robust EOS model for generating spinodal and binodal curves for pure systems was developed and was successfully applied to LNG constituents including nitrogen (for complex LNG systems).

The spinodals for the LNG mixture, however, did not show a very close match between PR-EOS and SRK-EOS. It is quite evident that the prediction by PR-EOS is more accurate.

Two EOS (PR-EOS and SRK-EOS) were used for both liquid and vapor phases and gave fairly close results.