WO2014027196A2 - Physical property modelling - Google Patents

Abstract

A method of determining K-values (partition coefficients) for a mixture comprises the following steps: (i) calculating K-values for the mixture explicitly at a number of pressures and at a temperature T; (ii) determining a projected critical pressure p_pc as the pressure above which a two-phase vapour-liquid phase equilibrium may no longer be estimated; and (iii) regressing the K-values as a function of p_r=p/p_pc.

Description

Physical property modelling The present invention relates to the field of physical property modelling, particularly, but not only, modelling involving simulations of petroleum fluids.

In chemical engineering, modern computer-based design tools rely on numerical models to estimate the thermodynamic and physical properties of fluids as the basis of their methodology. In the oil industry, many computationally intensive applications such as reservoir simulation, well performance studies, transient pipeline simulation, dynamic facilities design and life-of-field studies are restricted in practice by the speed at which a computer can solve the physical property models.

The conventional approach involves performing rigorous phase equilibrium calculations. In the majority of cases, simulation of the properties of a fluid requires as a first step the calculation of its phase equilibrium. The simplest example of a phase equilibrium problem is one where the pressure, temperature and overall composition (termed the feed) is defined, and the task is to predict how many phases the fluid will separate into, and what the compositions and relative amounts of the phases will be. Thermodynamically, such an equilibrium corresponds to the point where the Gibbs energy is at a minimum, and the necessary thermodynamic properties are calculated from a rigorous model such as an equation of state.

Considering the most common phase equilibrium, which involves combinations of a vapour (or gas-like) phase and a liquid phase, the first step is to establish whether the feed is stable and remains as a single phase, or whether it will split into two phases. The established technique which is used to detect possible phase splitting is the phase stability algorithm proposed by Michelsen (see M. L Michelsen: Fluid Phase Equilibria, 9, 1 , 1982). If it transpires from this algorithm that the feed will split, it is then necessary to calculate the actual phase equilibrium, i.e. what the composition and relative amount of each phase will be.

The most common method of calculating phase equilibria is the successive substitution method (see M. L Michelsen: Fluid Phase Equilibria, 9, 21 , 1982). The initial step in this method is to estimate the likely compositions y, x of the vapour and liquid phases, where y is a vector in composition hyperspace of the vapour mole fractions y,, and x is a vector of the liquid mole fractions x,. Subscript i denotes the i^th component which corresponds to the i^th dimension of the hyperspace.

A benefit of the Michelsen phase stability test is that it provides estimates for y and x in the case where it predicts that the feed will split into two phases.

The partition coefficient, which is also referred to as the K-value, is defined as the ratio of vapour to liquid mole fractions K,= y x, for each component i.

It follows from standard thermodynamic theory that the K-values can be evaluated from an equation of state, so

— = K, { p , T , x , y) where p is the pressure and T is the temperature.

Once the K-values are calculated, the next step is to obtain the vapour phase fraction F_v from the Rachford-Rice equation

where vector z represents the feed mole fractions z,.

Equation 2 is not explicit in F_v so it has to be solved iteratively using a suitable numerical technique such as Newton's method.

To complete the procedure, new values of the vapour and liquid mole fractions are obtained from the relationships y' F_v K_t+ 1- F_v (³) and x,=

F_VK,+ l - F, (4)

^{1 1} V

The complete rigorous algorithm to calculate the composition and relative amounts of each phase in the vapour-liquid phase equilibrium is therefore:

1. estimate the compositions of the vapour and liquid phases (e.g. from the Michelsen phase stability test);

2. determine the K-values from the thermodynamic model/equation of state;

3. solve Equation 2 iteratively for the vapour phase fraction F_v;

4. if F_v has not changed significantly exit as the problem is solved; 5. otherwise, update the values of the vapour and liquid mole fractions y, and x,with Equations 3 and 4; 6. go back to step 2.

The algorithm is computationally expensive because solving the equation of state to obtain the K-values (step 2) takes significant time even for relatively simple engineering models such as the Peng-Robinson (PR) or Soave-Redlich-Kwong (SRK) equations of state. Some more recent models such as the GERG or SAFT equations take much longer.

In addition, because the K-values are only implicitly defined by the rigorous model, the whole procedure (steps 1-6) is iterative so the K-values have to be evaluated (step 2) many times during one phase equilibrium calculation. In some cases, the procedure will not converge to a solution and alternative, more computationally intensive methods have to be used in order to find the equilibrium. In addition, the phase equilibrium calculation must often be preceded by the phase stability analysis procedure that takes comparable, if not more, computing time.

Fig. 1 shows a plot on a pressure-temperature plane of typical vapour-liquid phase behaviour exhibited by petroleum fluids, which are mixtures of predominantly non-polar hydrocarbon molecules. There is a region of two-phase behaviour bounded by the bubble-point and dew point lines. At higher pressures or temperatures the fluid is a single phase, although the identity of this phase gradually changes from being liquid-like at low temperatures to vapour-like at high temperatures. The bubble-point and dew-point lines meet at the critical point, the point where the vapour and liquid phases become just indistinguishable. It follows that at the critical point, the K-value of every component is unity.

Attempts have been made to speed up the conventional approach described above. In particular, the calculations can be made to proceed much faster if the required properties are expressed explicitly. Such a method starts from the observation that for most petroleum fluids the K-values are strongly dependent on temperature and pressure but almost independent of composition. The K-values may therefore be represented for a defined feed as explicit functions of pressure and temperature only. Such a method is referred to as an explicit K-value method.

The problem of finding the composition and relative amounts of each phase in the phase equilibrium then reduces to solving Equation 2 iteratively only once for F_v, and then using Equations 3 and 4 to obtain the compositions x and y for the vapour and liquid phases. If the feed is in the single phase region, the value of F_v obtained from Equation 2 lies outside the physically meaningful range of 0≤F_V≤1. When this happens, the solution of Equation 2 is abandoned as the fluid is known to be a single phase of composition z.

In order to set up such an explicit K-value method, it is necessary to know what the K-values should be. In principle the K-values could be regressed to match experimental data, but this is rarely possible as the amount of available data at the relevant conditions is normally very limited. A practical approach is to use a rigorous method (e.g. based on an accepted thermodynamic model) to calculate the phase equilibrium and hence the K-values at a small number of points, and then to regress the K-values explicitly for use in Equations 2, 3 and 4. Put another way, the method consists of interpolating between a small number of rigorous phase equilibrium calculations.

Such a method only takes a fraction of the computing time of the rigorous method, which only defines K-values implicitly. As it is a compositional method, it can be applied for any feed that is not too different from the specified feed z, the assumption being that the K-values remain approximately the same functions of pressure and temperature.

As the K-values at the critical point are all unity, it follows that the explicit K- value method places the critical point at a fixed temperature and pressure no matter what the feed composition.

If we consider an isotherm at the critical temperature T=T_C, the pressure dependence of the K-values for a typical petroleum fluid will appear as in Fig. 2. Each line in Fig. 2 represents the natural logarithm of the K-value as a function of pressure, at temperature T_c, for a different component of the feed. The more volatile the component, the higher its K-values (and ln(K)). The important point to note is that the K-values all converge to unity (zero in logarithmic form) at the critical pressure p_c. Above the critical pressure, the fluid is one phase and so there can be no phase equilibrium and hence there are no K-values.

For an isotherm below the critical temperature, T<T_C, the vapour fraction decreases as the pressure rises until it reaches the bubble point where F_v=0. From Equation 2, it can be seen that if the feed composition z is changed from the specified value, F_v is no longer necessarily zero, i.e. the bubble point will change with the feed composition. In order to apply the explicit function method for varying feed compositions, it is therefore necessary to obtain K-values from the rigorous model at pressures above the bubble point for the specified feed so that they can be regressed to form the explicit K-value functions at such pressures. However, as the pressure increases, there comes a point where it becomes impossible to find K- values, and the fluid is always one phase. This limit is marked by the dashed line in Fig. 1.

For temperatures T above the critical temperature but below the

cricondentherm, T_c <T<T_cri, the roles of the vapour and liquid phases are reversed. (The cricondentherm is the highest temperature at which a vapour-liquid equilibrium occurs for the mixture in question.) As the pressure rises above the lower dew point, the liquid fraction initially increases but then falls again to reach zero at the retrograde or upper dew point where F_v=1. As before, the dew point will move if the feed composition is altered so there is a region above the dew point for the specified feed where a two-phase solution is possible for different feed

compositions. However, there is ultimately a limiting pressure above which only one phase can exist. At temperatures above the cricondentherm, T>T_cri, only one phase exists for the specified feed, although the cricondentherm will also move if the feed composition is altered.

It is necessary to calculate the K-values at pressures above the bubble-point (or dew-point) rigorously (i.e. from a thermodynamic model/equation of state) so that they can be correlated for use in the explicit K-value method. The K-values have to be obtained using the specified feed with the intention of then applying them to different although similar feed compositions.

One technique which has been tried involves systematically altering the feed composition z in such a way as to raise the bubble-point or dew-point pressure so that a phase equilibrium can be computed at higher pressures and K-values obtained. The disadvantage of this and other similar methods is that the K-values will show a discontinuity with pressure at the bubble point because the methodology for finding them has changed.

An alternative approach to calculating K-values at pressures above the bubble point or dew point involves using the rigorous K-value method set out above but relaxing the requirement at step 2 that F_v must lie within the physically meaningful range of 0≤F_V≤1. In this approach, it is simply required that F_v≤1 (i.e. allow F_v to be negative) for pressures above the bubble point or that F_v≥0 (i.e. allow F_v to be greater than 1) for pressures above the dew point. This approach is more consistent than the previous approach and can avoid a discontinuity in K-values at the bubble or dew point. However, both approaches have the intrinsic limitation that they cannot calculate K-values at pressures where there is no phase equilibrium, i.e. above the limiting pressure for the given mixture.

Thus, one objective of the present invention is to provide a method which can be used to calculate K-values at pressures above the bubble-point or dew-point pressure (i.e. at pressures where there is a single phase for a given mixture).

Further objectives are to provide improved methods of calculating K-values for mixtures, for use in simulations.

According to a first aspect of the invention, there is provided a method of determining K-values for a mixture, the method comprising the following steps:

(i) calculating K-values for the mixture explicitly at a number of pressures and at a temperature T;

(ii) determining a projected critical pressure p_pc as the pressure above which a two-phase vapour-liquid phase equilibrium may no longer be estimated; and

(iii) regressing the K-values as a function of p_pc.

