A Direct Eulerian GRP Scheme for the Prediction of Gas-Liquid Two-Phase Flow in HTHP Transient Wells

Abstract

A coupled system model of partial differential equations is presented in this paper, which concerns the variation of the pressure and temperature, velocity, and density at different times and depths in high temperature-high pressure (HTHP) gas-liquid two-phase flow wells. A new dimensional splitting technique with Eulerian generalized riemann problem (GRP) scheme is applied to solve this set of conservation equations, where Riemann invariants are introduced as the main ingredient to resolve the generalized Riemann problem. The basic data of “X well” (HTHP well), 7100 m deep, located in Southwest China, is used for the case history calculations. Curve graphs of pressures and temperatures along the depth of the well are plotted at different times. The comparison with the results of Lax Friedrichs (LxF) method shows that the calculating results are more fitting to the values of real measurement and the new method is of high accuracy.

1. Introduction

The prediction of pressure and temperature of transient gas-liquid flow in a wellbore is important but difficult for well completion test because they are characterized by the dependence of pressure, density, velocity, and other flow parameters on both time and space. As for pressure prediction research, there exist empirical formulas, such as those given by Beggs and Brill [1], Mukherjee and Brill [2, 3], and so on. Different researchers such as Hurlburt and Hanratty [4] and Cazarez-Candia and Vásquez-Cruz [5] have proposed mechanistic models, assuming that flow is under steady-state conditions; other researchers such as Taitel et al. [6]; Ouyang and Aziz [7] have proposed unsteady-state gas-liquid two-phase flow models. Fontanilla and Aziz [8] and Ali [9] presented two simultaneous ordinary differential equations for estimating the steam pressure and quality and solved these equations by using the fourth-order Runge-Kutta method. However, those models can only predict the pressure profiles but not the temperature profiles and ignored their interdependence.

Concerning both pressure and temperature in HTHP wells, Wu et al. have presented a coupled system model of differential equations in [10], but this model only considered the single phase flow statement. In this paper, we build a set of coupled partial differential equations of pressure, temperature, density, and velocity in HTHP gas-liquid two-phase flow wells on the base of the model which was build by Xu et al. in [11]. The numerical model, which accords with the actual situation of the well, allows for the change of oblique angle,different heat transfer medium in annular and the depth of the physical properties of the formation.

We found an algorithm solving model with generalized Riemann problem (GRP) scheme, which is an analytic extension of the Godunov scheme in [12] and originally designed by Li and Chen in [13] for the shallow water equations. A direct and simple derivation of the Eulerian generalized Riemann problem scheme is presented to get the integration in time of the conservation laws. Riemann invariants are applied in order to resolve the singularity at the jump discontinuity. The approach has the advantage that the contact discontinuity in each local wave pattern is always fixed with speed zero, while the rarefaction and the shock waves are located on either side. Since the extension of this scheme to multidimensional cases is obtained using the dimensional splitting technique, getting the integration in time of the conservation laws is more direct and simple.

In this paper, we use GRP method for solving this problem and get more accurate prediction of pressure and temperature compared with those obtained from the existing correlations such as LxF method in [11]. The basic data for the calculation are from X well, 7100 m of depth in Southwest China. The curves of the gas pressure and temperature along the depth of the well are plotted. The results can provide a technical reliance for the process of designing well tests in HTHP gas-liquid wells and a dynamic analysis of production from wells.

2. Model Formulation

Considering the two-phase flow system shown in Figure 1, the mixture density and velocity are related to the in situ liquid volume fraction (holdup), 𝐻, as follows: 𝜌_𝑚=𝜌_𝑙𝐻+𝜌_𝑔()𝑣1−𝐻_𝑚=𝑣_𝑙𝐻+𝑣_𝑔()𝑢1−𝐻_𝑚=𝑢_𝑙𝐻+𝑢_𝑔(1−𝐻). (1)

Figure 1 

The two-phase flow.

