Accurate prediction of critical liquid-carrying velocity is essential for liquid loading judgment, keeping stable production of water-bearing gas well, and exploiting the potential of underground coal gasification^[1,2,3]. At the same time, it is of positive significance for promoting the rapid scale development of tight oil and gas and speed up the realization of “energy independence” of China^[4,5]. More and more experimental data suggests that inclined gas wells have the lowest liquid-carrying capacity, and accurately predicting the critical liquid-carrying velocity in an inclined gas well is also the most difficult and complicated^[6,7]. An in-depth study on the prediction of critical liquid-carrying velocity in inclined gas wells is necessary. Stratified flow is a common flow pattern of gas-liquid two-phase flow in inclined pipes. Usually, the circumferential distribution of two phases of gas and liquid in the wellbore has a crucial influence on the flow characteristics, especially in the case with a small flow rate of liquid phase and large flow rate of gas phase^[8,9]. There is no report on study of quantitative investigation of gas-liquid distribution and computation model of interfacial friction factor at critical conditions in inclined pipe. Based on this, the structural parameter of the gas-liquid distribution at critical conditions was determined according to the principle of minimum energy by investigating the variations in potential energy, gas-phase kinetic energy, and surface free energy of the two-phase system per unit length of tube. Meanwhile, the experimental data was used to modify the calculation model of the interfacial friction factor to close the prediction model. Finally, the prediction model was used to predict the critical liquid-carrying velocities in inclined gas wells to verify the effect of the model.

To simulate the dynamic process of liquid-carrying in gas wells, the visualization experiment was designed as shown in Fig. 1.

Water and glycerol-water solution (with the volume ratio of 1:1) were used to simulate the bottom-hole liquid loading, and air was used to simulate the natural gas. Under the experimental conditions, the air had a density ρ_a of 1.184 kg/m³ and a dynamic viscosity μ_a of 1.849×10^-5 Pa•s; water had a density ρ_w of 997.05 kg/m³ and a dynamic viscosity μ_w of 8.900 8×10^-4 Pa•s; the glycerol-water solution had a density ρ_gw of 1139.56 kg/m³ and a dynamic viscosity μ_gw of 6.882×10^-3 Pa•s; the surface tension of water and air σ_aw was 0.071 97 N/m and the surface tension of air and glycerol-water solution σ_agw was 0.067 11 N/m.

The experiments were carried out respectively in Plexiglas tubes with outer diameters of 50 mm and 70 mm and wall thickness of 5 mm. The experiments were performed with the following steps. (1) A certain volume of liquid was filled into the low-lying of the tube, and compressed air provided by fans was injected into the tube, and pressure and flow rate were measured continuously. (2) The frequency converter was used to continuously adjust the gas flow rate, and until the liquid accumulated in the low-lying spread in the upward inclined section, which remained at rest relative to the tube (critical condition was reached).

At critical gas flow rate of liquid-carrying, the flow pattern in the inclined tube is stratified wavy flow, and the liquid film is at equilibrium state under the joint action of gravity, interfacial shear stress and wall shear stress. In this case, the liquid film near the gas phase is transported downstream by the interfacial shear stress, while the liquid film near the wall flows back under the action of gravity, which results in a circulating flow in the liquid film. Force analysis of the liquid film per unit length in the inclined tube is shown in Fig. 2.

The above flow process can be simplified into a steady-state one-dimensional flow. Without consideration of the acceleration pressure drop, momentum equations of gas and liquid phases can be expressed as follows:

The average pressure gradients of gas and liquid phases can be assumed approximately equal during the flow process, and the governing equation can be obtained by combining Eqs. (1) and (2):

In the above formula, the geometric parameters are related to the gas-liquid distribution, and the stress parameters are related to the liquid holdup, velocity of each phase, and friction factor. The governing equation can be closed by closing geometrical and shear stress parameters.

According to the study by Brauner et al.^[10], the cases of circumferential distribution of liquid film are shown in Fig. 3. In Fig. 3, O and O₁ indicate the centers of the cross-section of the wellbore and the virtual circle, respectively. The geometric parameters in the governing equation can be calculated by the following formulas:

The calculation formulas of geometric parameters under the assumption of different gas-liquid distribution are listed in Table 1.

Gravity, shear stress, and surface tension will cause variations in the potential energy, kinetic energy, and surface free energy of the system in the process of interface shape changing from plane to a curved interface. The total energy variation per unit length of the gas-liquid two-phase system at critical conditions can be expressed as:

Based on the principle of minimum energy, the gas- liquid interface will be the most stable when the total energy per unit length of the gas-liquid two-phase system reaches the minimum value. This relationship is as follows:

(1) The variation per unit length in the gravitational potential energy can be written as follows^[10] :

where

(2) The variation in the surface free energy per unit length can be expressed as:

where

(3) The variation in the kinetic energy per unit length of the gas phase.

According to the one-seventh power law, the gas velocity distribution in turbulent flow can be expressed as:

where

The variation in the kinetic energy per unit length of the gas phase can be expressed as:

According to Taitel et al.^[11], the relationships of relevant stresses are given below:

Based on experimental investigation, Banafi et al.^[12] pointed out that the following relationship proposed by Taitel could predict gas-wall friction factor accurately of gas-liquid two-phase flow at low liquid holdup:

where

C_g=16, m=1 for laminar flow and C_g=0.046, m=0.2 for turbulent flow.

Biberg^[13] pointed out that the relation between the average liquid-to-wall shear stress and the average interfacial shear stress under smooth stratified flow in inclined pipes can be expressed as:

In the above formula, $f\left( {{\varphi }_{0\text{P}}} \right)$is an integral item that is difficult to obtain an analytical solution, so the rational approximation was given by Biberg^[13] as below:

The coefficients a_n and b_n used in Eq. (26) are listed in Table 2^[13]. The maximum absolute error in the approximation is less than 2×10^-6.

Accurate calculation of the interfacial friction factor is the key to the prediction of interfacial shear stress. Interface waves and liquid holdup have a great influence on the interfacial friction factor. The characteristics of interface waves are primarily determined by the rate of transfer of mechanical energy from the gas phase to the liquid phase, and the rate of viscous dissipation within the liquid, and these processes are affected by the inclination angle^[14]. As the inclination angle increases, the restriction of gravity on the growth of interface waves weakens, making it easier to observe waves in inclined pipes. This paper suggests that gravity has a more significant effect on wave fluctuations at critical conditions in the inclined pipe. Thus, on the one hand, the interface waves were equivalent to the relative roughness of the wall, and on the other hand, the inclination angle was introduced to modify the interfacial friction factor^[14]. The specific calculation process is as follows:

In which, the dimensionless interfacial wetted perimeter is defined as:

Then the modified interfacial friction factor considering the influence of the inclinations can be expressed as:

where A and B are related to fluid properties and pipe conditions^[15].

The governing equation is regarded as an function of liquid holdup:$f\left( \varepsilon \right)=0$. Zero is always one of the solu-tions of$f\left( \varepsilon \right)=0$. The function, $f{{\left( \varepsilon \right)}^{\frac{1}{3}}}$ is graphically depicted in Fig. 4. When u_sg=u_sg,cr, the equation $f\left( \varepsilon \right)=0$ has only one non-zero solution in the interval (0, 0.1), and the solution procedure is shown in Fig. 5.

The dimensionless governing equation is given as:

where

From the above formulas, the gas-liquid distribution at critical conditions is related to the Bond number, Weber number of the gas phase, liquid holdup, and inclination angle. Based on the principle of minimum energy, the variation of the total energy of the two-phase system with the interfacial curvature was quantified, and taking air and water as an example, the effects of the above four parameters on the gas-liquid distribution were investigated respectively. The results are shown in Fig. 6 (the red data points in the figures represent the minimum energy under the experimental conditions). From the figure: (1) As the pipe diameter increases, the interfacial curvature corresponding to the most stable gas-liquid distribution (corresponding to the minimum energy) gradually changes from 125° (slightly convex interface) to 180° (flat interface). Conversely, the smaller the pipe diameter, the smaller the interfacial curvature under the influence of wettability and surface tension is, resulting in a convex interface shape (Fig. 6a). (2) With the increase of the deviation angle, the offset effect of gravity to wettability and surface tension becomes more significant, and the interfacial curvature corresponding to the most stable gas-liquid distribution changes from 120° (slightly convex interface shape) to 180° (Fig. 6b). (3) Fig. 6c-d shows that as the gas velocity and density increase, the kinetic energy of the gas phase also increases, and the interfacial curvature decreases (more convex) under the influence of wettability and surface tension, which can reduce the surface free energy of the system. (4) With the increase of liquid holdup, the influence of gravitational potential energy of the liquid phase increases, which is good for keeping the original flat interface. At the same time, the convex interface can also reduce the increase of surface free energy of the system caused by the increase of the liquid holdup. But in a low-pressure system where the liquid density is much greater than that of the gas phase, the gravitation potential energy of the liquid phase plays a dominant role (Fig. 6e), and the interfacial curvature corresponding to the most stable gas-liquid distribution gradually changes from 105° to 165° (near horizontal interface) with the increase of the liquid holdup. In conclusion, the interface is slightly convex at low liquid-holdup under the influence of surface tension and wettability, which is consistent with the experimental result shown in Fig. 7.

