Along with the progressive development of high-temperature, high-pressure and high gas rate well, gas recovery becomes more and more difficult, higher requirements have been set on wellbore integrity^[1,2,3,4,5]. In the course of development of high gas rate well with open-flow capacity greater than 120×10⁴ m³/d in particular, the production fluctuates along with time constantly due to the impacts of geological and tubing conditions, and well opening and shut-in are frequently required^[6,7]. In view of the high production and high flow velocity of gas well, the opening and shut-in speed of high production well pose heavier impacts on the pressure distribution in the tubing. The faster operation of well opening and shut-in would lead to large pressure fluctuation, and trigger severe wellbore integrity issue^[8]. On the contrary, the slower the operation of well opening and shut-in, the longer time will consume, which will cause troubles for field operations. On October 4, 1981, a gas well blowout took place in Texas, and the reason was that the water hammering effect caused by fast well opening and shut-in gave rise to high pressure^[9]. During initial stage of production test of Well Heba 1 in the northeast of Sichuan Basin, fast shut-in impaired the normal production of the gas well and even caused natural gas explosion venting^[10]. It is crucial for maintaining the wellbore integrity to study the instant flow process in the tubing during well opening and shut-in and optimize the opening and shut-in time based on actual operating conditions.

The flow inside tubing can be deemed an unsteady pressure flow, under this flow state, any instant change of flow velocity would cause a drastic pressure fluctuation inside the tubing^[11]. During gas well production, when the valve opening starts to change, the flow velocity in the tubing will change accordingly and generate surge pressure. The change velocity of the valve opening will decide the intensity of the surge pressure. When the valve opens or closes in a sudden, the surge pressure will probably be several or even tens of times higher than the normal pressure, and the repeated variation of pressure will make oil tubing and equipment vibrate, and even cause damage to oil tubing and wellhead equipment^[12] and incident, destroying the wellbore integrity and safe production.

The surge pressure produced during instant closing of valve is also known as water hammer. The water hammer manifests the compressibility of fluid. It is a transient flow issue which can be interpreted as the changes of all dynamic physical variables (speed, pressure and etc.) with time and position in flow domain in the course of flow. At present, the water hammer issue has drawn attention from numerous experts and has become one of the major research focuses in the field of hydrodynamics. In 1993, Jarine et al.^[13] proposed the "hard shut-in" and "soft shut-in" and pointed out the strengths and weaknesses of the hard and soft shut-ins. Erika et al.^[14] predicted transient phenomena in the course of filling of fluid inside gas well. Jalal et al.^[15] studied the flow of a kind of micro-compressed Bingham fluid, established a new analytic equation based on horizontal well pressure and derivative of pressure, and made calculation and validation under six different operating conditions. The study of fluid transient in tubing is also reflected in the drilling process. He et al.^[16] made finite element simulation study on surge pressure occurring during overflow shut-in by using ADINA software. Peng et al.^[17] discussed the surge pressure in laminar and turbulent states and established drill stem motion-based steady-state surge pressure model on the basis of slot flow model and the actual speed distribution of fluid during tripping. Chen et al.^[18] presented the approach of drawing pressure wave velocity standard chart and drafted the pressure wave velocity chart to check the pressure wave velocities of gas and liquid flows under different temperatures and pressures and water contents inside tubing.

In spite of the fruitful results achieved in water hammer studies, the impacts of valve opening coefficient and liquid holdup etc. on tubing pressure during gas well production have not been discussed in previous studies of transient flow, there are few reports about fluid transient flow at tubing build-up section. Thus, it is necessary to establish transient flow model applicable to complex gas well working condition and make separate modeling of tubing build-up section to analyze the impacts of valve opening and liquid holdup on the transient flow, in the hope to provide technical support for studying the actual flow process and gas well pressure distribution^[19,20].

The transient flow process in tubing is similar to that in ordinary pipes. One round of pressure transfer can be divided into four stages (or a cycle). Assuming the well depth is L, pressure-wave speed is c_m, wellhead flow speed is u₀ and wellhead pressure is p₀, the four stages will be (1) 0-L/c_m, (2) L/c_m-2L/c_m, (3) 2L/c_m-3L/c_m and (4) 3L/c_m-4L/c_m.

