Pump System Model Parameter Identification Based on Experimental and Simulation Data

Wang, Sheldon; Gao, Dalong; Wester, Alexandria; Beaver, Kalyb; Edwards, Shanae; Taylor, Carrie Anne

doi:10.3390/fluids9060136

Open AccessArticle

Pump System Model Parameter Identification Based on Experimental and Simulation Data

by

Sheldon Wang

^1,*

,

Dalong Gao

²,

Alexandria Wester

¹,

Kalyb Beaver

¹,

Shanae Edwards

¹ and

Carrie Anne Taylor

³

¹

McCoy College of Science, Mathematics & Engineering, Midwestern State University, a Member of the Texas Tech University System, Wichita Falls, TX 76308, USA

²

Materials & Manufacturing Systems Research Laboratory, GM R&D, 30470 Harley Earl Blvd., Warren, MI 48092, USA

³

Echometer Company, 5001 Ditto Ln, Wichita Falls, TX 76302, USA

^*

Author to whom correspondence should be addressed.

Fluids 2024, 9(6), 136; https://doi.org/10.3390/fluids9060136

Submission received: 17 April 2024 / Revised: 21 May 2024 / Accepted: 30 May 2024 / Published: 4 June 2024

Download

Browse Figures

Versions Notes

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.

Keywords:

inverse; optimization; viscoelastic; computation

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

δ

, namely

C = 2 δ = 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 (t) = u_{o} H (t)

where

H (t)

is the Heaviside function. The total force function

f (t)

can be expressed as the combination of the location force

f_{o} (t) = K_{o} u (t)

and

f_{1} (t)

for the stiffness and dashpot couple

(K_{1}, C_{1})

with the relaxation time