The K-values, or partition coefficients, are defined as the ratio of vapour to liquid mole fractions K,= y for each component i.

Thus, by using this method K-values can be determined above the bubble- point or dew-point pressure for a given mixture. Such K-values may then be used to simulate a mixture with a slightly different composition, for example.

The method involves expressing the K-values as explicit functions of a small number of variables with the result that phase equilibrium calculations take only a fraction of the time required for a rigorous solution of the phase equilibrium with an equation of state. The procedure is well suited to larger-scale computer simulations where the time taken to solve the phase equilibria rigorously would render them impractical.

Besides K-values, other thermodynamic and transport properties needed in simulations may be correlated as simple functions using a similar and

complementary methodology. This is described in more detail below.

The method may be applied to any rigorous thermodynamic and transport property models. The resulting correlated expressions can be executed just as fast regardless of the time the rigorous methods take to execute. Therefore, the more computationally intensive the rigorous method, the greater the saving in computing time.

The method can provide a close approximation to the use of rigorous methods. In most engineering simulations the uncertainties in the data and the errors in the whole model will greatly exceed any error introduced by the explicit function method of the present invention. Provided regression points are close enough and sufficient coefficients are used to represent the K-values and physical properties as functions of pressure, etc., any reasonable degree of accuracy in the representation of the rigorous methods can be achieved. In this context, it should be remembered that the rigorous methods themselves are only approximations to physical reality.

Preferably, step (iii) comprises determining a reduced pressure p_r, where p_r = p/ppc and p is the actual pressure, and preferably then regressing the K-values as a function of p_r. p_r is an advantageous variable with which to regress the K-values.

The K-values may be regressed as a function of p_r using the equation

^{p_rK, ) = (c_l0 + c _Pr + c_ilP] + ... + c_imPr ^m y

for example, where i represents the components of the mixture, c₀-c_m are coefficients to be determined by the regression, q=1-p_r and φ is a constant for all components which is less than 1. Other, similar equations may also be used.

The method may further comprise regressing a variable other than the K- values as a function of p_r. The other variable could be the total volume, enthalpy, heat capacity, entropy, viscosity, thermal conductivity or interfacial tension of the mixture, for example. p_r is a convenient variable with which to regress these other variables, which may then be used in a simulation of the mixture.

In an example method, p_pc is determined by:

(a) rix Q, the elements of which are:

where A is the Helmholtz energy of the mixture, n is the mole number and subscripts i and j denote the components in the mixture;

(b) determining the composition n of the mixture and total volume V whilst requiring that the minimum eigenvalue of matrix Q, A_min=0 and the third derivative of A is zero in the direction of the associated eigenvector; and

(c) calculating p_pc from n and V.

At step (c) Ppc may be determined from an equation of state.

At step (b), n may be found by following the iterative procedure:

(b-1) estimate the mole number vector n of the projected critical fluid;

(b-2) form the matrix Q; (b-3) solve Q ^u to find A_min and its normalised eigenvector u; (b-4) calculate a hyperplane given by

_ _UT (Q(Aa)- Q(- Aa)) _ (Q(Aa)- Q(- Aa))

2Aa 2Aa

(b-5) extend a feed composition vector z in the direction u to intersect the hyperplane to determine the new mole number vector n;

(b-6) if n is not changed significantly, exit as the problem is solved; else (b-7) go back to step (b-2).

At step (b), V may found by following the iterative procedure:

(b-8) estimate V and n from previous projected phase equilibrium calculations;

(b-9) using the iterative procedure of claim 7, find the mole number vector n;

(b-10) calculate C, where C = u^TVA_mm = u^T fe(^Ag)- g(- )_u .

2Aa

(b-1 1) increment V;

(b-12) using the iterative procedure of claim 7, find the mole number vector n;

(b-13) re-calculate C;

(b-14) using the two most recently found values of V and C, extrapolate to a new value of V to satisfy C=0 by linear extrapolation;

(b-15) using the iterative procedure of claim 7, find the mole number vector n;

(b16) if V has not significantly changed from its most recently formed value, exit as calculation is complete; else

(b-17) re-calculate C; and

(b-18) go to step (b-14).

On convergence, the line from composition n in the direction of u is termed the critical tie line.

Other mathematical methods that satisfy the conditions A_min=0 and C=0, whilst imposing the condition that the feed composition z lies on the critical tie line may also be used to determined p_pc.

The K-values are preferably regressed on an isotherm at temperature T. The method preferably further comprises determining a vapour phase fraction F_v from each determined set of K-values. F_v may be determined iteratively, for example using the equation:

where i represents the components of the mixture, K, is the K-value for component i and is the amount of component i present in the mixture.

The method preferably further comprises determining vapour and liquid mole fractions from each K-value and F_v. This may be done using Equations 3 and 4.

Preferably, the K-values are treated as functions of pressure and a limited number of other variables but are independent of the composition of the mixture. In one example, the K-values are treated as a function of temperature and pressure only. In another example, the K-values are treated as a function of temperature, pressure and time only. Other variables could also be used depending on the particular simulation.

The K-values for the mixture are preferably determined using an equation of state. For example, the PR, SRK, GERG or SAFT equation of state may be used.

Preferably, the method is repeated at a different temperature. In this way, K-values covering a range of representative pressures and temperatures may be determined.

Preferably, the K-values are then splined between the different

temperatures. Splining may be performed using a cubic splining method, for example. The splining may be performed using a conventional cubic splining method such as the one described below.

However, preferably the splining is performed using a new splining method such as the one according to the tenth aspect of the invention described below. This is a less complicated method which avoids the need for calculating higher order derivatives.

The mixture may comprise a petroleum fluid, for example, and the K-values are then preferably used for a simulation of a petroleum fluid. This method is particularly useful when applied to petroleum fluids.

Preferably, the method further comprises using the K-values in a separator, distillation column, absorber, liquid-liquid separator, other unit operation, pipeline, reservoir or well simulation.

The K-values may be determined adaptively for pressures and temperatures as required in a simulation. This can help to avoid unnecessary computation. A benefit of this method is that it will make many design studies or control applications feasible which otherwise would not be carried out because of unacceptably long computing times.

The method may further comprise estimating K-values for a solid phase. In this case, the method could then further comprise modelling wax precipitation or asphaltenes, for example.

According to a second aspect of the invention, there is provided a computer program for performing the method of the first aspect.

According to a third aspect of the invention, there is provided a computer readable medium comprising a computer program according to the second aspect.

According to a fourth aspect of the invention, there is provided a method of determining K-values for a mixture explicitly at a number of pressures and at a temperature T, wherein the K-values are determined as a function of time, and the determined K-values are then used in simulation of a pipe or pipeline. This method has the benefit that the same explicit function can be used in the simulation all along the pipe or pipeline because it is not necessary to use the position on the pipeline explicitly in the expressions for the K-values and other physical properties. The method may be applied to other situations in which time can be used as a variable in place of one of more position variables.

In this case, preferably the method of the first aspect is used to determine the K-values.

The K-values may be determined as a function of time based on the flow into the pipe as a function of time, which may be a given piece of information. For example, in a real-time application the flow into the pipe as a function of time could be deduced from monitoring. In a simulation of past performance, the flow into the pipe as a function of time could come from historical records. If the pipe is connected to an oil well, the predicted flow as a function of time could come from a life-of-field reservoir model. In summary, the flow into the pipe as a function of time could be measured in real time or it could be estimated from some other plausible model.

According to a fifth aspect of the invention, there is provided a computer program for performing the method of the fourth aspect.

According to a sixth aspect of the invention, there is provided a computer readable medium comprising a computer program according to the fifth aspect According to a seventh aspect of the invention, there is provided a method of determining K-values for a mixture comprising vapour, liquid and solid phases, wherein first solid-liquid K-values are determined for the mixture and then the components of the solid phase are subsumed into the liquid phase to form an effective liquid phase before estimating resulting effective liquid-vapour K-values. By first determining K-values between the solid and liquid phases and then subsuming the solid phase components, determined from these K-values, into the liquid phase before estimating the resulting liquid-vapour K-values, the method is simplified and requires less computing time.

The method may further comprise calculating a phase equilibrium between the vapour and effective liquid phases, and preferably then splitting the effective liquid phase into oil and solid phases. This is a simple and effective way to simulate such a composition. The solid phase may be an asphaltene phase, for example.

Preferably, the effective liquid-vapour K-values are estimated using a rigorous calculation, for example using the method of the first aspect. However, since the solid phase does not display critical behaviour, any suitable method could be used to estimate the effective liquid-solid K-values.

The method may further comprise using the K-values for a separator, distillation column, absorber, liquid-liquid separator, other unit operation, pipeline, reservoir or well simulation.

According to an eighth aspect of the invention, there is provided a computer program for performing the method of the seventh aspect.

According to a ninth aspect of the invention, there is provided a computer readable medium comprising a computer program according to the eighth aspect.

As mentioned above, the method may involve a new splining method. This is considered to be inventive in its own right, regardless of whether or not it is used in combination with the method of the first aspect of the invention.

Thus, according to a tenth aspect of the invention, there is provided a method of splining K-value functions between splining points in two or more splining variables, wherein the splined K-value function is a function of switch functions and the K-value functions and their first order derivatives only.

Since only first order derivatives, and no higher order derivatives, are required, this method is less complicated than conventional cubic splining methods when applied to two or more splining variables. Preferably, the splining variable or variables is or comprises temperature. Preferably, the switch functions are polynomials, preferably cubic polynomials, of the splining variable.

Preferably, four switch functions are used and preferably these are given by:

S _{x)= (2-3x+ x³) 4 (I" X- X³)

T_{x)--^y '-

_T+(x)= (-i-*÷*²⁺*³) where x is a transformed splining variable which runs from -1 to +1 between the splining points.

There may only be one splining variable and in this case the splined function is preferably given by

f(x)= S_(x)f(-1)+ S₊(x)f(l)+ T_{x)f_x{-\)+ ₊(^χ)/_¾(1)

where f_x is df/dx.

If there are two splining variables, represented by the transformed splining variables x and y, each running from -1 to +1 between the splining points, the splined function is preferably given by

f(x,y)=S(x,y)+ S_x(x,y)+S_y(x,y)

where

S(x,y)= S_(x)S_(y)f(-l -1)+ S_(x)S₊(y)f (- 1,1)

+ S₊(x)S_(y)f(l,- 1)+ s₊(x)s₊(y)f(i,i)

S_x(x,y)= T_{x)S_{y)f_x{- \ -1)+ T_(x)S₊(y)f_x(- 1,1 )

+ T₊(x)S_(y)f_x(l - 1)+ T₊(x)S₊(y)f_x( l)

s_y(x,y)= s_{x)T_{y)f_y{- l - i)⁺ s_(x)T₊(y)fy(- U)

+ S₊(x)T_(y)f_y(l - l)+ S₊{x)T₊{y)f_y{\,\)

and f_x is df/dx and f_y is df/dy.

If there are three splining variables, represented by the transformed splining variables x, y and z, each running from -1 to +1 between the splining points, the splined function is preferably given by

fix, y,z) = S (x, y, z) + S_x (x, y, z) + S (x, y, z) + S_z (x, y, z) where

S(x,y,z)=S_(x)S_(y)S_(z)f(- 1 - 1,- 1)+ S_(x)S_(j/)S₊(z)/(- 1 - 1,1) + S_(x)S₊(y)S_(z)f(- U,- 1)+ S_(x)S₊(y)S₊(z)f(- 1,1,1)

+ S₊(x)S_(y)S_(z)f(l - 1 - 1)+ S₊(x)S_(j^;)S₊(z)/(l - 1,1)

+ S₊(x)S₊(^)S_(z)/(l,l - 1)+ S₊(x)S₊(y)S₊(z)f(l,l,l)

,y,z)= T_(x)S_(y)S_(z)f_x(- 1 - 1,- 1)+ r_(x)S_(j)S₊(z)/_x(- 1,- + Γ (x)S₊(y)S_(z)f_x(- LI.^" 1)⁺ Γ

1,1,1)

+ r₊(x)S_(^)S_(z)/_x(l - 1,- 1)+ r₊(x)S_(^)S₊(z)/_x(l - 1,1)

+ T x)S₊(y)S_(z)f_x(l,l,- l)+T₊(x)S₊(y)S₊(z)f_x(l,l,l)

,y,z)=S_{x)S_{y)T_ {z)f_z(-\ -\ -\)+S_{x)S_{y)T₊{z)f_z{-\ -\,\)

+S_(x)S₊(y)T_(z)f_z(-l,l -l)+S_(x)S₊(y)T₊(z)f_z(-lXl)

+S₊(x)Siy)T_(z)f_z(l -l -l)+S₊(x)S_(y)T₊(z)f_z(l -l,l)

+s₊(x)s₊(y)T-(z)f_z(^\ -i)+s₊(x)s₊(y)T₊(z)f₂(w) and f_x is df/dx, f_y is ef and f_z is δί/δζ.

The scheme can be extended in an analogous way to handle any number of splining variables.

The first order derivatives may conveniently be determined by a finite difference method.

According to an eleventh aspect of the invention, there is provided a computer program for performing the method of the tenth aspect.

According to a twelfth aspect of the invention, there is provided a computer readable medium comprising a computer program according to the eleventh aspect.

Preferably, the methods of the present invention are implemented using a computer system. Thus, the invention also extends to a computer system arranged for performing one or more of the methods of the present invention described above.

Key applications of the present invention involve modelling petroleum fluids which are predominantly composed of hydrocarbons. However, the method is not restricted to petroleum fluids and may be used to model other fluids. Preferred embodiments of the invention will now be described by way of example only and with reference to the accompanying drawings in which:

Fig. 1 is a plot on the pressure-temperature plane of typical vapour-liquid phase behaviour exhibited by petroleum fluids;

Fig. 2 is a plot of the pressure dependence of the K-values for a typical petroleum fluid for an isotherm at the critical temperature T=T_C;

Fig. 3 is a plot illustrating the procedure to extend a splining method to a function of two variables, f(x,y); and

Fig. 4 is a typical phase diagram plot for an asphaltenic crude oil.

As explained above, the explicit K-value method starts from the observation that for most petroleum fluids the K-values are strongly dependent on temperature and pressure, but almost independent of composition. The K-values are therefore represented for a defined feed as explicit functions of pressure and temperature only.

The problem of finding the phase equilibrium then reduces to solving

Equation 2 iteratively only once for F_v, and then using Equations 3 and 4 to obtain the compositions x and y for the vapour and liquid phases. If the feed is in the single phase region, the value of F_v obtained from Equation 2 lies outside the physically meaningful range of 0≤F_V≤1. When this happens, the solution of Equation 2 is abandoned as the fluid is known to be a single phase of composition z.

In order to set up the method, it is necessary to know what the K-values should be. This involves using a rigorous method to calculate the phase equilibrium and hence the K-values at a small number of points, and then regressing the K- values explicitly for use in Equations 2, 3 and 4.

In order to calculate the K-values rigorously at a small number of points, in a preferred embodiment the following method is used:

An equation of state defines the pressure of an isolated system of a fluid mixture as a function of the temperature, total volume and the mole numbers of all components present in the mixture, p=p(V,T,n). The equation of state also incorporates a number of model parameters. In most practical situations, the pressure and not the volume is known, so the first step in evaluating the K-values is to solve the equation of state for volume V for each coexisting phase. The residual pressure p* is defined as the difference between the pressure exerted by a real gas and that exerted by an ideal gas at the same volume, temperature and composition

p * (V, T, n) = p(V, T, n)- p^id (V, T, n) (5) The equation of state for an ideal gas is

n,_d- NRT

^P (6) where R is the gas constant and N is the total number of moles in the system

N= y n

The residual Helmholtz energy follows from the thermodynamic relationship p*=-dA*/dV. So

The fugacity coefficient of component i is defined as

1 dA^* , ( pV Λ

lnd = In — 9

' RT dn_t [NRT J

When two phases are in equilibrium, the components partition between the phases such that the free energy reaches a minimum. From the laws of thermodynamics, it can be shown that this is equivalent to requiring that

(j>_Viy_i = (j>_Lix_i for all components i.

Here, the equilibrium is expressed as being between a vapour and liquid phase: φ_ν,, φυ, are the fugacity coefficients for the vapour and liquid phases and y_h X, are the mole fractions for the vapour and liquid phases. It follows that the partition coefficient or K-value is iven by

The Soave-Redlich-Kwong (SRK) equation of state is used below as an example. Alternatively, other equations of state such as the Peng-Robinson (PR), GERG or SAFT equations could be used.

The SRK equation gives the pressure as

^{NRT a}

P = ? \

V - b V(V + b)

where the equation of state parameters a and b are defined by the mixing rules

' ' (13)

The subscripted values of a and b refer to their pure-component values. Those without subscripts refer to the averaged values used to describe a fluid mixture in Equation 1 1. k is the binary interaction parameter between components i and j.

In order to proceed further, Equation 1 1 must be solved for volume for the specified pressure, temperature and composition. In general, because an equation of state is not explicit in volume, the solution can only be found numerically by a suitable iterative technique. However, the SRK equation is a so-called cubic equation of state which can be rearranged into a cubic polynomial in volume

From the standard method for solving cubic equations, V can therefore be found analytically. Using Equations 6 and 8, the residual Helmholtz energy A* for the SRK equation is found to be

A^* = NRTln\—— | + -lnf—— | (15)

\y - b) b y+ b)

where the constant of integration is set to zero, which corresponds to setting the residual Helmholtz energy to zero at zero pressure.

It follows from Equation 9 that the fugacity coefficient for a given phase is

(16) where Sa/Sn, and Sb/Sn, are found from the mixing rules, Equations 12 and 13.

To summarise, the procedure for rigorously evaluating K-values from the SRK equation from the temperature, pressure and compositions of the vapour and liquid phases is as follows:

1. from the vapour composition calculate a and b for the vapour phase from Equations 12 and 13;

2. solve Equation 14 for the vapour volume; 3. substitute a and b and the volume in Equation 16 to get the vapour fugacity coefficients;

4. repeat steps 1-3 for the liquid phase;

5. use Equation 10 to obtain the K-values.

An alternative method for rigorously calculating K-values involves regressing them to match physical data. However, this is not usually possible because the amount of available data at the relevant conditions is normally very limited.

As discussed above, it is also necessary to calculate the K-values above the bubble-point (or dew-point) rigorously so that they can be correlated for use in the explicit K-value method. The K-values have to be obtained using the specified feed with the intention of then applying them to different although similar feed

compositions.

In order to obtain a smooth extrapolation of K-values above the bubble point, the same formal procedure as that used below the bubble point is followed. This is achieved by allowing F_v in Equations 2, 3 and 4 to become negative (but still less than or equal to one), but otherwise following the same iterative procedure to find the vapour-liquid phase equilibrium with K-values rigorously evaluated from an equation of state.

Similarly, in order to obtain a smooth extrapolation of K-values above the dew point, the same formal procedure as that used below the dew point is followed. This is achieved by allowing F_v in Equations 2, 3 and 4 to be greater than one (but not negative), but otherwise following the same iterative procedure to find the vapour-liquid phase equilibrium with K-values rigorously evaluated from an equation of state.

In a normal phase equilibrium calculation, the point in composition space corresponding to the feed vector z lies on a straight line joining the vapour and liquid composition vectors y and x. z also lies somewhere between y and x. As the pressure rises, x moves towards z until at the bubble point x and z coincide. Above the bubble-point pressure, where F_v<0, z still lies on a line joining y and x but now z no longer lies between x and y. This situation is termed a projected phase equilibrium, also called a 'negative flash'. The solution is not physically realistic as a negative vapour fraction is not directly meaningful, but it does correspond to an actual phase equilibrium that could result from a different feed composition. The projected phase equilibrium can therefore be used to generate K-values above the bubble point (or, similarly, above the dew point) as inputs for an explicit K-value method.

As the pressure continues to rise, the vapour and liquid phases in the projected phase equilibrium will become increasingly similar until a pressure is reached where the two phases will just become identical. This point is referred to as the projected critical point. It is a critical point for the fluid mixture arising from the projected phase equilibrium calculation (which is not the same as the feed composition). In other words, the projected critical point may be thought of as the critical point for a mixture with a different ratio of components to that in feed z.

The line joining composition vectors y and x is termed a tie line. As the critical point is approached the points y and x move together and the tie line becomes a critical tie line. The critical tie line is the limiting tie line as the separation between y and x tends to zero. The condition for a projected critical point is that the feed composition vector z must lie on the critical tie line.

Above the pressure of the projected critical point, p_pc, no further phase equilibria are possible as the K-values all converge to unity at the projected critical point. Consequently, the dashed line in Fig. 1 corresponds to the line of projected critical points. The actual critical point is the point where the critical mixture and the feed are the same. Below the critical temperature, the critical mixture is more liquid- like than the feed, and above the critical temperature it is more vapour-like.

As the pressure rises towards p_pc, the algorithm to calculate phase equilibrium from a rigorous thermodynamic model becomes difficult to apply because the vapour and liquid compositions become increasingly similar. At some point the algorithm will fail numerically because the compositions y and x will be assigned nearly identical values, and the procedure will give a trivial solution, a solution which shows that a phase can always be in equilibrium with itself.

However, because of the importance of the projected critical point in defining the K- values as functions of pressure (see below), it is very desirable to use a direct and reliable method to determine p_pc.

A critical point can be defined as a property of a single phase. The first property that it must exhibit is that the mixture is at the limit of stability, i.e. that it is about to split into two phases. Matrix Q is defined, the elements of which are

where A is the Helmholtz energy of the mixture given by an equation of state ,n is the mole number and subscripts i and j denote the components in the mixture.

A is a function of temperature T, total volume V and the vector of mole numbers n. In Equation 17, the total volume is constrained so the mole numbers and hence the mole number derivatives are all unconstrained and linearly independent. The mole numbers do not necessarily sum to unity, so they are not mole fractions.

The matrix Q is real and symmetric so it will have real eigenvalues. Using standard algorithms in linear algebra, the minimum eigenvalue A_min and associated eigenvector u are determined which satisfy

Q »= ^λ, , ¹¹ (18) The eigenvector is a direction vector and for convenience it is normalised so that u^r u= 1 (19) If A_min is positive, the phase is stable, and if A_min is negative the phase is unstable. The limit of stability occurs when

A_m„= u^rQ u= 0 ₂₀j

When all the critical conditions are satisfied, u will give the direction of the critical tie line.

The first step therefore in determining the projected critical point is to find a mole number vector n for which Equation 20 is satisfied, and for which the feed composition vector z lies on the critical tie line η+βιι, where β can be any value.

In order to do this we first need to find the derivatives of A_min with respect to the mole numbers by differentiating Equation 20. However, matrix Q in Equation 20 is already formed of second derivatives of the Helmholtz energy, and generally in engineering applications higher-order derivatives of the Helmholtz energy are avoided. Noting that the value of A_min in Equation 20 is equal to the curvature of A in the direction of the eigenvector u, the curvature can be approximated by using a finite difference approach.

Defining the Helmholtz energy as a function of displacement along the eigenvector n+au

_ _T (VA(Aa) - VA(- Aa))

m^m 2Aa ^{{ V} where V is the vector operator the i^th component of which is δ/δη,. Δα is chosen to be a suitably small value so that Equation 21 is a realistic approximation to

Equation 20.

Differentiating Equation 21 with respect to mole numbers then gives

_ T (Q -Q(- Aa)) _ (Q(Aa)-Q(- Aa))

In order to obtain an improved estimate of the composition n of the projected critical point, Newton's method can be applied to any component i as follows:

new ^ min However, as Newton's method essentially amounts to linearising the equation, it is equally possible to choose any linear combination of the N equations 23, where N is the total number of components in the mixture. These linear combinations define a hyperplane given by the equation

-_ d K mm / new \ \

L. _d—(ⁿ - - ⁿ i) - - ^K _mm (24) Moving to any point on the hyperplane is therefore equivalent to making one step of Newton's method to solve Equation 20. Making the further assumption that the eigenvector u is approximately independent of composition, which is analogous to assuming that K-values are approximately independent of composition, it follows that the new estimate of the vector n is that point on the hyperplane, Equation 24, that intersects a straight line through the feed composition in the direction u, i.e. on ζ+βιι where β can assume any value.

n can be therefore be found by the following iterative procedure:

1. estimate the mole number vector n of the projected critical fluid;

2. form the matrix Q;

3. use a standard procedure to solve Equation 18 to find the minimum eigenvalue A_min and its eigenvector u;

4. calculate the derivatives, Equation 20, and hence the hyperplane, Equation 22;

5. extend the feed composition z in the direction u to intersect the hyperplane to determine the new mole number vector n;

6. if n is not changed significantly, exit as the problem is solved;

7. go back to step 2. A critical point is defined by two criteria. The first is that the fluid is on the limit of stability which corresponds to A_min=0; the second is that in the direction of the eigenvector u, the derivative of the curvature of A is zero, or in other words, the third derivative of A is zero.

This is based on assuming that the volume of the system is fixed. An alternative approach could be based on assuming that the total number of moles is constant.

From Equation 20, A_min represents the curvature of A in direction u. Equation 22 therefore defines the mole number derivatives of the curvature, hence the derivative of the curvature in direction u is

c _{= u}r_{vx = u}r fe - g(- A_g))

2A

The second critical condition is therefore C=0. To satisfy C=0 the total volume V is varied iteratively until the condition is satisfied. The procedure is:

1. Estimate V and n from previous successful (projected) phase equilibrium calculations;

2. using the iterative procedure above, solve for the mole numbers n that satisfy the conditions of projected limiting stability, A_min=0;

3. use Equation 25 to calculate C;

4. increment V;

5. using the iterative procedure above, solve for mole numbers n that satisfy the conditions of projected limiting stability, A_min=0;

6. use Equation 25 to calculate new value of C;

7. using the two most recently formed values of V,C extrapolate to a new value of V to satisfy C=0 by linear extrapolation;

8. using the iterative procedure above, solve for mole numbers n that satisfy the conditions of projected limiting stability, A_min=0;

9. if V has not significantly changed from its most recently formed value, exit as calculation is complete;

10. use Equation 25 to calculate C;

11. go to step 7.

The procedure completes the determination of the projected critical point. The mole numbers n and total volume V can be substituted in the equation of state being used to give the projected critical pressure p_pc, and any other desired physical properties. The outlined iterative procedures are for illustrative purposes to demonstrate how the equations may be solved, but other methodologies are possible. Any procedure can be used that satisfies the conditions A_min=0 and C=0 while also imposing the constraint that the feed composition z lies on the critical tie line.

The next step is to regress the K-values as explicit functions of pressure on a given isotherm.

Since all the K-values on a defined isotherm converge to unity at the projected critical pressure p_pc, the pressure always needs to be related

appropriately to p_pc which is dependent on temperature. To achieve this, all relations are expressed in terms of a reduced pressure p_r defined as

^Pr p (26)

Conventionally, the reducing factor used for pressure is the normal critical pressure. Equation 26 is therefore a novel method for defining reduced pressure which is particularly advantageous for correlating K-values in a way that can extrapolate reliably between different temperatures.

At low pressures, liquid solution behaviour can be expressed by the standard relationship

Py, = Po, V , x, (27) where p is the pressure, p_0i is the pure-component saturated vapour pressure, γ, is the activity coefficient, x, is the liquid mole fraction and y, is the vapour mole fraction. Subscript i denotes the component number.

It follows that the K-values are defined by

K -— (28) or, since the activity coefficients of liquids are almost independent of pressure, at low pressures we may write

\n(p_rK_i ) = \η(ρ_ϋι γ. I p_pc ) = const (29) for a fixed temperature.

To extend the relation to higher pressures, Equation 29 can be converted to a polynomial in p_r, or some similar algebraic expression

^ {p_rK_t)= c_i0+ c_tlp + c_i2p²+ ... + c_im p^m _r (30) which has the correct limiting behaviour that

\n ( p_r K_t)= c_{0 (31 )} as p_r→0.

As the pressure approaches p_pc or p_r→1 , all the K-values tend to unity, so ln(p_rK|)→0 for all components.

To impose this condition, a switch function can be introduced into Equation 30. Defining q=1-p_r

p_rK_t ) = (c_m + c_aPr + c_i2p_r ² + ... ₊ c_imp )q* (32) or some other appropriate function. The power φ is set, for example by trial and error, to a suitable value appropriate for the equation of state used to generate the rigorous values. Provided φ<1 , Equation 32 shows the correct property that the derivatives of the K-values with respect to pressure are all infinite at the projected critical point.

The method as set out assumes that φ is fixed for all components and all temperatures. Typically a value of around 0.7 is found to be suitable for the SRK and PR equations.

The whole procedure for determining explicit K-values for a given isotherm can be summarised as follows:

Using a suitable equation of state, rigorous phase equilibrium calculations are performed to determine the K-values at a representative number of pressures lying within the two-phase region. The calculations are extended to higher pressures using the projected phase equilibrium methodology. When the pressure gets sufficiently high for the projected phase equilibrium algorithm to fail, the projected critical point algorithm is used to determine the projected critical pressure, p_pc. All pressures are then reduced in relation to p_pc. The K-values for each component are then regressed as a function of reduced pressure using Equation 32 or some similar expression. The number of coefficients Cii-c_im is kept to the minimum necessary to represent the K-value to the desired accuracy. The K- values for the isotherm are then explicitly represented by the coefficients c and the projected critical pressure p_pc.

In some situations it will be desirable or necessary to extrapolate the K- values above the retrograde dew-point pressure.

For temperatures greater than the critical temperature, as the pressure rises the vapour-liquid region of the phase diagram ends at a dew point not a bubble point. The situation is entirely analogous to the case of extrapolation above a bubble point except that the role of the vapour and liquid phases is reversed; the liquid fraction goes to zero. A projected phase equilibrium calculation can be performed as described above except that here F_v>1. The projected critical pressure can be calculated using the same algorithm as for extrapolation above a bubble point.

As the temperature rises, the normal lower dew point increases in pressure giving rise to a second problem, namely a lack of phase equilibria at low pressures. It is possible to provide K-values for the one-phase region below the lower dew point pressure also by using projected phase equilibria. As the pressure is reduced, the liquid fraction at the lower dew point goes to zero and then F_v>1.

The resulting correlations of the K-values as functions of pressure have the same functional form as for the bubble point case, and can be treated in an entirely equivalent way in what follows.

In some situations, it will be desirable or necessary to include an aqueous phase to the vapour-liquid equilibrium calculations.

In the oil industry, water and a number of associated components are frequently present in oil and gas streams. The possibility of the formation of an aqueous phase should therefore also be considered in many cases. Because of the very low mutual solubilities of water and hydrocarbons, many situations can be successfully approximated by assuming that water and its associated components, termed aqueous components, are not present in a liquid hydrocarbon phase, and hydrocarbons and associated components, termed non-aqueous components, are not present in the aqueous phase.

To support these calculations, if an aqueous phase forms during rigorous calculations with an equation of state, the K-values that need to be regressed for the aqueous components are the K-values between the vapour phase and the aqueous liquid phase. The K-values for the aqueous components between the vapour phase and the liquid hydrocarbon phase are ignored. Because the water phase will not show critical behaviour at any pressure or temperature of interest to the oil industry, the K-values for the aqueous components can be regressed using Equation 30 and there is no need to determine the projected critical point for the aqueous phase.

