Comparison of Phase Identification Methods Used ... - ACS Publications

Article pubs.acs.org/EF

Comparison of Phase Identification Methods Used in Oil Industry Flow Simulations Jim Bennett* and Kurt A. G. Schmidt Abingdon Technology Center, Schlumberger, Lambourn Court, Wyndyke Furlong, Abingdon OX14 1UJ, United Kingdom ABSTRACT: The determination of whether a phase is a liquid or a vapor is somewhat arbitrary, particularly when the fluid is in the single-phase region. However, the calculation of multiphase flow in oil industry simulators is dependent upon the labeling of the phases. Key properties in flow modeling, such as relative permeabilities and capillary pressures, are dependent upon which phases are present. Correct labeling of the fluid phases is therefore critical to the success of simulations. Many methods have been proposed for identifying phases. These are reviewed here. Five methods have been selected for detailed comparison, including a new method, based on the thermal expansion coefficient. The results of the comparison and the relative advantages of the methods are discussed.

■

One important aspect of multiphase flow is controlled by the end point of the relative permeability, which is the saturation where the relative permeability is zero and the phase does not flow. Typically, for liquids this occurs when the phase saturation is non-zero, representing connate water or unrecoverable oil, while for gas, it is usually zero or close to zero, meaning that all (or nearly all) gas will flow. Therefore, if oil is mislabeled, some will flow when it should not, increasing the recovery. If gas is mislabeled, some would become immobile, which would help to maintain pressure in the reservoir and would reduce the gas production. In both cases, the whole pattern of flow in the reservoir would be altered, invalidating the simulation. Correct phase labeling is also required in pipe flow simulators, because different flow models can be used for oil and gas phases.5 Phase labeling also causes difficulties when simulating surface equipment. If a fluid labeled gas is passed to a pump model or a fluid labeled oil is passed to a compressor model, the code has to decide whether to ignore the phase label and assume the phase type is correct or to issue an error message saying that the user has selected the wrong equipment type. Therefore, all of these flow simulators require a reliable and consistent way of labeling phases. In reservoir simulations, the phase identification method must be fast. This rules out the rigorous calculation of the critical point6 or drawing of the phase envelope. Therefore, a variety of alternative methods have been used to determine the phase type, discussed in the next section. Some methods are based on calculating a pseudo-critical point; others use the volume, the isothermal compressibility k, or derivatives of k. In this work, a new method based on the thermal expansion coefficient α, obtained from an equation of state (EoS), is proposed and compared to some of the other methods.

INTRODUCTION The determination of whether a single phase is a liquid or a vapor is somewhat arbitrary. Differences in liquid- and gas-like behavior in the supercritical region have been linked to the Widom lines for pure fluids,1,2 although these lines define a region where the properties change rather than a clear division between liquid and gas. Pragmatically, any line can be drawn on a PT phase envelope that divides the single-phase region into two. For example, Pedersen et al.3 suggest a line through the critical point parallel to the pressure axis. The line is arbitrary and artificial; fluid property models (such as density and viscosity) do not depend upon the phase label; therefore, in the supercritical region, two points on opposite sides of the line but close to the line will have the same properties, even though they will have different phase labels. This contrasts with two points on either side of the saturation line of a pure fluid, where there is a step change in properties across the line. The arbitrariness of the liquid−vapor line in the supercritical region could be overlooked if it has no effect. However, in flow simulation, the phase label turns out to be very important. In the petroleum industry, numerical simulators are used to model multiphase flow in hydrocarbon reservoirs, pipelines, and networks. In reservoir simulation, the flow of a hydrocarbon phase in the presence of a water phase depends upon the relative permeability kr of the hydrocarbon phase and the capillary pressure Pcw between the two phases.4 A simple model of the relative permeability would make it a function of the proportion of the phases (represented by the water saturation Sw), the surface tension between the phases σ, and the wettability of the rock w: k r = k r(Sw , σ , w)

(1)

However, rather than attempting to model this functionality explicitly, it is standard practice4 to use relative permeability functions that depend only upon the proportion of the phase and to use different relative permeability functions for the oil phase kro(Sw) and the gas phase krg(Sw). The same applies to the capillary pressure, which is specified by different functions for the oil phase Pcow(Sw) and the gas phase Pcgw(Sw). © XXXX American Chemical Society

Special Issue: 17th International Conference on Petroleum Phase Behavior and Fouling Received: September 13, 2016 Revised: November 11, 2016

A

DOI: 10.1021/acs.energyfuels.6b02316 Energy Fuels XXXX, XXX, XXX−XXX

Article

Energy & Fuels

Figure 1. PT plots of (a) log(V), (b) log(k), and (c) log(α) for pure methane. Low values are shown in blue, and high values are shown in red, although the scale is different in each plot. Pressure is in bar, and temperature is in kelvin. The saturation line is shown in black.