τ_{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{matrix} f_{o} & = & K_{o} u (t), \end{matrix}

(1)

\begin{matrix} \dot{u} (t) & = & \frac{f_{1}}{C_{1}} + \frac{{\dot{f}}_{1}}{K_{1}} . \end{matrix}

(2)

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

e^{t / τ_{1}} f_{1} (t) - 0 = \int_{- \infty}^{t} K_{1} e^{t / τ_{1}} \dot{u} (t) d t,

(3)

namely using integration by parts,

e^{t / τ_{1}} f_{1} (t) = K_{1} e^{t / τ_{1}} u (t) - \int_{- \infty}^{t} K_{1} / τ_{1} u (t) e^{t / τ_{1}} d t = K_{1} e^{t / τ_{1}} u_{o} - \int_{0}^{t} K_{1} / τ_{1} u_{o} e^{t / τ_{1}} d t .

(4)

Finally, we have

f_{1} (t) = K_{1} u_{o} e^{- t / τ_{1}}

and the total force

f (t) = f_{1} (t) + f_{o} (t) = (K_{o} + K_{1} e^{- t / τ_{1}}) u_{o}

for the Maxwell model. Likewise, denote the sudden load as

f (t) = f_{o} H (t)

where

H (t)

is a Heaviside function, the total displacement function

u (t)

can be expressed as the combination of the location displacement

u_{o} (t) = f (t) / K_{o}

and

u_{1} (t)

for the stiffness and dashpot couple

(K_{1}, C_{1})

with the relaxation time

τ_{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{matrix} f (t) & = & K_{o} u_{o} (t), \end{matrix}

(5)

\begin{matrix} f (t) & = & C_{1} {\dot{u}}_{1} + K_{1} u_{1}, 1 \leq i \leq n . \end{matrix}

(6)

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

e^{t / τ_{1}} u_{1} (t) - 0 = \int_{- \infty}^{t} e^{t / τ_{1}} f / C_{1} d t = \int_{0}^{t} e^{t / τ_{1}} f_{o} / C_{1} d t .

(7)

Finally, we have the local displacement

u_{1} (t) = f_{o} / K_{1} (1 - e^{- t / τ_{1}})

and the total displacement function can be expressed as

u (t) = u_{o} + u_{1} = f_{o} / K_{o} + f_{o} / K_{1} (1 - e^{- t / τ_{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 (t)

. 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} (t)

, whereas the displacement of the stiffness

k_{o}

is denoted as

u_{1} (t)

. Since the mass m is connected with the stiffness

k_{o}

in series, the total displacement of the mass

u (t)

is a combination of the two displacements

u_{1} (t)

and

u_{2} (t)

. In general, the external load

F (t)

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

m g

.

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

\begin{matrix} k_{o} u_{1} (t) & = & k_{1} u_{2} + c_{1} {\dot{u}}_{2}, \\ k_{o} u_{1} (t) & = & F (t) - m \ddot{u} . \end{matrix}

(8)

Using the kinematic relationship

u (t) = u_{1} (t) + u_{2} (t)

, we obtain the following third-order governing equation for

u_{2} (t)

,

\frac{m c_{1}}{k_{o}} {\overset{⃛}{u}}_{2} + \frac{m}{k_{o}} (k_{o} + k_{1}) {\ddot{u}}_{2} + c_{1} {\dot{u}}_{2} + k_{1} u_{2} = F (t) .

(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 (t)

is introduced to derive the closest

u_{2} (t)

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 (t) .

(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 (t) = u_{2} (t)

. Finally, we recover the familiar spring-mass-dashpot vibration system

m \ddot{u} + c_{1} \dot{u} + k_{1} u = F (t) .

Finally, in order to facilitate the solution, introduce the state variable

y = < u_{2}, {\dot{u}}_{2}, {\ddot{u}}_{2} >

, we can rewrite the third-order viscoelastic vibration system with a dynamical system format,

A = [\begin{matrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ - \frac{k_{o} k_{1}}{m c_{1}} & - \frac{k_{o}}{m} & - \frac{k_{o} + k_{1}}{c_{1}} \end{matrix}] and f = (\begin{matrix} 0 \\ 0 \\ \frac{k_{o} F}{m c_{1}} \end{matrix}) .

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 (t)

, whereas the displacement of the stiffness

k_{1}

is denoted as

u_{2} (t)

and the dashpot

c_{1}

as

u_{1} (t)

. Please note that the displacement of the stiffness

k_{o}

is

u (t)

. Thus, the total displacement of the mass

u (t)

is a combination of the two displacements

u_{1} (t)

and

u_{2} (t)

within one of the parallel branches. In general, the external load

F (t)

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

m g

. 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{matrix} F - m \ddot{u} & = & k_{o} u + F_{o}, \\ \dot{u} & = & \frac{{\dot{F}}_{o}}{k_{1}} + \frac{F_{o}}{c_{1}} . \end{matrix}

(11)

Using the kinematic relationship

u (t) = u_{1} (t) + u_{2} (t)

as well as

\dot{u} (t) = {\dot{u}}_{1} (t) + {\dot{u}}_{2} (t)

, we obtain the following third-order governing equation for

u (t)

,

\frac{m}{k_{1}} {\overset{⃛}{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} + k u = F (t) .

Likewise, as the infinitely stiff spring

k_{1}

, namely

k_{1} \to + \infty

, consequently,

u_{2} \to 0

and

u (t) = u_{1} (t),

we recover the familiar spring-mass-dashpot vibration system

m \ddot{u} + c_{1} \dot{u} + k_{o} u = F (t) .

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

\dot{y} = A y + f,

(13)

with the state variable

y = < u_{2}, {\dot{u}}_{2}, {\ddot{u}}_{2} >

, and

A = [\begin{matrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ - \frac{k_{o} k_{1}}{m c_{1}} & - \frac{k_{o} + k_{1}}{m} & - \frac{k_{1}}{c_{1}} \end{matrix}] and f = (\begin{matrix} 0 \\ 0 \\ \frac{\dot{F}}{m} + \frac{k_{1} F}{m c_{1}} \end{matrix}) .

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 k g

,

c_{1} = 0.2 N s / m

,

k_{1} = 4 N / m

, and

W_{o} = 10 N

. The external force

F (t)

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 (t)

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 (t)

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

E = \sum_{n = 1}^{N} \frac{1}{2} {(u_{n} - {\bar{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

{\bar{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

\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{matrix} f_{1} & = & \frac{\partial E}{\partial x_{1}} & = & \frac{\partial E}{\partial k_{o}} & = & \sum_{n = 1}^{N} (u_{n} - {\bar{u}}_{n}) \frac{\partial u_{n}}{\partial k_{o}} & = & 0, \\ f_{2} & = & \frac{\partial E}{\partial x_{2}} & = & \frac{\partial E}{\partial m} & = & \sum_{n = 1}^{N} (u_{n} - {\bar{u}}_{n}) \frac{\partial u_{n}}{\partial m} & = & 0, \\ f_{3} & = & \frac{\partial E}{\partial x_{3}} & = & \frac{\partial E}{\partial c_{1}} & = & \sum_{n = 1}^{N} (u_{n} - {\bar{u}}_{n}) \frac{\partial u_{n}}{\partial c_{1}} & = & 0, \\ f_{4} & = & \frac{\partial E}{\partial x_{4}} & = & \frac{\partial E}{\partial k_{1}} & = & \sum_{n = 1}^{N} (u_{n} - {\bar{u}}_{n}) \frac{\partial u_{n}}{\partial k_{1}} & = & 0 . \end{matrix}

(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

x = < k_{o}, m, c_{1}, k_{1} >

, can be rewritten as

f (x) = R,

(17)

where the given right-hand side vector

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

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

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

x^{k}

and the corresponding

f (x^{k - 1})

as well as the so-called Jacobian matrix

J (x^{k - 1})

with all entities

J_{i j}

defined as

\begin{matrix} J_{11} & = & \frac{\partial f_{1}}{\partial x_{1}} & = & \sum_{n = 1}^{N} (u_{n} - {\bar{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}}, \\ J_{22} & = & \frac{\partial f_{2}}{\partial x_{2}} & = & \sum_{n = 1}^{N} (u_{n} - {\bar{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}, \\ J_{33} & = & \frac{\partial f_{3}}{\partial x_{3}} & = & \sum_{n = 1}^{N} (u_{n} - {\bar{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}}, \\ J_{44} & = & \frac{\partial f_{4}}{\partial x_{4}} & = & \sum_{n = 1}^{N} (u_{n} - {\bar{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}} . \end{matrix}

Furthermore, the off-diagonal terms for the Jacobian matrix with

i \neq j

can also be elaborated as

J_{i j} = \frac{\partial f_{i}}{\partial x_{j}} = \sum_{n = 1}^{N} (u_{n} - {\bar{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

Δ x^{k}

J (x^{k - 1}) Δ x^{k} = - f (x^{k - 1}),

(19)

with the following update,

x^{k} = x^{k - 1} + Δ x^{k} .

(20)

The iteration identified in Equation (20) will stop with the relative incremental error

ϵ

smaller than prescribed small number

ϵ_{o}

,

\frac{∥ Δ x^{k} ∥}{∥ x^{o} ∥} = ϵ \leq ϵ_{o} .

(21)

In the actual implementation, since we do not have the analytical expressions for

\frac{\partial u_{n}}{\partial x_{i}}

,

\frac{\partial u_{n}}{\partial x_{i}} \frac{\partial u_{n}}{\partial x_{j}}

,

\frac{\partial^{2} u_{n}}{\partial x_{i}^{2}}

, and

\frac{\partial^{2} u_{n}}{\partial x_{i} \partial x_{j}}

, a central-difference-based numerical scheme is introduced. Assuming the increment

Δ x

, we will complete the numerical integration of the dynamic response of the viscoelastic model for three sets of parameters, namely

x - Δ x

,

x

,

x + Δ x

, and obtain the numerically the approximations as follows

\begin{matrix} \frac{\partial u_{n}}{\partial x_{i}} & = & \frac{u_{n} (x_{i} + Δ x_{i}) - u_{n} (x_{i} - Δ x_{i})}{2 Δ x_{i}}, \\ \frac{\partial u_{n}}{\partial x_{i}} \frac{\partial u_{n}}{\partial x_{j}} & = & \frac{u_{n} (x_{i} + Δ x_{i}) - u_{n} (x_{i} - Δ x_{i})}{2 Δ x_{i}} \frac{u_{n} (x_{j} + Δ x_{j}) - u_{n} (x_{j} - Δ x_{j})}{2 Δ x_{j}}, \\ \frac{\partial^{2} u_{n}}{\partial x_{i}^{2}} & = & \frac{u_{n} (x_{i} + Δ x_{i}) - 2 u_{n} (x_{i}) + u_{n} (x_{i} - Δ x_{i})}{Δ x_{i}^{2}} . \end{matrix}

(22)

As the most complicated term, the second partial derivative with respect to

x_{i}

and

x_{j}

with

i \neq j

, namely

\frac{\partial^{2} g}{\partial x_{i} x_{j}}

is evaluated as follows

\frac{g (x_{i} + Δ x_{i}, x_{j} + Δ x_{j}) + g (x_{i} - Δ x_{i}, x_{j} - Δ x_{j}) - g (x_{i} - Δ x_{i}, x_{j} + Δ x_{j}) - g (x_{i} + Δ x_{i}, x_{j} - Δ x_{j})}{4 Δ x_{i} Δ x_{j}}

In the implementation, we can always use a sufficiently small increment

Δ x

for the partial derivatives of

u_{n}

as a function of

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

τ

, expressed as

\frac{d u}{d t} + \frac{u}{τ} = f,

(23)

where

f (t)

is a time-dependent inhomogeneous term.

Instead of the actual experimental measurement, in the validation process, we employ the actual analytical solution

\bar{u}

, expressed as

\bar{u} (t) = u_{o} e^{- t / τ} + \int_{0}^{t} e^{- (t - s) / τ} f (s) d s,

(24)

with the initial solution

u (0) = u_{o}

and the convolution term with the so-called Green’s function for impulse

e^{- t / τ} .

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

{\bar{u}}_{n}

with

n = 1, \dots, 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

τ

and the scaling factor C. The analytical solution is based on the relaxation time

τ = 5

. As shown in Table 1, the only system parameter, namely the relaxation time

τ

, 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

\bar{E} = \sum_{n = 1}^{N} \frac{1}{2} {(D + ln u_{n} - ln {\bar{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

τ

,

{\bar{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 = ln C

.

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

\frac{\partial \bar{E}}{\partial D} = 0, \frac{\partial \bar{E}}{\partial τ} = 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

τ

, 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

x = < τ, D >

, with the nonlinear set of equations expressed as

\begin{matrix} {\bar{f}}_{1} & = & \frac{\partial \bar{E}}{\partial x_{1}} & = & \frac{\partial \bar{E}}{\partial τ} & = & \sum_{n = 1}^{N} (D + ln u_{n} - ln {\bar{u}}_{n}) \frac{\partial ln u_{n}}{\partial τ} & = & 0, \\ {\bar{f}}_{2} & = & \frac{\partial \bar{E}}{\partial x_{2}} & = & \frac{\partial \bar{E}}{\partial D} & = & \sum_{n = 1}^{N} (D + ln u_{n} - ln {\bar{u}}_{n}) = 0 . \end{matrix}

(27)

To effectively carry out the Newton–Raphson iterative procedures for the nonlinear and implicit governing equation about the unknown vector

x

, the so-called Jacobian matrix

J

must be properly evaluated,

\begin{matrix} {\bar{J}}_{11} & = & \frac{\partial {\bar{f}}_{1}}{\partial x_{1}} & = & \sum_{n = 1}^{N} (D + ln u_{n} - ln {\bar{u}}_{n}) \frac{\partial^{2} ln u_{n}}{\partial τ^{2}} + \sum_{n = 1}^{N} \frac{\partial ln u_{n}}{\partial τ} \frac{\partial ln u_{n}}{\partial τ}, \\ {\bar{J}}_{22} & = & \frac{\partial {\bar{f}}_{2}}{\partial x_{2}} & = & \sum_{n = 1}^{N} 1 . \end{matrix}

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

{\bar{J}}_{12} = \frac{\partial {\bar{f}}_{1}}{\partial x_{2}} = \sum_{n = 1}^{N} \frac{\partial ln u_{n}}{\partial τ} = \frac{\partial {\bar{f}}_{2}}{\partial x_{1}} = {\bar{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

\frac{\partial ln u_{n}}{\partial τ}

,

\frac{\partial ln u_{n}}{\partial τ} \frac{\partial ln u_{n}}{\partial τ}

, and

\frac{\partial^{2} ln u_{n}}{\partial τ^{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

τ

is much wider for the Newton–Raphson iteration procedures with the logarithmic modification. In Table 2, the converged solutions for the parameter

τ = 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

m \frac{d^{2} u}{d t^{2}} + c \frac{d u}{d t} + k u = f,

(28)

where

f (t)

is again a time-dependent inhomogeneous term.

Denote the natural frequency

ω_{o}

with

ω_{o}^{2} = k / m

, the damping ratio

ζ

with

c / m = 2 ζ ω_{o},

and the damped natural frequency

ω_{d} = \sqrt{1 - ζ^{2}} ω_{o} .

Instead of the actual experimental measurement, in the validation process, we employ the actual analytical solution

\bar{u}

, expressed as

\bar{u} (t) = e^{- ζ ω_{o} t} (u_{o} cos ω_{d} t + \frac{v_{o} + ζ ω_{o} u_{o}}{ω_{d}} sin ω_{d} t) + \int_{0}^{t} e^{- ζ ω_{o} (t - s)} \frac{f (s)}{m ω_{d}} sin ω_{d} (t - s) d s,

(29)

with the initial displacement

u (0) = u_{o}

and the initial velocity

v (0) = v_{o}

along with the convolution term with the so-called Green’s function for impulse

e^{- ζ ω t} \frac{1}{m ω_{d}} sin ω_{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, 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 lbf / in

,

m = 3.3063 lbf

·

s / {in}^{2}

,

c_{1} = 2.4585 lbf

·

s / in

, and

k_{1} = 3.4181 lbf / 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 lbf / in

,

m = 3.7094 lbf

·

s / {in}^{2}

,

c_{1} = 2.8966 lbf

·

s / in

, and

k_{1} = 3.9594 lbf / in

, and the scaling factor

C = 10.8569

.

The convergence information is listed in Table 4. 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 (10^{- 3}),

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.

Author Contributions

Conceptualization, S.W. and D.G.; methodology, S.W.; software, S.W., D.G., A.W. and K.B.; validation, S.W. and D.G.; formal analysis, S.W. and D.G.; investigation, S.W., D.G., A.W., K.B., S.E. and C.A.T.; resources, S.W. and D.G.; data curation, S.W. and D.G.; writing—original draft preparation, S.W. and D.G.; writing—review and editing, S.W., D.G., A.W., K.B., S.E. and C.A.T.; visualization, S.W.; supervision, S.W. and D.G.; project administration, S.W. and D.G.; funding acquisition, S.W. and D.G. All authors have read and agreed to the published version of the manuscript.

Funding

This research was partially supported by the GM Research Grant.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Acknowledgments

The financial support from the Midwestern State University, a member of Texas Tech University System and Manufacturing Systems Research Lab, GM R&D, is greatly appreciated.

Conflicts of Interest

Author Dr. Dalong Gao is employed by General Motors (GM). Author Carrie Anne Taylor is employed by Echometer Company. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Nomenclature

$D_{i}$	plunger outer diameter
$D_{o}$	barrel inner diameter
$H (t)$	Heaviside function
C	plunger and barrel clearance
$δ$	plunger and barrel gap with $C = 2 δ$
$L_{p}$	plunger length
$U_{p}$	plunger displacement
U	bridle displacement
m	lumped mass attached to the viscoelastic model
$F (t)$	bridle force measured at the surface
$p_{h}$	top plunger pressure
$p_{l}$	bottom plunger pressure
$k_{1}$	stiffness in viscoelastic stiffness and dashpot couple
$c_{1}$	dashpot in viscoelastic stiffness and dashpot couple
$u_{o}$	imposed displacement
$f_{o}$	imposed force

References

Romero, O.J.; Almeida, P. Numerical simulation of the sucker-rod pumping system. Ingeniería e Investigación 2014, 34, 4–11. [Google Scholar] [CrossRef]
Karhan, M.K.; Nandi, S.; Jadhav, P.B. Design and Optimization of Sucker Rod Pump Using Prosper. Int. J. Interdiscip. Res. Innov. 2015, 3, 108–122. [Google Scholar]
Takacs, G. Improved designs reduce sucker-rod pumping costs. Oil Gas J. 1996, 94, 5. [Google Scholar]
Yu, Y.; Chang, Z.; Qi, Y.; Xue, X.; Zhao, J. Study of a new hydraulic pumping unit based on the offshore platform. Energy Sci. Eng. 2016, 4, 352–360. [Google Scholar] [CrossRef]
Anderson, G.; Liang, B.; Liang, B. The successful application of new technology of oil production in offshore. Foreign Oilfield Eng. 2001, 17, 28–30. [Google Scholar]
Takacs, G. Modern Sucker-Rod Pumping; PennWell Books: Tulsa, OK, USA, 1993. [Google Scholar]
Takacs, G. Sucker-Rod Pumping Handbook—Production Engineering Fundamentals and Long-Stroke Rod Pumping; Elsevier: Amsterdam, The Netherlands, 2015. [Google Scholar]
Wang, S.; Rowlan, L.; Cook, D.; Conrady, C.; King, R.; Taylor, C. Dynamics of pump jacks with theories and experiments. Upstream Oil Gas Technol. 2023, 11, 100097. [Google Scholar] [CrossRef]
Lea, J.F. What’s New in Artificial Lift? Part I. World Oil 2012, 79–85. [Google Scholar]
Lea, J.F. What’s New in Artificial Lift? Part II. World Oil 2012, 85–94. [Google Scholar]
Bommer, P.M.; Podio, A.L. The Beam Lift Handbook; PETEX: Houston, TX, USA, 2012. [Google Scholar]
Bommer, P.M.; Podio, A.L.; Carroll, G. The Measurment of Down Stroke Force in Rod Pumps. In Proceedings of the Sixty Third Annual Meeting of the Southwestern Petroleum Short Course, Lubbock, TX, USA, 20–21 April 2016. [Google Scholar]
Podio, A. Artificial Lift. In Encyclopedia of Life Support Systems; UNESCO-EOLSS: Johannesburg, South Africa, 2013; pp. 1–9. [Google Scholar]
Takacs, G. Sucker Rod Pumping Manual. In Tulsa: PennWell Corporation; PennWell Books: Tulsa, OK, USA, 2003. [Google Scholar]
Dong, Z.; Zhang, M.; Zhang, X.; Pang, X. Study on reasonable choice of electric submersible pump. Acta Pet. Sin. 2008, 29, 128–131. [Google Scholar]
Karpuz-Pickell, P.; Roderick, R. From Failure to success: A Metallurgical Story on Sucker Rod Pump Barrels. In Proceedings of the Sixty Second Annual Meeting of the Southwestern Petroleum Short Course, Lubbock, TX, USA, 22–23 April 2015. [Google Scholar]
Nampoothiri, M.P. Evaluation of the Effectiveness of API-Modified Goodman Diagram in Sucker Rod Fatigue Analysis. Master’s Thesis, Texas Tech University, Lubbock, TX, USA, 2001. [Google Scholar]
Li, Z.; Song, J.; Huang, Y. Design and analysis for a new energy-saving hydraulic pumping unit. Proc. Inst. Mech. Eng. Part C: J. Mech. Eng. Sci. 2017, 232, 2119–2131. [Google Scholar] [CrossRef]
Rowlan, O.L.; McCoy, J.N.; Lea, J.F.; Use of the Pump Slippage Equation to Design Pump Clearances. Priv. Commun. 2012. Available online: https://www.echometer.com/ (accessed on 16 April 2024).
Wang, S.; Rowlan, L.; Elsharafi, M.; Ermila, M.; Grejtak, T.; Taylor, C. On Leakage Issues of Sucker Rod Pumping Systems. ASME J. Fluids Eng. 2019, 141, 111201. [Google Scholar] [CrossRef]
Nouri, J.M.; Whitelaw, J.H. Flow of Newtonian and Non-Newtonian Fluids in a Concentric Annulus with Rotation of the Inner Cylinder. ASME J. Fluids Eng. 1994, 116, 821–827. [Google Scholar] [CrossRef]
Fakher, S.; Khlaifat, A.; Hossain, M.; Nameer, H. A comprehensive review of sucker rod pumps’ components, diagnostics, mathematical models, and common failures and mitigations. J. Pet. Explor. Prod. Technol. 2021, 11, 3815–3839. [Google Scholar] [CrossRef]
Patterson, J.; Chambliss, K.; Rowlan, L.; Curfew, J. Fluid Slippage in Down-Hole Rod-Drawn Oil Well Pumps. In Proceedings of the 54th Southwestern Petroleum Short Course, Lubbock, TX, USA, 25–26 April 2007; Volume 11, pp. 45–59. [Google Scholar]
Patterson, J.; Williams, B. Fluid Slippage in Down-Hole Rod-Drawn Oil Well Pumps. In Proceedings of the 45th Southwestern Petroleum Short Course, Lubbock, TX, USA, 7–9 April 1998; Volume 11, pp. 180–191. [Google Scholar]
Zhao, R.; Zhang, X.; Liu, M.; Shi, J.; Su, L.; Shan, H.; Sun, C.; Miao, G.; Wang, Y.; Shi, L.; et al. Production Optimization and Application of Combined Artificial-Lift Systems in Deep Oil Wells. In SPE Middle East Artificial Lift Conference; SPE-184222-MS; SPE: Kuala Lumpur, Malaysia, 2016. [Google Scholar]
Patterson, J.; Dittman, J.; Curfew, J.; Hill, J.; Brauten, D.; Williams, B. Fluid Slippage in Down-Hole Rod-Drawn Oil Well Pumps. In Proceedings of the 46th Southwestern Petroleum Short Course, Lubbock, TX, USA, 21–22 April 1999; Volume 11, pp. 96–106. [Google Scholar]
Patterson, J.; Curfew, J.; Brock, M.; Brauten, D.; Williams, B. Fluid Slippage in Down-Hole Rod-Drawn Oil Well Pumps. In Proceedings of the 47th Southwestern Petroleum Short Course, Lubbock, TX, USA, 12–13 April 2000; Volume 11, pp. 117–136. [Google Scholar]
Pons, V. Modified Everitt-Jennings: A Complete Methodology for Production Optimization of Sucker Rod Pumped Wells. In Proceedings of the Sixty Second Annual Meeting of the Southwestern Petroleum Short Course, Lubbock, TX, USA, 20–23 April 2015. [Google Scholar]
Copeland, C. Fluid Extractor. In Proceedings of the Sixty Second Annual Meeting of the Southwestern Petroleum Short Course, Lubbock, TX, USA, 20–23 April 2015. [Google Scholar]
Ermila, M. Critical Evaluation of Sucker-Rod String Design Practices in the Hamada Field Libya. Master’s Thesis, University of Miskolc, Miskolc, Hungary, 1999. [Google Scholar]
Wang, S.; Yang, Y.; Wu, T. Model Studies of Fluid-Structure Interaction Problems. Comput. Model. Eng. Sci. 2019, 119, 5–34. [Google Scholar]
Wang, S.; Grejtak, T.; Moody, L. Structural Designs with Considerations of Both Material and Structural Failure. ASCE Pract. Period. Struct. Des. Constr. 2017, 22, 1–8. [Google Scholar] [CrossRef]
Wang, S. Scaling, complexity, and design aspects in computational fluid dynamics. Fluids 2021, 6, 362. [Google Scholar] [CrossRef]
Wang, S.; Rowlan, L.; Henderson, A.; Aleman, S.T.; Creacy, T.; Taylor, C.A. Viscoelastic representation of the operation of sucker rod pumps. Fluids 2021, 7, 70. [Google Scholar] [CrossRef]
Wang, S. Essentials of Mathematical Tools for Engineers; Sentia Publishing: Lakeway, TX, USA, 2022. [Google Scholar]
Chan, T.; Mahmood, R.; Zhu, I.Y. Inverse Optimization: Theory and Applications. Oper. Res. 2021, 1526–5463, 1–29. [Google Scholar] [CrossRef]
Yang, C.; Zhang, J. Two general methods for inverse optimization problems. Appl. Math. Lett. 1999, 12, 69–72. [Google Scholar] [CrossRef]
Lynch, K.M.; Park, F.C. Modern Robotics; Cambridge University Press: Cambridge, UK, 2017. [Google Scholar]
Wang, S. Fundamentals of Fluid-Solid Interactions-Analytical and Computational Approaches; Elsevier Science: Amsterdam, The Netherlands, 2008. [Google Scholar]
Bathe, K. Finite Element Procedures; Prentice Hall: New York, NY, USA, 1996. [Google Scholar]

Figure 1. Surface and downhole sucker rod-pump system models.

Figure 2. Typical viscoelastic models.

Figure 3. Typical responses for Maxwell and Kelvin models.

Figure 4. Kelvin viscoelastic vibration system response. (a) Phase Diagram; (b) Displacement; (c) Velocity; (d) Acceleration.

Figure 5. Maxwell viscoelastic vibration system response. (a) Phase Diagram; (b) Displacement; (c) Velocity; (d) Acceleration.

Figure 6. A Kelvin viscoelastic vibration system with the converged parameters.

Figure 7. V11 well dynamometer measurement at 7:14:19 a.m. on 7 August 2019 Stroke 17 and 18.

Figure 8. Kelvin model analytical solutions with converged parameters in comparison with experimental measurements for Stroke 17.

Figure 9. Kelvin model analytical solutions with converged parameters in comparison with experimental measurements for Stroke 18.

Table 1. Newton–Raphson iteration convergence of the relaxation time

τ

, the scaling factor C, and the relative incremental error

ϵ

as defined in Equation (21) with the logarithmic modification.

Table 1. Newton–Raphson iteration convergence of the relaxation time

τ

, the scaling factor C, and the relative incremental error

ϵ

as defined in Equation (21) with the logarithmic modification.

No.	$τ$	C	$ϵ$
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

Table 2. Newton–Raphson iteration convergence of the relaxation time

τ

, the scaling factor D, and the relative incremental error

ϵ

as defined in Equation (21) with the logarithmic modification.

Table 2. Newton–Raphson iteration convergence of the relaxation time

τ

, the scaling factor D, and the relative incremental error

ϵ

as defined in Equation (21) with the logarithmic modification.

No.	$τ$	D	$ϵ$
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

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).

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	$ϵ$
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

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).

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.	$k_{o}$	m	$c_{1}$	$k_{1}$	C	$ϵ$
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

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.

© 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/).

Share and Cite

MDPI and ACS Style

Wang, S.; Gao, D.; Wester, A.; Beaver, K.; Edwards, S.; Taylor, C.A. Pump System Model Parameter Identification Based on Experimental and Simulation Data. Fluids 2024, 9, 136. https://doi.org/10.3390/fluids9060136

AMA Style

Wang S, Gao D, Wester A, Beaver K, Edwards S, Taylor CA. Pump System Model Parameter Identification Based on Experimental and Simulation Data. Fluids. 2024; 9(6):136. https://doi.org/10.3390/fluids9060136

Chicago/Turabian Style

Wang, Sheldon, Dalong Gao, Alexandria Wester, Kalyb Beaver, Shanae Edwards, and Carrie Anne Taylor. 2024. "Pump System Model Parameter Identification Based on Experimental and Simulation Data" Fluids 9, no. 6: 136. https://doi.org/10.3390/fluids9060136

Article Menu

Pump System Model Parameter Identification Based on Experimental and Simulation Data

Abstract

1. Introduction

2. Theory and Modeling

3. Inverse Optimization

4. Implementation and Improvement

5. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Acknowledgments

Conflicts of Interest

Nomenclature

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI