Abstract

Translate

Simulating surface conditions by solving the wave equation of a sucker-rod string is the theoretical basis of a sucker-rod pumping system. To overcome the shortcomings of the conventional finite difference method and analytical solution, this work describes a novel hybrid method that combines the analytical solution with the finite difference method. In this method, an analytical solution of the tapered rod wave equation with a recursive matrix form based on the Fourier series is proposed, a unified pumping condition model is established, a modified finite difference method is given, a hybrid strategy is established, and a convergence calculation method is proposed. Based on two different types of oil wells, the analytical solutions are verified by comparing different methods. The hybrid method is verified by using the finite difference method simulated data and measured oil data. The pumping speed sensitivity and convergence of the hybrid method are studied. The results show that the proposed analytical solution has high accuracy, with a maximum relative error relative to that of the classical finite difference method of 0.062%. The proposed hybrid method has a high simulation accuracy, with a maximum relative area error relative to that of the finite difference method of 0.09% and a maximum relative area error relative to measured data of 1.89%. Even at higher pumping speeds, the hybrid method still has accuracy. The hybrid method in this paper is convergent. The introduction of the finite difference method allows the hybrid method to more easily converge. The novelty of this work is that it combines the advantages of the finite difference method and the analytical solution, and it provides a convergence calculation method to provide guidance for its application. The hybrid method presented in this paper provides an alternative scheme for predicting the behavior of sucker-rod pumping systems and a new approach for solving wave equations with complex boundary conditions.

Full text

Turn on search term navigation

Translate

1. Introduction

Petroleum is a strategic resource that plays a key role in the industrial field and manufacturing industry and still accounts for a relatively large proportion of global energy consumption [1]. In the petroleum industry, the sucker-rod pumping system has the merits of low comprehensive costs, simple equipment, and convenient operation [2] and has the highest number of installations worldwide [3].

A typical sucker-rod pumping system consists of a surface unit, a sucker-rod string, a tubing string, and a subsurface pump, which has a standing valve at the bottom of the well and a traveling valve attached to a rod [4,5]. A sucker-rod pumping system is generally unstable and can suffer various failures due to operation in open air and long-term operation, so monitoring its operating conditions is very important for oil production [6]. However, a pumping system’s operating conditions are difficult to measure directly because it operates in a small-diameter oil tube thousands of meters underground. In production practice, the pumping condition is usually identified by analyzing the surface dynamometer card [7], i.e., a closed curve of polished rod displacement and load, which is called diagnostic technology.

Many advanced methods have been applied in diagnostic technology based on fault dynamometer cards, such as the continuous hidden Markov model [8], supervised dictionary-based transfer subspace learning [9], optimized density peak clustering [10], oil-Net 1-D/2-D identification models from time-series and computer vision perspectives [11], four-segment time-frequency signature matrices, and deep learning [12].

In a real oil well, it is impossible to have every type of surface dynamometer card. Computer simulation of fault surface dynamometer cards is an economical and effective method. In this method, the card is simulated based on the solutions of the sucker-rod string wave equation with the description of the operation of the sucker-rod pump, which is called prediction technology [13]. Moreover, the behavior prediction of a sucker-rod pumping system is the theoretical basis for system design [14], energy consumption analysis [15], production parameter optimization [16], system control [4], and so on.

In 1963, Gibbs [17] established a one-dimensional wave equation to describe the vibration of a sucker-rod string and proposed the finite difference method to predict the behaviors of a sucker-rod pumping system. In the finite difference method, the upper boundary of the wave equation is the polished rod displacement, and the lower boundary is a “self-sensing” pumping condition model. It is a mixed function of the pump displacement and the load determined by the time of valve opening and closing, which permits automatic “sensing” of the valve operating times by computer tests.

Gibbs’s work was pioneering [18] and has since been used as a basis for research [19]. As a standard solution method [4], one of the merits of the finite difference method is that it can easily simulate the pump card (i.e., the downhole dynamometer card, which is a closed curve of pump displacement and load). However, the finite difference method should satisfy the stability criterion [17]. This results in a long computer simulation time, especially for tapered rod strings [20].

In addition to the finite difference method, there are other methods for solving the sucker-rod wave equation in prediction technology, such as the traveling wave method [21], the mode superposition method [22,23], the spring-mass-damper method (i.e., abbreviated as the SMD method, the sucker-rod string is dispersed into spring-mass-damper systems) [24,25], and the analytical solution.

This paper focuses on the analytical solution due to its advantages of being able to calculate the displacement and force at any position on the sucker-rod string [26] and not being limited by stability conditions [20]. As early as 1966 and 1967, Gibbs and Neely [27] and Gibbs [28] proposed an analytical solution using complex theory and separation of variables in diagnostic technology. The boundary conditions of the wave equation are the polished rod displacement and load, which are known and can be expanded into a Fourier series, so the pump displacement and load can be easily calculated [29].

Scientists have tried to apply it to prediction techniques but have no effective scheme. Li and Li [30] proposed an analytical solution using complex theory and separation of variables in prediction and noted that it required a periodic trail pump load-time function that should be modified according to the given conditions. However, modified details of the pumping condition model are lacking. Chen et al. [31] provided a unified expression for diagnostic and predictive analysis based on the analytical solution using complex theory and separation of variables, taking the displacement and force at any point of the sucker-rod string as known conditions. However, its boundary conditions are different from those in real predictive analyses [13,17]. Yin et al. [20] proposed a matrix solution based on the analytical solution. However, the pump load is a periodical pump load-time function determined by polished rod displacement, which is suitable only for the case of low pumping speeds. The analytical solution is essentially a frequency domain approach that requires a periodic pump load-time function, making it difficult to address the traditional “self-sensing” pumping condition model. Therefore, it is necessary to study how to introduce the finite difference method into the analytical wave solution of the prediction method.

This work proposes a hybrid method based on an analytical solution and the finite difference method. First, an analytical solution of the wave equation of the tapered sucker- rod based on a recursive matrix form is presented. Then, as part of the preparation for the hybrid method, a unified pumping condition model is established. As another part of the preparation for the hybrid method, a modified finite difference method for a single rod is introduced. Then, the hybrid method of fusing the analytical solution and the finite difference method is proposed, and the convergence of the hybrid method is given. Finally, the hybrid method is verified by simulated and measured cards based on oil field data.

The remainder of the paper is organized as follows: Section 2 describes the theory and algorithm of this work. Section 3 presents the validation and discussion. Section 4 presents the conclusion of the paper.

2. Theory and Algorithm

2.1. Fundamental Assumptions

A real well pumping system is complex, especially when considering the flow law of borehole fluids. To establish the model conveniently, the following basic assumptions are made: (1) The lateral vibration of the sucker-rod string is ignored. (2) The discharge pressure, intake pressure, and temperature at the pump are constant. (3) The valves are opened and closed instantaneously. (4) The hydraulic loss of the valves and the liquid friction of the pump plunger are ignored. (5) The gas-to-liquid ratio does not change at the same pump pressure and temperature. Only the compressibility of the gas phase is considered. Gas compression and expansion are isothermal polytropic processes. Notably, most of the assumptions are conventional assumptions [16,32,33], which do not deviate much from field applications.

2.2. Analytical Solution of the Wave Equation Based on a Recursive Matrix Form

2.2.1. One-Dimensional Wave Equation

The longitudinal vibration of the single rod string can be described by the one-dimensional wave equation given by Gibbs in a vertical well [17]. It is assumed that the total string consists of the K stage taper section and the wave equation of the i-th taper section is [20]:

(1) $\frac{\partial^{2} u_{i}}{\partial t^{2}} - c_{i}^{2} \frac{\partial^{2} u_{i}}{\partial x^{2}} + v_{i} \frac{\partial u_{i}}{\partial t} = 0 l - L_{i} < x < l$

where

c_{i} = \sqrt{E_{r i} g / ρ_{i}}

l = \sum_{k = 1}^{i} L_{k}

E_{r i}

is the Young’s modulus of the rod material,

g

is the gravitational acceleration,

ρ_{i}

is the rod density,

L_{k}

is the length of the k-stage taper section of the rod, and

v_{i}

is the viscous damping factor.

In prediction technology, the boundary conditions and continuity conditions are

(2) $\begin{array}{l} u (0, t) = - u_{a} (t), \\ E_{r K} A_{r K} \frac{\partial u (x, t)}{\partial x} |_{x = L} = P_{p} (t), \\ u_{i} (l, t) = u_{i + 1} (l, t), \\ E_{r i} A_{r i} \frac{\partial u_{i} (x, t)}{\partial x} |_{x = l} = E_{r i + 1} A_{r i + 1} \frac{\partial u_{i + 1} (x, t)}{\partial x} |_{x = l} \end{array}$

where

L = \sum_{k = 1}^{K} L_{k}

and

P_{p} (t) = A_{p} [p_{d} - p (t)]

A_{p}

is the plunger area,

p_{d}

is the discharge pressure,

p (t)

is the pump pressure, and

u_{a} (t)

is the polished rod displacement.