2.1. Mass Balance

Consider the flow model shown in Figure 2. According to the fluid moves through the fixed control volume depicted by John and Anderson in [14], we have _𝑆𝜌_𝑚𝑣_𝑚𝜕𝑑𝑠=−𝜕𝑡_𝜈𝜌_𝑚𝑑𝑣_𝑚. (2)

Figure 2 

Control volume 1.

Under transient conditions, applied to the control volume in Figure 3, in the limit as 𝑑𝑧 becomes very small, the volume and surface integral in (2) becomes 𝜕𝜕𝑡_𝜈𝜌_𝑚𝜕𝑑𝑣=𝜌𝜕𝑡^𝑚𝐴𝑑𝑧_𝑆𝜌_𝑚𝑣_𝑚𝑑𝑠=𝜌_𝑚𝑣_𝑚𝑑𝑣_𝑚+𝜌_𝑚𝐴𝑑𝑣_𝑚+𝐴𝑣_𝑚𝑑𝜌_𝑚=𝑑𝜌_𝑚𝑣_𝑚𝐴. (3)

Figure 3 

Control volume 2.

Substituting (3) into (2), we get the mass balance equation: 𝜕𝜌_𝑚+𝜕𝜌𝜕𝑡^𝑚𝑣_𝑚𝜕𝑧=0. (4)

2.2. Momentum Balance

As shown in Figure 4, the integral form of the 𝑧 component the momentum equation can be written as follows with the external forces: 𝜕𝜕𝑡_𝜈𝜌_𝑚𝑢_𝑚𝑑𝑉+_𝑆𝜌_𝑚𝑢_𝑚𝑣_𝑚⋅𝑑𝑆=−_𝑆(𝑃𝑑𝑆)𝑑𝑧−𝜌_𝑚𝑔cos𝜃𝐴𝑑𝑧−𝜆𝜌_𝑚𝑣²_𝑚2𝑑𝐴𝑑𝑧, (5) where𝜌_𝑚𝑔cos𝜃𝐴𝑑𝑧 is the force of gravity,(𝜆𝜌_𝑚𝑣²𝐴/2)𝑑𝑧is the shear stress, and 𝜕𝜕𝑡_𝜈𝜌_𝑚𝑢_𝑚𝜕𝑑𝑉=𝜌𝜕𝑡^𝑚𝑣_𝑚𝐴𝑑𝑧_𝑆𝜌_𝑚𝑢_𝑚𝑣_𝑚⋅𝑑𝑆=−𝜌_𝑚𝑣²_𝑚𝐴+𝜌_𝑚+𝑑𝜌_𝑚𝑣_𝑚+𝑑𝑣_𝑚²(𝐴+𝑑𝐴)−_𝑆(𝑃𝑑𝑆)_𝑧=−𝑃𝐴+(𝑃+𝑑𝑃)(𝐴+𝑑𝐴)−2𝑃𝑑𝐴2. (6)

Figure 4 

Control volume 3.

Substituting (6) into (5), we obtain momentum balance equation: 𝜕𝜌𝜕𝑡^𝑚𝑣_𝑚+𝜕𝜕𝑧𝑃+𝜌_𝑚𝑣²_𝑚=−𝜌_𝑚𝑔cos𝜃−𝜆𝜌_𝑚𝑣²_𝑚.2𝑑 (7)

2.3. Energy Balance

For the transient flow, it leads directly to the energy equation in terms of temperature. As shown in Figure 5, we should consider the heat transmission within wellbore and from wellbore to formation as transient.

Figure 5 

The radial transfer of heat.

According to the energy balance law, the heat variation flowing on control volume that is equal to the combination heat of inflow and outflow, and the heat transferring to the second dimension, we get the energy balance equation of transient flow: 𝑤𝐶_𝑝𝑇𝑧−𝑤𝐶_𝑝𝑇(𝑧+𝑑𝑧)−2𝜋𝑟_𝑡𝑜𝑈_𝑡𝑜𝑇−𝑇_𝑟=𝜕𝜌𝑑𝑧_𝑚𝐶_𝑝𝑇𝜕𝑡𝐴𝑑𝑧, (8) where 𝑇_𝑟=(𝐾_𝑒𝑇_𝑒𝑖+𝑟_𝑡𝑜𝑈_𝑡𝑜𝑇_𝑒𝑇_wbD)/(𝐾_𝑒+𝑟_𝑡𝑜𝑈_𝑡𝑜𝑇_wbD) and 𝑤=𝜌_𝑚𝑣_𝑚𝐴. Equation (8) equals the following equation: 𝜕𝜌_𝑚𝑣_𝑚𝑇+𝜕𝜌𝜕𝑧^𝑚𝑇=𝜕𝑡2𝜋𝐾_𝑒𝑟_𝑡𝑜𝑈_𝑡𝑜𝑇_𝑒𝑖−𝑇𝐶_𝑃𝐾_𝑒+𝑟_𝑡𝑜𝑈_𝑡𝑜𝑇_wbD. (9)

Finally, we obtain the coupled system model of partial differential equations: 𝜕𝜌_𝑚+𝜕𝜌𝜕𝑡^𝑚𝑣_𝑚𝜕𝜌𝜕𝑧=0,_𝑚𝑣_𝑚+𝜕𝜕𝑡𝑃+𝜌_𝑚𝑣²_𝑚𝜕𝑧=−𝜌_𝑚𝜆𝑔cos𝜃−,𝜕𝜌2𝑑_𝑚𝑇+𝜕𝜌𝜕𝑡^𝑚𝑣_𝑚𝑇=𝜕𝑧2𝜋𝐾_𝑒𝑟_𝑡𝑜𝑈_𝑡𝑜𝑇_𝑒𝑖−𝑇𝐶_𝑃𝐾_𝑒+𝑟_𝑡𝑜𝑈_𝑡𝑜𝑇_wbD,𝜌_𝑚=𝑀𝑃𝛾_𝑔,𝑍𝑅𝑇𝑃,𝑇,𝜌_𝑚,𝑣_𝑚=𝑃₀,𝑇₀,𝜌_𝑚0,𝑣_𝑚0𝑡=0,𝑧=0. (10)

3. Format Construction

We unify the conservation equations (4), (7), and (9) which are also included in (10) into the following formation: 𝜕𝐴_𝑚+𝜕𝑡𝜕𝐵_𝑚𝜕𝑧=𝐶_𝑚𝐴,𝑚=1,2,3,𝑈=¹=𝜌_𝑚𝐴₂=𝜌_𝑚𝑣_𝑚𝐴₃𝐹(𝑈)=𝐵=0,¹=𝜌_𝑚𝑣_𝑚𝐵₂=𝑃+𝜌_𝑚𝑣²_𝑚𝐵₃=𝜌_𝑚𝑣_𝑚𝑆(𝑈)=𝐶𝑇,₁𝐶=0₂=−𝜌_𝑚𝜆𝑔cos𝜃−𝐶2𝑑₃=2𝜋𝐾_𝑒𝑟_𝑡𝑜𝑈_𝑡𝑜𝑇_𝑒𝑖−𝑇𝐶_𝑃𝐾_𝑒+𝑟_𝑡𝑜𝑈_𝑡𝑜𝑇_wbD. (11)

We define the equally spaced grid points, the interface points, and the cells as 𝑧_𝑗=𝑗Δ𝑧,𝑧_𝑗+1/2=𝑧_𝑗+𝑧_𝑗+12,𝐶_𝑗=𝑧_𝑗−1/2,𝑧_𝑗+1/2. (12)

We assume that the data at time𝑡=𝑡_𝑛are piecewise linear with a slope𝜎^𝑛_𝑗and we have𝑈(𝑧,𝑡_𝑛)=𝑈^𝑛_𝑗+𝜎^𝑛_𝑗(𝑧−𝑧_𝑗),𝑧∈(𝑧_𝑗−1/2,𝑧_𝑗+1/2).

The second-order Godunov scheme for (11) takes the following form: 𝑈_𝑗^𝑛+1=𝑈^𝑛_𝑗−(Δ𝑡/Δ𝑥)(𝐹(𝑈^𝑛+1/2_𝑗+1/2)−𝐹(𝑈^𝑛+1/2_𝑗−1/2)), where 𝑈^𝑛+1/2_𝑗+1/2is the midpoint value or the value of 𝑈 at the cell interface(𝑧_𝑗+1/2,𝑡_𝑛)with accuracy of second order. More specifically, the mid-point value𝑈^𝑛+1/2_𝑗+1/2is computed with the formulas 𝑈^𝑛+1_𝑗+1/2=𝑈^𝑛_𝑗+1/2+(Δ𝑡/2)(𝜕𝑈/𝜕𝑡)^𝑛_𝑗+1/2and 𝑈^𝑛_𝑗+1/2=𝑅^𝐴(0;𝑈^𝑛_𝑗+1/2,−,𝑈^𝑛_𝑗+1/2,+). Also,  𝑅^𝐴((𝑧−𝑧_𝑗+1/2)/(𝑡−𝑡_𝑛);𝑈^𝑛_𝑗+1/2,−,𝑈^𝑛_𝑗+1/2,+)is the solution of the Riemann problem centered at(𝑧_𝑗+1/2,𝑡_𝑛). Moreover, 𝑈^𝑛_𝑗+1/2,− and 𝑈^𝑛_𝑗+1/2,+ are the limiting values of initial data𝑈(𝑧,𝑡_𝑛) on both sides of (𝑧_𝑗+1/2,𝑡_𝑛). We present a direct and simple derivation of the Eulerian geralized Riemann problem (GRP) scheme and apply Riemann invariants in order to resolve the singularity at the jump discontinuity.

The local wave configuration is usually piecewise smooth and consists of rarefaction waves, shocks, and contact discontinuities. As the general rarefaction waves are considered, the initial data can be regarded as a perturbation of the Riemann initial data𝑈_𝐿 and 𝑈_𝑅.

The GRP scheme assumes piecewise linear data for the flow variables, which leads to the generalized Riemann problem for (11) subject to the initial data: 𝑈()=𝑈𝑧,0^𝐿+𝑧𝑈^′_𝐿𝑈,𝑧<0_𝑅+𝑧𝑈^′_𝑅,𝑧>0, (13) where𝑈_𝐿,𝑈_𝑅,𝑈^′_𝐿, and𝑈^′_𝑅are constant vectors.

The initial structure of the solution is determined by the associated Riemann solution, denoted by lim_𝑡→0𝑈(𝜆𝑡,𝑡)=𝑅^𝐴(𝜆;𝑈_𝐿,𝑈_𝑅),𝜆=𝑥/𝑡.

4. Solving Process

Step 1. Set the step length. In this paper, ℎ=1(m),𝜏=60(s). (14)

Step 2. Obtain each point’s inclination:𝜃_𝑗=𝜃_𝑗−1+(𝜃_𝑗′−𝜃_𝑗′−1)ℎ/Δ𝑠_𝑗′−1.

Step 3. The in situ liquid volume fraction (holdup) in (1) can be calculated from 𝐻𝜃_𝑗=0.98𝐸_𝑙^0.4846𝐹𝑟^0.0868×1+1−𝐸_𝑙ln4.7𝑁^0.1244_𝑣𝑙𝐸_𝑙^0.8692𝐹𝑟^0.5056×sin1.8𝜃_𝑗−13sin³1.8𝜃_𝑗. (15)