At the moment the wellhead valve is closed, the flow speed of micro-unit with the height of Δz close to the valve will turn into 0 instantly and the flow satisfies the Bernoulli’s equation of actual stream line. When a section of flow unit is defined from section A₁ to A₂ (Fig. 1), the Bernoulli’s equation of actual main stream line can be expressed as equation (1) after integration.

This equation shows the universal rule of mutual conversion between potential energy and kinetic energy, flow speed and intensity of pressure in actual main stream line. If the micro-unit fluid with height Δz is taken as the study object, the positional head and friction along tubing can be neglected, in other words, the sum of velocity head and pressure head is a constant value. If the speed falls, the pressure will rise.

Fig. 2 shows the four stages of the formation of water hammer. After the wellhead valve is closed, the pressure wave propagates to the bottom of the well at the wave velocity of c_m, and then reaches the bottom after t= L/c_m. The fluid pressure in the tubing is p₀+Δp (where Δp=Δp₁+ Δp₂+…+Δp_i). The pressure in the tubing is higher than that at the bottom of the well at this time. The fluid near the well bottom flows back to the bottom at the velocity of u_m,0, and the fluid above flows into the bottom in turn. In the range of L/c_m-2L/c_m, the pressure difference gradually disappears and returns to normal pressure condition. When t=2L/c_m, the fluid still flows to the well bottom at the speed of u_m,0 under the effect of gravity, forming negative pressure in the tubing, and a decompression wave surface which propagates to the bottom of the well at the wave speed of c_m until the pressure drops to p₀-Δp. When t=3L/c_m, the bottom hole pressure is higher than the pressure in tubing, the fluid moves upward at the velocity of u_m,0, and the bottom hole pressure returns to p₀. This imbalance propagates to the wellhead at wave velocity of c_m in turn and reaches the wellhead when t=4L/c_m. To this point, the water hammer pressure produced by the well shut-in completed a cycle of propagation.

The classic transient flow model is for horizontal pipeline single phase flow and most of them only consider liquid phase. In view of complicate down-hole conditions, multiple factors which haven’t been considered in classic models should be taken into account. For instance, the classic models cannot accurately describe the impacts of gravity in vertical or nearly vertical tubing; in build-up section, it is necessary to consider the counter-elastic expansion force produced when fluid impacts with elbow and local friction caused in such section; the transmission media in the tubing is multi-phase fluid^[21,22], which is quite different in parameters like pressure wave speed and fluid speed from single liquid phase. Either for drilling or production, the single-phase flow models are unable to satisfy the computing needs.

The multi-phase flow involves interactions between fluid media in the tubing, and complicated stream pattern conversion, transfer processes of inter-phase mass, momentum and energy. The physical properties of the media are dependent on the volumetric ratio and density ratio of the media. Therefore, simplification is needed during modeling, and the following assumptions are made: (1) The fluid media and tubing are linear elastomer with constant elasticity modulus and the tubing has no longitudinal elastic deformation. (2) The fluid media in the tubing flow is one-dimension and the multi- phase fluid on the same section of the tubing is even in distribution and same in flow speed. (3) The tubing constant flow formula can still be used to calculate the friction loss in transient flow.

The basis for studying the transient flow process in tubing is the transient flow basic differential equation system which consists of continuity equation and equation of motion. As shown in Fig. 3, the fluid in a section of tubing with finite length (radian) each is taken from the vertical well section and build-up section respectively as the object of study, which has a length (radian) of dz, fluid inlet node of j and fluid outlet node of j+1.

Since the vertical well section and build-up section will not affect the continuous medium hypothesis of fluid, the same continuity equation can be used for two kinds of well section conditions. According to the conservation of mass of finite element control volume within the time dt:

Which is simplified as:

When considering the multi-phase flow factor, Darcy formula isn’t applicable, and the frictional resistance needs to be calculated with Fanning equation:

Correspondingly, the water hammer pressure can be expressed as:

c_m hereby stands for the pressure wave speed under multi-phase flow condition (solid state isn’t considered here), and is defined as:

With given tubing material and size, if the change of elastic modulus of gas-liquid fluid in the tubing is ignored, the change of liquid holdup will directly affect the pressure wave speed and in turn the cross-sectional water hammer pressure. According to the equation above, decrease of liquid holdup will lead to increase of pressure wave speed, thereby resulting in the rise of water hammer pressure.

Thus, the continuity equation describing the gas well transient flow is:

Similarly, the fluid in a section of tubing with finite length (radian) from the vertical well section and build-up section each are taken respectively as the object of study, which have a length (radian) of dz, fluid inlet node of j and fluid outlet node of j+1, as shown in Fig. 4.

When fluid mixture flows through vertical well section, according to the stress analysis above, the resultant force on this section is:

When fluid mixture flows through the build-up section of tubing, the fluid will impact the elbow section due to inertia. If the tubing is made of linearly elastic material, an elastic expansion force will be given back to the fluid:

The resultant force on the build-up section is:

Thereby, the equation of motion meeting the gas well transient flow process can be obtained; so the equations for Newton's second law of motion for the vertical well section and build-up section can be written.

Equation for vertical well section:

Equation for build-up section:

Therefore, the gas well transient flow process can be calculated in three separate well sections, as shown in Fig. 5.

On the basis of classic gas well transient flow model, the effects of gravity of vertical tubing, multi-phase flow and build-up section on fluid are considered. This model can better reflect the complicated down-hole environment and give a precise mathematical description of the tubing fluid transient motion.

To seek solution of transient flow mathematical model, analytic method and numerical method are commonly used. The transient flow model is a pair of quasilinear hyperbolic partial differential equations, and it is very difficult to obtain its analytical solution. In contrast, the numerical method does not have high requirement on the equation and does not need it to be further simplified. By replacing model experiment with computer simulation, the numerical method makes the study of engineering issues easier. Therefore, the most common characteristic method in the numerical methods is used in this study.

The vertical well section and build-up section need to be described separately. Taking the well section I. (vertical) as an example, the eigenvalue equation system is set up and solved; the solution steps for the build-up section and vertical well section II (inclined) are similar to those of the vertical well section I.

By introducing the undetermined coefficient ω, the basic differential equation of simplified transient flow for the vertical well section can be linearized:

The equation linearized is converted into ordinary differential equation and appropriate ω is selected to obtain the characteristic equations to describe the transient flow process. The downwave characteristic line (C⁺) and the backward wave characteristic line (C^-) are as followings:

The tubing needs to be divided into three sections to discretely analyze the well depth and time with finite difference discrete method, i.e., vertical well section I, build-up section and vertical well section II, each with a length of Δz. The time is equally divided into M sections at the step of Δt. Taking the vertical well section I as an example, the feature mesh L-t planar graph is made as shown in Fig. 6.

If the well depths of points P₁ and P₂ are L₁ and L₂, the cross-sectional pressures on the corresponding tubing nodes are $p_{{{P}_{1}},i}^{{}}$ and $p_{{{P}_{2}},i}^{{}}$ at time i,, and the flow speeds are ${{u}_{\text{m}{{P}_{1}}\text{,}i}}$ and ${{u}_{\text{m}{{P}_{2}}\text{,}i}}$, then the pressure and flow speed of point P at time i+1 can be calculated with equations (13) and (14). By repeating these steps, the pressures and flow speeds at all nodes of the tubing at subsequent time nodes can be obtained. The solution steps of the build-up section and vertical well section II are same as above.

By differentiating the equations (13) and (14) along the characteristic lines C⁺ and C^-, pressure and flow speed at cross section of any node j and the time i can be obtained:

In the solution of gas well surge pressure and flow speed, the initial values of p and u_m of each node on the tubing are already known, so p and u_m of each node at time Δt are calculated first, then the p and u_m of the next moment, until the moment needs to be calculated is reached.

Manual calculation is not recommended due to huge amount of calculation of joint discrete equations of vertical well section and build-up section; instead, computer program is used for the gas well transient flow simulated calculation. The example well is divided into three sections for the computation: vertical well section I is divided into 500 segments; the build-up section into 150 segments and the vertical well section II into 350 segments. Assuming the wellhead valve opening of K_v satisfies the following function:

where, β is valve opening coefficient, when β=1.0, the valve is linearly close.

When converting the gas production into flow rate in tubing under standard condition, in consideration of the natural gas compressibility, the actual gas state equation will be used for the conversion. The natural gas compressibility factor Z is a very important parameter. In reference to the formula SGERG-88 prescribed in national standard GB/T 17747^[23], gross calorific value, relative density and CO₂ content are taken as the input variables. According to the standard GB/T 17747, the formula for calculating natural gas compressibility factor with physical property values is as below:

The natural gas compressibility factor of micro-element in tubing can be calculated with aforementioned formulas. In formula (18), B and C stand for the second and third virial coefficients of natural gas respectively, which can be calculated according to the standard GB/T 17747.

Taking Well XX6 as an example, the wellhead pressure fluctuation after valve closing was simulated to analyze the impacts of valve opening coefficient, sectional liquid holdup and valve closing time on the wellhead surge pressure. The basic data of gas well is shown in Table 1.

By using the multi-phase flow gas well transient flow analysis model presented in this paper, wellhead pressure fluctuation of Well XX6 was calculated (Fig. 7). Fig. 7 shows that at the moment of valve closing, the wellhead fluid flow speed changes to zero suddenly, and the pressure increases sharply and reaches 32.65 MPa at 3.28 s, after which it drops in different gradient and hits the lowest point of 26.48 MPa at 9.06 s, then the pressure wave starts another round of transmission. Since the natural gas in tubing is easy to be compressed, and under the action of frictional resistance, surge pressure decays quickly and after a short while, it reaches the balanced pressure.

The valve opening coefficient (β) is an important parameter to simulate transient flow process. It directly influences the peak value and magnitude of the surge pressure. Different β values were taken to simulate the wellhead pressure variations in the process of transient flow (Fig. 8a) and the part in red rectangle in Fig. 8a is amplified (Fig. 8b). According to Fig. 8a, when β=1.5, wellhead surge pressure reaches its peak value of 33.07 MPa, and the greater the β, the higher the wellhead surge pressure peak value will be. As shown in Fig. 8b, the higher the β, the smaller the wellhead surge pressure variation, and the smaller the pressure change range during the pressure transmission will be; the smaller the β (taking β=0.5 for example), the greater the change of wellhead surge pressure and the larger the pressure change range will be, the peak value hit 32.79 MPa, higher than the surge pressure reaching wellhead earliest. Although surge pressure is greater when β value is larger, lower β value would cause repeated fluid pressure disturbance in tubing, resulting in fluid-solid coupling vibration. This is a very tricky issue at present. Therefore, a higher valve opening coefficient shall be adopted on the basis of not exceeding the maximum shut-in pressure.

Sectional liquid holdup (H_L) is another important parameter in multi-phase flow study. It refers to the ratio of flow cross section area of liquid phase to the total flow cross section area in the course of gas-liquid flow. By selecting the fluid layer close to the wellhead as the object of study, the pressure fluctuations of this section under different H_L conditions were calculated (Fig. 9a), and the part in the red rectangle in Fig. 9a is amplified (Fig. 9b). The calculation results show the higher the H_L, the higher the pressure wave transmission speed and the shorter the transmission cycle will be. For example, when H_L increases from 0.03 to 0.10, pressure inflection points appear at 3.47, 3.28, 3.08 and 2.79 s respectively (See the red rectangle part in Fig. 9b).

A higher H_L corresponds to a greater magnitude of surge pressure variation, which will lead to pressure disturbance in tubing. Meanwhile, a higher H_L also corresponds to a higher pressure value; and the change of the H_L will lead to an obvious change of surge pressure. In actual production, the liquid holdup in tubing should be precisely calculated and an appropriate H_L value should be set by tuning production parameters to control the range and magnitude of surge pressure and mitigate water hammering pressure.

Different valve closing time also greatly affects the wellhead surge pressure peak value. Fig. 10a shows the simulation results of wellhead surge pressure at three different valve closing time, the three dotted line represent the pressure peak values at three valve closing time. It can be seen that the increase of valve closing time will make the wellhead maximum pressure value decrease and the occurring moment of peak value delay. Fig. 10b shows the wellhead fluctuating pressure peak values and the corresponding moments at the valve closing time of 0.5-10.0 s. It can be seen that the pressure wave speed after quick shut-in is faster than that after slow shut-in.

Well XX6 is a shallow gas storage well with high production. A larger pressure wave shock and higher wellhead pressure will be generated when the valve is closed instantaneously. By increasing the valve closing time, the water hammering pressure would reduce and reach a balanced pressure faster. In deep wells or ultra-deep wells, due to the greater frictional resistance and higher compressibility of the gas itself, the pressure wave decays even faster and poses less impact on the wellhead.