Equations (1) and (2) show that when $u_{a} (t)$ and $P_{p} (t)$ or $p (t)$ are given, the rod displacement function $u (x, t)$ can be calculated after solving Equation (1), then the dynamic properties of the rod can be analyzed.

2.2.2. Block Matrix Expression of the Analytical Solution

By separating variables and assuming complex periodic solutions $T (t) = e^{i n ω t}$ , Gibbs and Neely [27] and Gibbs [28] proposed an analytical solution in diagnostic technology. Equation (1) of the wave equation of the i-th taper section is defined as follows:

Definition 1

([27]). The generic solution to Equation (1) is as follows:

(3) $\begin{array}{l} u_{i} (x, t) = ξ_{i} + η_{i} x + \sum_{n = 1}^{N} [Q_{i n} (x) \cos n ω t + P_{i n} (x) \sin n ω t], \\ D_{i} (x, t) = E_{r i} A_{r i} η_{i} + E_{r i} A_{r i} \sum_{n = 1}^{N} [Q_{i n} ’ (x) \cos n ω t + P_{i n} ’ (x) \sin n ω t] \end{array}$

where

(4) $\begin{array}{l} Q_{i n} (x) = ε_{i n} c h c (x) + π_{i n} s h s (x) + κ_{i n} c h s (x) + μ_{i n} s h c (x), \\ P_{i n} (x) = - ε_{i n} s h s (x) + π_{i n} c h c (x) + κ_{i n} s h c (x) - μ_{i n} c h s (x), \\ Q_{i n} (x) = \frac{\partial Q_{i n} (x)}{\partial x}, P_{i n} (x) = \frac{\partial P_{i n} (x)}{\partial x}, \\ c h c (x) = \cosh β_{i n} x \cos α_{i n} x, s h s (x) = \sinh β_{i n} x \sin α_{i n} x, \\ c h s (x) = \cosh β_{i n} x \sin α_{i n} x, s h c (x) = \sinh β_{i n} x \cos α_{i n} x, \\ α_{i n} = \frac{n ω}{c_{i} \sqrt{2}} \sqrt{1 + \sqrt{1 + {(\frac{v_{i}}{n ω})}^{2}}}, β_{i n} = \frac{n ω}{c_{i} \sqrt{2}} \sqrt{- 1 + \sqrt{1 + {(\frac{v_{i}}{n ω})}^{2}}} \end{array}$

where

ω = 2 π / T

is the angular frequency of pumping and

T

is the pumping cycle period.

Definition 2.

$Q_{i n} (x)$ , $P_{i n} (x)$ , $Q_{i n}^{'} (x)$ , and $P_{i n}^{'} (x)$ in Equation (3) are expressed by a block matrix as follows:

(5) ${[\begin{matrix} Q_{i n} (x) & P_{i n} (x) & Q_{i n}^{'} (x) & P_{i n}^{'} (x) \end{matrix}]}^{T} = {[\begin{matrix} M_{i, 1}^{x} & F_{i, 1}^{x} \end{matrix}]}^{T} {[ε π_{n}]}_{i} + {[\begin{matrix} M_{i, 2}^{x} & F_{i, 2}^{x} \end{matrix}]}^{T} {[κ μ_{n}]}_{i}$

where

(6) $\begin{array}{l} {[ε π_{n}]}_{i} = {[\begin{matrix} ε_{i n} & π_{i n} \end{matrix}]}^{T}, {[κ μ_{n}]}_{i} = {[\begin{matrix} κ_{i n} & μ_{i n} \end{matrix}]}^{T}, \\ M_{i, 1}^{x} = {[\begin{matrix} c h c (x) & s h s (x) \\ - s h s (x) & c h c (x) \end{matrix}]}_{i}^{x}, M_{i, 2}^{x} = {[\begin{matrix} c h s (x) & s h c (x) \\ s h c (x) & - c h s (x) \end{matrix}]}_{i}^{x}, \\ F_{i, 1}^{x} = \frac{\partial M_{i, 1}^{x}}{\partial x}, F_{i, 2}^{x} = \frac{\partial M_{i, 2}^{x}}{\partial x} \end{array}$

2.2.3. Recursive Matrix Solving Method

According to Definitions 1–2, the solution of the definite solution problem, i.e., Equations (1) and (2), is transformed into the solution of the coefficients ${[ε π_{n}]}_{i}$ and ${[κ μ_{n}]}_{i}$ .

According to the displacement expression and the load expression of Equation (3), the continuity conditions of Equation (2) can be rewritten as follows:

(7) $[\begin{matrix} M_{i, 1}^{l} & M_{i, 2}^{l} \\ E_{r i} A_{i} F_{i, 1}^{l} & E_{r i} A_{r i} F_{i, 2}^{l} \end{matrix}] {[\begin{matrix} π ε_{n} \\ κ μ_{n} \end{matrix}]}_{i} = [\begin{matrix} M_{i + 1, 1}^{l} & M_{i + 1, 2}^{l} \\ E_{r i + 1} A_{r i + 1} F_{i + 1, 1}^{l} & E_{r i + 1} A_{r i + 1} F_{i + 1, 2}^{l} \end{matrix}] {[\begin{matrix} π ε_{n} \\ κ μ_{n} \end{matrix}]}_{i + 1}$

Resolving Equation (7) yields the recursive matrix

(8) ${[\begin{matrix} π ε_{n} \\ κ μ_{n} \end{matrix}]}_{i} = {[\begin{matrix} M_{i, 1}^{l} & M_{i, 2}^{l} \\ E_{r i} A_{i} F_{i, 1}^{l} & E_{r i} A_{r i} F_{i, 2}^{l} \end{matrix}]}^{- 1} [\begin{matrix} M_{i + 1, 1}^{l} & M_{i + 1, 2}^{l} \\ E_{r i + 1} A_{r i + 1} F_{i + 1, 1}^{l} & E_{r i + 1} A_{r i + 1} F_{i + 1, 2}^{l} \end{matrix}] {[\begin{matrix} π ε_{n} \\ κ μ_{n} \end{matrix}]}_{i + 1}$

The coefficients ${[ε π_{n}]}_{K}$ and ${[κ μ_{n}]}_{K}$ are the key to the solution of the definite solution problem.

Proposition 1.

The coefficients ${[ε π_{n}]}_{K}$ and ${[κ μ_{n}]}_{K}$ of the last section of the tapered rod are determined from the rod’s characteristic matrix $[M F_{n}]$ and the boundary condition coefficient matrix $[δ v_{n}]$ and $[σ τ_{n}]$ as follows:

(9) ${[\begin{matrix} κ μ_{n} \\ π ε_{n} \end{matrix}]}_{K} = [M F_{n}] [\begin{matrix} δ v_{n} \\ σ τ_{n} \end{matrix}]$

The poof Proposition 1 is shown in Appendix A.1.

According to Equations (A1) and (A3), the coefficients ${[π ε_{n}]}_{1}$ and ${[κ μ_{n}]}_{1}$ can be resolved as follows:

(10) $\begin{array}{l} {[π ε_{n}]}_{1} = [v δ_{n}], \\ {[κ μ_{n}]}_{1} = [F_{t 1}] {[κ μ_{n}]}_{N} + [F_{t 2}] {[π ε_{n}]}_{N} \end{array}$

The other coefficients ${[\begin{matrix} κ μ_{n} \\ π ε_{n} \end{matrix}]}_{i}$ can be calculated according to the recursive matrix of Equation (8).

Thus, the pump displacement and pump load can be calculated as follows:

(11) $\begin{array}{l} u_{p} (t) = - u_{K} (L, t) = \frac{ν_{0}}{2} - \frac{σ_{0}}{2} \sum_{k = 1}^{K} \frac{L_{k}}{E_{r k} A_{r k}} - \sum_{n = 1}^{N} [Q_{K n} (L) \cos n ω t + P_{K n} (L) \sin n ω t], \\ P R L (t) = W_{f} + D_{1} (0, t) = W_{f} + \frac{σ_{0}}{2} + E_{r 1} A_{r 1} \sum_{n = 1}^{N} [Q'_{1 n} (0) \cos n ω t + P'_{1 n} (0) \sin n ω t] \end{array}$

where

(12) $\begin{array}{l} {[\begin{matrix} Q_{K n} (L) & P_{K n} (L) \end{matrix}]}^{T} = [M_{K, 1}^{L}] {[ε π_{n}]}_{K} + [M_{K, 2}^{L}] {[κ μ_{n}]}_{K}, \\ {[\begin{matrix} Q_{1 n}^{'} (0) & P_{1 n}^{'} (0) \end{matrix}]}^{T} = [F_{1, 1}^{0}] {[ε π_{n}]}_{1} + [F_{1, 2}^{0}] {[κ μ_{n}]}_{1} \end{array}$

$W_{f}$ is the weight of the sucker-rod in the fluid.

As shown in Equation (11), $u_{p} (t)$ and $P R L (t)$ can be calculated after $u_{a} (t)$ and $P_{p} (t)$ are approximated by the truncated Fourier series. Then, the key cards of this paper, i.e., the surface dynamometer card $[u_{a} (t), P R L (t)]$ and the pump card $[u_{p} (t), P_{p} (t)]$ , can be drawn.

The flowchart of the solution procedure is shown in Figure 1.