Fig. 8 shows the experimental results of the critical liquid-carrying velocity and pressure gradient. It can be seen that the critical liquid-carrying velocity and pressure gradient increase non-linearly with the increase of the inclination angle. At the same inclination angle, the greater the density or viscosity of the liquid phase, the greater the critical liquid-carrying velocity will be, which is consistent with the experimental results of Birvalski et al.^[16]. At the same time, the increase in pipe diameter makes the critical liquid-carrying velocity increase. This is because the change of the interface shape is not as significant as the influence of Bo and We_sg on the critical liquid-carrying velocity at the same liquid holdup. Therefore, to make the governing equation still hold at critical conditions, the gas velocity must be increased to offset the increment of Bo/We_sg ratio caused by the increase in pipe diameter.

The correction factor for the interfacial friction factor at each inclination angle was back-calculated based on the flat interface shape by using both the experimental data obtained by Birvalski et al.^[16] and current study, and the correction factor ($f\left( \alpha \right)$) was obtained by fitting of data. The results are listed in Table 3. After modified with the flow correction factor listed in Table 3, the calculation model of interfacial friction factor (Eq. (31)) gave results with high accuracy and high Pearson correlation coefficient.

From Table 3, the slope of $f\left( \alpha \right)$ is positive, which means that the interfacial friction factor increases with the increase of inclination angle on the whole. This is mainly because the increase of inclination angle makes the gravity component constraining the interface fluctuation decrease, resulting in more intense and obvious interface fluctuation. At the same time, $f\left( \alpha \right)$ is less than 1, indicating that the interfacial friction factor in slightly inclined tubes is smaller than that in horizontal tubes at the same gas velocity at low liquid holdup, and the main reason is that the influence of the inclination angle on the interface fluctuation is not as significant as that of liquid holdup in near-horizontal pipes. In this paper, based on the experimental data and dimensional analysis, the equation of correction factor of interfacial friction factor is given as:

The experimental results were compared with the calculated results from the model, as shown in Fig. 9. The critical liquid-carrying velocity predicted by the model increases with the increase of the pipe diameter, the liquid density and viscosity (the critical liquid-carrying velocity is greater for the glycerol-solution than that for the water at the same tube diameter), which also tallies with the results of Birvalski et al.^[16] and Rastogi et al.^[17]. The predicted pressure gradient decreases with the increase of the pipe diameter, which agrees with the experimental results. The increase of hydraulic diameter of liquid phase caused by tube diameter increase is an important reason for pressure gradient drop. Compared with the experimental data of this study, the critical liquid-carrying velocity and critical pressure gradient calculated by the model have an average relative error of 1.19% and 3.02%, respectively. Meanwhile, as the inclination angle increases, the critical liquid-carrying velocity first increases and then decreases, and the inclination angle corresponding to the maximum critical liquid-carrying velocity increases with the increase of pipe diameter and is also affected by the fluid properties. In the case of natural gas carrying mineralized water, the inclination angles, calculated based on the model under standard engineering condition, corresponding to the maximum critical liquid-carrying velocity are 30° (D=50.67 mm), 37° (D=62.00 mm), 40° (D=75.9 mm), and 45° (D=100.53 mm) respectively.

The predicted results from the model proposed in this paper need to be converted to the high-pressure condition. From the energy perspective, it is considered that the mechanical energy transferred between the gas and liquid phases of the unit volume at critical conditions in pipe systems with high pressure and low pressure are equal:

For easy comparison, the critical liquid-carrying velocity is converted to critical liquid-carrying flowrate under standard engineering conditions:

When the natural gas velocity is greater than the maximum critical liquid-carrying velocity in the inclined gas well, the gas produced from the well can carry liquid continuously, and the maximum critical liquid-carrying velocity in field can be worked out in conjunction with logging data. A total of 42 gas wells in Puguang gas field and Yanchang gas field were selected, and the field data were compared with the predicted results (Fig. 10). It can be seen that the model proposed in this paper misjudges only 1 gas well, with an error rate of 2.38%, indicating that the predicted results of the model have high credibility and can predict liquid loading in inclined gas wells effectively.

The interface is more likely to curve under the joint effect of wettability and surface tension at low liquid holdup. The smaller the pipe diameter, the smaller the liquid holdup, the smaller the deviation angle, the greater gas velocity, the greater the gas density, the more curved the interface will be.

The critical liquid-carrying velocity increases nonlinearly and then decreases with the increase of deviation angle. The deviation angle corresponding to the maximum critical liquid-carrying velocity decreases with the increase of the diameter of the wellbore, and it is also affected by the properties of the gas phase and liquid phase.

The critical gas velocity and critical pressure gradient predicted by the model in this paper have a mean relative error of 1.19% and 3.02% respectively, and the misjudgement rate is 2.38% in the field examples, implying the new model can provide valid judgment on the situation of liquid loading in inclined gas wells.