Compositional fluids used in simulation are usually modeled using either the Peng−Robinson7 (PR) or Soave−Redlich− Kwong8 (SRK) cubic EoS with the van der Waals mixing rule. In reservoir simulators, a two-phase flash is typically used to determine the distribution of the hydrocarbon components within the non-aqueous phases, while the aqueous phase contains only aqueous components, thermodynamically isolated from the hydrocarbons.4 In pipeline and network simulation, three-phase liquid−vapor−aqueous flash calculations are used, with either a cubic EoS combined with the Huron−Vidal mixing rule or a more advanced EoS, such as cubic plus association (CPA). Multiphase flashes are also often used as a screening tool, particularly to determine when a solid phase might appear. In this investigation, five methods of labeling phases have been selected and compared by performing flashes on a variety of reservoir fluids ranging from natural gases to heavy oils. The fluids consist of one-, two-, or three-phase systems as determined by a multiphase flash calculation.

Nc

Tc = r1

wj = xjVcj

⎛ b ⎞ Vc = Zc⎜ ⎟ ⎝ Ωb 0 ⎠

(4)

(5)

The phase is labeled vapor for temperatures above the estimated critical temperature. The value of r1 can be used to tune the phase identification, so that it gives the correct results at known conditions. The Multiflash package16 labels a phase gas if VT2 > VcTc2, where V c and Tc are the pseudo-critical volume and temperature, respectively. The exact formulas used for the phase critical properties are not specified; therefore, this method has not been tested here. However, in practice, this method works well, with phase identification comparable to the best results found in this study. Volume Methods. The condition V > Vc together with eq 3 can be written as:

REVIEW OF PHASE IDENTIFICATION TECHNIQUES Phase identification can be used in two ways during the flash calculation. First, it can be used to determine the root to use9 in a cubic EoS, which determines whether the phase has a liquid or vapor density. Second, it can be used at the end of the flash to identify the final state of the fluid.10 Some of the methods in the literature are presented below. Critical Point Methods. Gosset et al.11 label a phase as vapor if either the temperature exceeds the phase critical temperature, T > Tc, or if the phase volume exceeds the phase critical volume, V > Vc. For the phase critical properties, they followAsselineau et al12: ⎛ 1 ⎞⎛ a ⎞⎛ Ω b 0 ⎞ ⎜ ⎟⎜ ⎟⎜ ⎟ ⎝ R ⎠⎝ b ⎠⎝ Ωa0 ⎠

Nc

∑ j = 1 wj

The weights wj are calculated from the component mole fraction and component critical volume:

■

Tc =

∑ j = 1 wT j cj

V > r2 b

(6)

where r2 is a constant. This can be regarded as a criteria based on volume rather than the critical point. The b coefficient is the hard sphere volume in the cubic EoS, i.e., the minimum molar volume of the phase that occurs at infinite pressure. Pedersen et al.3 recommend r2 = 1.75 to determine the phase label, although it can also be used to tune the method. Figure 1a shows a contour plot of log(V) for pure methane. It can be seen that there is a jump in log(V) at the saturation line, but the size of the discontinuity becomes smaller near the critical point. This confirms that, away from the critical point, using a relative volume is a sensible way of distinguishing phases. Negative Flash Method. Perschke17 described a technique based on the negative flash18 with constant K values. The twophase Rachford−Rice equation19 can be solved for the vapor fraction β, given a feed composition {xi} and K values for each component:

(2)

(3)

The terms in these equations are defined in the Nomenclature. A constant value is used for the critical Z factor, Zc, depending upon the EoS. This method is used in the reservoir simulator IPARS.13 The ECLIPSE compositional reservoir simulator14 uses the simple condition T > Tc. The phase critical temperature is estimated with a mixing rule from Li.15 This mixing rule uses the weighted average of the component critical temperatures:

⎫ ⎪ ⎬ ⎪ ⎪ 1 ( 1) + − K β ⎭ j j=1 ⎩ Nc

G (β ) =

⎧ (Kj − 1)xj ⎪

∑⎨

(7)

In the traditional flash technique, eq 7 will not have a solution for a stable phase in the range β ∈ (0, 1).20 However, in the B


Article

Energy & Fuels negative flash technique, eq 7 can be used to identify the phase. In this approach, if eq 7 is solved and the solution results in a value of β > 1, the phase is vapor, and if eq 7 is solved and the solution results in a value of β < 0, the phase is liquid. In fact eq 7 does not need to be solved, because it is a monotonically decreasing function. Equation 7 can simply be evaluated at β = 0.5, and if G(0.5) > 0, the solution is β > 0.5 and the phase can be identified as vapor. Although K values are only defined in the two-phase region, the Wilson K value correlation,21 which is defined everywhere in PT space, can be used to obtain a solution in the single- and two-phase regions. Isothermal Compressibility Methods. The isothermal compressibility k is given by: ⎛ ∂ log V ⎞ 1 ⎛ ∂V ⎞ ⎟ =− ⎜ ⎟ k = −⎜ ⎝ ∂P ⎠T V ⎝ ∂P ⎠T