The polished rod displacement $u_{p} (t)$ is easily determined, and the next key is to determine the pump load $P_{p} (t)$ according to the operating conditions of the pump, i.e., the pump condition model.

2.3. Uniform Pumping Condition Model

Under common pumping conditions, the model in this paper considers only the ideal pumping conditions and the gas interference pumping conditions.

When the pumping speed is low, the vibration of the sucker-rod string in Equation (1) can be neglected, and the sucker-rod string is in a state of elastic deformation and is equivalent to a spring. The fundamental force balance equation follows according to Equation (1) to Equation (2):

(13) $u_{a} (t) - u_{p g} (t) - A_{p} k_{t} [p_{d} - p (t)] = A_{p} k_{r} [p_{d} - p (t)]$

where

u_{p g} (t)

is the pump displacement caused by gas interference and

A_{p} k_{t} [p_{d} - p (t)]

is the pump displacement caused by the elastic deformation of the tubing string.

k_{r} = \sum_{i = 1}^{K} \frac{L_{i}}{E_{r i} A_{r i}}

and

k_{t} = \sum_{i = 1}^{K} \frac{L_{t i}}{E_{t i} A_{t i}}

L_{t i}

and

A_{t i}

are the length and the cross-sectional area of the i-th tubing string, respectively.

E_{t i}

is the Young’s modulus of the i-th tubing string.

According to the gas expansion and gas compression [23], the pump displacement can be rewritten and differentiated. Differentiating Equation (13) yields

(14) $v_{a} (t) - v_{p g} (t) \frac{d p (t)}{d t} = - A_{p} k_{e} \frac{d p (t)}{d t}$

where

$\begin{array}{l} k_{e} = k_{r} + δ_{t} k_{t}, \\ v_{p g} (t) = - [r_{p} \frac{l_{o g} \sqrt[n]{p_{d}}}{n} + (1 - r_{p}) \frac{l_{g} \sqrt[n]{p_{s}}}{n}] {[p (t)]}^{- \frac{1}{n} - 1}, \\ l_{o g} = R_{s} l_{0}, n = 1.1, R_{s} = \frac{R}{1 + R}, l_{g} = R_{d} (l_{0} + l_{p}), \\ R_{d} = \frac{l_{o g} {[\frac{p (t_{m})}{p_{d}}]}^{- \frac{1}{n}}}{l_{0} - l_{o g} + l_{o g} {[\frac{p (t_{m})}{p_{d}}]}^{- \frac{1}{n}}} \end{array}$

$R$ is the gas-to-liquid ratio at the pump suction. In the loading portion, $r_{p} = 1$ , while in the unloading portion, $r_{p} = 0$ . $δ_{t} = 0$ indicates that the tubing string is anchored, while $δ_{t} = 1$ indicates that the tubing string is unanchored. $p_{s}$ is the intake pressure, and $t_{m}$ is the downstroke start time. $l_{0}$ is the clearance column length of the pump and $l_{p}$ is the length of the pump stroke.

To calculate the pump pressure, the forward difference is used, and the pump pressure is calculated in the loading and unloading portions according to Equation (14) as follows:

(15) $p (t_{j}) = \frac{v_{p g} (t_{j - 1}) p (t_{j - 1}) + A_{p} k_{e} p (t_{j - 1}) - v_{a} (t_{j}) Δ t}{v_{p g} (t_{j - 1}) + A_{p} k_{e}}, j = 1, 2, 3, \dots, J$

where

Δ t

is the time increment and

T = (J - 1) Δ t

Obviously, $p (t) = f (δ_{t}, R_{s})$ . The model naturally includes the conditions of faulty pumping and tubing anchoring.

As indicated by Equation (15), the pump pressure is a function of the polished rod displacement, and it only applies to static cases. Therefore, it can be called a static model.

For the dynamic model, the pump pressure is

(16) $p (t_{j}) = \frac{v_{p g} (t_{j - 1}) p (t_{j - 1}) + A_{p} k_{t} p (t_{j - 1}) - v_{p} (t_{j}) Δ t}{v_{p g} (t_{j - 1}) + A_{p} k_{t}}, j = 1, 2, 3, \dots, J$

It is worth noting that the model can be easily extended to include other failure conditions, such as leak conditions, including the clearance leakage of the pump plunger and pump barrel, the leakage of the traveling valve, and the leakage of the standing valve. A more comprehensive model is under investigation.

The initial value of the pump pressure calculation expression is

(17) $\begin{array}{l} v_{a} (0) = 0, \\ v_{p} (0) = 0, \\ p (0) = p_{d} \end{array}$

The judgment condition is

(18) $\begin{array}{l} p (t) \leq p_{s} \to p (t) = p_{s}, \\ p (t) \geq p_{s} \to p (t) = p_{d} \end{array}$

As shown in Equation (16), the pump pressure $p (t)$ is determined by the pump speed $v_{p} (t)$ , while Equation (18) shows that it is also determined by itself, so it is a kind of “self-sensing” model. It is difficult to simulate the behavior of the analytical solution because a periodical pump load-time function is needed for Fourier series expansion. However, finite difference methods can easily handle this challenge.

2.4. Finite Difference Method

In this section, the classical finite difference method for a single rod is employed in the next hybrid method, while the finite difference method for a tapered rod string is used for comparison with the hybrid method.

The classical finite difference solution equation of a single rod is adapted below [13,18]:

(19) $u_{m, j + 1} = γ_{s 1} γ_{s 2} (u_{m + 1, j} + u_{m - 1, j}) + γ_{s 1} γ_{s 3} u_{m, j} - γ_{s 1} u_{m, j - 1}; m = 2, 3, \dots, M; j = 2, 3, \dots, J$

where

(20) $\begin{array}{l} γ_{s 1} = \frac{1}{1 + v Δ t}, \\ γ_{s 2} = {(\frac{Δ t c}{Δ x})}^{2}, \\ γ_{s 3} = (γ_{s 1}^{- 1} - 2 γ_{s 2} + 1) \end{array}$

To apply the finite difference method in the hybrid method, the upper boundary condition of the finite difference solution should be modified; that is, the displacement $u_{u} (t)$ and the load $P_{u} (t)$ should be modified. Here, they are called the middle displacement and middle load, respectively. Considering Equation (19) and the finite difference form of the pump load [17], the following modified expression is obtained:

(21) $\begin{array}{l} u_{1, j} = - u_{u, j}, \\ u_{2, j + 1} = γ_{s 1} (γ_{s 2} + 0.25 γ_{s 3}) u_{3, j} + γ_{s 1} (γ_{s 2} + 0.75 γ_{s 3}) u_{1, j} - γ_{s 1} u_{2, j - 1} + \frac{γ_{s 1} γ_{s 3} Δ x}{2 E_{r K} A_{r K}} P_{u, j} \end{array}$

The lower boundary condition is [17]

(22) $(1 - δ_{t}) u (L, t) + δ_{t} \frac{\partial u (L, t)}{\partial x} = (1 - δ_{t}) u_{c} + δ_{t} \frac{A_{p}}{E_{r K} A_{r K}} [p_{d} - p (t)]$

According to the dynamic model of Equation (16), Equation (22) can be rewritten as

(23) $u_{M, j + 1} = (1 - δ_{t}) u_{c} + δ_{t} \frac{2 u_{M - 1, j + 1} - 0.5 u_{M - 2, j + 1} + \frac{Δ x}{E_{r K} A_{r K}} F_{j}}{1.5 + \frac{Δ x}{E_{r K} A_{r K}} A_{p} p_{r e 2, j}}$

where

(24) $\begin{array}{l} F_{j} = A_{p} (p_{d} - p_{r e 1, j} p_{r e 2, j} + p_{r e 2, j} u_{M, j}), \\ p_{r e 1, j} = v_{p g, j} p_{j} + A_{p} k_{t} p_{j}, \\ p_{r e 2, j} = \frac{1}{v_{p g, j} + A_{p} k_{t}} \end{array}$

Thus, the pump load, pump pressure, and pump displacement are as follows:

(25) $\begin{array}{l} P_{p, j + 1} = \frac{E_{r K} A_{r K}}{Δ x} (1.5 u_{M, j + 1} - 2 u_{M - 1, j + 1} + 0.5 u_{M - 2, j + 1}), \\ p_{j + 1} = p_{d} - P_{p, j + 1} / A_{p}, \\ u_{p, j + 1} = - u_{M, j + 1} \end{array}$

As shown in Equations (23)–(25), the pump load $P_{p, j + 1}$ can be calculated according to $u_{M, j}$ , so the information of the j-th time step is used to recursively infer the information of the j + 1-th time step, which can easily simulate the pump operating conditions. This strategy was first proposed by Gibbs [17].

The finite difference method for tapered rod strings was proposed by Zhang et al. [34], Takacs [13], and Yin [20].

The polished rod displacement $u_{a} (t)$ is calculated according to the equation from the kinematic analysis of the movement of different types of pumping units, i.e., the long-stroke pumping unit (Rotaflex) [35] or the beam pumping unit [36].

2.5. Hybrid Method of the Tapered Rod

2.5.1. Design of the Hybrid Method

