Pump System Model Parameter Identification Based on Experimental and Simulation Data

Abstract

In this paper, the entire downhole fluid-sucker rod-pump system is replaced with a viscoelastic vibration model, namely a third-order differential equation with an inhomogeneous forcing term. Both Kelvin’s and Maxwell’s viscoelastic models can be implemented along with the dynamic behaviors of a mass point attached to the viscoelastic model. By employing the time-dependent polished rod force measured with a dynamometer as the input to the viscoelastic dynamic model, we have obtained the displacement responses, which match closely with the experimental measurements in actual operations, through an iterative process. The key discovery of this work is the feasibility of the so-called inverse optimization procedure, which can be utilized to identify the equivalent scaling factor and viscoelastic system parameters. The proposed Newton–Raphson iterative method, with some terms in the Jacobian matrix expressed with averaged rates of changes based on perturbations of up to two independent parameters, provides a feasible tool for optimization issues related to complex engineering problems with mere information of input and output data from either experiments or comprehensive simulations. The same inverse optimization procedure is also implemented to model the entire fluid delivery system of a very viscous non-Newtonian polymer modeled as a first-order ordinary differential equation (ODE) system similar to the transient entrance developing flow. The convergent parameter reproduces transient solutions that match very well with those from fully fledged computational fluid dynamics models with the required inlet volume flow rate and outlet pressure conditions.

1. Introduction

Towards the end of an oil and gas reservoir’s productive life, natural lift forces due to fluid pressure gradually decay and diminish. Thus, it is common to use artificial lift methods to transfer the oil and gas from the underground formation to the surface [1,2,3]. In practice, sucker rod pumping units, as illustrated in Figure 1, are the most popular artificial lift means and are still widely utilized along with hydraulic and electric submersible pumping systems [4,5,6,7]. Teaming with such pumping systems, various four-bar linkage-based horse-head pump jacks have also been invented to support and drive the entire artificial lift system [8]. As one of the earliest inventions for oil industries, these sucker rod-pump systems have been proven to be efficient and adaptable [9,10,11,12,13]. Over $80\%$ of the wells worldwide still depend on this traditional and most extensively used mechanism [14,15]. Nevertheless, contemporary computational tools have been introduced in recent research efforts on leakage, safety, and reliability and shed light on these popular artificial lift systems [16,17,18,19].

The fluid mechanics aspect of the sucker rod-pump unit, such as the leakage, has been discussed and presented in refs. [20,21]. The so-called relaxation time, which can be derived numerically based on the Bessel functions for cylindrical coordinates and Fourier series for Cartesian coordinates is calculated and confirmed to be less than a tenth of the typical sampling period, which corresponds to a sampling frequency, normally ranging from 30 to 60 samples per second, currently used in oil industries with comparable physical dimensions [22,23,24]. The issues related to the leakage also exposed an imminent desire and need in engineering practice that sophisticated perturbation methods and analytical approaches continue to be essential to derive information from fully fledged computational simulations and, in some cases, to help the establishment of viable and more insightful simulation setups [25,26,27,28].

As the traveling unit, the plunger consists of a so-called traveling valve fixed at its bottom and is moving within a barrel or a chamber with a so-called standing valve fixed at the bottom [29,30]. Overall, the sucker rod pumping unit involves very complicated four strokes with respect to the open and closed positions of the traveling valve and the standing valve, as illustrated in Figure 1. In fact, through the exercises of applying fully fledged computational fluid dynamics (CFD) and finite element method (FEM) to sucker rod-pump systems, it is becoming clear that even for a simple developing viscous fluid within sucker rod-pump systems, engineers must learn how to handle extreme aspect ratios. For example, the so-called clearance C, namely the difference between the plunger outer diameter $D_o$ and the barrel inner diameter $D_i$, or the gap between the plunger and the barrel $\delta$, namely $C=2\delta =D_i-D_o$, is often measured in mills, one-thousandth of an inch, whereas, the plunger system length $L_p$ is often measured in feet [31,32]. This extreme aspect ratio renders the fluid mechanics studies of sucker rod-pump systems very challenging. It is more so when we must consider the difference between concentric and eccentric situations [33]. However, through our previous works [34], we have identified the possibilities of using a simple viscoelastic model to represent the entire complicated pumping strokes. The measured data at the surface are the sucker rod displacement as well as the dynamometer force measurements. We are encouraged by the very first confirmation of the Inversed Optimization method. This success might be a new beginning in using measured data to help the construction of mathematical and physical models when the engineering fluid system is very complicated. The structural and mechanisms aspects of this work have already been published in refs. [8,32].

To follow up with the ideas reported in ref. [34], the entire downhole sucker rod-pump systems are modeled with Kelvin or Maxwell viscoelastic systems, as illustrated in Figure 2. In this paper, we focus more on the iterative procedure to identify the optimal choice of needed parameters in these models with a so-called Inverse Optimization Method. The preliminary results do confirm that the parameters derived with the Inverse Optimization Method yield virtually identical results in comparison with the dynamometer measurements of the polished rod load and the displacement in the field [12]. In fact, the same inverse optimization procedure is also implemented to model the entire fluid delivery system of a very viscous non-Newtonian polymer as a mere first-order ordinary differential system. The convergent parameter matches very well with that identified with the fully fledged computational fluid dynamics model, with the required inlet and outlet pressure drop as the input and the flow rate as the output.

2. Theory and Modeling