determining the phase. Nevertheless, results from this investigation suggest that it has merits for reservoir fluids. It does however have the disadvantage that the second temperature derivative of volume (∂2V/∂T2) is needed, which requires some effort and care to calculate. A Widom line can be defined as a line where a given property reaches an extrema or inflection point.2 Therefore, the lines (∂k∂T)P = 0 and (∂α/∂T)P = 0, which demark the liquid and vapor regions in the isothermal compressibility method and the thermal expansion coefficient method are both Widom lines. The isothermal compressibility and the thermal expansion coefficient are also related by: ⎛ ∂k ⎞ ⎛ ∂α ⎞ ⎟ ⎜ ⎟ = −⎜ ⎝ ∂P ⎠T ⎝ ∂T ⎠P

(8)

Therefore, the method of Pasad and Venkatarathnam,22 which is defined as measuring the change in k in Figure 1b in the horizontal T direction, is equivalent to measuring the change in α in Figure 1c in the vertical P direction. The proposed method described in this section is equivalent to measuring the changes in α in Figure 1c in the horizontal T direction. Other Methods. This study focuses on the correct labeling of the vapor phase, and no attempt is made to distinguish between multiple liquid phases. However, this will be important when simulating flow of more than two hydrocarbon phases. Perschke et al.26 used the fraction of the heaviest hydrocarbon to identify an oil phase and the density of the other phases to identify a second liquid phase and a vapor phase. The phases were then tracked during simulation. Recently, Xu and Okuno27 described an improvement on the method of Perschke et al.26 If three phases are present, their method is similar but they do not use the phase densities; instead, they identify the vapor phase as the phase with the highest methane concentration. However, when two phases are present, they use a very different and interesting approach based on critical end points and tie triangles. The method has not been tested here because it is more complicated to implement than the other selected methods. Use of key components (such as methane) to identify a phase is problematic when the component is not present in the fluid. This is unlikely to be the case in a reservoir simulation but could occur in other flow simulations, with modeling injection of non-reservoir fluids (e.g., a CO2 injection well).

Figure 1b shows contour plots of log(k) for pure methane calculated using the PR EoS. At low temperatures, the fluid is a liquid and k increases with the temperature. There is a discontinuity with a sharp jump in k as the temperature crosses the saturation line and the fluid becomes vapor. The compressibility decreases as the temperature increases further, although there is a stronger pressure dependence at high temperatures. At low pressures in the vapor region, k decreases relatively slowly with the temperature. Poling et al.9 use the fact that the isothermal compressibility is much smaller for liquids than for vapors. Their criteria use the dimensionless term kP, with kP < 0.005 for a liquid phase and 0.9 < kP < 3 for a vapor phase. Pasad and Venkatarathnam22 recommend use of the temperature derivative of the isothermal compressibility. For values of (∂k∂T)P > 0, the fluid is liquid (or a liquid-like vapor). In some circumstances, this derivative changes sign at high temperatures and k increases again, even though the fluid is still a vapor. Venkatarathnam and Oellrich10 therefore recommend checking whether the point (∂k∂T)P = 0 occurs at a temperature below the actual temperature. The temperature derivative of the isothermal compressibility has been recently used to identify solid phases.23 Thermal Expansion Coefficient Method. Williams and Reza24 looked at derivatives to determine the phase type when non-cubic EoSs were used. In their work, they found that (∂2V/ ∂T2)P was a good indicator of a vapor phase whenever the derivative was negative. This approach has been modified here to use the temperature derivative of the thermal expansion coefficient α, given by: ⎛ ∂ log V ⎞ 1 ⎛ ∂V ⎞ ⎟ = ⎜ ⎟ α=⎜ ⎝ ∂T ⎠P V ⎝ ∂T ⎠P

(10)

■

METHOD OF COMPARISON

For this study, phase identification methods were selected that could be easily added to multiphase flash codes. In particular, methods that could calculate a score Q depending upon the phase type were used to rank the phases. Thus, if a scheme identified more than one phase as a vapor, the phase with the highest Q could be chosen as the vapor phase. The same method was used for each phase, regardless of whether the system had one, two, or more phases present. Five different methods were selected for comparison: (a) the pseudo-critical temperature method,14 (b) the volume method of Pedersen et al.,3 (c) the negative flash method, (d) the k derivative method,22 and (e) the new α derivative method. These five methods are all fast; therefore, the time to calculate the phase label is negligible compared to the time taken to perform the flash. The first two methods are tunable but have been used here with the recommended values r1 = 1 and r2 = 1.75.3 These methods are summarized in Table 1, where a phase is deemed to be a vapor if Q < 0 and a liquid if Q > 0. The methods were tested on over a 100 reservoir fluids taken from the literature and real reservoir studies. The tests were performed by flashing each fluid on a 100 × 100 pressure−

(9)

