Algorithm for Constructing Complete Asphaltene ... - ACS Publications

Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

pubs.acs.org/IECR

Algorithm for Constructing Complete Asphaltene PT and Px Phase Diagrams Christian S. Agger and Henrik Sørensen* Calsep A/S, Parallelvej 12, DK-2800 Kgs. Lyngby, Denmark S Supporting Information *

ABSTRACT: An algorithm is presented for constructing complete asphaltene PT phase diagrams. The algorithm works with any equation of state capable of providing consistent fugacity coefficients and appropriate derivatives for gas, oil, and asphaltene-rich phases. SRK, PCSAFT, and CPA were successfully applied. The algorithm has been validated against multiphase PT flash calculations, and the results are further shown to agree with asphaltene simulation results in the open literature. Construction of the asphaltene PT phase diagram is fast and robust even in near-critical regions. Examples are shown illustrating the multitude of possible appearances of an asphaltene PT phase diagram. The algorithm has been modified to take an injection gas amount as specification, which allows for construction of asphaltene Px phase diagrams showing the response of phase boundaries to gas injection. Usage of the algorithm to get a full overview of the behavior of the asphaltene precipitation of the fluid is demonstrated.

■

INTRODUCTION Asphaltenes are large molecules present in hydrocarbon reservoir fluids and are generally thought to be aromatic and to some extent polar due to heteroatoms like oxygen, sulfur, nitrogen, and metals.1 Asphaltenes can be kept in solution by the smaller similar molecules referred to as resins also present in hydrocarbon reservoir fluids, whereas the presence or addition of paraffinic molecules or smaller alkanes such as light gas molecules increases the tendency of precipitation. In situations where asphaltenes cannot be kept in solution, they will drop out along with other components and form a highly viscous separate phase rich in asphaltene components. For simplicity and briefness, such a phase rich in asphaltene is in this paper denoted as an “asphaltene phase” or simply “asphaltene”. The pressure−temperature conditions also influence the ability of the hydrocarbon reservoir fluid to keep asphaltenes in solution. Precipitation can occur due to pressure depletion of the reservoir or due to the drop in pressure and temperature occurring when the hydrocarbon reservoir fluid enters and flows up the well. Asphaltene precipitation is a concern during oil production. The open literature contains numerous examples of production impairment following buildup of asphaltene deposition over time in the near well-bore reservoir and in the well tubing.2−5 It is important to get an overview of the risk of asphaltene precipitation from a hydrocarbon reservoir fluid when exposed to the pressure−temperature conditions experienced during production. The pressure and temperature decrease as the reservoir fluid enters the well and may decrease further in a pipeline and process plant. In the refinery the temperature may, © XXXX American Chemical Society

however, be much higher than in the reservoir. This calls for constructing a complete asphaltene PT phase diagram covering large spans in pressure and temperature. In a PT phase diagram, asphaltene precipitation typically appears as a band on both sides of the saturation pressure line. This band may be continuously open over the whole temperature range or it may collapse onto the saturation pressure line. In the latter case, it may or may not open up again. Gas injection may be applied during oil production and will increase the risk of asphaltene precipitation.6,7 Published asphaltene data and simulation results for petroleum reservoir fluids only cover a limited temperature span, which has led to discussions on whether it can be correct that some data sources report asphaltene onset pressures increasing with increasing temperature, whereas other data sources find the onset pressures to decrease with increasing temperature. Being able to generate the full asphaltene PT phase diagram may help in the investigation of this question and improve the understanding of the mechanisms behind asphaltene precipitation. The full asphaltene PT phase diagram is also an important tool for determining appropriate asphaltene model parameters, in case the initial fluid model does not match experimental asphaltene data. A full asphaltene PT phase diagram reveals how Received: October 12, 2017 Revised: December 11, 2017 Accepted: December 13, 2017

A

DOI: 10.1021/acs.iecr.7b04246 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research far the initial fluid model is from matching the data and is in some cases important for initiating parameter estimation. Lindeloff and Michelsen8 have presented an algorithm for construction of PT phase diagrams for mixtures of hydrocarbons and aqueous components. This paper will look into whether the same calculation principles can be used to generate PT gas−oil− asphaltene phase diagrams. It is further to be investigated whether similar calculation principles can be used to generate asphaltene Px-phase diagrams showing the effect of gas injection on the gas−oil−asphaltene phase boundaries. To the extent needed, the simulation principles will be modified to account for differences between handling an aqueous phase and an ashaltene phase as the third phase, in addition to gas and oil. The paper will accept the general conception that the asphaltene phase is a heavy liquid phase dominated by heavy aromatics.9−11 The same equation of state is used for all phases (vapor, liquid, and asphaltene). Examples are presented using the Soave−Redlich−Kwong (SRK),12 the PC-SAFT,13 and the CPA14 equations of state.

equilibrium point on the two-phase PT diagram commonly referred to as a phase envelope.16

■

ASPHALTENE PT PHASE DIAGRAM ALGORITHM The algorithm presented in this paper is based on similar principles as outlined by Lindeloff and Michelsen8 for construction of PT phase diagrams for mixtures of hydrocarbons and water. However, even though the formulation of the governing equations is similar, and parts of the solution procedure are readily reused, important differences exist as outlined in the following. In particular, stability analysis require different trial phase compositions in the search for an asphaltene-rich phase rather than an aqueous phase. Likewise, the search for an isolated three-phase area is revised. These differences are attributed to the higher mutual solubility between asphaltene and hydrocarbon than between water and hydrocarbon. Figure 1 shows a generic example of an asphaltene PT phase diagram with three-phase points, constructed using the algorithm.

■

TWO-PHASE EQUILIBRIUM To introduce the concepts used to simulate the complex threephase equilibrium (gas−oil−asphaltene), the equations governing two-phase equilibrium are briefly summarized.15,16 All systems spontaneously seek a global minimization of the Gibbs free energy. At this equilibrium state, the chemical potentials (μ) of the N chemical species in the mixture are the same in all phases. For a two-phase equilibrium between a vapor (v) and a liquid (l) phase, μvi = μlι leads to ln K i + ln φi v − ln φi l = 0

i = 1, ..., N

(1)

where Ki = xi/yi is the equilibrium K-factor of component i. xi and yi are the liquid and vapor phase component mole fractions, and φli and φνi , the liquid and vapor phase component fugacity coefficients. Introducing and rearranging the material balance equations, the component mole fractions, xi and yi, are found as zi xi = and yi = K ixi 1 − β + βK i (2)

Figure 1. Generic asphaltene PT phase diagram with three-phase points.

where zi is the mole fraction of the ith component in the overall fluid, z, and β is the vapor phase mole fraction defined from the material balance, z = βy + (1−β)x. Conservation of mass requires that

∑ yi = ∑ xi = 1 ⇒ ∑ (yi − xi) = 0 i

i

i

Initial Two-Phase Line. The algorithm starts by identifying a low pressure two-phase dew point, noted point A in Figure 1. To locate this point, a pressure of P0 = 1 bar is specified (i.e., S = ln P0, and is = N + 2). Solving eqs 1, 3, and 4 (equivalent to eq 5) will result in a high temperature dew point, where an asphaltene-rich phase is incipient from the vapor phase. The initial estimate for starting the solution process is found from the Wilson K-factor approximation.17 Following this, successive substitution18 (SS) and Newton−Raphson18 (NR) iterations are used to obtain tight convergence. To track the dew point line, one option could be to repeatedly specify an appropriately small step in either temperature or pressure and retrace the steps mentioned above. It turns out, however, that, even knowing and using the previous u-vector (eq 6) as an initial estimate for the next point, such an approach is infeasible at best and impossible at worst. This is because converging discrete points of the phase diagram is troublesome in the region around a critical point.8,16 Instead, it is advisible to design an extrapolation scheme where few known points on the phase line are used to extrapolate to the next point.15,16 A scheme, for the extrapolated solution vector, ue, can be a polynomial expansion

(3)

A specification equation uiS − S = 0

(4)

is added where S is the specified value (temperature, pressure, etc.) and iS is the index of the specified value type in the vector of independent variables, u (see below). Equations 1, 3, and 4 complete a set of N + 2 equations, conveniently written as g (u ) = 0

(5)

to be solved for the vector u = (u1 , u 2 , ..., uN + 2) = (ln K1, ln K 2 , ..., ln KN , α , ln P) (6)

of N + 2 unknown variables, where P is the absolute pressure. If the independent variable α is set to be ln T, where T is absolute temperature, u holds all information about a two-phase B


Article

Industrial & Engineering Chemistry Research J

uie(S) =

picked to be a potential candidate for dominating a separate phase. The component is also a potential candidate if the dominating component in W1 is identical to the one resulting in the largest Ki = φi(z)/φi(W0) factor. This simple screening typically returns between three and six trial phases. In these trial phases, the dominating components are typically either the lightest gaseous components or the heaviest components including asphaltene components. The strategy for selecting trial phases is the same for two-phase points as for three-phase points. If full stability analysis using the candidate trial phases reveals the presence of a third stable type of phase y, at (T0, P0), then the extrapolation has overshot a three-phase point, and a search for the exact three-phase PT point, (T3, P3), is initiated. As indicated in Figure 1, (T0, P0) is not located on an equilibrium line, but on a metastable continuation of the twophase line. The search for the exact three-phase PT point (T3, P3) uses the (T0, P0) point, the phase from stability analysis, y, and the phases investigated for instability, x and z. Particular to a three-phase point on a two-phase line is that two out of the three phases are incipient, and consequently, the nonincipient phase is the same as the overall fluid, z. When any phase, x, is in equilibrium with z, the reduced tangent plane distance15 becomes 0, i.e.,

∑ CijΔS j (7)

j=0

where ΔS is the extrapolation step.15,16 By experience we have found it practical and robust to use J = 3 (i.e., a third degree polynomial) for extrapolating to the next point. This leads to high quality estimates, only requiring few NR iterations for tight convergence. The NR iteration count used for convergence of the current point, is used to adjust the step-size, ΔS, used for extrapolating to the next point. The step-size adjustment from point n to n + 1 is performed so that ΔSn+1 = m·ΔSn, where m is the multiplier given in Table 1 below. The value of the multiplier is determined from practical experience. Table 1. Step-Size Multiplier for Different NR Iteration Counts NR iterations for convergence multiplier, m (two-phases) multiplier, m (three-phase)

1

2

3

4

5

6 or more

2 2

1.5 1.5

1.25 1.2

1 1.2

0.5 1

0.5 0.65

Upon closing in on a critical or a semicritical point, the number of NR iterations increases because the problem becomes increasingly ill-conditioned.8 A semicritical point is a point where three phases exist out of which two of the phases are mutually critical, and as a consequence, identical in composition. Furthermore, estimates from extrapolation in the critical or semicritical regions with high-order polynomials may lead to substantial error. For these two reasons, we have found it practical, when crossing a critical or semicritical point, to reduce the order of the polynomial expansion from J = 3 to J = 1. This need for change in polynomial expansion is not present when PT phase diagrams for mixtures of hydrocarbons and water are constructed. To find the coefficient matrix Cij in the case of J = 3, the previous (marked with (−1)) and the current (marked with (0)) point and corresponding derivatives, i.e., the phase line tangent (or sensitivity) vectors, u(−1)

∂u(−1) ∂S

u(0)

∂u(0) ∂S

N

tpdx =

∑ xi(ln xi + ln φi(x) − ln zi − ln φi(z)) = 0 i

(9)

The x and y phases are assumed to represent the incipient phases well, which is true if the distance from (T3, P3) to (T0, P0) is not too large. It is now possible to define the approximating polynomial of degree 1 for tpdx, tpdx(T ,P) ≈ P1x(T ,P) = tpdx(T0 ,P0) + ΔT

∂tpdx ∂T

+ ΔP P ,x

∂tpdx ∂P

T ,x

(10)

where ΔT = T − T0 and ΔP = P − P0 is the difference between a temperature and pressure point in the vicinity of (T0, P0) and (T0, P0) itself. An analogous approximating polynomial, Py1(T,P), can be written using the other incipient phase. To get an approximation for the three-phase point location (T3, P3), one can now solve Px1(T3, P3) = Py1(T3, P3) = 0. This leads to a 2 × 2 system of equations ÄÅ ÉÑ ÅÅÅ ∂tpdx ÑÑÑ ∂tpdx ÅÅ ÑÑ É ÄÅ ÅÅ ∂T ∂P T , x ÑÑÑÑ ÅÅ−tpdx ÑÑÑ ÅÅ P , x ÑÑ ÅÅ ÅÅ ÑÑ ÑÑ ÑÑ Aγ = b A = ÅÅÅ b = ÅÅÅ ÑÑ ÑÑ Å ∂tpdy tpd − ÅÅ ∂tpdy y ÑÑÖ Å Ñ ÅÇ ÅÅ ÑÑ ÅÅ ÑÑ ∂P ÅÅ ∂T ÑÑ P ,y T ,y Ñ ÅÅÇ ÑÖ

(8)

are needed. After each extrapolation to a new point, SS and NR iterations are applied to obtain tight convergence, as was the case for the initial point. Furthermore, upon passage of critical or semicritical points or other points of interest in the diagram (see section Critical and Semicritical Point Handling), the same scheme can be used for interpolation to locate such points. Stability Analysis. While generating points on the high temperature dew point branch in the direction from point A to B in Figure 1 by using extrapolation and subsequent tightening of convergence, we performed a stability analysis at each calculated point to detect whether it is energetically favorable for a third type of phase to form. Stability is investigated using modified tangent plane distance analysis.15,19 For this purpose, trial phases composed of molar amounts, W, must be defined. To obtain such trial phases, we make a simple screening of artificial mixtures, one for each component, being 95 parts pure, with the remaining 5 parts evenly distributed over the remaining components. A single SS step with ln W1i = ln φi(z) + ln zi − ln φi(W0), will do a simple update of the composition of the trial phase, W0 → W1. If the dominating component in W0 is still dominating the updated trial phase, W1, the component is

(11)

for the temperature and pressure differences, γ = [ΔTΔP], between the three-phase point (T3, P3) and (T0, P0). A stability analysis at the approximate (T3, P3) is performed on the known equilibrium phase, z, using x and y as trial phases, to update the incipient compositions x and y at (T3, P3). This procedure is repeated in an SS iterative manner, until convergence is obtained. In Figure 1 the first three-phase point found in this manner is marked as point B and as (T3, P3). Three-Phase Lines. Introduction of a third phase of mole fraction composition w requires the equation system (5) be T

C


Article

Industrial & Engineering Chemistry Research expanded to accommodate the extra phase. The equilibrium Kfactors are now written for two pairs of phases as y x K i1 = i K i2 = i wi wi (12) Because w is incipient (there will always be at least one incipient phase on a phase boundary), it still holds that z = βy + (1 − β)x. With the definitions above zi xi = wK wi = i i1 , yi = wK i i2 βK i2 + (1 − β)K i1 (13) Equation 3 is equally valid in the three-phase case, and now an additional mass balance equation is added

∑ wi − 1 = 0

(14)

i

Figure 2. Generic asphaltene PT phase diagram with no three-phase point on the two-phase line.

For a three-phase system, the two-phase equilibrium relations in eq 1 are expanded to ln K i1 + ln φi(x) − ln φi(w) = 0 ln K i2 + ln φi(y) − ln φi(w) = 0

and (15)

a bubble-point line. At the bubble-point, asphaltene precipitation peaks and the remaining fluid has a low asphaltene content. Thus, temporarily removing all asphaltene components from the mixture will allow a two-phase saturation pressure calculation at a low temperature, T3, marked with A in Figure 2. This intermediate calculation will give a reasonable initial estimate for the composition of the incipient vapor phase, y, and a good estimate of the pressure P3 on the three-phase line for the temperature T3. The temporary, asphaltene-free, mixture is subsequently discarded and a temporary pure asphaltene phase is defined instead, and used as a trial phase for inspecting the stability of the overall fluid z at T3, P3. An asphaltene-rich phase, w, is known to exist for pressures below the upper asphaltene onset pressure (AOP) (point B in Figure 2). Analyzing the stability of the overall fluid will give a good estimate of this asphaltene-rich phase. Again, the stability analysis is made using the modified tangent plane distance. Finally, the Rachford-Rice (RR) equation is solved20 to find an estimate for the composition of a liquid phase x ≠ z in two-phase equilibrium with the asphaltene phase, w. At this point, an estimate has been found for the liquid phase, x, the vapor phase, y, the asphaltene phase w, and the threephase pressure P3 at the temperature T3. All this information can now be used to start the solution process of converged points on the three-phase line, and tracking of the three-phase line using extrapolation as described above. Critical and Semicritical Point Handling. While traversing the phase lines, we find several points on the diagram are of special interest. These could be critical points on the two-phase line, semicritical points on the three-phase lines where two of three phases are identical, the cricondenbar, cricondentherm, or specified tabulation points for temperatures or pressures of specific interest. All of the mentioned points are conveniently determined by using the polynomial expansion in eq 7, but for interpolating rather than extrapolating. Particularly for the critical and semicritical points, converging to a solution using SS or NR iterations is difficult,8,16 and interpolation must be preferred. Using eq 7 for interpolation means that the current and previous points, needed to define the coefficient matrix, Cij, must be on either side of the particular point of interest, and it will be necessary to supervise when the calculation passes any such point. In the case of a critical or semicritical point, this

Equations 3, 4, 14, and 15 constitute the three-phase analogue of eq 5. The vector of independent variables u from eq 5 must be expanded to accommodate the extra set of equilibrium K-factors, and one-phase fraction, giving u → r = (ln K11, ln K 21, ..., ln KN1, α , ln P , β , ln K12 , ln K 22 , ..., ln KN 2)

(16)

The solution process is started using the phases found at the three-phase point. The overall solution strategy remains the same as for the two-phase case, but as the equation system for points on the three-phase line is larger than in the two-phase case, it is more complex to find the Jacobian and solve it using the NR method. Similarly, the extrapolation scheme of eq 7 remains unaltered but must be restarted, because no two-phase point can aid the extrapolation to a three-phase point. As for the two-phase case, the step-size for extrapolation may be adjusted during calculation, and the multiplier is found from Table 1. In the three-phase case we have found it practical to further limit step-size selection to a maximum value of ΔS = 0.2, if the specified variable is temperature, and ΔS = 0.6 if it is pressure, to prevent too large steps in pressure and temperature. In Figure 1, the three-phase line starting at point B returns to point B. The next line to construct is a two-phase line between points B and C in Figure 1. The two-phase line from point A to B had an asphaltene-rich phase as the incipient phase. On the twophase line from point B to C the incipient phase is a nonasphaltene phase. This nonasphaltene incipient phase was found when point B was located, and this starts the calculation of the two-phase line from B to C. The two-phase line is constructed in the same way as the line from A to B, and it ends at three-phase point C, from which a new set of three-phase and two-phase lines are constructed as described above. No Three-Phase Points on the Initial Two-Phase Line. An asphaltene PT phase diagram may not have three-phase instability at any point on the two-phase line, in which case the equivalent of points B and C in Figure 1 does not exist. Three phases may still be present at some P and T. Figure 2 shows an asphaltene PT phase diagram representing such case. In this case, no a-priori knowledge on how to locate the inner isolated three-phase (dashed) line is available at any temperature. The low temperature part of the dashed line in Figure 2 is D


Article

Industrial & Engineering Chemistry Research

In both cases, the calculation is started by performing a twophase saturation pressure calculation with α = 0. This is to establish a would-be saturation pressure, Psat, assuming no asphaltene precipitation occurs, between an incipient vapor phase y and a liquid phase x = z = z̃. The calculation is initiated from the Wilson K-factor approximation,17 and SS iterations18 subsequently tighten convergence. Once this would-be saturation pressure is found, the nonincipient phase is investigated for instability using a pure asphaltene phase as a trial phase to establish whether an asphaltene-rich phase is present. If the stability analysis confirms the presence of an asphaltene-rich phase, then the situation corresponds to the second case mentioned above, else to the first case. Case 1: When no asphaltene precipitation occurs at α = 0, the asphaltene Px phase diagram can appear as in Figure 3.

supervision is readily managed by inspecting whether the equilibrium ln Ki-factors in the subvectors of u and r ũ = (ln K1, ln K 2 , ..., ln KN ) rj̃ = (ln K1j , ln K 2j , ..., ln KNj)

j = 1, 2

(17)

change sign reflecting a change in incipient phase. For illustration purposes, consider passing a true critical point (0) (0) on a two-phase line. In this case u(0) = (ln K(0) 1 , ln K2 , ..., ln KN , (−1) are known, and it is also known that they are α, ln P) and u located on either side of ũ = 0. Interpolating to the critical point ũ = 0, is a matter of selecting the specification variable S to be the lth ln K-factor, taking ΔS = −ln K(0) l , and evaluating eq 7. Figure 2 shows two semicritical points, one at approximately 397 °C and 211 bar, and one at approximately 449 °C and 1.91 bar. At the point where the saturation pressure line (dashed line) meets the lower AOP line (blue line), the incipient asphaltene phase composition on the lower AOP line and the incipient hydrocarbon phase composition on the saturation pressure line are identical and form a semicritical point. This cusp requires special care to pass. As the algorithm approaches a semicritical point, the step-size, ΔS, between consecutive points decreases because it is adjusted according to the number of NR iterations required for convergence, as discussed earlier (Table 1). The more iterations, the smaller the step. Furthermore, in this semicritical region, the step-size is also bounded by a limit decided by the largest ln Kij-factor in the jth phase-pair closest to being critical. Finally, as mentioned earlier, it has been found most robust to pass the semicritical points by extrapolation with a J = 1 polynomial (eq 7) only relying on the current point (0) rather than a J = 3 polynomial relying on both the current and previous points. After extrapolation, the point is converged using NR iterations, and the calculation of the remaining phase line can continue in the usual manner. The location of the semicritical point is calculated using J = 3 interpolation between the two converged points on either side of the semicritical point. A sign change of ∂ ln T/∂S or ∂ ln P/∂S reveals that a cricondentherm or cricondenbar, respectively, has been passed. The point can be found by solving either of the second degree polynomials ∂ueN+1/∂S = 0 or ∂ueN+2/∂S = 0, respectively, for ΔS, and interpolating appropriately once ΔS is found.

Figure 3. Generic asphaltene Px phase diagram with no asphaltene precipitation for the original fluid.

The calculation traces the two-phase line (dashed line in Figure 3) from α = 0 and on, as described in the section Initial Two-Phase Line. Tracking continues until stability analysis reveals the extrapolation scheme in eq 7 has overshot a threephase point (α3, P3). An additional result of this stability analysis is an estimate of the asphaltene phase composition, w, separating from z̃ at the overshot point (α0, P0). Apart from the estimate of the composition of the asphaltene phase, the stability analysis gives the minimum in the tangent plane distance to z̃. This distance is the same for all components, so the difference in fugacities between the asphaltene phase and nonasphaltene liquid phase x = z̃ is15

■

ASPHALTENE Px PHASE DIAGRAM ALGORITHM The above algorithm can be modified to handle asphaltene precipitation in hydrocarbon fluids subject to gas injection under isothermal conditions. If the independent variable α in eq 6 or 16 is the phase fraction of a known injection gas with mole fraction composition, zinjection, then the composition of the mixture of the original and the injected fluid is given by z̃ = α· zinjection + (1 − α)·z. The technical details in terms of converging and extrapolating points when α is the phase fraction of an injection fluid, remain similar to when α is ln T. The Jacobians associated with eq 5 and the three-phase equivalent are not identical though. However, because the technical details are similar, the details of setting up the Jacobians are not covered here, only the strategy and approach to solve the problem. There can be two possible appearances of the Px diagram. In the first case, the original fluid does not show asphaltene precipitation, i.e., for α = 0 the asphaltene Px diagram only has the saturation pressure. After some addition of injection gas (α increasing from 0), asphaltene precipitation starts. Alternatively, the original fluid shows asphaltene precipitation in which case the asphaltene precipitation band is open at α = 0 and onward.

Δ = ln fi ̂ (w) − ln fi ̂ (z)̃ = −ln ∑ Wi = −ln WT i

(18)

where W is the absolute (mole number) composition of the asphaltene phase. The fugacity difference, Δ, can be incorporated into the three-phase equilibrium equation system, by modifying eq 15 to get ln K i1 + ln φi(z)̃ − ln φi(w) + Δ = 0 ln K i2 + ln φi(z)̃ − ln φi(y) = 0

and (19)

The corresponding Δ between the overall main phase z̃ and its equilibrium vapor phase, y, is zero from the would-be saturation pressure calculation and it will remain zero as expressed in the last part of eq 19. The formulation of the equations in our implementation further ensures that the vapor phase is kept incipient, i.e., the vapor phase fraction is zero. The asphaltene E


Article


now by increasing the asphaltene phase fraction. Once found, all phases are in equilibrium and the vapor phase is incipient, but the asphaltene phase has a nonzero phase fraction. From here, curve tracking for increasing α establishes the remaining part of the saturation pressure line. Essentially, the procedure is repeated to find the lower AOP line. However, here the asphaltene phase must be kept incipient instead of the vapor phase. It has turned out to be more practical to search for a point Δ = 0 where Δ is the fugacity difference between both asphaltene and liquid phase and between asphaltene and vapor phase, rather than for just a single phase pair. This means both sets of K-factors in eq 19 are modified with a Δ in this case. From here, curve tracking for increasing α establishes the remaining part of the lower AOP line.

phase fraction is also kept constant at zero to ensure the asphaltene phase is the second incipient phase at the three-phase point. Δ is substituted for the asphaltene phase fraction in the vector of independent variables r, and the curve tracking scheme can now locate the three-phase point (α3, P3) where Δ = 0. Once the three-phase point at α3 and the associated incipient phases are found, the algorithm is ready to calculate the remaining line sections. Using the overall fluid composition z̃(α3) and the incipient asphaltene phase composition w(α3) as the starting point, a two-phase line tracking starts from the threephase point to calculate the upper AOP line (solid black line in Figure 3). The saturation pressure line (dashed black line in Figure 3) and the lower AOP line (solid blue line in Figure 3) are both three-phase equilibrium lines for α ≥ α3. As is the case of the asphaltene PT phase diagram, a vapor phase is incipient along the saturation pressure line, while the asphaltene phase is incipient along the lower AOP line. In both cases, all information needed to start the tracking of the lines is available once the three-phase point has been found. Case 2: When the stability analysis finds that an asphaltene phase w is present at α = 0, the asphaltene Px phase diagram can appear as in Figure 4.

■

VALIDATION OF ASPHALTENE PHASE DIAGRAM ALGORITHM The location of the lines generated by the asphaltene PT phase diagram algorithm has been validated by carrying out multiphase PT flash calculations18,19 in a PT grid covering the same temperature and pressure range as the PT phase diagram. The phases present have been mapped in each PT point. Three validation examples are presented. Example 1 uses crude A from ref 21. Three upper AOPs and four saturation pressures have been measured and a PC-SAFT equation of state model has been developed to match the data in the mentioned reference. Figure 5 shows a mapping of the phases generated

Figure 4. Generic asphaltene Px phase diagram with asphaltene precipitation for the original fluid.

In this case, there is no three-phase point with two incipient phases for any positive value of α, but the stability analysis at α = 0 provides an estimate of the composition of the asphaltene phase and its fugacity difference, Δ, to the overall fluid z̃ = z. The algorithm starts by detecting the upper AOP at α = 0. The upper AOP line is a two-phase equilibrium line, where an asphaltene-rich phase is incipient and in equilibrium with the overall fluid z̃. The asphaltene phase found using stability analysis can be tracked to a point where Δ = 0 by modifying the two-phase equi-fugacity equation as was done to get from eq 15 to eq 19, while keeping the asphaltene phase incipient. Once this point is found, the incipient w is in equilbirum with z̃ and the upper AOP is established at α = 0. From here, curve tracking for increasing α establishes the remaining part of the line. To establish the true pressure at the saturation pressure line, the would-be saturation pressure is a good starting point, along with the would-be vapor and liquid phases, y and z̃ = z, and the estimated asphaltene phase, w. By increasing the asphaltene phase fraction, and maintaining a vapor fraction β of the value 0, we again use curve tracking to search for a point where Δ between the asphaltene phase and the liquid phase is zero, but

Figure 5. Phase mapping and lines from asphaltene PT phase diagram algorithm for example 1.