Let us consider a lumped mass linked with either Maxwell or Kelvin models as illustrated in Figure 2. A sudden displacement is applied to the Maxwell model, whereas a sudden load is introduced to the Kelvin model [34]. Consequently, the Maxwell model will yield a stiffness function or relaxation function to relate the responding force with the sudden imposed displacement $u_o$ as shown in Figure 3. Likewise, the Kelvin model will produce a compliance function or creep function to relate the resulting displacement with the suddenly imposed force $f_o$, also shown in Figure 3. For the Maxwell model, we also introduce the displacement $u\left(t\right)=u_{o}H\left(t\right)$ where $H\left(t\right)$ is the Heaviside function. The total force function $f\left(t\right)$ can be expressed as the combination of the location force $f_o\left(t\right)=K_{o}u\left(t\right)$ and $f_1\left(t\right)$ for the stiffness and dashpot couple $(K_1,C_1)$ with the relaxation time ${\tau }_1=C_1/K_1$. Furthermore, using the same force or stress for the series component and the same displacement or strain for the parallel component, we obtain the following governing equations for the Maxwell model:

$$\begin{array}{ccc}\hfill f_o& =& K_{o}u\left(t\right),\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \dot{u}\left(t\right)& =& {\displaystyle \frac{f_1}{C_1}+\frac{{\dot{f}}_1}{K_1}.}\hfill \end{array}$$

(2)

Hence, using the concepts we have been discussing for the first-order ODE, we obtain:

$${\displaystyle e^{t/{\tau }_1}{f}_1\left(t\right)-0={\int }_{-\infty }^t{K}_1{e}^{t/{\tau }_1}\dot{u}\left(t\right)dt,}$$

(3)

namely using integration by parts,

$${\displaystyle e^{t/{\tau }_1}{f}_1\left(t\right)=K_1{e}^{t/{\tau }_1}u\left(t\right)-{\int }_{-\infty }^t{K}_1/{\tau }_{1}u\left(t\right)e^{t/{\tau }_1}dt=K_1{e}^{t/{\tau }_1}{u}_o-{\int }_0^t{K}_1/{\tau }_1{u}_o{e}^{t/{\tau }_1}dt.}$$

(4)

Finally, we have $f_1\left(t\right)=K_1{u}_o{e}^{-t/{\tau }_1}$ and the total force $f\left(t\right)=f_1\left(t\right)+f_o\left(t\right)=(K_o+K_1{e}^{-t/{\tau }_1})u_o$ for the Maxwell model. Likewise, denote the sudden load as $f\left(t\right)=f_{o}H\left(t\right)$ where $H\left(t\right)$ is a Heaviside function, the total displacement function $u\left(t\right)$ can be expressed as the combination of the location displacement $u_o\left(t\right)=f\left(t\right)/K_o$ and $u_1\left(t\right)$ for the stiffness and dashpot couple $(K_1,C_1)$ with the relaxation time ${\tau }_1=C_1/K_1$. Using the same force or stress for the series component and the same displacement or strain for the parallel component, we obtain the following governing equations for the Kelvin model:

$$\begin{array}{ccc}\hfill f\left(t\right)& =& K_o{u}_o\left(t\right),\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill f\left(t\right)& =& C_1{\dot{u}}_1+K_1{u}_1,1\le i\le n.\hfill \end{array}$$

(6)

Thus, using the concepts we have been discussing for the first-order ODE, we obtain:

$${\displaystyle e^{t/{\tau }_1}{u}_1\left(t\right)-0={\int }_{-\infty }^t{e}^{t/{\tau }_1}f/C_{1}dt={\int }_0^t{e}^{t/{\tau }_1}{f}_o/C_{1}dt.}$$

(7)

Finally, we have the local displacement ${\displaystyle u_1\left(t\right)=f_o/K_1(1-e^{-t/{\tau }_1})}$ and the total displacement function can be expressed as ${\displaystyle u\left(t\right)=u_o+u_1=f_o/K_o+f_o/K_1(1-e^{-t/{\tau }_1})}$ for the Kelvin model. To further the discussion with relaxation and creep functions, we can also introduce a general Kelvin viscoelastic model as elaborated in ref. [35]. Similar approaches are also implemented for the viscoelastic vibration model with a mass m and a concentrated load $F\left(t\right)$. We first introduce the Kelvin viscoelastic model in a dynamic case as shown in Figure 2, essentially a typical Kelvin viscoelastic setup for creep test combined in series with a spring with a stiffness $k_o$ and a mass m. The displacement of the parallel section of the stiffness $k_1$ and the dashpot $c_1$ share the same displacement $u_2\left(t\right)$, whereas the displacement of the stiffness $k_o$ is denoted as $u_1\left(t\right)$. Since the mass m is connected with the stiffness $k_o$ in series, the total displacement of the mass $u\left(t\right)$ is a combination of the two displacements $u_1\left(t\right)$ and $u_2\left(t\right)$. In general, the external load $F\left(t\right)$ is applied to the mass m. We can imagine that in the creep test, we can add a dead weight $W_o$ in addition to the weight of the mass $mg$.

Using the same procedure and considering each section in series and the consequent continuity of axial forces, we have

$$\begin{array}{ccc}\hfill k_o{u}_1\left(t\right)& =& k_1{u}_2+c_1{\dot{u}}_2,\hfill \\ \hfill k_o{u}_1\left(t\right)& =& F\left(t\right)-m\ddot{u}.\hfill \end{array}$$

(8)

Using the kinematic relationship $u\left(t\right)=u_1\left(t\right)+u_2\left(t\right)$, we obtain the following third-order governing equation for $u_2\left(t\right)$,

$${\displaystyle \frac{m{c}_1}{k_o}{\stackrel{⃛}{u}}_2+\frac{m}{k_o}(k_o+k_1){\ddot{u}}_2+c_1{\dot{u}}_2+k_1{u}_2=F\left(t\right).}$$

(9)