Step 4. Calculate the following parameters by Liao and Feng in [15]: 𝑈_𝑡𝑜=𝑟_𝑡𝑜𝑟_𝑡𝑖ℎ_𝑟+𝑟_𝑡𝑜ln𝑟_𝑡𝑜/𝑟_𝑡𝑖𝑟_𝑡𝑖ℎ_𝑟+1ℎ_𝑐+ℎ_𝑟+𝑟_𝑡𝑜ln𝑟_𝑐0/𝑟_𝑐𝑖𝑟_cas+𝑟_𝑡𝑜ln𝑟_ℎ/𝑟_𝑐0𝑟_cem⁻¹,𝑇_wbD=2𝜋𝐾_𝑒𝑇_𝑒𝑖−𝑇_wb∑^𝑚_𝑗=1𝑄_𝑗−1−𝑄_𝑗. (16)

Step 5. For piecewise given initial data𝑈^𝑛(𝑧)=𝑈^𝑛_𝑗+𝜎^𝑛_𝑗(𝑧−𝑧_𝑗),𝑧∈(𝑧_𝑗−1/2,𝑧_𝑗+1/2), we solve the Riemann problem for (11) to define the Riemann solution𝑈^𝑛_𝑗+1/2=𝑅(0;𝑈^𝑛_𝑗+(Δ𝑧/2)𝜎^𝑛_𝑗,𝑈^𝑛_𝑗+1−(Δ𝑧/2)𝜎^𝑛_𝑗+1), which is the same as the classical Godunov scheme and the Riemann solver in [16] is used in the solution.

Step 6. Determine(𝜕𝑈/𝜕𝑡)^𝑛_𝑗+1/2and evaluate the new cell averages𝑈_𝑗^𝑛+1. We apply monotonic algorithm slope limiters to suppress the local oscillations near discontinuities. We use parameter𝛼=1.9in𝜎_𝑗^𝑛+1=minmod(𝛼(𝑈_𝑗^𝑛+1−𝑈^𝑛+1_𝑗−1)/Δ𝑧,𝜎_𝑗^𝑛+1,−,𝛼(𝑈^𝑛+1_𝑗+1−𝑈_𝑗^𝑛+1)/Δ𝑧), where 𝑈^𝑛+1,−_𝑗+1/2=𝑈^𝑛_𝑗+1/2+Δ𝑧(𝜕𝑈/𝜕𝑧)^𝑛_𝑗+1/2 and 𝜎_𝑗^𝑛+1,−=(1/Δ𝑧)(Δ𝑈)_𝑗^𝑛+1,−=(1/Δ𝑧)(𝑈^𝑛+1,−_𝑗+1/2−𝑈^𝑛+1,−_𝑗−1/2).

5. Results and Discussion

In this simulation, we study a pipe in X well located in Sichuan Basin, Southwest China. All the needed parameters are given in [17] as follows: fluid density is 1000 kg/m³; depth of the well is 7100 m; friction coefficient is 1.2; ground temperature is 160°C; ground thermal conductivity parameter is 2.06; ground temperature gradient is 0.0218°C/m. Parameters of pipes are given in Table 1. Inclination, azimuth, and vertical depth are given in Table 2.

Table 1 

Parameters of pipes.

Table 2 

Parameters of azimuth, inclination, and vertical depth.

Through the simulation, we use GRP method to calculate the prediction of pressure and temperature of the oil in the pipe and draw a sensitive analysis for the results. We compare the results of pressure and temperature calculated for the well head at 1200 s by GRP and LxF scheme with the measurement results, which also shows that GRP scheme is more accurate in the real calculation. We obtain series of results contained in tables and figures and analyze these results as follows.

When the bottom pressure is 70 MPa, temperatures are plotted in Figure 6 at different depths and shown in detail in Table 3. When the output keeps constant, the temperature increases with the increasing depth of the well and when the depth fixed, the temperature increases with the increasing time. In addition, it can be seen from the figure that the temperature changes quickly in the early stage but stabilizes over time, especially after 1200 s.

Table 3 