using multiphase PT flash calculations overlaid with lines (white dashed) generated by the asphaltene PT phase diagram algorithm. The flash and phase diagram results are fully consistent. Figure 6 shows the asphaltene PT phase diagram for the same fluid together with the data from ref 21. Because the applied PCSAFT model was developed to match the upper AOP and saturation pressure data, the simulation results should agree closely with the experimental data and that is seen to be the case. Example 2 uses the tuned CPA equation model from ref 22 developed for fluid C1 from ref 23 and the experimental upper AOP and saturation pressure data for the same fluid. Figure 7 shows mapping of the phases using multiphase PT flash F


Article


Figure 6. Asphaltene PT phase diagram and experimental data for example 1.




calculations overlaid with the lines generated by the asphaltene PT phase diagram algorithm. The results are fully consistent. Figure 8 shows the asphaltene PT phase diagram for the same fluid together with the data from ref 23. The applied CPA model was developed to match the upper AOP and saturation pressure data and the simulation results in Figure 8 are seen to agree closely with the experimental data. Example 3 is an SRK equation model developed in this work for fluid 1 from ref 24 tuned to the experimental upper AOP and saturation pressure data for the same fluid. The characterized fluid compositions and model parameters are provided in Tables S1 and S2 of the Supporting Information. Figure 9 shows mapping of the phases using multiphase PT flash calculations overlaid with the lines generated by the asphaltene PT phase diagram algorithm. The results are fully consistent. Figure 10 shows the asphaltene PT phase diagram for the same fluid together with the data from ref 24. The applied SRK model was developed to match the upper AOP and saturation pressure data, and the simulation results in Figure 10 agree closely with the experimental data. At temperatures of several hundred °C, asphaltene and other heavy molecules will not be stable.25 A two-phase PT phase envelope for a reservoir oil, which considers gas−oil phase


equilibria, may therefore not reflect the real phase behavior at high temperatures. Some of the heavier molecules may have broken down and that will not be reflected by the phase G


Article


(2) Haskett, C. E.; Tartera, M. A Practical Solution to the Problem of Asphaltene Deposits-Hassi Messaoud Field, Algeria. JPT, J. Pet. Technol. 1965, 17 (4), 387−391. (3) Kabir, C. S.; Jamaluddin, A. K. M. Asphaltene Characterization and Mitigation in South Kuwait’s Marrat Reservoir. Paper SPE 53155 presented at 1999 SPE Middle East Oil Show, Bahrain, 20−23 February 1999; SPE, 1999. (4) von Albrecht, C.; Diaz, B.; Salathiel, W. M.; Nierode, D. E. Stimulation of Asphaltic Deep Wells and Shallow Wells in Lake Maracaibo, Venezuela. Advances in Methods of Increasing Well Productivity or Injectivity 1979, No. 1, 55−62. (5) Bouts, M. N.; Wiersma, R. J.; Muijs, H. M.; Samuel, A. J. An Evaluation of New Asphaltene Inhibitors; Laboratory Study and Field Testing. JPT, J. Pet. Technol. 1995, 47 (9), 782−787. (6) Nagel, R. G.; Hunter, B. E.; Peggs, J. K.; Fong, D. K.; Mazzocchi, E. Tertiary Application of a Hydrocarbon Miscible Flood: Rainbow Keg River ″B″ Pool. SPE Reservoir Eng. 1990, 5 (3), 301−308. (7) Srivastava, R. K.; Huang, S. S.; Dong, M. Asphaltene deposition during CO2 flooding. SPE Prod. Facil. 1999, 14 (4), 235−245. (8) Lindeloff, N.; Michelsen, M. L. Phase Envelope Calculations for Hydrocarbon-Water Mixtures. Paper 77770 SPE Journal 2003, 8 (3), 298−303. (9) Rydahl, A.; Pedersen, K. S.; Hjermstad, H. P. Modeling of Live Oil Asphaltene Precipitation. AIChE Spring National Meeting March 9−13, Houston, TX, USA. 1997. (10) Ting, P. D.; Hirasaki, G. J.; Chapman, W. G. Modeling of asphaltene phase behavior with the SAFT equation of state. Pet. Sci. Technol. 2003, 21 (3−4), 647−661. (11) Zhang, X.; Pedrosa, N.; Moorwood, T. Modeling Asphaltene Phase Behavior: Comparison of Methods for Flow Assurance Studies. Energy Fuels 2012, 26 (5), 2611−2620. (12) Soave, G. Equilibrium Constants From a Modified RedlichKwong Equation of State. Chem. Eng. Sci. 1972, 27, 1197−1203. (13) Gross, J.; Sadowski, G. Perturbed chain SAFT: An equation of state based on perturbation theory for chain molecules. Ind. Eng. Chem. Res. 2001, 40, 1244−1260. (14) Kontogeorgis, G. M.; Voutsas, E. C.; Yakoumis, I. V.; Tassios, D. P. An equation of state for associating fluids. Ind. Eng. Chem. Res. 1996, 35 (11), 4310−4318. (15) Michelsen, M. L.; Moellerup, J. Thermodynamic Models: Fundamentals & Computational Aspects, 2nd ed.; Tie-Line Publications: Holte, 2007. (16) Michelsen, M. L. Calculation of phase envelopes and critical points for multicomponent mixtures. Fluid Phase Equilib. 1980, 4 (1), 1−10. (17) Wilson, G. M. A Modified Redlich-Kwong Equation of State, Application to General Physical Data Calculation. Paper No. 15C presented at the 1969 AIChE 65th National Meeting, Cleveland, Ohio; AIChE, 1969. (18) Michelsen, M. L. The isothermal flash problem. Part II. Phasesplit calculation. Fluid Phase Equilib. 1982, 9, 21−40. (19) Michelsen, M. L. The isothermal flash problem. Part I. Stability. Fluid Phase Equilib. 1982, 9, 1−19. (20) Rachford, H. H.; Rice, J. D. Procedure for use of electronic digital computers in calculating flash vaporization hydrocarbon equilibrium. JPT, J. Pet. Technol. 1952, 4, 327−328. (21) Panuganti, S. R.; Vargas, F. M.; Gonzalez, D. L.; Kurup, A. S.; Chapman, W. G. PC-SAFT characterization of crude oils and modeling of asphaltene phase behavior. Fuel 2012, 93, 658−669. (22) Arya, A.; von Solms, N.; Kontogeorgis, G. M. Determination of asphaltene onset conditions using the cubic plus association equation of state. Fluid Phase Equilib. 2015, 400, 8−19. (23) Buenrostro-Gonzales, E.; Lira-Galeana, C.; Gil-Villegas, A. Asphaltene Precipitation in Crude Oils: Theory and Experiments. AIChE J. 2004, 50 (10), 2552−2570. (24) Gonzalez, D. L.; Mahmoodaghdam, E.; Lim, F.; Joshi, N. Effect of Gas Additions to Deepwater Gulf of Mexico Reservoir Oil: Experimental Investigation of Asphaltene Precipitation and Deposition. Paper SPE 159098 presented at 2012 SPE Annual Technical