In this paper, for simplicity, in the implementation of the inverse optimization approaches, Equation (9) is directly applied with adjusted parameters $k_o$, m, $c_1$, and $k_1$. The measured force $F\left(t\right)$ is introduced to derive the closest $u_2\left(t\right)$ which matches with the measured displacement. As validated in ref. [34], with no dashpot, namely $c_1=0$, Equation (8) yields

$$k_1{u}_2=k_o{u}_1,$$

hence

$$u=u_1+u_2=\frac{k_o+k_1}{k_o}{u}_2.$$

Finally, the governing Equation (9) yields the familiar spring-mass vibration system

$$m\ddot{u}+\frac{k_o{k}_1}{k_o+k_1}u=F\left(t\right).$$

(10)

Notice the effective stiffness of two springs in series. Moreover, the infinitely stiff spring $k_o$, namely $k_o\to +\infty$, consequently, $u_2\to 0$ and $u\left(t\right)=u_2\left(t\right)$. Finally, we recover the familiar spring-mass-dashpot vibration system

$$m\ddot{u}+c_1\dot{u}+k_{1}u=F\left(t\right).$$

Finally, in order to facilitate the solution, introduce the state variable $\mathbf{y}=<u_2,{\dot{u}}_2,{\ddot{u}}_2>$, we can rewrite the third-order viscoelastic vibration system with a dynamical system format,

$$\mathbf{A}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ {\displaystyle -\frac{k_o{k}_1}{m{c}_1}}& {\displaystyle -\frac{k_o}{m}}& {\displaystyle -\frac{k_o+k_1}{c_1}}\end{array}\right]\mathrm{and}\mathbf{f}=\left(\begin{array}{c}0\\ 0\\ {\displaystyle \frac{k_{o}F}{m{c}_1}}\end{array}\right).$$

Consider the Maxwell viscoelastic model in a dynamic case as shown in Figure 2, essentially a Maxwell viscoelastic setup for relaxation test in combination with a mass m connecting in series with two parallel branches, one with a stiffness $k_o$ and another one with a stiffness $k_1$ and a dashpot $c_1$ in series. The displacement of the parallel section shares the same displacement $u\left(t\right)$, whereas the displacement of the stiffness $k_1$ is denoted as $u_2\left(t\right)$ and the dashpot $c_1$ as $u_1\left(t\right)$. Please note that the displacement of the stiffness $k_o$ is $u\left(t\right)$. Thus, the total displacement of the mass $u\left(t\right)$ is a combination of the two displacements $u_1\left(t\right)$ and $u_2\left(t\right)$ within one of the parallel branches. In general, the external load $F\left(t\right)$ is applied to the mass m. In a relaxation test, we can simply add a dead weight $W_o$ in addition to the weight of the mass $mg$. Similarly, assuming the force in the dashpot and stiffness series as $F_o$, and considering each section in series and the consequent continuity of axial forces, we have

$$\begin{array}{ccc}\hfill F-m\ddot{u}& =& k_{o}u+F_o,\hfill \\ \hfill \dot{u}& =& {\displaystyle \frac{{\dot{F}}_o}{k_1}+\frac{F_o}{c_1}.}\hfill \end{array}$$

(11)

Using the kinematic relationship $u\left(t\right)=u_1\left(t\right)+u_2\left(t\right)$ as well as $\dot{u}\left(t\right)={\dot{u}}_1\left(t\right)+{\dot{u}}_2\left(t\right)$, we obtain the following third-order governing equation for $u\left(t\right)$,

$${\displaystyle \frac{m}{k_1}{\stackrel{⃛}{u}}_2+\frac{m}{c_1}{\ddot{u}}_2+\frac{k_o+k_1}{k_1}{\dot{u}}_2+\frac{k_o}{c_1}u=\frac{\dot{F}}{k_1}+\frac{F}{c_1}.}$$

(12)

Again, with no dashpot, namely $c_1=0$. Equation (8) yields $F_o=0$, hence the governing Equation (12) yields the familiar spring-mass vibration system

$$m\ddot{u}+ku=F\left(t\right).$$

Likewise, as the infinitely stiff spring $k_1$, namely $k_1\to +\infty$, consequently, $u_2\to 0$ and $u\left(t\right)=u_1\left(t\right),$ we recover the familiar spring-mass-dashpot vibration system

$$m\ddot{u}+c_1\dot{u}+k_{o}u=F\left(t\right).$$

As presented in ref. [34], in numerical solutions, we rewrite the third-order viscoelastic vibration system with a dynamical system format,

$$\dot{\mathbf{y}}=\mathbf{A}\mathbf{y}+\mathbf{f},$$

(13)

with the state variable $\mathbf{y}=<u_2,{\dot{u}}_2,{\ddot{u}}_2>$, and

$$\mathbf{A}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ {\displaystyle -\frac{k_o{k}_1}{m{c}_1}}& {\displaystyle -\frac{k_o+k_1}{m}}& {\displaystyle -\frac{k_1}{c_1}}\end{array}\right]\mathrm{and}\mathbf{f}=\left(\begin{array}{c}0\\ 0\\ {\displaystyle \frac{\dot{F}}{m}+\frac{k_{1}F}{m{c}_1}}\end{array}\right).$$