Making the assumptions set out above, the phase equilibrium relation for the explicit K-value method, Equation 2, when aqueous components are present becomes

^ ,+ 1- „- £ ^- /· A', ++ F /·_w

where subscript i denotes non-aqueous components and subscript j denotes aqueous components. F_w is the aqueous phase fraction. The phase compositions for the vapour phase are then, y' F_vK_t+ 1- F_v- F_w (³⁴)

_ κ,^ζ,

y' F_vK_j+F_w

for the liquid hydrocarbon phase

z_t

X -

F_rK,+ 1- F_r- F_w (³⁵) and for the a ueous phase

In Equations 33-36, F_v, F_w and 1-F_V -F_w may not be negative. If a solution occurs for which F_v<0, then there is no vapour phase and Equations 33-36 reduce to

X,-

(l-F_w) (³⁸)

If a solution arises for which 1-F_V-F_W<0, there is no liquid hydrocarbon phase and Equations 33-36 reduce to ∑TTJ-⁺∑_{l-_Fw)Kj+Fw-» (40) y'⁼^F^~ _w (⁴ ) y' {\-F_w)K_]+F_w (⁴²) x ,=

{ \- F_w)K_{J +} F_w (⁴³)

If a solution arises for which F_w<0, there is no aqueous liquid phase and Equations 33-36 revert to Equations 2-4.

If two of the quantities F_v, F_w and 1-F_V -F_w are negative, there is only one phase present which must have a composition corresponding to the feed z.

The calculation of phase equilibria is the most important and time- consuming part of physical property simulation. However, once the phase equilibrium has been computed there is usually a need to calculate a number of thermodynamic and transport properties of the coexisting phases. The procedure for this can be treated as an extension of the method of correlating the K-values with the proviso that suitable correlating functions should be used to reflect the correct limiting behaviour of each property.

For any extensive property such as total volume V, it is always the case for a mixture with mole numbers n that y-∑ ⁿ' ^_t (44) or setting the total amount of substance to one mole

^{V =}∑ ' ^V' (45) where v denotes molar volume and x mole fractions, v, = δν/δη, is the partial molar volume of component i, and, like the K-values, it is only weakly dependent on composition. The volumetric properties of a fluid can therefore be correlated by storing the partial molar volumes together with the K-values.

At low pressures all vapours obey the ideal gas law which predicts that the partial molar volume of all components is RT/p; for liquids the partial molar volume is finite at zero pressure and shows a small dependency on pressure.

At the projected critical point, the vapour and liquid partial molar volumes become the same because the liquid and vapour phases become identical at this point. When the algorithm to find the projected critical point converges, the mole number vector n of the critical fluid is determined. It is therefore possible to use this composition in the equation of state model to calculate the partial molar

thermodynamic quantities at the projected critical point. The projected critical partial molar volume of component i is denoted by v_pc,i. These considerations suggest functions similar to the following would be suitable for partial molar volumes.

For the vapour phase

and for the liquid phase

^V _l = (^CW + ^CuPr+^C _l2Pr²+-^+C _lmP™)q +^^~<1 )^V _Pc,_t (47) ζ is set to an appropriate value which is kept constant for all partial molar volumes from a particular equation of state. Provided that ζ<1, the pressure derivatives of the partial molar volumes are correctly predicted to be infinite at the projected critical point. As with the K-values, the coefficients c have to be fitted to rigorously calculated values of each component's partial molar volume.

In addition to the vapour and liquid molar volumes, the molar volumes should be correlated for cases where the phase equilibrium calculation predicts that there is a single phase. In order to do this, the partial molar volumes of the feed as a single phase must be regressed. As there is no critical phenomenon arising in this case, the correlating equations should not contain the switch function. So, for the vapour region (at low pressures)

v_t={RT I p){\ + c_i2p²+ c_i3p + ...+ c_imp^m _r) (48) and for the liquid region (at higher pressures)

V'⁼ (^C'0^{+ c}ilPr⁺ C_i2p² _r+ ...+ C_imp^m _r) (₄g)

Equation 49 is also suitable for the aqueous liquid phase.

The molar volume for the mixtures is then found from the partial molar volumes using Equation 45.

Like volume, thermal properties such as enthalpy and entropy are extensive thermodynamic properties and can be correlated with a similar treatment.

In the ideal gas state, partial molar enthalpy is independent of pressure, so a suitab our phase might be

-- ^{+ c} _imP™)^cl ⁺{^l-^ciⁱ)^h _Pc,₁ (50) while for the liquid phase ^={^cw^+cnPr^+caPl+---^+c _imP™)^cl ⁺{^l-^ciⁱ)^h _Pc,₁ (51)

For the single-phase vapour case

h,= {c_i0 ⁺ c_i2p²+ c_i3p + ...+ c_imp^m _r) (52) while for the single-phase liquid case including the aqueous liquid phase h_t= (c_i0+ c_ap_r+ c_i2p² _r+ ...+ c_imp^m _r) (53) Finally, the molar enthalpy h of a phase is given by

h= L. x < A < (54) Molar heat capacity C_p can in principle be obtained by differentiation of the enthalpy with respect to temperature. However, differentiating a correlated representation of a property can be unsatisfactory because the systematic deviations introduced by the correlation method will distort the derivatives. It is therefore preferable to correlate the rigorously calculated values of C_p directly. The same functional forms used for the enthalpy can also be used for C_p.

Entropy is more complicated as it shows a marked pressure and composition dependence. The partial molar entropy in the ideal gas state has the form

_S;^d=f(T)-RH^) ₍₅₅₎ where p₀ is the entropy datum reference pressure.

In order to obtain a quantity that is only weakly dependent on pressure and suitable for correlation, we define

S=s+Rln{^) ₍₅₆)