In the design of the hybrid method, the last section of the tapered rod near the pump is discretized into at least six units, and the whole tapered rod is divided into two calculation domains, as shown in Figure 2. Domain 2 is the remaining five units of the last section of the rod. In domain 1, the analytical solution is used to solve the wave equation, while in domain 2, the modified finite difference method proposed in this paper is used. Thus, it is possible to conveniently simulate the pump card using the classical “self-sensing” pumping condition model.

The calculation scheme of the hybrid method is shown in Figure 3.

The details are listed as follows:

Setp1, in domain1, calculate the static pump loads according to Equation (15), i.e., $P_{p s} (t) = A_{p} {p_{d} - p [v_{a} (t)]}$ . According to the boundary conditions of $P_{p s} (t)$ and $u_{a} (t)$ , the middle load and displacement can be calculated using the analytical solution, i.e., $P_{u} (t) = W_{f} + D_{K} (L - 5 Δ x, t)$ , $u_{u 1} (t) = - u_{K} (L - 5 Δ x, t)$ .

Setp2, in domain 2, takes the middle load $P_{u} (t)$ and displacement $u_{u 1} (t)$ as the upper boundary conditions of the finite difference method. After discretization, the pump load $P_{p} (t)$ and displacement $u_{p} (t)$ are calculated according to the finite difference method with the pumping condition model of Equation (23).

Setp3, in domain 1, takes the new load $P_{p} (t)$ and polished rod displacement $u_{a} (t)$ as the boundary, and the new middle load $P_{u} (t)$ , displacement $u_{u 2} (t)$ and polished rod load $P R L (t)$ can be calculated according to the analytical solution.

Setp4, calculate $ε = \frac{\max | u_{u 1} (t) - u_{u 2} (t) |}{u_{u 1 m}}$ , where $u_{u 1 m}$ is the maximum value of $u_{u 1} (t)$ . If $ε > ε_{0}$ , update the displacement $u_{u 1} (t)$ and return to step 2; otherwise, output the results, i.e., $P R L (t)$ , $u_{p} (t)$ , $P_{p} (t)$ , and end the program. Thus, the surface dynamometer card, which is $[u_{a} (t), P R L (t)]$ , and the pump card, which is $[u_{p} (t), P_{p} (t)]$ , can be drawn.

The algorithm’s pseudocode is shown in Algorithm 1.

Algorithm 1: Hybrid method

Calculate:

P_{p s} (t)

according to Equation (15)Calculate:

P_{u} (t)

and u_{u 1} (t)

according to the analytical solution taking P_{p s} (t)

and u_{a} (t)

as the boundary conditionsSet: ε₀ = 0.20%For i = 1:300 Calculate:

P_{p} (t)

and u_{p} (t)

according to the finite difference method, taking P_{u} (t)

and u_{u 1} (t)

as the boundary conditions with the pumping condition model of Equation (23) Calculate:

P_{u} (t)

, u_{u 2} (t)

, and P R L (t)

according to the analytical solution, taking P_{p} (t)

and u_{a} (t)

as the boundary conditions Calculate:

ε = \max | u_{u 1} (t) - u_{u 2} (t) | / u_{u 1 m}

Update:

u_{u 1} (t) = 0.5 [u_{u 1} (t) - u_{u 2} (t)]

If ε<=ε₀ break endend Output:

u_{a} (t)

, P R L (t)

, u_{p} (t)

and P_{p} (t)

2.5.2. Convergence of the Hybrid Method

As shown in the calculation scheme of the hybrid method, there is an iterative process.

Proposition 2.

The iterative coefficient matrix $[M_{i t}]$ of the hybrid method is the product of the iterative coefficient matrix $[A_{M}]$ of the analytical solution method and the iterative coefficient matrix $[F_{M}]$ of the finite difference method and is given as follows:

(26) $[M_{i t}] = [F_{M}] [A_{M}]$

The poof Proposition 2 is shown in Appendix A.2.

In domain 1, according to Equation (3), the deviation of the middle displacement $Δ u_{u 1} (t)$ is

(27) $\begin{array}{l} Δ u_{u 1} (t) = - \sum_{n = 1}^{N} [Q_{K n} (L - 5 Δ x) \cos n ω t + P_{K n} (L - 5 Δ x) \sin n ω t], \\ [\begin{array}{l} Q_{K n} (L - 5 Δ x) \\ P_{K n} (L - 5 Δ x) \end{array}] = [\begin{matrix} M_{K, 1}^{L - 5 Δ x} & M_{K, 2}^{L - 5 Δ x} \end{matrix}] {[\begin{array}{l} ε π_{n} \\ κ μ_{n} \end{array}]}_{K} \end{array}$

According to Equation (A6), the transfer coefficient matrix of the analytical solution is given as follows:

(28) $[A_{M, n}] = [\begin{matrix} M_{K, 1}^{L - 5 Δ x} & M_{K, 2}^{L - 5 Δ x} \end{matrix}] [M F_{n}]$

In domain 2, the finite difference method has the famous stability criterion, i.e., $\frac{Δ t c}{Δ x} \leq 1$ [17,18]. Combined with Equation (19), the transfer coefficient matrix of the finite difference method is set to $[F_{M}] = γ_{s 1} γ_{s 2}$ .

Thus, according to Equation (26), the iterative matrix can be given as follows:

(29) $[M_{i t, n}] = \frac{k_{r}}{k_{e}} \frac{σ_{0}}{2 E_{r K} A_{r K}} γ_{s 1} γ_{s 2} [\begin{matrix} M_{K, 1}^{L - 5 Δ x} & M_{K, 2}^{L - 5 Δ x} \end{matrix}] [M F_{n}]$

where,

σ_{0} / 2 E_{r K} A_{r K}

adjusts the order of magnitude and

k_{r} / k_{e}

is the effect coefficient of the anchored state of the tubing [37].

The spectral radius of the iteration matrix, i.e., $ρ (M_{i t, n})$ , can be calculated. Only when the spectral radius of the iteration matrix is less than 1 will the iteration converge.

3. Validation and Discussions

To ensure the credibility of the hybrid method, measured oilfield data and published data from a different oil well are used. The measured data from well 1 (number X2) are from the Shengli oilfield of the Oilfield Company of China. The published measured data and simulated data, which consider hydraulic loss and Newtonian fluid leakage, are from well 2, as proposed by Xing [38]. The basic parameters of the wells are listed in Table 1.

In this paper, a comparison method is used to validate the proposed method. The validation is divided into two parts: one is the validation of the analytical solution compared with the finite difference method, and the SMD method is employed. The finite difference method is a well-known and proven solution for the predictive analysis of sucker-rod pumping systems [37]. Samuel and Anatolii [39] argued that the SMD method is more applicable for simulating the behaviors of sucker-rod pumping systems. The SMD method can also be extended to failure analysis and prevention of pipe strings in horizontal wellbores [40]. The details of the SMD method are given in [24,25].

The other is the validation of the hybrid method. The surface dynamometer card and the pump card are simulated using the hybrid method and the classical finite difference method according to the same oil well parameters. The results of the finite difference method are taken as the standard results and compared. At the same time, the simulated surface dynamometer card is also compared with the measured card.

The computing platform used in this paper is the MATLAB R2021a. In the finite difference method and the SMD method, the length increment is $L_{i m} / 6$ , and the time increment is $0.95 Δ x / c_{i m}$ , where $L_{i n}$ is the minimum length and $c_{i m}$ is the maximum sound velocity of the tapered rod. In the analytical solution, the Fourier series number $N$ is 200, so its time increment is $60 / (400 n_{p})$ , where $n_{p}$ is the pumping speed. Since the time increments of the two methods are different, the interpolation function of MATLAB, which is “interp1”, is used to unify the vector lengths of the two methods. In well 2, there is a gas interference pumping condition with the parameter $R = 0.2$ and $l_{0} = 0.1 l_{p}$ . In the hybrid method and finite difference method, $ε_{0} = 0.20 %$ .

3.1. Validation of Analytical Solution by the Different Methods

To validate the proposed analytical solution, the upper boundary and the lower boundary are given. The upper boundary is the polished rod displacement, and the lower boundary is the static pump load calculated according to Equation (15).

Validation 1: For well 1, in the SMD method, the numbers of elements are 15, 51, and 6. The mass $m_{n}$ , elastic modulus $k_{n}$ and friction coefficient $b_{n}$ of the i-th element are [99.9305 77.3510 123.2582] kg, [3.9749 3.0796 5.0709] MPa·m, and [64.5298 49.5274 82.4773] Pa·s·m, respectively. The given upper boundary and the lower boundary are shown in Figure 4a,b, and the pump displacements calculated by the analytical solution, the finite difference method, and the SMD method are shown in Figure 4c. As shown in Figure 4c, the calculated results of the three solutions are consistent, indicating that all three methods can be used to solve the wave equation. However, as seen from the magnified portion of the curves, the curves of the analytical solution coincide with those of the finite difference method, while the curves of the SMD method deviate somewhat. The relative error of the curves between the analytical solution and the finite difference method is 0.013%, while that between the SMD method and the finite difference method is 0.40%, indicating that the proposed analytical solution has the same accuracy as the finite difference method. Here, the relative error is defined as $\max | u_{p}^{'} (t) - u_{p d} (t) | / u_{p d m} \times 100 %$ , where $u_{p}^{'} (t)$ is the pump displacement curve calculated by the analytical solution or the SMD method, $u_{p d} (t)$ is the pump displacement curve calculated by the finite difference method, and $u_{p d m}$ is the maximum value of $u_{p d} (t)$ .