This modification has two advantages: first, the derivative of the thermal expansion coefficient has more physical meaning than the volume derivative (∂2V/∂T2)P, and second, the modification gave better results for the fluids tested in this study. Figure 1c shows relative values of log(α) for pure methane. As seen, log(α) shows a ridge, with high values along the saturation line (the largest being at the critical point). The contours run roughly parallel to the saturation line at low and high temperatures. For values of ∂α/∂T > 0, the fluid is liquid. Troncoso et al.25 report that α can exhibit both positive and negative temperature dependence, depending upon the (pure) fluid, which suggests that its temperature derivative, ∂α/∂T, would not be a good choice in all circumstances for C


Article

Energy & Fuels

pressure line at low pressures but diverges from it at high pressures. The other three methods all follow the saturation line. As shown in Figure 2d, the isothermal compressibility method identifies vapors as liquids at high temperatures and pressures, which confirms the problems identified by Venkatarathnam and Oellrich.10 At high temperatures, α at first increases with pressure but then starts to decrease (Figure 1c), which would mean a change of phase if the isothermal compressibility method is used. The problem does not arise if the thermal expansion coefficient method is used (Figure 2e), because α decreases with the temperature at all pressures in the vapor region. The improvement recommended by Venkatarathnam and Oellrich10 for the isothermal compressibility method was not explored further because this method makes the isothermal compressibility method more complicated. For the fluids used in the remainder of this study, the correction was not expected to affect the phase identification, because the pressures and temperatures used for the study are outside the superheated region. Only a few pure component fluids are used in reservoir flow simulations. Water is often modeled as a single component, and pure methane and pure carbon dioxide may be injected into a reservoir. However, water is often modeled as an inert phase for simplicity, and reservoir temperatures are always above the critical temperature of methane. Therefore, it is only when pure CO2 is modeled in the injection well that phase identification of a pure component fluid needs to be considered.

Table 1. Phase Identification Methods Compared in This Study a b c d e

method

Q

pseudo-critical temperature volume negative flash isothermal compressibility derivative thermal expansion coefficient derivative

T − Tc r2 − V/b −G(0.5) (∂k∂T)P (∂α/∂T)P

temperature grid of cells. Therefore, a million flashes were performed for each of the five phase identification methods. The results were then plotted and compared visually.

■

COMPARISON OF RESULTS Five representative fluids were selected for illustrative purposes and subsequent discussions. In each case, the tests were re-ran with a 500 × 500 grid to show more detail. Each grid cell was colored according to the fluid state, and the following color scheme was used in all of the plots in this section: red, singlephase vapor (V); green, single-phase liquid (L); cyan, liquid− liquid (LL); yellow, liquid−vapor (LV); and black, three phase (both LLL and LLV). Pure Methane. PT phase envelopes are shown in Figure 2 for pure methane. The pseudo-critical temperature method (Figure 2a) stands out because it does not attempt to follow the saturation line. Instead, it labels the fluid depending upon whether the temperature is greater or less than the critical temperature. The volume method (Figure 2b) does follow the saturation

Figure 2. PT plots for pure methane, with pressure in bar and temperature in kelvin. The saturation line is shown in black. The phases are determined by the five methods: (a) pseudo-critical temperature, (b) volume, (c) negative flash, (d) isothermal compressibility derivative, and (e) thermal expansion coefficient derivative. D


Article

Energy & Fuels

Figure 3. PX plots for a binary fluid CO2 and CH4, with pressure in bar and the x axis representing the fraction of methane in the mixture. The phases are determined by the five methods: (a) pseudo-critical temperature, (b) volume, (c) negative flash, (d) isothermal compressibility derivative, and (e) thermal expansion coefficient derivative.

Figure 4. PX plots for a gas injection study from Li and Firoozabadi,28 with pressure in bar and the x axis representing the fraction of injected gas in the mixture. The phases are determined by the five methods: (a) pseudo-critical temperature, (b) volume, (c) negative flash, (d) isothermal compressibility derivative, and (e) thermal expansion coefficient derivative.

Binary Fluid. PX plots are shown at 180 K in Figure 3 for a binary CH4/CO2 mixture following Michelsen and Møllerup.20 In this figure, the x axis represents the fraction of methane in the feed. The figure shows a two-phase region split by a threephase line at a constant pressure.

The pseudo-critical temperature (Figure 3a) method identifies all phases as liquids. The other four methods do differentiate between liquid and vapor. They show the twophase region splits from LL to LV, a single-phase vapor at pressures below the two-phase region and a single-phase liquid above the two-phase region. However, the volume method E


Article

Energy & Fuels

Figure 5. PT envelopes for the synthetic fluid described in Table 2, with pressure in bar and temperature in kelvin. The true critical point is shown as +, and measured dew points are shown as ×. The phases are determined by the five methods: (a) pseudo-critical temperature, (b) volume, (c) negative flash, (d) isothermal compressibility derivative, and (e) thermal expansion coefficient derivative.

There is a single-phase liquid region at high pressures near x = 1. However, this phase is classified as a vapor in panels a and b of Figure 4. Synthetic Fluid. PT plots for a five-component synthetic fluid taken from Gozalpour et al.29 are shown in Figure 5. The fluid was modeled using the PR EoS.30 The composition and binary interaction coefficients are given in Table 2.