Temperature at different depths on 300 s, 900 s, 1200 s, and 3600 s.

Figure 6 

Temperature distribution at different depths.

It is established that, when depth is constant, the pressure shown in Figure 7 and Table 4 increased with an increase of the time. When the output keeps constant, the pressure increased with the increasing depth of the well. This is because, with time increasing, the flow increases and then the frictional heat leads to an increase in the pressure. It can also be seen that the pressure changes quickly in the early stage but stabilizes over time.

Table 4 

Pressure at different depths on 300 s, 900 s, 1200 s, and 3600 s.

Figure 7 

Pressure distribution at different depths.

As shown in Table 5, for the comparative results of the well head temperature at 1200 s, the relative error between the calculation results and the measurement results of GRP scheme method is 5.12% and by LxF method is 6.70%, while the relative error between the results in pressure predition at the same time calculated by GRP scheme method is 8.81% and by LxF method is 9.73%, which shows that the distribution prediction of the two-phase flow is more accurate in actual calculation by GRP scheme method.

Table 5 

Comparative results of the well head at 1200 s.

6. Conclusion

In this paper, considering the variation of pressure, temperature, velocity; and density at different times and depths in gas-liquid two-phase flow, we present a system model of partial differential equations according to mass, momentum, and energy. We establish an algorithm solving model with a new difference method with a direct Eulerian GRP scheme which is proven to be efficient for the numerical implementation in this paper. The basic data of the X well (HTHP well), 7100 m deep in Sichuan Basin, Southwest China, was used for case history calculations, and a sensitivity analysis is completed for the model. The gas-liquid’s pressure and temperature curves along the depth of the well are plotted, and the curves intuitively reflect the flow law and the characteristics of heat transfer in formation. The results can provide the technical reliance for the process of designing well tests in high temperature-high pressure gas-liquid two-phase flow wells and dynamic analysis of production. Furthermore, the works in this paper can raise safety and reliability of deep completion test and will yield notable economic and social benefits and avoid or lessen accidents caused by improper technical design.

Nomenclature

𝐴:	A total length of conduit (m2)
𝐶𝐽:	Joule-Thompson coefficient (K/pa)
𝐶𝑝:	Heat capacity (J/kg·K)
𝐷:	A hydraulic diameter (m)
𝐺:	Acceleration constant of gravity (m/s2)
𝐾𝑒:	Formation conductivity (J/m·K)
𝑃:	Pressure (KPa)
𝑟𝐷:	Dimensionless radius
𝑟𝑡𝑜:	Outer radius of conduit (m)
𝑇:	Temperature (K)
𝑡𝐷:	Dimensionless time
𝑇𝑒:	Temperature of the stratum (K)
𝑇wb:	Wellbore temperature (K)
𝑇wbD:	Dimensionless wellbore temperature (K)
𝑇𝑟:	Temperature of the second surface (K)
𝑇𝑒𝑖:	Initial temperature of formation (K)
𝑈𝑡𝑜:	Overall-heat-transfer coefficient (W/m·K)
𝑉:	Velocity (m/s)
𝑍:	A total length of conduit (m)
𝑧:	The distance coordinate in the direction along the conduit
ℎ𝑐:	Heat transfer coefficient for natural convection based on outside tubing surface and the temperature difference between outside tubing and inside casing surface
ℎ𝑟:	Heat transfer coefficient for radiation based on the outside tubing surface and the temperature difference between the outside tubing and inside casing surface
𝐾cas:	Thermal conductivity of the casing material at the average casing temperature
𝐾cem:	Thermal conductivity of the cement at the average cement temperature and pressure
𝜆:	The friction coefficient, dimensionless
𝛾𝑔:	Euler constant 1.781
𝜌:	Density (kg/m3)
𝜃:	Inclination angle flow conduit.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research was supported by the Key Program of NSFC (Grant no. 70831005) and the Key Project of China Petroleum and Chemical Corporation (Grant no. GJ-73-0706).