For the vapour case, S, can be correlated with a function of the type s₁=(^c ₁o^+c ₁2Pl^+c ₁3Pl⁺-^+c _imP^h ⁺(^l-^)^S _Pc,₁ (57) for the liquid case

for the single-phase vapour

S,= C_i0+ C_i2p²+ C_i3p + ...+ C_imp^m _r (₅g) and for the single phase liquid and aqueous liquid

S,= c_i0+ c_up_r+ c_i2p²+ ...+ c_imp^m _r (50)

The molar entropy of the whole phase is then given by

5=∑ x,(_?-/to(-¾) ₍₆₁)

Unlike the thermodynamic properties above, which are all extensive, the transport properties are all intensive. osity η, for example, it follows that

A different approach is therefore needed to represent transport properties in a way that can account for variations in feed composition. Taking the simpler single- phase case first, we can introduce the quantity

where η₀ is the viscosity corresponding to the defined feed composition with mole fractions denoted by vector z₀. The viscosity of a single phase with mole fractions z can then be evaluated from

n= n.⁺∑ ( ¾) ₍₆₄₎

In the two-phase vapour and liquid cases, the situation is more complicated because the composition of each phase is also changing.

Taking the vapour phase, η₀ now denotes the vapour viscosity at composition y₀ corresponding to the defined feed composition z₀. Equation 64 becomes

n= n₀ ⁺∑ n, ( ,- o,) ₍₆₅₎

The vapour phase fraction F_v will have been found in the phase equilibrium procedure by solving Equation 2, so y, -y₀i can be approximated from Equation 3 to give

K, (z - z_0i)

yr y₀,= F _rVKK,,++ I1"- FFy- (66)

Alternatively, rigorously calculated values of y, and y_0i can be used.

For the liquid phase, a similar procedure may be used. η₀ now denotes the liquid viscosity at composition x₀ corresponding to the defined feed composition z₀. Equation 64 becomes now

n= n₀ ⁺∑ n, ( , - ^,) ₍₆₇₎ and from Equation 4 To correlate viscosity compositionally as a function of pressure, it is therefore necessary to regress both η₀ and the η, as functions of pressure. A suitable correlating function for liquid and vapour cases is

n = (^C ₁0 ^{+ C} ₁l Pr ^{+ C} ₁2 Pl ⁺ - - - ^{+ C} _{l m} P h^{e +}(^l -<i^e _Pc , ₁ (69)

which is also suitable for η₀. Θ is given a constant value for all viscosity correlations.

For the single phase case where critical effects are immaterial, Equation 69 can be simplified to

Π_/= ( ^cw⁺ c_tlp_r+ c_i2 p²+ ...+ c_imp^m _r ) (70)

Thermal conductivity may be treated in an analogous way to viscosity.

Interfacial tension (I FT) is different because it is a single quantity that only exists when two phases are present in equilibrium, so only the two-phase case needs to be considered. The IFTs can be evaluated from

where o, represents the variation in I FT that arises from variations in the feed composition σ - - (72) Finally, σ, can be correlated with pressure using an equation such as, ai ⁼ (^CW ^{+ C}il Pr ^+Ci2 Pr^{2 +} - ^{+ C} _im P^{" 9} (73)

since the I FT at the projected critical point must be zero; and likewise for o₀.

The explicit function method of the present invention may be extended to include additional factors.

As set out above, the explicit function method is based on using simple expressions for the K-values to accelerate phase equilibrium calculations. The K- values are expressed as functions of pressure at fixed temperature that are regressed to the results of rigorous phase equilibrium calculations performed by an equation of state. The K-value functions can then be interpolated between temperatures by cubic splining, for example. The method is valid so long as the feed does not change to such an extent that the K-value functions are no longer an acceptable approximation to the actual K-values for the altered feed. However, the methodology can be extended by introducing other factors as splining variables besides temperature. These extra factors have the role of accounting for changes in the K-values arising from changes in the feed composition, and likewise in all the other physical property correlations.

The conventional approach to cubic splining, which can be used in the present invention, will now be described.

In order to simplify the expressions in what follows, it is assumed that all variables in the splining process are subject to a linear transformation. The splining points are always at the points where the transformed variables are equal to -1 and +1. To spline a function in one variable, f(x), the procedure is well understood: the function is represented by a cubic polynomial

f (x )= a₀+ a^ x + a₂ x²+ a₃ x³ (74) where x is the transformed variable that runs from -1 to +1 in the splining interval. The coefficients a₀, a a₂, a₃ are determined by fitting them to the values of the function and its x-derivatives at the splining points, i.e. to f(-1 ), f(1), f_x (-1) and f_x (1), where f_x is an abbreviation for df/dx.

Fig. 3 illustrates the procedure to extend the splining to a function of two variables, f(x,y).

In order to define the splining function at point P, we first need to find values of the function and its x-derivatives at points A and B, i.e. we require values for f(- 1.y). f ( , y) , fx(-1 ,y) and f_x(1 ,y). These four quantities can in turn be found by performing a cubic spline in the y-direction. f(-1 ,y), for example, can be found from a cubic polynomial in y fitted to f(-1 ,-1), f (- , 1), f_y(- ,-1) and f _y(- 1 , 1); an analogous expression gives f(1 ,y). However, to set up a cubic polynomial in y to determine f_x(-1 ,y), it is necessary to know the values of f_x(-1 ,-1 ), f_x(-1 , 1), ^(-1 ,-1) and fxy(-1 , 1) where fxy is the second derivative of f with respect to x and y; determining f_x(1 ,y) is similar. Note that if the role of the variables x and y are interchanged and the function at point P is found by splining the values and y-derivatives as points X and Y, the result is the same as before.

So, in order to extend the cubic splining method to two variables, it is necessary to know not only the values of the function and its first derivatives at the four splining points, (-1 ,-1 ), (1 ,-1), (-1 , 1) and (1 , 1), but also the second cross- derivatives as well.

It follows that, to extend the splining procedure to a function of three variables, third derivatives would be required, etc. Consequently, the need for higher-order derivatives limits the usefulness of conventional splining methodology in cases when several variables are involved.

As explained above, the conventional splining method requires higher-order derivatives at the splining points as more variables are introduced. This can make the process complicated. As an alternative to conventional splining methods, a new method of splining may be used instead that only requires first derivatives regardless of how many variables are used. This will now be described.

It can shown that the values of the coefficients a₀, a^ a₂, a₃ for single- variable splining are

(2/(l)+2/(- l)- /,(!)+ /,(-!))

ao= 4 (75)

(3/(l)-3/(- l)- f_x(l)- /,(-!))

,=

¹ 4

(/,(!)- /,(-!))

a2= 4

(/(-i)-/(i)⁺ (ΐ)⁺ (-ΐ))

³ 4

in order to make Equation 74 satisfy the values f(-1), f(1), f_x (-1) and f_x (1).

Combining Equations 74 and 75 gives

f(x)= S_(x)f(-1)+ S₊(x)f(l)+ T_{x)f_x{-\)+ ₊(^χ)/_¾(1) ₍₇₆₎ where we introduce four switch functions

e _{/ \}- (2^_ 3 x+ x³)

SAx)-- ^K (77)

T ^rr I (x \)-= ¹ (1^{~ x~ χ2+ χ3})

- ' 4 τ_Λχ)= (- ι-*÷*²⁺*³)

Equation 76 can be generalised to a splining function in two variables as follows f(x,y)=S(x, y)+ S_x(x,y)+ S_y(x, y) ₍₇₈₎

S(x,y)= S_(x)S_(y)f(-l -1)+ S_(x)S₊(y)f (- 1,1)

+ S₊(x)S_(y)f(l - 1)+ s₊(x)s₊(y)f(i,i)

S_x(x,y)= T_(x)S_(y)f_x(- l -1)+ Γ_(χ)5₊( )/_χ(- 1,1 )

+ ( )^-(^)Λ(ΐ - i)⁺ T₊(x)s₊(y)f_x( i) S_y(x,y)= S_(x)T_(y)f_y{- 1 - 1)⁺ S_(x)T₊(y)f_y(- 1,1)

+ S₊(x)T_(y)f_y(l - l)+ S₊(x)T₊(y)f _y(l,l)

Equation 78 has the following properties:

1. the function f(x,y) is smooth throughout the domain -1≤x≤1, -1≤y≤1;

2. if the same method of splining is used for adjoining domains in the variables, the splining function joins the splining functions of the adjoining domains smoothly;

3. the function reproduces the values of the function and its first derivatives at the splining points, i.e. f(-1-1), f(-1,1), f(1-1), f(1,1), f_x (-1,-1), f_x (-1,1), f_x (1,-1), f_x (1,1),f_y(-1,-1),f_y(-1,1),f_y(1,-1)andf_y(1,1);

4. the method does not require the value of any function derivatives above first order at the splining points.

Here, the term "smooth" means continuous in the value of the function and all first derivatives.

The method can be generalised further albeit with mounting complexity. For three variables the equations become

f{x,y,z)=S{x,y,z)+S_x{x,y, z)+S_y (x , y , z)+S_z (x , y , z) ₍₇₉₎

S(x,y,z)=S_(x)S_(y)S_(z)f(- 1 - 1,- 1)+ S_(x)S_(j/)S₊(z)/(- 1 - 1,1) + S_(x)S₊(y)S_(z)f(- U,- 1)+ S_(x)S₊(y)S₊(z)f(- 1,1,1)

+ S₊(x)S_(y)S_(z)f(l,- 1 - 1)+ S₊(x)S_(y)S₊(z)f(l - 1,1)

+ S₊(x)S₊(y)S_(z)f (1,1,- 1)+ S₊(x)S₊(y)S₊(z)f(l,l,l)

,y,z)= T_(x)S_(y)S_(z)f_x(- 1 - 1,- 1)+ T_(x)S_(y)S₊(z)f_x(- ^~

+ T _(x)S ₊(y)S _(z) f _x(- 1,1 - 1)+ T_(x)S₊(y)S₊(z)f _x(- 1,1,1)

+ T₊(x)S_(y)S_(z)f_x(l,- 1 - 1)+ T₊(x)S_(y)S₊(z)f_x(l - 1,1)

+ T x)S₊(y)S_(z)f_x(l,l - l)+T₊(x)S₊(y)S₊(z)f_x(l,l,l)

+S_{x)S₊{y)T_{z)f_z{-l,l -l)+S_{x)S₊{y)T₊{z)f_z{-l,l,l)

+S₊(x)S_(y)T_(z)f_z(l -l -l)+S₊(x)S_(y)T₊(z)f_z(l -l,l)

+s₊(x)s₊(y)T-(z)f_z(i,i -i)+s₊(x)s₊(y)T₊(z)f_z(w) The scheme can be extended in an analogous way to handle any number of variables.

In what follows below, it is assumed that the new splining procedure is used, whenever there is more than one variable.

The K-values may be extended to different temperatures as described below.

In order to apply the explicit K-value method of the present invention efficiently at any arbitrary temperature, it is necessary to establish a grid of temperatures at which the K-values are regressed as explicit functions of pressure, as described above. The regressions may be performed adaptively in the course of a simulation as needed to satisfy the requested phase equilibrium calculations. In other words, the regressions need only be performed at the p, T values where the simulation requires a K-value.

The method of calculation at a given temperature T involves extrapolating the K-values from the highest temperature on the grid of regressed K-values for which Τ>Τ_Ϊ and from the lowest temperature T₂ for which T<T₂. A good method of interpolation is cubic splining as this preserves smoothness.

Starting with p_pc as an illustration, the property is expressed as a cubic function of temperature of the form

P_Pc {T)= a₀+ α_λ Τ+ a₂ T²+ a₃ T³ (_Q0)

The four coefficients a₀, a^ a₂, a₃ are found by fitting Equation 80 to four values at the splining points, namely _pc(T^, dp_pc(J^/dJ, p_pc(T₂) and <3p_pc(T₂)/aT.

The temperature derivatives are a necessary part of this cubic splining method, and they may be estimated by finite difference as follows. When the regression method for the K-values described above is applied at any grid temperature T_n, the whole procedure is repeated at a similar temperature T_n +ΔΤ_Η where ΔΤ_Η is a small value. The temperature derivative is then

d Ppc _ Ρρ ^Τ η^{+ T}n )~ Ρρ ^Τη)

d T Δ Γ„ (⁸¹ ) where ΔΤ_Η is chosen so that Equation 81 provides a good approximation to the temperature derivative.

The method described for splining p_pc is also applied to each of the coefficients c for the K-values in Equation 32. All the splining coefficients can be stored in computer memory for the temperature interval T T₂ so that no further regression calculations are required for this interval during a simulation. As a result, each K-value in interval T T₂ is determined as an explicit function of pressure and temperature, Kj=Kj(p,T).

In order to calculate the phase equilibrium, Equations 2, 3 and 4 need to be applied only once to determine the phase fractions and the compositions of the vapour and liquid phases. All the other explicit physical property functions can be extended to cover temperature by splining their coefficients in the same way as the K-values.

The explicit K-value method may also or alternatively be extended to account for variation in feed composition as described below.

In a sequential modular simulation, each stream of a flow sheet may be represented by a set of explicit K-values and associated physical properties defined as described above. In a situation where significant changes in feed composition occur for a stream, a new set of explicit functions may be set up. The K-values and other properties will be subject to a discontinuous change but that may be acceptable provided it is not too large and the flow sheet continues to converge. The important consideration is to set the criterion for regressing new explicit functions so that it does not happen too often but convergence still occurs.

For simulations where smoothness is required, for example in an equation orientated simulation or for process optimisation, periodic recalculation of the explicit functions for the K-values, etc., may not be acceptable at it introduces discontinuities. Situations where significant changes in feed composition occur can be incorporated into the method by introducing one or more additional variables besides temperature that can be used as splining variables.

A simple example in process simulation is modelling the phase behaviour in a unit operation which has two feeds of different composition where the flow rates may vary. The blend fraction between the two feeds may be introduced as a second variable to complement temperature when regressing K-values. The K-values may be expressed as explicit functions of pressure on a grid of points at different temperatures and blend fractions. As before, the regressions can be performed adaptively as required to allow the K-values to be interpolated to the required temperature and blend fraction. In this case, four splining points are needed at lower and higher temperatures and lower and higher blend fractions. At each splining point the blend fraction derivatives are estimated by finite difference in addition to the temperature derivatives. Equation 78 can then be used to interpolate the projected critical pressure p_pc and each of the coefficients c used to represent the K-values.

The methodology may be extended to unit operations with more feeds. If there are three feeds, for example, then there needs to be two blend fractions and splining must be done using Equation 79 to interpolate the three variables:

temperature and two blend fractions.

The method of the present invention may be applied in simulations of various situations such as a distillation column, a pipeline, or a reservoir or well.

A distillation column is a complex group of interconnected phase equilibria.

However, the fluid mixture on each tray usually varies smoothly up the column. A rigorous steady-state simulation of the column, for example, could provide typical feed compositions on each tray. These may be denoted by height up the column (or tray number). So whereas the previous example used blend fraction as a second variable, the K-values in a column could be related to height or tray number. If a dynamic simulation is to be performed, the compositions of the column are likely to be perturbed from the steady-state values. The explicit K-value method of the present invention can therefore provide an approach to performing a rapid but reasonably accurate dynamic simulation of a column.

Simulating a pipeline has its own particular requirements that reflect the essential geometrical fact that the length hugely exceeds all the other dimensions. As a consequence, the residence time of fluid in the pipeline is high. In practice, many pipelines are part of a network, but the methodology described here may be applied to each single un-branched section of pipe.

Initially, the phase behaviour at the inlet of the pipe may be regressed to form explicit K-values. It is assumed that the composition of the inlet fluid is known over time either from direct measurement or else from simulation of the supply pipes to the inlet point in question. The second variable is ideally time. A grid of temperature and time points may be set up and the K-values regressed as functions of pressure at each grid point. The temperature and time derivatives may be found by finite difference as before. The time intervals are ideally selected to reflect changes in inlet composition; if the composition is changing slowly, the time intervals may be large whereas if the composition is changing rapidly, the time intervals should preferably be short. The K-values at any other temperature or time can then be found by extrapolation using Equation 78. The model assumes that the explicit K-value (and other) functions travel along the pipe at the same velocity as the average velocity of the fluids. Except for steady-state flow, there will be some phase intermingling and so overall

compositions will alter as the fluid flows through the pipe. However, the effect is normally small enough for the assumption to be made that the K-value functions will remain good approximations.

The usefulness of the method is that the same K-value functions may be applied to any point on the pipe by replacing the time variable with a displaced time variable. At any defined point along the pipe at a given time, the displaced time is the time corresponding to when the same fluid was at the inlet. The time difference is calculated by finding from the simulation the average fluid velocity from the inlet to the point in question; the time difference is then the length of the pipe from the inlet to the point divided by the average fluid velocity. The displaced time is defined as the actual time minus the time difference.

Reservoir simulation and well-performance simulations are good potential applications of the explicit K-value method as they notably involve intensive computations. Many reservoirs show quite small variations in temperature, so it is relatively easy to spline the K-values and physical properties at a small number of temperatures to cover the entire simulation. In cases where different zones contain fluid with distinctly different compositions, the approximation method may be extended by representing each fluid composition with its own explicit functions and splining them by introducing blend fractions, just like the case of units operations with multiple fluids. If necessary, gas drives may also be handled by treating the injection gas as another fluid for which blend fractions can be introduced.

Another benefit of the method is that it can be used in conjunction with rigorous methods. For example, in a multilateral well it may be desirable to represent different zones of the reservoir with different approximate functions, but the commingled fluid in the well bore, pipelines or processing facilities could be simulated using rigorous thermodynamic methods.

For some applications a compositional treatment is not required. The proposed methodology may be reduced to handle this case. For example, for a fluid at a defined temperature, the projected critical pressure p_pc can be calculated from a rigorous model as before. All pressures can still be reduced relative to p_pc, as that is still advantageous in providing an appropriate pressure variable for use in the correlating functions. Now, instead of correlating individual K-values, the vapour fraction F_v, and if necessary the water fraction F_w, are themselves expressed as explicit functions of p_r.

With regard to physical properties, only the total physical properties then need to be considered. Taking molar volume as an example, for the vapour and liquid phases in the two-phase region, equations such as 46 and 47 may be used, but now for the molar volume of the whole fluid, not the partial molar volumes. An equation such as 48 can be used for the single-phase molar volumes at low pressure (vapour-like) and 49 for the single-phase molar volumes at high pressure (liquid-like) or for the water phase. An analogous approach may be adopted for the thermal properties. For the transport properties, only the total quantities need to be regressed as in the compositional treatment; the derivatives of the transport properties with respect to composition are not required.

Up to this point, the K-values and physical properties have all been regressed as functions of pressure using simple algebraic functions. However, it is also possible to treat reduced pressure as a splining variable like temperature and any other variables that may be introduced to account for changes in composition such as displaced time in a pipeline, for example.

Taking K-values as an example, consider the case where variations in temperature and pressure must be modelled. A grid of temperatures and reduced pressures should be established. Consider the case of interpolation between temperatures and T₂ and reduced pressures p_r1 and p_a. As before, the projected critical pressure must be calculated using the algorithm above so that the pressures can be converted to reduced pressures. It remains a key advantage of the method that projected critical pressure at any temperature is used as the reducing property. The K-values and their derivatives with respect to temperature and reduced pressure are evaluated using a rigorous equation-of-state model at the four splining points (Ti ,p_ri), ΟΊ ,ρ^), (T₂,p_ri) and (T₂, Pr2)- Whereas before the projected critical pressure p_pc and the regression coefficients Cii-c_im were interpolated, with this approach the K-values themselves are directly interpolated from their values at the splining points using Equation 78 to give values at a specified pressure and temperature. As before, the rigorous values at the splining points can be calculated adaptively as required.

The method can be extended to include modelling solids formation.

For oil industry applications, phase equilibrium problems can include the formation of certain solids such as waxes and asphaltenes. The formation of these solids is operationally very undesirable and should be avoided if at all possible. There are a number of computationally intensive applications available which model the formation of waxes and asphaltenes. However, these take a long time to run. Waxes can deposit in sub-sea pipelines and applying the explicit K-value method of the present invention would reduce the time of these calculations, especially for pipeline networks or for transient simulations. Asphaltenes tend to deposit in oil reservoirs and wells, and these also require large amounts of computing time with the presently available applications.

Extending Equations 33-43 to include a solid phase increases their complexity considerably, and increases computing time. In order to minimise the impact of including a solid phase on computing time, a sequential approximation may be introduced. The approximation is valid for any situation where the components of the solid phase are only found in one of the fluid phases in any appreciable concentration. The concept is illustrated below for the case of wax precipitation, but the method is general and may be used in other situations where solids form.

Petroleum wax is a solid solution predominantly of n-paraffins (i.e. straight- chain alkanes) that occur naturally in petroleum fluids. In nearly all cases, the wax- forming components are only found in the hydrocarbon liquid (or oil) phase. When applying a rigorous model to calculate the phase equilibrium in order to form the explicit fluid-phase K-values, the presence of a wax phase may be predicted. If that happens the solid-liquid K-values

w_ilx_i (₈₂) are stored, where w_h x, denote the mole fractions of component i in the solid (i.e. wax) phase and liquid (i.e. oil) phase, respectively. The solid-liquid K-values can be regressed using an expression such as Equation 30, since a solid phase will not exhibit critical behaviour.

Once the solid-liquid K-values are stored, the sequential approximation proceeds by subsuming the wax components into the liquid hydrocarbon phase, and treating the resulting system as a fluid phase equilibrium involving vapour, liquid hydrocarbon and water phases exactly as before. The effective mole fractions of the components of the hydrocarbon liquid are given by where F_s is the solid phase fraction and F_L is the hydrocarbon liquid phase fraction. The effective hydrocarbon liquid fraction F_L ^eff=F_s+F_L.

Once the explicit K-values have been correlated for the given temperature, etc., the sequential approximation involves first calculating the fluid phase equilibrium as usual with Equations 33-43. The next stage then involves applying Equation 2 to the effective liquid hydrocarbon phase to split it into a solid and liquid hydrocarbon phase. In this, Equation 2 becomes

the resulting real liquid phase fraction would be

the solid phase mole fractions would be

and the hydrocarbon liquid mole fractions

Equation 84 is solved by an iterative numerical procedure, like Equation 2. If during the solution F_s becomes less than zero, the calculation is terminated as no wax will form. In the unlikely event that F_s becomes greater than F_L ^eff, it is concluded that the hydrocarbon liquid phase is to be interpreted entirely as a wax phase.

Like wax, asphaltenes are solid phases that form from petroleum fluids and can be handled using the same solid-liquid formalism. However, one aspect of asphaltene behaviour is different from waxes. Whereas the wax K-values are only mildly dependent on composition, for asphaltenes they are very sensitive to the presence of light components.

Fig. 4 shows a typical phase diagram for an asphaltenic crude oil. Above the bubble point line, the oil phase may or may not form an asphaltene phase. If asphaltene does form, the solid-liquid K-values can be treated in an analogous way to the wax case. However, below the bubble point line, as the light gases leave the liquid phase, the solubility of asphaltene increases rapidly because light

components are destabilising to asphaltenes. The scheme for correlating the solid- liquid K-values to form explicit functions therefore needs to be elaborated for asphaltenes.

An effective procedure is first to calculate the phase behaviour above the bubble point. Any points where asphaltene forms are used to obtain the solid-liquid K-values, and the asphaltene is subsumed into the oil phase exactly as for the wax case. The solid-liquid K-values are regressed as explicit functions of pressure using an expression such as Equation 30. The points below the bubble point line are then processed and the solid-liquid K-values are obtained. As before, the asphaltene phase is subsumed into the oil phase. The new solid-liquid K-values are now split into two parts

Κ?= ΚΓ+ Δ Κ?" ₍₈₈₎ where Kf°" is the correlated function for the K-values already obtained by regressing the points above the bubble point. AKj^SL represents the effect of compositional change and is found from Equation 88.

AK|^SL is then regressed as a simple function of some parameter representing the presence of light components in the oil phase, such as oil molecular weight, or methane mole fraction, etc.

In principle, A ^SL could also be made dependent on reduced pressure p_r, but in most cases that would not be necessary as the lower asphaltene envelope is normally quite close to the bubble point line.

In summary, for asphaltenes, the phase equilibrium procedure consists of solving Equations 33-43 to find the fluid phase equilibrium. If there is no gas phase present, i.e. the pressure is above the bubble point, the solid-liquid K-values are evaluated from a function of reduced pressure such as Equation 30. If there is a gas present, the solid-liquid K-values are formed from Equation 88. The effective oil phase is then split into oil and asphaltene phases using Equations 84-87 in exactly the same way as for the wax case.

Claims

Claims

1. A method of determining K-values for a mixture, the method comprising the following steps:

(i) calculating K-values for the mixture explicitly at a number of pressures and at a temperature T;

(ii) determining a projected critical pressure p_pc as the pressure above which a two-phase vapour-liquid phase equilibrium may no longer be estimated; and

(iii) regressing the K-values as a function of p_pc.

2. A method as claimed in claim 1 , wherein step (iii) comprises determining a reduced pressure p_r, where p_r = p/p_pc and p is the actual pressure, and preferably then regressing the K-values as a function of p_r.

3. A method as claimed in claim 2, wherein the K-values are regressed as a functio

where i represents the components of the mixture, c₀-c_m are coefficients to be determined by the regression, q=1-p_r and φ is a constant for all components which is less than 1.

4. A method as claimed in claim 2 or 3, further comprising regressing a variable other than the K-values as a function of p_r, preferably wherein the other variable is the total volume, enthalpy, heat capacity, entropy, viscosity, thermal conductivity or interfacial tension of the mixture.

5. A method as claimed in any preceding claim, wherein p_pc is determined by:

(a) defining a matrix Q, the elements of which are:

where A is the Helmholtz energy of the mixture, n is the mole number and subscripts i and j denote the components in the mixture; (b) determining the composition n of the mixture and total volume V whilst requiring that the minimum eigenvalue of matrix Q, A_min=0 and the third derivative of A is zero in the direction of the associated eigenvector; and

(c) calculating p_pc from n and V.

6. A method as claimed in claim 5, wherein at step (c) p_pc is determined from an equation of state.

7. A method as claimed in claim 5 or 6, wherein at step (b), n is found by following the iterative procedure:

(b-1) estimate the mole number vector n of the projected critical fluid; (b-2) form the matrix Q;

(b-3) solve Q ^u~ ^™^{« u} to find A_min and its normalised eigenvector u; (b-4) calculate a new hyperplane given by

,_ d h

d n

where

_ _UT (QM- Q(- Aa)) _ (Q(Aa)- Q(- Aa))

2Aa 2Aa

(b-5) extend a feed composition vector z in the direction u to intersect the hyperplane to determine the new mole number vector n;

(b-6) if n is not changed significantly, exit as the problem is solved; else

(b-7) go back to step (b-2).

8. A method as claimed in claim 7, wherein at step (b), V is found by following the iterative procedure:

(b-8) estimate V and n from previous projected phase equilibrium calculations;

(b-9) using the iterative procedure of claim 7, find the mole number vector n;

(b-10) calculate C, where C = u^TVA_mm = u^T (g(^{AQ ~} g(- ^AQ )_u .

2Aa

(b-1 1) increment V;

(b-12) using the iterative procedure of claim 7, find the mole number vector n;

(b-13) re-calculate C;

(b-14) using the two most recently found values of V and C, extrapolate to a new value of V to satisfy C=0 by linear extrapolation;

(b-15) using the iterative procedure of claim 7, find the mole number vector n;

(b16) if V has not significantly changed from its most recently formed value, exit as calculation is complete; else

(b-17) re-calculate C; and

(b-18) go to step (b-14).

9. A method as claimed in any preceding claim, wherein the K-values are regressed on an isotherm at temperature T.

10. A method as claimed in any preceding claim, further comprising determining a vapour phase fraction F_v from each determined K-value.

11. A method as claimed in claim 10, wherein F_v is determined iteratively, preferably using the following equation:

^ ^{K'^{~ l )z}' -= 0

F_vK,+ \- F

where i represents the components of the mixture, K, is the K-value for component and z, is the amount of component i present in the mixture.

12. A method as claimed in claim 10 or 1 1 , further comprising determining vapour and liquid mole fractions from each K-value and F_v.

13. A method as claimed in any preceding claim, wherein the K-values are treated as functions of pressure and temperature but are independent of the composition of the mixture.

14. A method as claimed in any preceding claim, wherein the K-values for the mixture are determined using an equation of state.

15. A method as claimed in any preceding claim, further comprising repeating that method at a different temperature.

16. A method as claimed in claim 15, wherein the K-values are splined between the different temperatures. 17. A method as claimed in claim 16, wherein splining is performed using a cubic splining method.

18. A method as claimed in any preceding claim, wherein the mixture comprises a petroleum fluid and the K-values are preferably used in a simulation of a petroleum fluid.

19. A method as claimed in any preceding claim, further comprising using the K- values in a separator, distillation column, absorber, liquid-liquid separator, other unit operation, pipeline, reservoir or well simulation.

20. A method as claimed in any preceding claim, wherein the K-values are determined adaptively for pressures and temperatures as required in a simulation.

21. A method as claimed in any preceding claim, further comprising estimating K-values for a solid phase.

22. A method as claimed in claim 21 , further comprising modelling wax precipitation or asphaltenes.

23. A computer program for performing the method of any preceding claim.

24. A computer readable medium comprising a computer program as claimed in claim 23.

25. A method of determining K-values for a mixture explicitly at a number of pressures and at a temperature T, wherein the K-values are determined as a function of time, and the determined K-values are then used in simulation of a pipe or pipeline.

26. A method as claimed in claim 25 comprising using the method of any of claims 1-22 to determine the K-values.

27. A method as claimed in claim 25 or 26, wherein the K-values are determined as a function of time based on the flow into the pipe as a function of time.

28. A method as claimed in claim 27, wherein the flow into the pipe as a function of time is measured in real time or estimated from a model.

29. A computer program for performing the method of any of claims 25-28.

30. A computer readable medium comprising a computer program as claimed in claim 29.

31. A method of determining K-values for a mixture comprising vapour, liquid and solid phases, wherein first solid-liquid K-values are determined for the mixture and then the components of the solid phase are subsumed into the liquid phase to form an effective liquid phase before estimating resulting effective liquid-vapour K- values.

32. A method as claimed in claim 31 , further comprising calculating a phase equilibrium between the vapour and effective liquid phases, and preferably then splitting the effective liquid phase into oil and solid phases.

33. A method as claimed in claim 31 or 32, wherein the effective liquid-vapour K-values are estimated using the method of any of claims 1-22.

34. A method as claimed in claim 31 , 32 or 33, further comprising using the K- values for a separator, distillation column, absorber, liquid-liquid separator, other unit operation, pipeline, reservoir or well simulation.

35. A method as claimed in any of claims 31-34, wherein the mixture comprises a petroleum fluid and the K-values are preferably used in a simulation of a petroleum fluid.

A computer program for performing the method of any of claims 31-35.

37. A computer readable medium comprising a computer program as claimed in claim 36.

38. A method of splining K-value functions between splining points in one or more splining variables, wherein the splined K-value function is a function of switch functions and the K-value functions and their first order derivatives only.

39. A method as claimed in claim 38, wherein the splining variable or variables is or comprises temperature.

40. A method as claimed in claim 38 or 39, wherein the switch functions are polynomials, preferably cubic polynomials, of the splining variable.

A method as claimed in claim 38, 39 or 40, wherein four switch functions are

42. A method as claimed in claim 41 , wherein the switch functions are given by:

(2-3 ₊ x³)

- ' 4

(1- x- x²+ x³)

T_{x)--^y '-

_T+{X)= (-i-*÷*²⁺*³) where x is a transformed splining variable which runs from -1 to +1 between the splining points.

43. A method as claimed in claim 42, wherein there is only one splining variable and the splined function is given by

f(x)= S_(x)f(- 1)+ S₊(x)f(l)+ T_{x)f_x{- \)+ ₊(χ)/_¾(1) where f_x is df/dx.

44. A method as claimed in claim 42, wherein there are two splining variables, represented by the transformed splining variables x and y, each running from -1 to +1 between the splining points, and the splined function is given by

f(x,y)=S(x,y)+ S_x(x,y)+S_y(x,y)

where

S(x,y)= S_(x)S_(y)f(-l -1)+ S_(x)S₊(y)f (^~ U)

+ S₊(x)S_(y)f(l,- 1)+ S₊(x)S₊(y)f (1,1)

S_x(x,y)= T_(x)S_(y)f_x(- l -1)+ T_(x)S₊(y)f_x(^~ U )

+ T₊(x)S_(y)f_x(l,- 1)+ T₊(x)S₊(y)f_x(l,l)

S_y(x,y)= S_{x)T_{y)f_y{- 1 - 1)+ 5_(χ)Γ₊(> /,(- 1,1)

+ 5₊(x)r (>>)/,(i - i)+ s₊(x)T₊(y)f _y(i,i)

and f_x is δί/δχ and f_y is 5f/5y.

45. A method as claimed in claim 42, wherein there are three splining variables, represented by the transformed splining variables x, y and z, each running from -1 to +1 between the splining points, and the splined function is given by

f{x,y, z)=S{x , y,z)+S_x{x , y, z)+S_y{x ,y, z)+S_z{x ,y, z)

where

S(x,y,z)=S_(x)S_(y)S_(z)f(- 1 - 1,- 1)+ S_(x) s_(y)s f (^~ i - ι,ΐ) + S_(x)S₊(y)S_(z)f(- 1,1,^" 1)⁺ S_(x)S₊(y)S₊(z)f(- 1,1,1)

+ S₊(x)S_(y)S_(z)f(l - 1 - 1)+ S₊(x)S_(y)S₊(z)f(l - 1,1)

+ S₊(x)S₊(y)S_(z)f (1,1,- 1)+ S₊(x)S₊(y)S₊(z)f(l,l,l)

S_x(x,y,z)= T_(x)S_(y)S_(z)f_x(- 1 - 1,- 1)+ T_(x)S_(y)S₊(z)f_x(- 1,- 1,1) + T_(x)S₊(y)S_(z)f_x(- 1,1 - 1)+ T_(x)S₊(y)S₊(z)f_x(- 1,1,1)

+ T₊(x)S_(y)S_(z)f_x(l,- 1 - 1)+ T₊(x)S_(y)S₊(z)f_x(l - 1,1)

+ T x)S₊(y)S_(z)f_x(l,l,- l)+ T₊(x)S₊(y)S₊(z)f _x(l,l,l)

S₂{x,y,z)=S_{x)S_{y)T_ [z] f _z{-\ -\ -\)+S _{x)S _{y)T ₊{z)f _z{-\ -\,\) +S_(x)S₊(y)T_(z)f_z(-l,l -l)+S_(x)S₊(y)T₊(z)f_z(-lXl)

+S₊(x)Siy)T_(z)f_z(l -l -l)+S₊(x)S_(y)T₊(z)f_z(l -l,l)

+s₊(x)s₊(y)T-(z)f_z(^\ -i)+s₊(x)s₊(y)T₊(z)f₂(w) and f_x is df/dx, f_y is df/dy and f_z is df/dz.

46. A method as claimed in any of claims 38 to 45, wherein the first order derivatives are determined by a finite difference method.

47. A method as claimed in any of claims 38 to 46, further comprising using the splined K-value functions in a simulation of a petroleum fluid, and/or in a separator, distillation column, absorber, liquid-liquid separator, other unit operation, pipeline, reservoir or well simulation.

48. A computer program for performing the method of any of claims 38 to 47.

49. A computer readable medium comprising a computer program as claimed in claim 48.