Validation 2: For well 2, in the SMD method, the numbers of elements are 40 and 45. The reason why the numbers are not 6 and 6, i.e., $L_{i m} / 6$ , is that the simulation error is too large. The mass $m_{n}$ , elastic modulus $k_{n}$ and friction coefficient $b_{n}$ of the i-th element are [59.6808 44.5140] kg, [3.9914 2.9771] MPa·m, [60.7456 46.1426] Pa·s·m, respectively. The given upper boundary and the lower boundary are shown in Figure 5a,b, and the pump displacements calculated by the analytical solution, the finite difference method, and the SMD method are shown in Figure 5c. As shown in Figure 5c, the calculated results of the three solutions are also the same. As shown in the magnified portion of the curve, the curve of the analytical solution coincides well with that of the finite difference method, while that of the SMD method has a certain deviation. The relative error of the curves between the analytical solution and the finite difference method is 0.062%, while that between the SMD method and the finite difference method is 0.43%.

Thus, it can be concluded that, under the same boundary conditions, the proposed analytical solution has the same accuracy as the finite difference method, while the SMD method has some errors, and the sources of error are under study. Therefore, the proposed analytical solution is feasible. In the next hybrid method validation, the finite difference method is taken as the standard and compared.

3.2. Validation of the Hybrid Method by the Simulated Results

Validation 1: For well 1, the surface dynamometer card and the pump card are simulated by the hybrid method and the finite difference method, as shown in Figure 6. As shown in Figure 6, the surface dynamometer card and the pump card simulated by both methods are perfectly matched, with relative area errors of only 0.06% and 0.09%, respectively. This indicates that the hybrid method gives the same simulation results as the finite difference method. Here, the relative area error is defined as the ratio of the absolute value of two area differences to one of the areas.

Validation 2: For well 2, the cards simulated by both methods are shown in Figure 7. As shown in Figure 7, the cards simulated by both methods are also the same, with relative area errors of 0.02% and 0.03%, respectively.

It can be concluded that the hybrid method gives the same simulation results as the finite difference method.

3.3. Validation of the Hybrid Method by the Measured Results

Validation 1: For well 1, the measured surface dynamometer card and the simulated card obtained by the hybrid method are shown in Figure 8. It should be noted that the buoyant rod weight in the simulated dynamometer card is adjusted by −1.5 kN. This error may be due to the calculation method of the buoyant rod weight considering the coupling of the sucker-rod string or the calibration of the polished rod transducers [35]. The relative area error excludes the influence of the weight of the buoyant rod. Therefore, it is used in this paper. The relative area error between the measured and simulated cards is 1.89%. Therefore, the results of the hybrid method in this paper are in good agreement with the field test results.

Validation 2: For well 2, the surface dynamometer card and pump card simulated by the hybrid method and by Xing [38] (defined as an old method), as well as the measured surface dynamometer card, are shown in Figure 9. As indicated in Figure 9, the three types of surface dynamometer cards exhibited the same trends. However, the area relative error between the surface dynamometer card simulated by the hybrid model and the measured card is 1.00%, while that between the surface dynamometer card simulated by the old model and the measured model is 23.12%. Therefore, the results of the proposed method are closer to the measured data. There are also the same trends for the pump card, and the area relative error between the pump card simulated by the hybrid method and those simulated by the old model is 2.16%. This is because the hybrid method does not consider hydraulic loss or Newtonian fluid leakage [38]. Validation by surface dynamometer card measurements from wells 1 and 2 shows that the hybrid method is feasible.

3.4. Sensitivity Analysis of the Pump Speed

To further illustrate the effectiveness of the hybrid method, a simulation is carried out with pumping speeds of $n_{p}$ , $2 n_{p}$ , and $3 n_{p}$ .

Sensitivity analysis 1: For well 1, the surface dynamometer cards and pump cards are simulated by the hybrid method, and the finite difference method is shown in Figure 10. After the first loop of the hybrid method, the simulated pump cards are shown in Figure 10a–c. Figure 10a shows that at a low pumping speed, there are slight differences in the loading and unloading portions of both cards. Figure 10b,c shows that at higher speeds, the results of the hybrid method obviously deviate from those of the finite difference method, especially in the loading and unloading portions of the cards. Figure 10c also shows that there is an obvious smoothing effect at the highest peak. In the first cycle of the hybrid method, the pump load $P_{p s} (t)$ is the static pump load calculated according to the static pumping condition model, i.e., calculated according to the polished rod velocity. Therefore, it can be concluded that the static model is only applicable to the case of low pumping speed, and a smoothing effect may occur at high pumping speed. Our previous study [20], which is the conclusion of the analytical solution, is sensitive to the viscous damping factor and may belong to this case. The sensitivity of the analytical solutions to the viscous damping factor is being further investigated.

When the hybrid method meets the accuracy requirements, the loop is finished, and the simulation results are shown in Figure 10d–f. As illustrated in Figure 10d–f, the cards simulated by the hybrid method perfectly match those simulated by the finite difference method. Therefore, even at higher pump speeds, the hybrid method can still eliminate the influence of the static model and has high accuracy.

Sensitivity analysis 2: For well 2, the surface dynamometer cards and pump cards are simulated by the hybrid method, and the finite difference method at different pumping speeds is shown in Figure 11. After the first loop of the hybrid method, the simulated pump cards are shown in Figure 11a–c. Figure 11a–c shows that the trend is the same as that of Figure 10a–c, i.e., the static model is suitable for the low pumping speed case, and the higher the pumping speed is, the greater the simulation error. When the loop is finished in the hybrid method, the simulation results are shown in Figure 11d–f. As illustrated in Figure 11d–f, the cards simulated by the hybrid method also perfectly match with those simulated by the finite difference method. As a result, the hybrid method eliminates the effects of the static pumping condition model and has high accuracy despite being applied to different pumping wells and higher pumping speeds.

3.5. Numerical Analysis of the Convergence

In this section, the curves between the spectral radius $ρ (M_{i t, n})$ and the Fourier series number $n$ , i.e., the convergence curves of the hybrid method are calculated according to Equation (29). The convergence curves in domain 1 are calculated by setting $γ_{s 1} γ_{s 2} = 1$ in Equation (29). The convergence curves in the whole domain are calculated by setting $γ_{s 1} γ_{s 2} = 1$ and $[\begin{matrix} M_{K, 1}^{L - 5 \times 0} & M_{K, 2}^{L - 5 \times 0} \end{matrix}]$ in Equation (29).

Convergence analysis 1: The convergence curves for oil well 1 at pump speeds of 1 n_p and 3 n_p are shown in Figure 12 and Figure 13. As shown in Figure 12, there are four peaks in the convergence curves, with the largest peak occurring at a Fourier series number of 26. The spectral radius of the largest peak in all three convergence curves is less than 1, particularly in the whole domain, indicating that the algorithm converges when using the analytical solution iteration method according to reference [37]. The spectral radius of domain 1 is slightly smaller than that of the whole domain. According to Equation (28), domain 2 shortens the length of the rod but accounts for a smaller percentage of 5 $∆ x$ /L = 6.84%. The maximum peak spectral radius of the hybrid method is the smallest, indicating that the introduction of the finite difference method allows the algorithm to more easily converge, i.e., it can improve the convergence of the algorithm.

As illustrated in Figure 13, the overall trend is the same as that of Figure 12, with the difference that the number of peaks in the convergence curve increases to 11, indicating an increase in the number of harmonics, and the largest peak of the spectral radius occurs at a Fourier series number of 9, which indicates that a frequency redshift occurs, since the Fourier series number corresponds to the frequency of the harmonics. The value of the spectral radius is smaller than that in Figure 12, indicating that the convergence of the algorithm is enhanced after increasing the pumping speed.

Convergence analysis 2: The convergence curves for oil well 2 at pump speeds of 1 n_p and 3 n_p are shown in Figure 14 and Figure 15. As illustrated in Figure 14, the convergence curves have seven peaks, with the largest peak occurring at a Fourier number of 15. The maximum peak spectral radius of domain 1 and the whole domain is greater than 1, indicating that the iterative algorithm of reference [37] does not converge. When the maximum peak spectral radius of the hybrid method is less than 1, the hybrid method converges, i.e., the introduction of the finite difference method makes the algorithm convergent. The spectral radius of domain 1 is significantly smaller than that of the whole domain, which is different from that of well 1 in Figure 12, and according to Equation (28), the reason is that the length of the rods in domain 2 accounts for a larger proportion of the total length of the rods, which is 5 $∆ x$ /L = 44.12%.