As shown in Figure 4 and Figure 5, for such a viscoelastic model representing the downhole sucker rod system and starting from the reference position with no initial velocity and acceleration, we employ $k_o=4 N/m$, $m=1 kg$, $c_1=0.2 Ns/m$, $k_1=4 N/m$, and $W_o=10 N$. The external force $F\left(t\right)$ is a deadweight $W_o$. Based on the Matlab simulation, the vibration eventually settles down to an equilibrium position due to the dashpot damping. Of course, the deadweight will be replaced with the actual dynamometer load measurements, and the actual polish rod displacements will be compared with the displacement solutions. As shown in Figure 6, even with a very preliminary Kelvin model, we can still recreate the hysteresis loop of load and displacement for the sucker rod-pump system. In the simulation, the magnification factor is 85 for both Maxwell and Kelvin models to match the actual motion of the sucker rod with the viscoelastic model in the simulation. Notice that the load $F\left(t\right)$ used in the viscoelastic model is the exact load measurement recorded in TAM software 1.6 examples. In general, by reducing the parameter m, the peak of the displacement will decrease and shift to the left, which corresponds to the increase of the natural frequency. In addition, the endpoint will also decrease accordingly. Moreover, by reducing the parameter $k_o$, the peak of the displacement will increase, and the endpoint will also be elevated accordingly. By increasing the parameter $c_1$, the viscous damping is increased, and both the displacement peak and the ending points will be increased. Finally, by increasing the stiffness parameter $k_1$, both the displacement peak and the end displacement tend to decrease. Notice that here, the magnification factor C, along with the parameters m, $c_1$, and $k_1$, are chosen with intuition through trial-and-error methods and human interventions, as reported in ref. [34]. However, in this paper, with the initial guesses of these parameters, we will implement a so-called Inversed Optimization Method to identify the optimal set of parameters with iterations in the steepest descents. With this Newton–Raphson iteration-based approach, these parameters will be searched and improved automatically. A subsequent Kelvin viscoelastic vibration system with adjusted parameters m, $c_1$, and $k_1$ demonstrates a much closer displacement response with the measured force $F\left(t\right)$ as shown in Figure 6, which confirms the potential of inverse optimization approaches with the Newton–Raphson iterations.

3. Inverse Optimization

In engineering practice, many complex systems are impossible to characterize with a simple physical and mathematical model; therefore, an implicit matrix-free iterative method is very useful in providing better guidance to the optimum operation conditions with only the availability of the input and output data [36]. In operation research, to accomplish an optimal route or desired outcome, so-called inverse optimization approaches are introduced to identify provided key parameters [37]. Similarly, in robotics, to reach a particular position, multiple joint angles can be identified with the same Newton–Raphson iterative schemes for the solution of nonlinear position equations governing the desirable positions of the robotic arm along with the unknown joint angles [38]. The inverse optimization is based on the Newton–Raphson iterations, which can be used to search more strategically for the optimal system parameters. Since engineers function within an approximated reality, this inverse optimization approach can be used to search without knowing the true relationship between the optimal material compositions to achieve the best mechanical properties or any type of property. Of course, similar acceleration methods commonly employed in the Newton–Raphson iterations can also be available for the design of more efficient and purposeful search processes.

In the inverse optimization approaches introduced in this paper, in order to identify these intrinsic viscoelastic properties by considering different physical units, we also introduce a scaling factor C. We set up an inverse engineering problem to minimize the difference measured by the following variational indicator

$${\displaystyle E=\sum _{n=1}^N\frac{1}{2}{(u_n-{\overline{u}}_n)}^2,}$$

(14)

where N represents the number of time stations, $u_n$ is the displacement evaluated with the viscoelastic model with the stiffness k, m, $c_o$, and $k_o$, whereas ${\overline{u}}_n$ is the experimental measurement of the displacement in question at the same time intervals. For the minimized error, cost function, or variational indicator E, we have

$${\displaystyle \frac{\partial E}{\partial k_o}=0,\frac{\partial E}{\partial m}=0,\frac{\partial E}{\partial c_1}=0,\frac{\partial E}{\partial k_1}=0.}$$

(15)

This estimate evaluated with finite difference schemes will be directly coupled with the system parameters, for example, k, m, $c_o$, and $k_o$,

$$\begin{array}{ccccccccc}\hfill {\displaystyle f_1}& =& {\displaystyle \frac{\partial E}{\partial x_1}}\hfill & =& {\displaystyle \frac{\partial E}{\partial k_o}}\hfill & =& {\displaystyle \sum _{n=1}^N(u_n-{\overline{u}}_n)\frac{\partial u_n}{\partial k_o}}\hfill & =& 0,\hfill \\ \hfill {\displaystyle f_2}& =& {\displaystyle \frac{\partial E}{\partial x_2}}\hfill & =& {\displaystyle \frac{\partial E}{\partial m}}\hfill & =& {\displaystyle \sum _{n=1}^N(u_n-{\overline{u}}_n)\frac{\partial u_n}{\partial m}}\hfill & =& 0,\hfill \\ \hfill {\displaystyle f_3}& =& {\displaystyle \frac{\partial E}{\partial x_3}}\hfill & =& {\displaystyle \frac{\partial E}{\partial c_1}}\hfill & =& {\displaystyle \sum _{n=1}^N(u_n-{\overline{u}}_n)\frac{\partial u_n}{\partial c_1}}\hfill & =& 0,\hfill \\ \hfill {\displaystyle f_4}& =& {\displaystyle \frac{\partial E}{\partial x_4}}\hfill & =& {\displaystyle \frac{\partial E}{\partial k_1}}\hfill & =& {\displaystyle \sum _{n=1}^N(u_n-{\overline{u}}_n)\frac{\partial u_n}{\partial k_1}}\hfill & =& 0.\hfill \end{array}$$

(16)

We must introduce the Newton–Raphson iterative procedures to obtain the solution of the nonlinear and implicit set of equations. In general, the Newton–Raphson iterative method should be used for this type of nonlinear set of equations. The nonlinear and implicit governing equation about the unknown vector $\mathbf{x}=<k_o,m,c_1,k_1>$, can be rewritten as