(Figure 3b) also shows vapor in the high-pressure single-phase region. For this binary mixture, the negative flash method (Figure 3c) and the two derivative methods (panels d and e of Figure 3) all give the same results, with natural transitions between the phases. These are consistent with the phase calculated by Michelsen and Mollerup.20 There is a three-phase line between the LL and LV regions, which continues all of the way across the two-phase region, splitting the LV region into two at high methane concentrations. Here, the three-phase line marks the transition between L1V and L2V states. As mentioned earlier, no attempt was made to distinguish between the liquid states in Figure 3. CO2 Injection Study. In this example, the North Ward Estes oil is taken from Li and Firoozabadi.28 A mixture of 95% CO2 and 5% CH4 was injected into the oil, and PX diagrams are shown in Figure 4. The x axis represents the mole fraction of injected gas in the feed. The thermal expansion coefficient method (Figure 4e) shows a good match between the phase identification and the threephase region. The three-phase region can be seen as the overlap of the LV region and the LL region. In addition, the LV−LL line, where the state changes without going through a threephase region, can be seen as a continuation of the three-phase region. This correspondence between the LV−LL line and the three-phase region is poor with the pseudo-critical temperature method (Figure 4a) and the volume method (Figure 4b), where the transition from LV to LL is predicted to occur at pressures greater than the pressure in the three-phase region. The correspondence between the line and the region is good with the negative flash method (Figure 4c) and the isothermal compressibility method (Figure 4d), although not quite as good as the thermal expansion coefficient method (Figure 4e), because the LV−LL lines do not quite meet the apex of the three-phase region.

Table 2. Composition of the Batch 2 Synthetic Fluid from Gozalpour et al.29 a

a

component

mole fraction

CH4 BICs

C3H8 BICs

CH4 C3H8 n-C5 n-C10 n-C16

0.8205 0.0895 0.0500 0.0199 0.0201

0.0140 0.0236 0.0501 0.0370

0.0100 0.0250 0.0250

Non-zero binary interaction coefficients (BICs) are shown.

None of the methods show a L−V line that is particularly close to the calculated critical point. The L−V line in the isothermal compressibility method (Figure 5d) stands out as being further away from the critical point than the others. This means that some (or all) of the measured dew points will be misidentified as bubble points. This could be corrected for the methods that allow for tuning, such as the pseudo-critical temperature method and the volume method. In Figure 5, the L−V lines continue into the two-phase regions, splitting them into LV and LL regions. This is a result of applying the phase identification methods to each phase. It is a sensible approach because it labels the low-temperature/highpressure two-phase region as liquid−liquid. However, it differs from the approach used in reservoir simulation, when a twoF


Article

Energy & Fuels phase flash is used. There, the two-phase region is always considered to be liquid−vapor. In three cases (panels c−e of Figure 5), the LL−LV line appears as an extension of the three-phase region. In these cases, the three-phase region can be thought of as the overlap between the LL and LV regions. However, in the pseudocritical temperature method (Figure 5a), there is a LL region below the three-phase region, and in the volume method (Figure 5b), there is a small LV region above the three-phase region. Reservoir Fluid. PT plots for a 10 component reservoir fluid (Table 3) are shown in Figure 6. The fluid is one of the

coefficient method (Figure 6e) all give similar phase identifications of this fluid. The fluid has two three-phase regions shown in black, one at low temperatures and the other at a higher temperature but only at a low pressure. This second three-phase LLV region corresponds with the boundary of the LL and LV regions in all methods, except the pseudo-critical temperature method (Figure 6a), where there is a LL region above and below the LLV region. On the other hand, the L−V line is closest to the critical point in the pseudo-critical temperature method (Figure 6a). Only the isothermal compressibility method (Figure 6d) shows continuity between the L−V line and the LL−LV line. However, this line is again not close to the critical point and is at a relatively low pressure compared to the other methods.

Table 3. Example Reservoir Fluid Composition Taken from the ECLIPSE Compositional Reservoir Simulator14 a component

mole fraction

molecular weight

specific gravity

N2 CO2 CH4 C2−C3 C4−C6 C7 C13 C18 C26 C43

0.0069 0.0314 0.5280 0.1515 0.0703 0.0867 0.0529 0.0340 0.0238 0.0145

35.88 67.98 110.14 173.11 248.85 361.77 600.98

0.9752 0.6236 0.7378 0.7679 0.8127 0.8194 0.8953

■

DISCUSSION Visual comparison of the phase identification results from each method has given insight into the strengths and weaknesses of the five phase identification methods. P−T grids used here are a powerful tool for the visualization of the phase behavior of the fluids tested and give a means of indicating where problems with the phase identification may occur. The range of the pressure and temperature in the P−T grids was originally chosen to focus on two-phase and three-phase regions, which is where problems in the flash will show up. To evaluate phase identification, it would be worth testing on a larger pressure and temperature area to check single-phase phase identification and ensure no issue exists with the overall flow simulations. When phases are identified, it would be desirable for the L− V line to pass through the critical point (determined by the rigorous method). Direct methods of calculating the critical point6 or drawing the phase envelope are slow and, therefore, cannot be used. Furthermore, not all (modeled) fluids have a