As shown in Figure 15, the overall trend is the same as that of Figure 14, with the difference that the number of peaks in the convergence curve increases, indicating an increase in the number of harmonics and the maximum peak of the spectral radius occurs at a Fourier number of 5, i.e., the frequency redshift. The value of the maximum peak spectral radius is reduced compared to that in Figure 14, indicating that the algorithm is more likely to converge after increasing the pump speed.

3.6. Spectral Analysis of the Polisher Rod Load and Pump Load

Based on frequency spectrum analysis, He et al. [41] estimated the typical characteristics of operating conditions and affecting factors of the dynamometer cards. Therefore, frequency spectrum analysis is a useful tool for the character analysis. To further understand the relationship between the polisher rod load and pump load, the frequency spectrum is analyzed based on the simulated polisher rod load and pump load.

Spectral analysis 1: The polisher rod load and pump load simulated by the hybrid method for well 1 are shown in Figure 16a, and their frequency spectrums are shown in Figure 16b. As shown in Figure 16b, the signal characteristics of the two spectrums are similar, and the amplitude of the main frequency is basically the same; only the amplitude increases in the frequency region of [0.3192–0.9044] Hz. After deleting the frequency region and reconstructing the signal, the modified frequency spectrum and the polisher rod load are obtained, as shown in Figure 16c,d. The vibration signal disappears in the modified polisher rod load. This indicates that the spectrum analysis method is helpful for separating the vibration signal. The maximum value of the frequency region is 0.67 Hz, as shown in Figure 16b, 0.67 Hz × T = 26, which is the position of the first peak in the convergence curves in Figure 12. Therefore, the resonant frequency of the sucker-rod for the frequency spectrum of the polisher rod load can be obtained.

Spectral analysis 2: The polisher rod load and pump load simulated by the hybrid method for well 2 are shown in Figure 17a, and their frequency spectrums are shown in Figure 17b. As shown in Figure 17b, the shape characteristics of the two spectrums are basically the same; only the amplitude increases in the frequency region of [0.6344–0.9516] Hz. After reconstructing the signal after deleting the frequency region, the modified frequency spectrums and the polisher rod load are obtained, as shown in Figure 17c,d. The vibration signal also disappears in the modified polisher rod load. The maximum value of the frequency region is 0.80 Hz, as shown in Figure 17b, 0.80 Hz × T = 15; it is also the position of the first peak in the convergence curves in Figure 14. This trend is the same as that of spectrum analysis 1, in which the low-frequency portion of the spectrum reflects the pump conditions, while the high-frequency portion corresponds to the sucker-rod vibration information.

4. Conclusions

This paper presents a new hybrid method for predicting the behaviors of sucker-rod pumping units. The hybrid method includes an analytical solution of the sucker-rod wave equation, a unified pumping condition model, and a modified finite difference method of the sucker-rod wave equation. The proposed analytical solution is the recursive matrix form based on the truncated Fourier series of the polished rod displacement and the pump load. The established unified pumping condition model includes both the static and dynamic models. In the modified finite difference method, the upper boundary is the displacement and load. In the hybrid method, the whole rod string is divided into two calculation domains. In the small domain near the pump, the sucker-rod string wave equation is solved by the modified finite difference method to handle the “self-sensing” pumping condition model. In the other large domain, it is solved by the analytical solution; thus, the disadvantage of not being able to handle the “self-sensing” pumping condition model is overcome. The convergence calculation method of the hybrid method is proposed.

Based on two different types of oil wells, the proposed analytical solution of the sucker-rod wave equation is verified by the finite difference method and the SMD method with the same given upper boundary and lower boundary. The hybrid method is verified by using simulated data from the classical finite difference method and measured data. The sensitivities of the pumping speed and the convergence of the hybrid method are investigated. The frequency spectrums of the polisher rod load and pump load simulated by the hybrid method are analyzed. The main results can be summarized as follows:

(1). The proposed analytical solution has the same accuracy as the finite difference method, while the SMD method has a lower accuracy. The maximum relative error between the analytical solution and the finite difference method is 0.062%, while that between the SMD method and the finite difference solution is 0.43%;
(2). The proposed hybrid method is feasible and has a high simulation accuracy. The maximum relative area error between the hybrid method and the finite difference method is 0.09%, and the maximum relative area error between the hybrid method and the measured data is 1.89%;
(3). Despite the increase in the pumping speed, the hybrid method can still eliminate the influence of the static pumping condition model and has the same accuracy as the finite difference method. The hybrid method converges in this paper because the maximum spectral radius of the iterative matrix is less than 1. The introduction of the finite difference method allows the algorithm of the hybrid method to more easily converge;
(4). The frequency spectrums of the polisher rod load and pump load simulated by the hybrid method are similar. The low-frequency portion of the spectrum can reflect the pump conditions, while the high-frequency portion corresponds to the sucker-rod vibration information.

The research presented in this paper can be applied to fault diagnosis, system design, energy consumption analysis, parameter optimization, and system control of sucker-rod pumping systems. It also provides a new way to solve the wave equation with complex boundary conditions.

Author Contributions

Conceptualization, J.Y.; methodology, J.Y.; software, J.Y.; validation, J.Y.; formal analysis, J.Y.; investigation, J.Y.; resources, J.Y.; data curation, J.Y.; writing—original draft preparation, J.Y.; writing—review and editing, J.Y.; visualization, J.Y.; supervision, J.Y.; project administration, H.M.; funding acquisition, H.M. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Acknowledgments

Comments and suggestions from the editor and reviewers are very much appreciated.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Table

Figure 1. Flowchart of the solving procedure.

Figure 2. Schematic of the calculation domain in the hybrid method.

Figure 3. Flowchart of the hybrid method procedure.

Figure 4. Calculated results of the analytical solution, the finite difference method, and the SMD method for well 1.

Figure 5. Calculated results of the analytical solution, the finite difference method, and the SMD method for well 2.

Figure 6. Simulation cards of the hybrid method and the finite difference method for well 1.

Figure 7. Simulation cards of the hybrid method and the finite difference method for well 2.

Figure 8. Measured surface dynamometer card and simulated card by the hybrid method for well 1.

Figure 9. Measured surface dynamometer card and simulated cards by the hybrid method and the old method for well 2.

View Image - Figure 10. Cards simulated by the hybrid method and the finite difference method with different pump speeds for well 1. (a) First loop at a speed of np. (b) First loop at a speed of 2 np. (c) First loop at a speed of 3 np. (d) Finished loop at a speed of np. (e) Finished loop at a speed of 2 np. (f) Finished loop at a speed of 3 np.

Figure 10. Cards simulated by the hybrid method and the finite difference method with different pump speeds for well 1. (a) First loop at a speed of np. (b) First loop at a speed of 2 np. (c) First loop at a speed of 3 np. (d) Finished loop at a speed of np. (e) Finished loop at a speed of 2 np. (f) Finished loop at a speed of 3 np.

View Image - Figure 11. Cards simulated by the hybrid method and the finite difference method with different pump speeds for well 2. (a) First loop at a speed of np. (b) First loop at a speed of 2 np. (c) First loop at a speed of 3 np. (d) Finished loop at a speed of np. (e) Finished loop at a speed of 2 np. (f) Finished loop at a speed of 3 np.

Figure 11. Cards simulated by the hybrid method and the finite difference method with different pump speeds for well 2. (a) First loop at a speed of np. (b) First loop at a speed of 2 np. (c) First loop at a speed of 3 np. (d) Finished loop at a speed of np. (e) Finished loop at a speed of 2 np. (f) Finished loop at a speed of 3 np.

Figure 12. Convergence curve of well 1 at a pump speed of 1 np.

Figure 13. Convergence curve of well 1 at a pump speed of 3 np.

Figure 14. Convergence curve of well 2 at a pump speed of 1 np.

Figure 15. Convergence curve of well 2 at a pump speed of 3 np.

Figure 16. Load curves and its frequency spectrums of well 1.

Figure 17. Load curves and its frequency spectrums of well 2.

Table 1

Basic parameters of the oil wells.

Item	Value	Value
Well number	1	2
Pumping unit	DX700 (Rotaflex)	CYJ14-6-73HB (Beam pumping)
Pump stroke, m	5.9	4.0
Pumping speed, min⁻¹	1.6	3.18
Sucker-rod string, mm × m	25 × 389 + 22 × 1322 + 28 × 153	22 × 800 + 19 × 900
Tubing string, mm × m	76 × 543 + 62 × 982, Unanchored	62 × 1400, Anchored
Pump diameter, mm	44	56
Pump depth, m	1904	1700
Fluid density, kg/m³	915.4	959
Dynamic liquid level, m	1554	1400
Oil pressure, MPa	1.2	0.6
Casing pressure, MPa	0	0.3
Fluid viscosity, mPa.s	329.84	430
Gas/oil ratio, m³/m³	0	20
Rod and tube’s density, kg/m³	7850 (Steel)	7850 (Steel)
Rod and tube’s Young’s modulus, GPa	210 (Steel)	210 (Steel)

Appendix A

Appendix A.1. Proof of Proposition 1

Proof.