$$\mathbf{f}\left(\mathbf{x}\right)=\mathbf{R},$$

(17)

where the given right-hand side vector $\mathbf{R}=<0,0,0,0>$.

It is very important to point out that the initial guess must be fairly close to the actual solution for the unknown $\mathbf{x}$ to ensure the convergence of the Newton–Raphson iteration scheme. In practice, we often start with a few tryouts and narrow down the true solution neighborhood. With an educated guess of the initial set of parameters ${\mathbf{x}}^o=<k_o^o,m^o,c_1^o,k_1^o>$, not too far from the converged solution, the Jacobian matrix can be defined and evaluated. Assume we have all the information before the kth iteration, namely ${\mathbf{x}}^k$ and the corresponding $\mathbf{f}\left({\mathbf{x}}^{k-1}\right)$ as well as the so-called Jacobian matrix $\mathbf{J}\left({\mathbf{x}}^{k-1}\right)$ with all entities $J_{ij}$ defined as

$$\begin{array}{ccccc}\hfill {\displaystyle J_{11}}& =& {\displaystyle \frac{\partial f_1}{\partial x_1}}\hfill & =& {\displaystyle \sum _{n=1}^N(u_n-{\overline{u}}_n)\frac{{\partial }^2{u}_n}{\partial k_o^2}+\sum _{n=1}^N\frac{\partial u_n}{\partial k_o}\frac{\partial u_n}{\partial k_o},}\hfill \\ \hfill {\displaystyle J_{22}}& =& {\displaystyle \frac{\partial f_2}{\partial x_2}}\hfill & =& {\displaystyle \sum _{n=1}^N(u_n-{\overline{u}}_n)\frac{{\partial }^2{u}_n}{\partial m^2}+\sum _{n=1}^N\frac{\partial u_n}{\partial m}\frac{\partial u_n}{\partial m},}\hfill \\ \hfill {\displaystyle J_{33}}& =& {\displaystyle \frac{\partial f_3}{\partial x_3}}\hfill & =& {\displaystyle \sum _{n=1}^N(u_n-{\overline{u}}_n)\frac{{\partial }^2{u}_n}{\partial c_1^2}+\sum _{n=1}^N\frac{\partial u_n}{\partial c_1}\frac{\partial u_n}{\partial c_1},}\hfill \\ \hfill {\displaystyle J_{44}}& =& {\displaystyle \frac{\partial f_4}{\partial x_4}}\hfill & =& {\displaystyle \sum _{n=1}^N(u_n-{\overline{u}}_n)\frac{{\partial }^2{u}_n}{\partial k_1^2}+\sum _{n=1}^N\frac{\partial u_n}{\partial k_1}\frac{\partial u_n}{\partial k_1}.}\hfill \end{array}$$

Furthermore, the off-diagonal terms for the Jacobian matrix with $i\ne j$ can also be elaborated as

$${\displaystyle J_{ij}=\frac{\partial f_i}{\partial x_j}=\sum _{n=1}^N(u_n-{\overline{u}}_n)\frac{{\partial }^2{u}_n}{\partial x_i\partial x_j}+\sum _{n=1}^N\frac{\partial u_n}{\partial x_i}\frac{\partial u_n}{\partial x_j}.}$$

(18)

It is straightforward to confirm that the Jacobian matrix is indeed symmetric. After the solution of the following incremental linear system of equations for the unknown $\Delta {\mathbf{x}}^k$

$$\mathbf{J}\left({\mathbf{x}}^{k-1}\right)\Delta {\mathbf{x}}^k=-\mathbf{f}\left({\mathbf{x}}^{k-1}\right),$$

(19)

with the following update,

$${\mathbf{x}}^k={\mathbf{x}}^{k-1}+\Delta {\mathbf{x}}^k.$$

(20)

The iteration identified in Equation (20) will stop with the relative incremental error $ϵ$ smaller than prescribed small number ${ϵ}_o$,

$${\displaystyle \frac{\parallel \Delta {\mathbf{x}}^k\parallel }{\parallel {\mathbf{x}}^o\parallel }=ϵ\le {ϵ}_o.}$$

(21)

In the actual implementation, since we do not have the analytical expressions for ${\displaystyle \frac{\partial u_n}{\partial x_i}}$, ${\displaystyle \frac{\partial u_n}{\partial x_i}\frac{\partial u_n}{\partial x_j}}$, ${\displaystyle \frac{{\partial }^2{u}_n}{\partial x_i^2}}$, and ${\displaystyle \frac{{\partial }^2{u}_n}{\partial x_i\partial x_j}}$, a central-difference-based numerical scheme is introduced. Assuming the increment $\Delta \mathbf{x}$, we will complete the numerical integration of the dynamic response of the viscoelastic model for three sets of parameters, namely $\mathbf{x}-\Delta \mathbf{x}$, $\mathbf{x}$, $\mathbf{x}+\Delta \mathbf{x}$, and obtain the numerically the approximations as follows

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial u_n}{\partial x_i}}& =& {\displaystyle \frac{u_n(x_i+\Delta x_i)-u_n(x_i-\Delta x_i)}{2\Delta x_i},}\hfill \\ \hfill {\displaystyle \frac{\partial u_n}{\partial x_i}\frac{\partial u_n}{\partial x_j}}& =& {\displaystyle \frac{u_n(x_i+\Delta x_i)-u_n(x_i-\Delta x_i)}{2\Delta x_i}\frac{u_n(x_j+\Delta x_j)-u_n(x_j-\Delta x_j)}{2\Delta x_j},}\hfill \\ \hfill {\displaystyle \frac{{\partial }^2{u}_n}{\partial x_i^2}}& =& {\displaystyle \frac{u_n(x_i+\Delta x_i)-2u_n\left(x_i\right)+u_n(x_i-\Delta x_i)}{\Delta x_i^2}.}\hfill \end{array}$$