a

The oil density is 38.0° API gravity/835 kg/m3, and the gas/oil ratio (GOR) is 1896 scf/sbbl/338 sm3/sm3.

example fluids that come with the ECLIPSE compositional reservoir simulator.14 The volume method (Figure 6b), the negative flash method (Figure 6c), and the thermal expansion

Figure 6. PT envelopes for the reservoir fluid described in Table 3, with pressure in bar and temperature in kelvin. The true critical point is shown as +. The phases are determined by the five methods: (a) pseudo-critical temperature, (b) volume, (c) negative flash, (d) isothermal compressibility derivative, and (e) thermal expansion coefficient derivative. G


Article

Energy & Fuels

technically more difficult than any of the other methods tested here, because it requires second-order derivatives with respect to the temperature. However, it has performed well with the fluids in this study, giving a reasonable match between the L−V line and the critical point and the LV−LL line and the threephase region. It should be mentioned that only reservoir fluids, i.e., fluids containing hydrocarbons and gas components (CO2, H2S, and N2), have been tested in this study. No water or aqueous components were constituents of these fluids. Additionally, only cubic EoS with the van der Waals mixing rule have been studied. With these caveats kept in mind, the method warrants further investigation as a practical method of identifying phases.

critical point; thus, a method that accurately determines the critical point would still not be reliable for these fluids. Several of the methods do show a L−V line passing through the critical point for a single-component mixture, but none of the methods studied here does this consistently for multicomponent fluids. During the evaluation of the phase identification techniques, particular attention was focused on whether the LL−LV line connects with the three-phase region. All of the methods did this for some fluids, but none of the methods did this for all of the 100 fluids tested. The negative flash method and the two derivative methods performed better than the other two methods, with the thermal expansion coefficient method performing the most consistently. Satisfaction of this criterion could be useful because it would allow for a fairly efficient search for a three-phase region; a bisection method between a LL state and a LV state would either hit a three-phase region or converge on a transition line between the two states. Pseudo-critical Temperature Method. In reservoir flow simulation, the flow depends upon whether a phase is labeled oil or gas. The problem with this is that a single-phase may change suddenly from oil to gas if the temperature, pressure, or composition changes and may cause problems with the flow simulator results. The pseudo-critical temperature method has the advantage that the single-phase identification does not change as the pressure changes. In oil reservoirs, the reservoir temperature is fairly constant and the pseudo-critical temperature method can be easily tuned to ensure that the pseudocritical temperature is always less than (or greater than) the reservoir temperature, ensuring consistent phase identification throughout the reservoir. However, the phase can still change from pure vapor to pure liquid as a result of a change in composition, and therefore, flow models need to be robust enough to handle this situation. Because the pseudo-critical temperature method looks only at the critical temperature and ignores the critical pressure, it is unable to identify phase changes that are due to pressure, e.g., in the PX diagram (Figure 2a). Therefore, it is not a good method for a priori phase identification, although it can be successfully used within reservoir simulation. Volume Method. The volume method does not have the same drawbacks as the pseudo-critical temperature method, although it does have limitations; for example, it fails to match the full length of the saturation curve for a single component (Figure 1b). However, as with the pseudo-critical temperature method, it is also tunable. Negative Flash Method. The negative flash method gives better results than the other two non-derivative methods tested here and is simpler to implement than the derivative methods. Isothermal Compressibility Derivative Method. This method showed better consistency between the L−V line and the LL−LV line than any other method. Nevertheless, the method cannot be recommended. In comparison to the other methods, the L−V line was often at a relatively low pressure (e.g., Figures 5d and 6d) and met the two-phase region further from the critical point. In these circumstances, the method will indicate a liquid when vapor is a more realistic choice. In addition, Venkatarathnam and Oellrich10 have pointed out that the method will give the wrong phase at high pressures and temperatures (e.g., Figure 1d). Their improvement has not been tested here because it is somewhat more complicated to implement. Thermal Expansion Coefficient Derivative Method. The thermal expansion coefficient derivative method is

■

CONCLUSION Sometimes different flash packages are used to model the same fluid in different parts of the workflow: characterization, reservoir simulation, and well and surface simulation. If the same fluid models are used in these packages, the numerical flash results are usually fairly similar. However, if different phase identification methods are used, the phase labeling will be inconsistent and may result in incorrect modeling of flow. Ideally, all packages would use the same method. Realistically, this is unlikely to happen, and perhaps, the best that can be hoped for is that the packages are flexible enough to allow the user to choose their own preferred method. The pseudo-critical temperature method works well in reservoir simulation but is a poor choice for a general thermodynamics program. The negative flash method gives consistent results but lacks the ability to tune the results, which may be necessary in some practical applications. The new method based on the thermal expansion coefficient has been shown to give a consistent identification of reservoir fluid phases but requires more effort to implement.