Considering all the number of tapers in the string yields (A1) ${[\begin{matrix} π ε_{n} \\ κ μ_{n} \end{matrix}]}_{1} = [\begin{matrix} M_{t 1} & M_{t 2} \\ F_{t 1} & F_{t 2} \end{matrix}] {[\begin{matrix} π ε_{n} \\ κ μ_{n} \end{matrix}]}_{K}$ where (A2) $[\begin{matrix} M_{t 1} & M_{t 2} \\ F_{t 1} & F_{t 2} \end{matrix}] = \prod_{i = 1}^{K - 1} {[\begin{matrix} M_{i, 1}^{l} & M_{i, 2}^{l} \\ E_{r i} A_{i} F_{i, 1}^{l} & E_{r i} A_{r i} F_{i, 2}^{l} \end{matrix}]}^{- 1} [\begin{matrix} M_{i + 1, 1}^{l} & M_{i + 1, 2}^{l} \\ E_{r i + 1} A_{r i + 1} F_{i + 1, 1}^{l} & E_{r i + 1} A_{r i + 1} F_{i + 1, 2}^{l} \end{matrix}]$

After the polished rod displacement $u_{a} (t)$ and the pump load $P_{p} (t)$ are approximated by the truncated Fourier series, the boundary conditions of Equations (2) can be expressed as follows:(A3) $\begin{array}{l} [M_{1, 1}^{0}] {[π ε_{n}]}_{1} + [M_{1, 2}^{0}] {[κ μ_{n}]}_{1} = [v δ_{n}], \\ E_{r K} A_{r K} ([F_{K, 1}^{L}] {[π ε_{n}]}_{K} + [F_{K, 2}^{L}] {[κ μ_{n}]}_{K}) = [σ τ_{n}] \end{array}$ where (A4) $\begin{array}{l} [v δ_{n}] = - {[\begin{matrix} v_{n} & δ_{n} \end{matrix}]}^{T}, \\ [σ τ_{n}] = {[\begin{matrix} σ_{n} & τ_{n} \end{matrix}]}^{T}, \\ u_{a} (t) = \frac{ν_{0}}{2} + \sum_{n = 1}^{N} (ν_{n} \cos n ω t + δ_{n} \sin n ω t), \\ P_{p} (t) = \frac{σ_{0}}{2} + \sum_{n = 1}^{N} (σ_{n} \cos n ω t + τ_{n} \sin n ω t) \end{array}$

$v_{n}$ and $δ_{n}$ are the Fourier coefficients of the displacement function, and $σ_{n}$ and $τ_{n}$ are the Fourier coefficients of the dynamic load function.

Considering $x = 0$ and according to Equations ((A1), (A3)), the following is obtained:(A5) $[\begin{matrix} M_{t 1} & M_{t 2} \\ E_{r K} A_{r K} F_{K, 1}^{L} & E_{r K} A_{r K} F_{K, 2}^{L} \end{matrix}] {[\begin{matrix} π ε_{n} \\ κ μ_{n} \end{matrix}]}_{K} = [\begin{matrix} δ v_{n} \\ σ τ_{n} \end{matrix}]$

Thus, Equation (A5) can be resolved as follows:(A6) ${[\begin{matrix} κ μ_{n} \\ π ε_{n} \end{matrix}]}_{K} = {[\begin{matrix} M_{t 1} & M_{t 2} \\ E_{r K} A_{r K} F_{K, 1}^{L} & E_{r K} A_{r K} F_{K, 2}^{L} \end{matrix}]}^{- 1} [\begin{matrix} δ v_{n} \\ σ τ_{n} \end{matrix}]$ where $[M F_{n}] = {[\begin{matrix} M_{t 1} & M_{t 2} \\ E_{r K} A_{r K} F_{K, 1}^{L} & E_{r K} A_{r K} F_{K, 2}^{L} \end{matrix}]}^{- 1}$ .

Proposition 1 is therefore proved. □

Appendix A.2. Proof of Proposition 2

Proof.

During the iterative process, the polished rod displacement does not change. The input of the analytical solution is the pump load $P_{p} (t)$ and its output is the middle displacement $u_{u 1} (t)$ , so $u_{u 1} (t) = [A_{M}] P_{p} (t)$ , where $[A_{M}]$ is the transfer coefficient matrix of the analytical solution. If $Δ P_{p 0} (t)$ is set to be the deviation between the initial value of the pump load $P_{p s} (t)$ and the true value of the pump load $P_{p r} (t)$ , thus, the deviation of the middle displacement $Δ u_{u 1} (t)$ can be calculated according to the analytical solution, that is, $Δ u_{u 1} (t) = [A_{M}] Δ P_{p 0} (t)$ . The new deviation of the pump load $Δ P_{p 1} (t)$ is calculated according to the finite difference method, that is, $Δ P_{p 1} (t) = [F_{M}] Δ u_{u 1} (t)$ , and $[F_{M}]$ is the transfer coefficient matrix of the finite difference method. It can be concluded that $Δ P_{p 1} (t) = [F_{M}] [A_{M}] Δ P_{p 0} (t)$ , that is, the iterative coefficient matrix is $[M_{i t}] = [F_{M}] [A_{M}]$ .

Proposition 2 is therefore proved. □

References

1. Wu, Y.; Liu, S.; Ma, W.; Ran, X.; Qu, B. Machine Learning Method for Predicting the Fatigue Life of Sucker Rods. Eng. Fract. Mech.; 2023; 282, 109161. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2023.109161]