By simulating gas well shut-in under transient state with gas well transient flow model, the valve opening coefficient and closing time can be optimized to minimize damages to wellhead equipment and tubing caused by water hammering and ensure the integrity of wellbore.

The larger the opening coefficient of the wellhead valve when it is closed, the higher the wellhead pressure peak value, the smoother the surge pressure variation, and the less prominent the pressure change range will be. Without exceeding the maximum shut-in pressure of tubing, using a higher opening coefficient can reduce the shock of pressure wave.

The higher the sectional liquid holdup, the faster the pressure wave speed and the shorter the transmission cycle, but the larger the change magnitude of surge pressure and the greater the pressure will be. In actual production, production parameters can be adjusted to get an appropriate liquid holdup to better control the range and magnitude of surge pressure and reduce the shock of water hammering.

The increase of valve closing time would make the wellhead surge pressure peak value decrease, the occurring time of the peak value delay, and the pressure change range disappear gradually. On the contrary, the shorter the valve closing time, the faster the pressure wave transmission speed will be.

The gas well transient flow model can be used to optimize valve opening coefficient and valve closing time to reduce damages to wellhead equipment and tubing caused by water hammering to ensure the wellbore integrity.

A, A₁, A₂—tubing cross sectional area, m²;

B—the second virial coefficient of natural gas, m³/kmol;

C—the third virial coefficient of natural gas, m³/kmol²;

C⁺—downwave characteristic line;

C^-—backward wave characteristic line;

c_m—pressure wave speed , m/s;

D—tubing outer diameter, m;

e—pipe wall thickness, m;

E_g—gas elastic modulus, Pa;

E_L—liquid elastic modulus, Pa;

E_p—pipe wall elastic modulus, Pa;

f—Fanning friction coefficient, dimensionless;

F_be—resultant of external forces on the tubing of build-up section, Pa;

F_c—elastic expansion force on the tubing of build-up section, Pa;

F_s—average lateral force caused by tubing expansion, Pa;

F_ve—resultant of external forces on the tubing of straight section, Pa;

g—acceleration of gravity, m/s²;

G—tubing gravity, N;

h_w,1,2—streamwise fiction from water head position 1 to 2, m;

H_L—sectional liquid holdup, %;

i—serial number of time node;

j—serial number of space node;

K_v—valve opening, dimensionless;

L—gas well depth, m;

M—time slot number;

p—pressure inside the tubing, Pa;

p₀—wellhead pressure, Pa;

p₁, p₂—pressures at water head position 1 and 2, Pa;

${{p}_{{{P}_{\text{1}}},i}}$, ${{p}_{{{P}_{\text{2}}},i}}$—sectional pressures at P₁ and P₂ tubing nodes at time i, Pa;

P, P₁, P₂—dots on feature mesh L-t;

p_in—the inlet flow pressure of micro-element, Pa;

p_out—the outlet flow pressure of micro-element, Pa;

R—universal gas constant, J/(mol·K);

t—time, s;

T—gas temperature, K;

T_s—valve closing time, s;

u₀—wellhead fluid flow speed, m/s;

u_m—gas-liquid mixed phase flow speed, m/s;

u_m,1, u_m,2—flow speed of gas-liquid mixed phase at water head position 1 and 2, m/s;

${{u}_{\text{m}{{P}_{1}}\text{,}i}}$, ${{u}_{\text{m}{{P}_{2}}\text{,}i}}$—sectional speed at P₁ and P₂ tubing nodes at time i, m/s;

z, z₁, z₂—water head position, m;

Z—natural gas compressibility factor, dimensionless;

α₁, α₂—kinetic energy correction factors of water head position 1 and 2, dimensionless;

β—valve opening coefficient, dimensionless;

Δp—tubing pressure increment, Pa;

Δt—time step, s;

Δz—well depth step, m;

ρ_m—density of mixed gas-liquid phase, kg/m³;

ρ_mol—molar density, kmol/m³;

τ—tubing wall shearing stress, Pa;

φ, φ₁, φ₂—included angles between pipe axle and horizontal position, (°);

Ψ, Ψ₁, Ψ₂—included angle between pipe side wall and axle, (°);

ω—fitting undetermined coefficient, m²•s/kg.