■

APPENDIX Two methods require derivatives of the volume. Because the cubic EoSs specify the pressure as a function of the volume and temperature, it is easier to convert the volume derivatives to pressure derivatives. The temperature derivative of a function F(V,T) keeping pressure constant is: ⎛ dF ⎞ ⎛ dV ⎞ ⎜ ⎟ = F + F ⎜ ⎟ T V ⎝ dT ⎠ P ⎝ dT ⎠ P

(A1)

where FT = (dF/dT)V and FV = (dF/dV)T. Applying eq A1 to F = P, noting (dP/dT)P = 0, and rearranging gives the cyclic rule: ⎛⎛ d P ⎞ ⎛ d P ⎞ ⎞ − P ⎛ dV ⎞ T ⎜ ⎟ = − ⎜⎜ ⎟ /⎜ ⎟ ⎟ = ⎝ dT ⎠ P PV ⎝⎝ d T ⎠ V ⎝ d V ⎠ T ⎠

(A2)

Therefore, eq A1 can be written as: ⎛ dF ⎞ ⎜ ⎟ = F − F (P / P ) T V T V ⎝ dT ⎠ P

(A3)

The temperature derivative of k is required for the isothermal compressibility method. Equation 8 can be written as k = −1/ (VPV). Applying eq A3 to F = k gives: ⎛ dk ⎞ ⎜ ⎟ = ⎝ dT ⎠ P H

{VP

VT

P

− (PV + VPVV ) PT V

2

(VPV )

} (A4)


Energy & Fuels

■

Because V, PT, and (VPV)2 are positive, Venkatarathnam and Oellrich10 use the following condition for a liquid:

PVT P 1 − VV > PT PV V

The temperature derivative of α, required for the thermal expansion coefficient method, is given by differentiating eq 9:

(A6)

Because V > 0, the following criteria can be used for a liquid: ⎛ ∂ 2V ⎞ ⎛ ∂V ⎞ V⎜ 2 ⎟ − ⎜ ⎟ > 0 ⎝ ∂T ⎠P ⎝ ∂T ⎠P

(A7) 2

2

where (dV/dT)P is given by eq A2 and (∂ V/∂T )P can be calculated by applying eq A3 to F = (dV/dT)P = −PT/PV:

{

■

2

( ddVT )P + PVV( ddVT )P }

PTT + 2PVT ⎛ ∂ 2V ⎞ ⎜ 2⎟ = − ⎝ ∂T ⎠P

PV

REFERENCES

(1) Corradini, D.; Rovere, M.; Gallo, P. The Widom line and dynamical crossover in supercritical water: Popular water models versus experiments. J. Chem. Phys. 2015, 143 (11), 114502. (2) Imre, A. R.; Ramboz, C.; Deiters, U. K.; Kraska, T. Anomalous fluid properties of carbon dioxide in the supercritical region: application to geological CO2 storage and related hazards. Environ. Earth Sci. 2015, 73 (8), 4373−4384. (3) Pedersen, K. S.; Christensen, P. L.; Shaikh, J. A. Phase Behavior of Petroleum Reservoir Fluids, 2nd ed.; CRC Press: Boca Raton, FL, 2014. (4) Mattox, C. C.; Dalton, R. L. Reservoir Simulation; Society of Petroleum Engineers: Richardson, TX, 1990. (5) Brill, J. P.; Mukherjee, H. Multiphase Flow in Wells; Society of Petroleum Engineers: Richardson, TX, 1999. (6) Michelsen, M. L. Calculation of Critical Points and Phase Boundaries in the Critical Region. Fluid Phase Equilib. 1984, 16 (1), 57−76. (7) Peng, D.-Y.; Robinson, D. B. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 1976, 15 (1), 59−64. (8) Soave, G. Equilibrium constants from a modified Redlich− Kwong equation of state. Chem. Eng. Sci. 1972, 27 (6), 1197−1203. (9) Poling, B. E.; Grens, E. A.; Prausnitz, J. M. Thermodynamic properties from a cubic equation of state: Avoiding trivial roots and spurious derivatives. Ind. Eng. Chem. Process Des. Dev. 1981, 20 (1), 127−130. (10) Venkatarathnam, G.; Oellrich, L. R. Identification of the phase of a fluid using partial derivatives of pressure, volume, and temperature without reference to saturation properties: Applications in phase equilibria calculations. Fluid Phase Equilib. 2011, 301 (2), 225−233. (11) Gosset, R.; Heyen, G.; Kalitventzeff, B. An efficient algorithm to solve cubic equations of state. Fluid Phase Equilib. 1986, 25 (1), 51− 64. (12) Asselineau, L.; Bogdanic, G.; Vidal, J. A versatile algorithm for calculating vapour−liquid equilibria. Fluid Phase Equilib. 1979, 3 (4), 273−290. (13) Pope, G. A.; Sepehrnoori, K.; Delshad, M. A New Generation Chemical Flooding Simulator; U.S. Department of Energy: Washington, D.C., 2005; Final Report DE-FC-26-00BC15314. (14) Schlumberger. ECLIPSE: Reservoir Simulation Software, Technical Description, Version 2015.1; Schlumberger: Abingdon, U.K., 2015. (15) Li, C. C. Critical temperature estimation for simple mixtures. Can. J. Chem. Eng. 1971, 49 (5), 709−710. (16) KBC Process Technology, Ltd. User Guide for Multiflash for Windows, Version 6.0; KBC Process Technology, Ltd.: London, U.K., 2015. (17) Perschke, D. R. Development and application of an equation of state compositional simulator. Ph.D. Dissertation, The University of Texas at Austin, Austin, TX, 1988. (18) Whitson, C. H.; Michelsen, M. L. The Negative Flash. Fluid Phase Equilib. 1989, 53, 51−71. (19) Rachford, H. H., Jr.; Rice, J. D. A comparative study of reducedvariables-based flash and conventional flash. Trans. Am. Inst. Min. Metall. Eng. 1952, 195, 327−328. (20) Michelsen, M. L.; Mollerup, J. Thermodynamic Models: Fundamentals & Computational Aspects, 2nd ed.; Tie-Line Publications: Holte, Denmark, 2007. (21) Wilson, G. A modified Redlich−Kwong equation of state applicable to general physical data calculations. Proceedings of the 65th Annual AIChE National Meeting; Cleveland, OH, May 4−7, 1968; Paper 15C. (22) Pasad, G. V.; Venkatarathnam, G. A method for avoiding trivial roots in isothermal flash calculations using cubic equations of state. Ind. Eng. Chem. Res. 1999, 38 (9), 3530−3534. (23) Jayanti, P. C.; Venkatarathnam, G. Identification of the phase of a substance from the derivatives of pressure, volume and temperature, without prior knowledge of saturation properties: Extension to solid phase. Fluid Phase Equilib. 2016, 425, 269−277. (24) Williams, M. J.; Reza, Z. Personal Communication, 2007.