a_n, b_n, A, B—constant coefficients, dimensionless;

A_T—cross-sectional area of tubing, m²;

A_g—tube cross-sectional area occupied by gas phase, m²;

ΔA_gl—the variation of the contact area per unit length between gas and liquid, m²/m;

A_l—tube cross-sectional area occupied by liquid phase, m²;

ΔA_lw—the variation of contact area per unit length between liquid and wall, m²/m;

Bo—Bond number, dimensionless;

C_g, m—constant coefficient related to the flow pattern, dimensionless;

D—inner diameter of the tubing, m;

D_H,g—hydraulic diameter of the gas phase, m;

D_l(ϕ_0P)—equivalent laminar diameter at different distribution angles of liquid phase, m;

f_g—friction factor between gas and wall, dimensionless;

f_i—interfacial friction factor, dimensionless;

f_l—friction factor between liquid and wall, dimensionless;

f(α)—correction factor of the interfacial friction factor, dimensionless;

f(ϕ_0P)—function related to liquid distribution, dimensionless;

F—intermediate variable, dimensionless;

g—acceleration of gravity, m/s²;

k—interfacial absolute roughness, m;

n—ordinal number of constant coefficient, dimensionless;

p—pressure, Pa;

Q_cr—critical liquid-carry flowrate, m³/d;

r—distance between wellbore axis and the radial position, m;

R—radius of the tubing, m;

Re_g—Reynolds number of gas phase, dimensionless;

R₁—radius of the virtual circle, m;

S_g—wetted perimeter of gas phase, m;

S_i—interfacial wetted perimeter, m;

s_i—dimensionless interfacial wetted perimeter, dimensionless;

S_l—wetted perimeter of liquid phase, m;

T—thermodynamic temperature, K;

u_g—actual gas velocity, m/s;

${{\bar{u}}_{\text{g}}}$—average gas velocity, m/s;

u_g,max—the maximum gas velocity, m/s;

u_i—actual interface velocity, m/s;

u_l—actual liquid velocity, m/s;

u_sg—superficial gas velocity, m/s;

Δu_sg—the step size for iteration of superficial gas velocity, m/s;

u_sg,cr—critical superficial gas velocity, m/s;

u_sg,cr,hp—critical superficial gas velocity in high pressure system, m/s;

We_sg—superficial gas Weber number;

x—coordinate in the flow direction, m;

Y_GPC—the center of gravity of the two phases with curved interface, m;

Y_GP—the center of gravity of the two phases with flat interface, m;

Z—compression factor;

α—tube inclination angle from horizontal, °;

β—deviation angle, °;

ΔE_k,g—the variation of gas kinetic energy per unit length, J/m;

ΔE_p—the variation of gravitational potential energy of gas-liquid two-phase per unit length, J/m;

ΔE_s—the variation of surface energy of gas-liquid two-phase per unit length, J/m;

ΔE_tot—the variation of total energy of gas-liquid two-phase per unit length, J/m;

δ_l—thickness of liquid film, m;

δ_l(ϕ)—thickness distribution of liquid film, m;

ε—liquid holdup;

ε₀—the initial value of liquid holdup for iteration;

Δε—the step size for iteration of liquid holdup;

θ—the wetting angle, rad;

μ_a—air dynamic viscosity, Pa•s;

μ_g—gas dynamic viscosity, Pa•s;

μ_gw—dynamic viscosity of glycerol-water solution, Pa•s;

μ_l—liquid dynamic viscosity, Pa•s;

μ_w—water dynamic viscosity, Pa•s;

ρ_a—air density, kg/m³;

ρ_g—gas density, kg/m³;

ρ_g,hp—gas density at high pressure, kg/m³;

ρ_gw—density of glycerol-water solution, kg/m³;

ρ_l—liquid density, kg/m³;

ρ_w—water density, kg/m³;

σ_aw—surface tension between air and water, N/m;

σ_agw—surface tension between air and glycerol-water solution, N/m;

σ_gl—surface tension between gas and liquid, N/m;

τ_i—interfacial shear stress, Pa;

${{\bar{\tau }}_{\text{i}}}$—average interfacial shear stress, Pa;

τ_ig—interfacial shear stress acting on the gas, Pa;

τ_gi—interfacial shear stress acting on the interface, Pa;

τ_wg—shear stress between gas and wall, Pa;

τ_wl—shear stress between liquid and wall, Pa;

$\overline{{{\tau }_{\text{wl}}}}$—average shear stress between liquid and wall, Pa;

ϕ—liquid-phase distribution angle, rad;

ϕ₀—liquid-phase distribution angle for the curved interface, rad;

ϕ_0P—liquid-phase distribution angle for the flat interface, rad;

ϕ_PA—interfacial curvature, rad.

g—gas phase;

l—liquid phase.