2. He, Y.; Zang, C.; Zeng, P.; Wang, M.; Dong, Q.; Wan, G.; Dong, X. Few-Shot Working Condition Recognition of a Sucker-Rod Pumping System Based on a 4-Dimensional Time-Frequency Signature and Meta-Learning Convolutional Shrinkage Neural Network. Petrol. Sci.; 2023; 20, pp. 1142-1154. [DOI: https://dx.doi.org/10.1016/j.petsci.2023.02.017]

3. Takács, G.A. Critical Analysis of Power Conditions in Sucker-Rod Pumping Systems. J. Petrol. Sci. Eng.; 2022; 210, 110061. [DOI: https://dx.doi.org/10.1016/j.petrol.2021.110061]

4. Langbauer, C.; Langbauer, T.; Fruhwirth, R.; Mastobaev, B. Sucker Rod Pump Frequency-Elastic Drive Mode Development—From the Numerical Model to the Field Test. Liq. Gas Energy Resour.; 2021; 1, pp. 64-85. [DOI: https://dx.doi.org/10.21595/lger.2021.22074]

5. Bangert, P.; Sharaf, S. Diagnosing and Predicting Problems with Rod Pumps Using Machine Learning. Proceedings of the SPE Middle East Oil and Gas Show and Conference; Manama, Bahrain, 18–21 March 2019; [DOI: https://dx.doi.org/10.2118/195845-ms]

6. Wang, X.; He, Y.; Li, F.; Dou, X.; Wang, Z.; Xu, H.; Fu, L. A Working Condition Diagnosis Model of Sucker Rod Pumping Wells Based on Big Data Deep Learning. Proceedings of the International Petroleum Technology Conference; Beijing, China, 26–28 March 2019; [DOI: https://dx.doi.org/10.2523/iptc-19242-ms]

7. Xu, J. A Method for Diagnosing the Performance of Sucker Rod String in Straight Inclined Wells. Proceedings of the SPE Latin America/Caribbean Petroleum Engineering Conference; Buenos Aires, Argentina, 27–29 April 1994; [DOI: https://dx.doi.org/10.2118/26970-ms]

8. Zheng, B.; Gao, X. Sucker Rod Pumping Diagnosis Using Valve Working Position and Parameter Optimal Continuous Hidden Markov Model. J. Process Control; 2017; 59, pp. 1-12. [DOI: https://dx.doi.org/10.1016/j.jprocont.2017.09.007]

9. Zheng, B.; Gao, X.; Li, X. Diagnosis of Sucker Rod Pump Based on Generating Dynamometer Cards. J. Process Control; 2019; 77, pp. 76-88. [DOI: https://dx.doi.org/10.1016/j.jprocont.2019.02.008]

10. Han, Y.; Li, K.; Ge, F.; Wang, Y.; Xu, W. Online Fault Diagnosis for Sucker Rod Pumping Well by Optimized Density Peak Clustering. ISA Trans.; 2021; 120, pp. 222-234. [DOI: https://dx.doi.org/10.1016/j.isatra.2021.03.022] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33810843]

11. Ma, R.; Tian, H.; Cheng, X.; Xiao, Y.; Xu, Q.; Yu, X. Oil-Net: A Learning-Based Framework for Working Conditions Diagnosis of Oil Well through Dynamometer Cards Identification. IEEE Sens. J.; 2023; 23, pp. 14406-14417. [DOI: https://dx.doi.org/10.1109/jsen.2023.3275933]

12. He, Y.-P.; Cheng, H.-B.; Zeng, P.; Zang, C.-Z.; Dong, Q.-W.; Wan, G.-X.; Dong, X.-T. Working Condition Recognition of Sucker Rod Pumping System Based on 4-Segment Time-Frequency Signature Matrix and Deep Learning. Pet. Sci.; 2024; 21, pp. 641-653. [DOI: https://dx.doi.org/10.1016/j.petsci.2023.08.031]

13. Takacs, G. Calculation of Operational Parameters. Sucker-Rod Pumping Handbook: Production Engineering Fundamentals and Long-Stroke Rod Pumping; Gulf Professional Publishing: Houston, TX, USA, 2015; Chapter 4 pp. 247-376. [DOI: https://dx.doi.org/10.1016/b978-0-12-417204-3.00004-2]

14. Lv, X.; Wang, H.; Zhang, X.; Liu, Y.; Chen, S. An Equivalent Vibration Model for Optimization Design of Carbon/Glass Hybrid Fiber Sucker Rod Pumping System. J. Petrol. Sci. Eng.; 2021; 207, 109148. [DOI: https://dx.doi.org/10.1016/j.petrol.2021.109148]

15. Zhao, H.; Wang, Y.; Zhan, Y.; Xu, G.; Cui, X.; Wang, J. Practical Model for Energy Consumption Analysis of Beam Pumping Motor Systems and Its Energy-Saving Applications. IEEE Trans. Ind. Appl.; 2018; 54, pp. 1006-1016. [DOI: https://dx.doi.org/10.1109/tia.2017.2779103]

16. Li, W.; Dong, S.; Xiu-rong, S. An Improved Sucker Rod Pumping System Model and Swabbing Parameters Optimized Design. Math. Probl. Eng.; 2018; 2018, 4746210. [DOI: https://dx.doi.org/10.1155/2018/4746210]

17. Gibbs, S.G. Predicting the Behavior of Sucker-Rod Pumping Systems. J. Pet. Technol.; 1963; 15, pp. 769-778. [DOI: https://dx.doi.org/10.2118/588-pa]

18. Schäfer, D.J.; Jennings, J.W. An Investigation of Analytical and Numerical Sucker Rod Pumping Mathematical Models. Proceedings of the SPE Annual Technical Conference and Exhibition; Dallas, TX, USA, 27–30 September 1987; [DOI: https://dx.doi.org/10.2118/16919-ms]

19. Eisner, P.; Langbauer, C.; Fruhwirth, R.K. Comparison of a Novel Finite Element Method for Sucker Rod Pump Downhole Dynamometer Card Determination Based on Real World Dynamometer Cards. Upstream Oil Gas Technol.; 2022; 9, 100078. [DOI: https://dx.doi.org/10.1016/j.upstre.2022.100078]

20. Yin, J.-J.; Sun, D.; Yang, Y. Predicting Multi-Tapered Sucker-Rod Pumping Systems with the Analytical Solution. J. Petrol. Sci. Eng.; 2021; 197, 108115. [DOI: https://dx.doi.org/10.1016/j.petrol.2020.108115]

21. Torgaeva, D.S.; Shalyapina, N.A.; Sukhorukov, M.P. Simulation of Load on a Polished Rod of Sucker Rod Pump for Oil Production. Proceedings of the 2019 International Multi-Conference on Engineering, Computer and Information Sciences; Novosibirsk, Russia, 21–27 October 2019; [DOI: https://dx.doi.org/10.1109/sibircon48586.2019.8958101]

22. Xing, M.; Dong, S. A New Simulation Model for a Beam-Pumping System Applied in Energy Saving and Resource-Consumption Reduction. SPE Prod. Oper.; 2015; 30, pp. 130-140. [DOI: https://dx.doi.org/10.2118/173190-pa]

23. Li, W.; Vaziri, V.; Aphale, S.S.; Dong, S.; Wiercigroch, M. Dynamics and Frequency and Voltage Control of Downhole Oil Pumping System. Mech. Syst. Signal Process.; 2020; 139, 106562. [DOI: https://dx.doi.org/10.1016/j.ymssp.2019.106562]

24. Jian, L.; Xing, M. Characteristic Analysis of Longitudinal Vibration Sucker Rod String with Variable Equivalent Stiffness. Proceedings of the 4th International Conference on Mechatronics and Mechanical Engineering (ICMME 2017); Kuala Lumpur, Malaysia, 28–30 November 2017; [DOI: https://dx.doi.org/10.1051/matecconf/201815306009]

25. Zhang, B.; Gao, X.; Li, X. Complete Simulation and Fault Diagnosis of Sucker-Rod Pumping. SPE Prod. Oper.; 2021; 36, pp. 277-290. [DOI: https://dx.doi.org/10.2118/204215-pa]

26. Yin, J.; Sun, D.; Yang, Y. A Novel Method for Diagnosis of Sucker-Rod Pumping Systems Based on the Polished-Rod Load Vibration in Vertical Wells. SPE J.; 2020; 25, pp. 2470-2481. [DOI: https://dx.doi.org/10.2118/201228-pa]

27. Gibbs, S.G.; Neely, A.B. Computer Diagnosis of Down-Hole Conditions in Sucker Rod Pumping Wells. J. Pet. Technol.; 1966; 18, pp. 91-98. [DOI: https://dx.doi.org/10.2118/1165-pa]

28. Gibb, S.G. Method of Determining Sucker Rod Pump Performance. U.S. Patent; 3,343,409A, 26 September 1967.

29. Lin, K.S. Mathematical Modelling and Solution Procedure for the Behavior of Sucker-Rod Pumping System. Proceedings of the Ninth International Conference On Computer Applications; Yangon, Myanmar, 5–6 May 2011.

30. Li, J.; Li, Z. Prediction and diagnosis of sucker-rod pumping systems in directional wells. Soc. Petrol. Eng.; 1999; 1–8, SPE-57014-MS.Available online: https://www.semanticscholar.org/paper/Prediction-and-Diagnosis-of-Sucker-Rod-Pumping-in-Li-Li/4ab285ceb388f1d43a83223cc0604e13306f2276 (accessed on 6 May 2024).

31. Chen, Z.; White, L.W.; Zhang, H. Predicting Sucker-Rod Pumping Systems with Fourier Series. SPE Prod. Eng.; 2018; 33, pp. 928-940. [DOI: https://dx.doi.org/10.2118/189991-pa]

32. Xing, M. Response Analysis of Longitudinal Vibration of Sucker Rod String Considering Rod Buckling. Adv. Eng. Softw.; 2016; 99, pp. 49-58. [DOI: https://dx.doi.org/10.1016/j.advengsoft.2016.05.004]

33. Lv, X.; Wang, H.; Liu, Y.; Chen, S.; Lan, W.; Sun, B. A Novel Method of Output Metering with Dynamometer Card for SRPS under Fault Conditions. J. Pet. Sci. Eng.; 2020; 192, 107098. [DOI: https://dx.doi.org/10.1016/j.petrol.2020.107098]

34. Zhang, Z.; Dong, S.; Chen, Z.; Feng, G.; Leng, L. Reliability Prediction Method of Fatigue Life for Rod String. Proceedings of the International Conference on Quality, Reliability, Risk, Maintenance, and Safety Engineering; Chengdu, China, 15–18 June 2012; [DOI: https://dx.doi.org/10.1109/ICQR2MSE.2012.6246370]

35. Takacs, G. Long-Stroke Sucker-Rod Pumping. Sucker-Rod Pumping Handbook: Production Engineering Fundamentals and Long-Stroke Rod Pumping; Gulf Professional Publishing: Houston, TX, USA, 2015; Chapter 7 pp. 505-551. [DOI: https://dx.doi.org/10.1016/B978-0-12-417204-3.00007-8]

36. Xing, M.; Dong, S.; Tong, Z.; Tian, R.; Chen, H. Dynamic Simulation and Efficiency Analysis of Beam Pumping System. J. Cent. South Univ.; 2015; 22, pp. 3367-3379. [DOI: https://dx.doi.org/10.1007/s11771-015-2877-z]

37. Yin, J.; Ma, H. A Novel Method for Predicting the Behavior of a Sucker Rod Pumping Unit Based on the Polished Rod Velocity. Mathematics; 2024; 12, 1318. [DOI: https://dx.doi.org/10.3390/math12091318]

38. Xing, M. An Improved Numerical Simulation Research for Plunger Pump in the Condition of Newtonian Fluid. J. Meas. Eng.; 2016; 4, pp. 32-42.

39. Samuel, I.T.; Anatolii, Z. Analysis of Motor Power Curve for Detecting Fault Conditions in Sucker Rod Pump. Proceedings of the 2020 27th International Workshop on Electric Drives: MPEI Department of Electric Drives 90th Anniversary (IWED); Moscow, Russia, 27–30 January 2020; [DOI: https://dx.doi.org/10.1109/IWED48848.2020.9069555]

40. Zhang, Q.; Wang, J.; Cui, W.; Xiao, Z.; Yue, Q. Post-Buckling Transition of Compressed Pipe Strings in Horizontal Wellbores. Ocean. Eng.; 2020; 197, 106880. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2019.106880]

41. He, Y.; Wu, X.; Han, G.; Xiao, W.; Li, W.; Yu, X. Frequency spectrum analysis method for recognition of dynamometer card. Acta Pet. Sin.; 2008; 29, pp. 619-624. [DOI: https://dx.doi.org/10.7623/syxb200804029]

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

A Hybrid Method for Solving the One-Dimensional Wave Equation of Tapered Sucker-Rod Strings

Content area

Abstract

Full text