envelope. Generation of phase envelopes is still considered useful to get a complete overview of the phase behavior predicted by the applied model. In this respect the asphaltene phase envelope algorithm is no different. It will at all conditions reflect the applied model. Correctness of the entire diagram from a model point of view ignoring alterations of molecules has been demonstrated through Figures 5, 7, and 9. The published data for the examples and the published simulation results only cover a limited temperature span. Figures 6, 8, and 10 show the full asphaltene PT phase diagrams and reveal significant differences between the modeling results for the fluids outside the temperature range with experimental data. The low temperature branches in Figures 6 and 10 suggest that the asphaltene and oil phases are immiscible at temperatures below approximately 50 °C. Figure 8, however, suggests that oil and asphaltene become more miscible at temperatures below 50 °C. Only experimental data can tell whether the models are right, but the differences revealed by the full phase diagram illustrate the potential of the developed algorithm to analyze whether an asphaltene model is physically sound or just capable of matching data points in a limited temperature interval.

■

CONCLUSION An algorithm for constructing complete asphaltene PT phase diagrams covering large spans in pressure and temperature has been presented. With appropriate modifications it can also be used for constructing asphaltene Px phase diagrams. Examples given using the PC-SAFT, CPA, and SRK equations show that the algorithm can work with any equation of state, providing consistent fugacity coefficients and appropriate derivatives for vapor, liquid, and asphaltene phases. Examples are shown illustrating the multitude of possible appearances of an asphaltene PT phase diagram. Construction of the asphaltene PT phase diagram is fast and robust even in the near-critical and semicritical regions. The algorithm can also handle cases with more than one three-phase point, multiple semicritical points, and the cusps, which can occur where lines in the asphaltene PT phase diagram meet. The algorithm provides results consistent with multiphase PT flash calculations and agree well with asphaltene simulation results reported in the open literature.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.7b04246.

■

Characterized fluid composition and model parameters for fluid 1 from ref 24 (PDF)

AUTHOR INFORMATION

Corresponding Author

*H. Sørensen. E-mail: [email protected]. ORCID

Henrik Sørensen: 0000-0002-8213-0882 Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Hirschberg, A.; DeJong, L. N. J.; Schipper, B. A.; Meijer, J. G. Influence of temperature and pressure on asphaltene flocculation. SPEJ, Soc. Pet. Eng. J. 1984, 24 (3), 283−293. H


Article

Industrial & Engineering Chemistry Research Conference and Exhibition, San Antonio, Tx, USA, 8−10 October 2012; SPE, 2012. (25) Martínez, M. T.; Benito, A. M.; Callejas, M. A. Thermal cracking of coal residue: kinetics of asphaltene decomposition. Fuel 1997, 76 (9), 871−877.

I


Algorithm for Constructing Complete Asphaltene ... - ACS Publications

Recommend Documents