(22)

As the most complicated term, the second partial derivative with respect to $x_i$ and $x_j$ with $i\ne j$, namely ${\displaystyle \frac{{\partial }^{2}g}{\partial x_i{x}_j}}$ is evaluated as follows

$${\displaystyle \frac{g(x_i+\Delta x_i,x_j+\Delta x_j)+g(x_i-\Delta x_i,x_j-\Delta x_j)-g(x_i-\Delta x_i,x_j+\Delta x_j)-g(x_i+\Delta x_i,x_j-\Delta x_j)}{4\Delta x_i\Delta x_j}}$$

In the implementation, we can always use a sufficiently small increment $\Delta \mathbf{x}$ for the partial derivatives of $u_n$ as a function of $\mathbf{x}$.

4. Implementation and Improvement

To verify this proposed inverse optimization procedure, we start with a simple first-order differential equation with one parameter, a so-called relaxation time $\tau$, expressed as

$${\displaystyle \frac{du}{dt}+\frac{u}{\tau }=f,}$$

(23)

where $f\left(t\right)$ is a time-dependent inhomogeneous term.

Instead of the actual experimental measurement, in the validation process, we employ the actual analytical solution $\overline{u}$, expressed as

$${\displaystyle \overline{u}\left(t\right)=u_o{e}^{-t/\tau }+{\int }_0^t{e}^{-(t-s)/\tau }f\left(s\right)ds,}$$

(24)

with the initial solution $u\left(0\right)=u_o$ and the convolution term with the so-called Green’s function for impulse ${\displaystyle e^{-t/\tau }.}$

Please note that although we have the analytical solution as a continuous function, in engineering practice, the experimental measurements often exist in discrete form. Therefore, a proper number of time steps N must be taken before we start the inverse optimization process for the search of the system parameter, which yields the solution ${\displaystyle {\overline{u}}_n}$ with $n=1,\cdots ,N.$ In the initial implementation, the maximum number of iterations is set to be 200, and the time station number is 2001. In addition, the finite difference approximation of some of the implicit terms of the Jacobian matrix is based on the $0.1\%$ perturbation of the parameter $\tau$ and the scaling factor C. The analytical solution is based on the relaxation time $\tau =5$. As shown in

Table 1. Newton–Raphson iteration convergence of the relaxation time $\tau$, the scaling factor C, and the relative incremental error $ϵ$ as defined in Equation (21) with the logarithmic modification.

No.	$\mathit{\tau }$	C	$\mathit{ϵ}$
24	4.113353577899	1.108643414669	0.8481099309715
25	4.808147443754	1.012180228335	0.3507291199584
26	4.991132754710	1.001324920436	0.0916535074825
27	5.004938782496	1.000004471816	0.0069345149036
28	5.005001665203	1.000000000128	0.0000315207506
29	5.005001666667	1.000000000000	0.0000000007343

, the only system parameter, namely the relaxation time $\tau$, does converge to a value very close to 5 within approximately $0.1\%$ range. Moreover, the scaling factor C is also simultaneously converged to 1. The quadratic convergence rate towards the end of the iteration processes is evident in Table 1.

It is well accepted that there is no guarantee for the Newton–Raphson iteration to converge, particularly when the initial guess is very far from the desired solution. Moreover, it is often possible that different initial solutions will lead to different branches of solutions [35,39]. Different acceleration procedures have been proposed to improve the convergence behaviors of these nonlinear iterative solution procedures, such as line search methods [40]. In this paper, based on the understanding of the possible solution key features, we can also start the optimization procedure from the very beginning. For example, for the first-order ordinary differential equation (ODE) exponential solutions, rather than connecting with the solution itself, we can choose to connect with the logarithmic of the solution. In this way, the nonlinear system can be less challenging, which can fundamentally change the basin of convergence in the Newton–Raphson iterative procedures. Hence, we set up an inverse engineering problem to minimize the difference measured by the following variational indicator

$${\displaystyle \overline{E}=\sum _{n=1}^N\frac{1}{2}{(D+ln{u}_n-ln{\overline{u}}_n)}^2,}$$

(25)

where N represents the number of time stations, $u_n$ is the displacement evaluated with the system parameter, in this case, the relaxation time $\tau$, ${\overline{u}}_n$ is the experimental measurement of the displacement in question at the same time intervals, and D stands for the new constant which is related to the initial scaling factor C with $D=lnC$.

Similarly, to the minimized error, cost function, or variational indicator E, we have

$${\displaystyle \frac{\partial \overline{E}}{\partial D}=0,\frac{\partial \overline{E}}{\partial \tau }=0.}$$

(26)

Please note that the scaling constant D still has an explicit expression. The derivative with respect to D will be evaluated directly, unlike the system parameter $\tau$, which must be evaluated with finite difference schemes. This estimate of the scaling factor D will be directly coupled with the system parameter and form the unknown vector with $\mathbf{x}=<\tau ,D>$, with the nonlinear set of equations expressed as

$$\begin{array}{ccccccccc}\hfill {\displaystyle {\overline{f}}_1}& =& {\displaystyle \frac{\partial \overline{E}}{\partial x_1}}\hfill & =& {\displaystyle \frac{\partial \overline{E}}{\partial \tau }}\hfill & =& {\displaystyle \sum _{n=1}^N(D+ln{u}_n-ln{\overline{u}}_n)\frac{\partial ln{u}_n}{\partial \tau }}\hfill & =& 0,\hfill \\ \hfill {\displaystyle {\overline{f}}_2}& =& {\displaystyle \frac{\partial \overline{E}}{\partial x_2}}\hfill & =& {\displaystyle \frac{\partial \overline{E}}{\partial D}}\hfill & =& {\displaystyle \sum _{n=1}^N(D+ln{u}_n-ln{\overline{u}}_n)=0.}\hfill \end{array}$$