(A5)

⎛ 1 ⎞⎛ d V ⎞ ⎛ dα ⎞ ⎛ 1 ⎞⎛ ∂ 2V ⎞ ⎟⎜ ⎜ ⎟ = −⎜ ⎟ + ⎜ ⎟⎜ ⎟ ⎝ dT ⎠ P ⎝ V 2 ⎠⎝ dT ⎠P ⎝ V ⎠⎝ ∂T 2 ⎠ P

Article

(A8)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jim Bennett: 0000-0002-6167-4899 Notes

The authors declare no competing financial interest.

■

NOMENCLATURE α = thermal expansion coefficient β = vapor-phase fraction σ = surface tension Ωa0 = cubic equation of state constant Ωb0 = cubic equation of state constant a = cubic equation of state a parameter, representing attractive forces b = cubic equation of state b parameter, representing the volume at infinite pressure G = Rachford−Rice equation k = isothermal compressibility kr, kro, and krg = relative permeabilities Kj = K value for component j P = pressure Pcw, Pcow, and Pcgw = capillary pressures Q = phase identification parameter r1 and r2 = tuning parameters R = gas constant Sw = water saturation (volume fraction) T = temperature Tc = phase pseudo-critical temperature Tcj = critical temperature of component j in the phase V = phase molar volume Vc = phase pseudo-critical molar volume Vcj = critical volume of component j in the phase w = rock wettability wj = weights used in the Li mixing rule xj = mole fraction of component j in the phase Zc = phase pseudo-critical Z factor I


Article

Energy & Fuels (25) Troncoso, J.; Navia, P.; Romaní, L.; Bessieres, D.; Lafitte, T. On the isobaric thermal expansivity of liquids. J. Chem. Phys. 2011, 134 (9), 094502. (26) Perschke, D. R.; Pope, G. A.; Sepehrnoori, K. Phase Identification during Compositional Simulation. Soc. Pet. Eng. 1990, SPE-19442-MS. (27) Xu, Z.; Okuno, R. Numerical Simulation of Three-Hydrocarbon-Phase Flow with Robust Phase Identification. Proceedings of the SPE Reservoir Simulation Symposium; Houston, TX, Feb 23−25, 2015; SPE-173202-MS, DOI: 10.2118/173202-MS. (28) Li, Z.; Firoozabadi, A. General strategy for stability testing and phase-split calculation in two and three phases. SPE J. 2012, 17 (4), 1096−1107. (29) Gozalpour, F.; Danesh, A.; Todd, A. C.; Tehrani, D.-H.; Tohidi, B. Vapour−liquid equilibrium volume and density measurements of a five-component gas condensate at 278.15−383.15 K. Fluid Phase Equilib. 2003, 206, 95−104. (30) Robinson, D. B.; Peng, D.-Y. The Characterization of the Heptanes and Heavier Fractions for the GPA Peng−Robinson Programs; Gas Processors Association (GPA): Tulsa, OK, 1978; Research Report RR-28.

J


Comparison of Phase Identification Methods Used ... - ACS Publications

Recommend Documents