(27)

To effectively carry out the Newton–Raphson iterative procedures for the nonlinear and implicit governing equation about the unknown vector $\mathbf{x}$, the so-called Jacobian matrix $\mathbf{J}$ must be properly evaluated,

$$\begin{array}{ccccc}\hfill {\displaystyle {\overline{J}}_{11}}& =& {\displaystyle \frac{\partial {\overline{f}}_1}{\partial x_1}}\hfill & =& {\displaystyle \sum _{n=1}^N(D+ln{u}_n-ln{\overline{u}}_n)\frac{{\partial }^{2}ln{u}_n}{\partial {\tau }^2}+\sum _{n=1}^N\frac{\partial ln{u}_n}{\partial \tau }\frac{\partial ln{u}_n}{\partial \tau },}\hfill \\ \hfill {\displaystyle {\overline{J}}_{22}}& =& {\displaystyle \frac{\partial {\overline{f}}_2}{\partial x_2}}\hfill & =& {\displaystyle \sum _{n=1}^{N}1.}\hfill \end{array}$$

Furthermore, the off-diagonal terms for the Jacobian matrix can also be elaborated as

$${\displaystyle {\overline{J}}_{12}=\frac{\partial {\overline{f}}_1}{\partial x_2}=\sum _{n=1}^N\frac{\partial ln{u}_n}{\partial \tau }=\frac{\partial {\overline{f}}_2}{\partial x_1}={\overline{J}}_{21}.}$$

Again, it is straightforward to confirm that the Jacobian matrix is indeed symmetric. In the actual implementation, since we do not the analytical expressions for ${\displaystyle \frac{\partial ln{u}_n}{\partial \tau }}$, ${\displaystyle \frac{\partial ln{u}_n}{\partial \tau }\frac{\partial ln{u}_n}{\partial \tau }}$, and ${\displaystyle \frac{{\partial }^{2}ln{u}_n}{\partial {\tau }^2}}$, central difference schemes similar to Equation (22) will be employed by replacing $u_n$ with $ln{u}_n$. It has been discovered that the range of the system parameter $\tau$ is much wider for the Newton–Raphson iteration procedures with the logarithmic modification. In

Table 2. Newton–Raphson iteration convergence of the relaxation time $\tau$, the scaling factor D, and the relative incremental error $ϵ$ as defined in Equation (21) with the logarithmic modification.

No.	$\mathit{\tau }$	D	$\mathit{ϵ}$
7	5.3502492856171	0.1709734105086	0.3327036938013
8	5.7906235191341	0.0500242875435	0.2283408394767
9	5.9764971771645	0.0062271928244	0.0954819384169
10	6.0044277938157	0.0001228805887	0.0142949464357
11	6.0050011491776	0.0000000513466	0.0002931822785
12	6.0050013888888	0.0000000000000	0.0000001225744

, the converged solutions for the parameter $\tau =6$ and the scaling factor $D=0$, which is equivalent to $C=1$ are listed with a similar quadratic rate near the solution and $0.1\%$ system error due to the finite difference approximations with the perturbation of about $0.1\%$ of the system parameter.

The same implementation has also been extended to a more complicated second-order differential equation with three parameters, namely the mass m, the damping c, and the stiffness k, expressed as

$${\displaystyle m\frac{d^{2}u}{d{t}^2}+c\frac{du}{dt}+ku=f,}$$

(28)

where $f\left(t\right)$ is again a time-dependent inhomogeneous term.

Denote the natural frequency ${\omega }_o$ with ${\displaystyle {\omega }_o^2=k/m}$, the damping ratio $\zeta$ with $c/m=2\zeta {\omega }_o,$ and the damped natural frequency ${\displaystyle {\omega }_d=\sqrt{1-{\zeta }^2}{\omega }_o.}$ Instead of the actual experimental measurement, in the validation process, we employ the actual analytical solution $\overline{u}$, expressed as

$${\displaystyle \overline{u}\left(t\right)=e^{-\zeta {\omega }_{o}t}(u_{o}cos{\omega }_{d}t+\frac{v_o+\zeta {\omega }_o{u}_o}{{\omega }_d}sin{\omega }_{d}t)+{\int }_0^t{e}^{-\zeta {\omega }_o(t-s)}\frac{f\left(s\right)}{m{\omega }_d}sin{\omega }_d(t-s)ds,}$$

(29)

with the initial displacement $u\left(0\right)=u_o$ and the initial velocity $v\left(0\right)=v_o$ along with the convolution term with the so-called Green’s function for impulse $e^{-\zeta \omega t}\frac{1}{m{\omega }_d}sin{\omega }_{d}t.$

In the second implementation, the maximum number of iterations is again set to be 200, and the time station number is 2001. In addition, the finite difference approximation of some of the implicit terms of the Jacobian matrix is based on the $0.1\%$ perturbation of the parameters m, c, and k along with the scaling factor C. The analytical solution is based on the set of parameters $m=9.88$, $c=4.94$, and $k=247$. As shown in

Table 3. Newton–Raphson iteration convergence of the system parameters m, c, and k, the scaling factor C, and the relative incremental error $ϵ$ as defined in Equation (21).

No.	m	c	k	C	$\mathit{ϵ}$
2	14.0088697014	38.8505562741	252.2630700	4.4695001143	7.5404751904
3	9.8962841404	14.3998130515	260.9769377	1.4358822096	0.8809455612
4	9.8402652289	7.0242940976	247.2775516	1.0914951778	0.5182224546
5	9.8739033983	5.2106554463	247.1593685	1.0123649473	0.0605886943
6	9.8798105710	4.9450642071	246.9998836	1.0002193760	0.0103258241
7	9.8799373442	4.9399706629	246.9979158	1.0000000949	0.0001820245
8	9.8799373906	4.9399686098	246.9979147	1.0000000016	0.0000000768

, the system parameters, namely m, c, and k, do converge to the set of values very close to the true system parameters, again within approximately $0.1\%$ range. Moreover, the scaling factor C is also simultaneously converged to 1. The quadratic convergence rate towards the end of the iteration processes, when the solutions are very close to the targets, is again evident in Table 3, which also indicates proper programming and implementation of the Newton–Raphson iteration procedures.

Finally, we use the experimentally measured data from V11, which is well documented in Echometer Co and utilized in ref. [34] as the analytical solution of the equivalent viscoelastic vibration system modeled with the Kelvin model. With the same inverse optimization procedures, we obtain the system parameters with $k_o=4.4668\mathrm{lbf}/\mathrm{in}$, $m=3.3063\mathrm{lbf}$·$s/{\mathrm{in}}^2$, $c_1=2.4585\mathrm{lbf}$·$s/\mathrm{in}$, and $k_1=3.4181\mathrm{lbf}/\mathrm{in}$, and the scaling factor $C=12.0113$. The displacement and force/displacement loop match very well with the experimental data documented by Echometer Co., as shown in Figure 6, Figure 7, Figure 8 and Figure 9.

To validate the inverse optimization approach further and to test our hypothesis that the entire downhole sucker rod-pump system can be approximated with the viscoelastic model, we employ the same procedures for another set of dynamometer measurements taken from the same V11 well at 7:14:19 a.m. on 7 August 2019 Stroke 17, the exact load measurement recorded in TAM software examples, and obtain a different but similar set of system parameters with $k_o=4.9006\mathrm{lbf}/\mathrm{in}$, $m=3.7094\mathrm{lbf}$·$s/{\mathrm{in}}^2$, $c_1=2.8966\mathrm{lbf}$·$s/\mathrm{in}$, and $k_1=3.9594\mathrm{lbf}/\mathrm{in}$, and the scaling factor $C=10.8569$.

The convergence information is listed in

Table 4. Newton–Raphson iteration convergence of the system parameters $k_o$, m, $c_1$, and $k_1$, the scaling factor C, and the relative incremental error $ϵ$ as defined in Equation (21).

No.	${\mathit{k}}_{\mathit{o}}$	m	${\mathit{c}}_1$	${\mathit{k}}_1$	C	$\mathit{ϵ}$
10	4.8145200	3.7227113	3.0371959	4.0113022	9.9459747	0.0583832
11	4.8157944	3.7250974	3.0323557	4.0186058	9.9727830	0.0336033
12	4.8769304	3.71661320	2.9368406	3.9800812	10.5630007	0.0518076
13	4.8762824	3.7159651	2.9388687	3.9795838	10.5762036	0.0011525
14	4.8948972	3.7115266	2.9074872	3.9649666	10.7825937	0.0180758
15	4.8951385	3.7107580	2.9057328	3.9636912	10.7948712	0.0010746
16	4.8989279	3.7098281	2.8993963	3.9606940	10.8378530	0.0037609
17	4.8999700	3.7095979	2.8977176	3.9599161	10.8495151	0.0010198
18	4.9004231	3.7094863	2.8969610	3.9595585	10.8546888	0.0004526
19	4.9005655	3.7094541	2.8967290	3.9594508	10.8562851	0.0001396
20	4.9006119	3.7094470	2.8966611	3.9594212	10.8567655	0.0000420
21	4.9006292	3.7094476	2.8966428	3.9594152	10.8569094	0.0000126

. The viscoelastic model does not match exactly with the sucker rod-pump system. In addition, unlike the previous mathematical models from which the numerical or analytical solutions are derived, the data from the experiments are heuristic. Thus, the quadratic convergence properties of the Newton–Raphson iterative procedures are replaced with an oscillatory convergence with a rate much slower than the typical quadratic one observed for the Newton–Raphson iteration when the iterative solutions are within the neighborhood. Again, since the numerical evaluation of some of the Jacobian matrix entities has a truncation error around $O\left({10}^{-3}\right),$ the converged system parameters will hover around that relative error range.

Again, the displacement and force/displacement loop match very well with the experimental data documented by Echometer Co, as shown in Figure 6. Interestingly, the same system parameters have been applied to Stroke 18 and demonstrate a remarkable match with the actual experimental measurement.

5. Conclusions

In the petroleum industry, it is very difficult to quantify the dynamical behaviors of the downhole sucker rod-pump systems. The proposed inverse optimization method and the viscoelastic vibration model shed light on the understanding of the reciprocal nature of the intricate relationship between the displacement of the polish rod and the total force exerted through the pump jack. The optimal parameters of the proposed viscoelastic model for the downhole sucker rod-pump system yield almost identical results in comparison with the experimental measurements in the oil field. With this confirmation, the same approach will be implemented for the delivery system of a very viscous non-Newtonian fluid in EV manufacturing, as well as other complex engineering systems. It is promising that even one or two parameters can be utilized to approximate the complicated relationship between the required pressure drop and the flow rate with non-Newtonian internal fluid modeled with a power law distribution and bypass fully fledge computational models, which normally require a lot of resources and experience.