Scaling, Complexity, and Design Aspects in Computational Fluid Dynamics

: With the availability of more and more efﬁcient and sophisticated Computational Fluid Dynamics (CFD) tools, engineering designs are also becoming more and more software driven. Yet, the insights in temporal and spatial scaling issues are still with us and very often imbedded in complexity and many design aspects. In this paper, with a revisit to a so-called leakage issue in sucker rod pumps prevalent in petroleum industries, the author would like to demonstrate the need to use perturbation approaches to circumvent the multi-scale challenges in CFD with extreme spatial aspect ratios and temporal scales. In this study, the gap size between the outer surface of the plunger and the inner surface of the barrel is measured with a mill (one thousandth of an inch) whereas the plunger axial length is measured with inches or even feet. The temporal scales, namely relaxation times, are estimated with both expansions in Bessel functions for the annulus ﬂow region and expansions in Fourier series when such a narrow circular ﬂow region is approximated with a rectangular one. These engineering insights derived from the perturbation approaches have been conﬁrmed with the use of full-ﬂedged CFD analyses with sophisticated computational tools as well as experimental measurements. With these conﬁrmations, new perturbation studies on the sucker rod leakage issue with eccentricities have been presented. The volume ﬂow rate or rather leakage due to the pressure difference is calculated as a quadratic function with respect to the eccentricity, which matches with the early prediction and publication with comprehensive CFD studies. In short, a healthy combination of ever more powerful modeling tools along with the physics, mathematics, and engineering insights with dimensionless numbers and classical perturbation approaches may provide a balanced and more ﬂexible and efﬁcient strategy in complex engineering designs with the consideration of parametric and phase spaces. The results derived from the full-ﬂedged three-dimensional computational ﬂuid dynamics (CFD) simulation based on the same transient (or rather quasi-static) pressure difference and plunger velocity which by itself is very challenging due to extreme spatial dimensions does match very well with the analytical solutions based on steady Poiseuille and Couette ﬂows for both concentric and eccentric annulus regions and their rectangular domain approximates. In short, a combination of analytical and computational approaches is presented with comparison to existing empirical measurements. This is in particular evident for the calculation of relaxation times with both cylindrical and Cartesian coordinates which is consistent with the petroleum industry practice in the ﬁeld. In conclusion, a healthy combination of ever more powerful modeling tools along with the physics, mathematics, and engineering insights with dimensionless numbers and classical perturbation approaches may provide a balanced and more ﬂexible and efﬁcient strategy in complex engineering designs with the consideration of parametric and phase spaces.


Introduction
In oil industry, artificial lift methods are often utilized to transfer the oil and gas from the formation to the surface, especially during the later phase of a petroleum reservoir's productive life due to the diminishing fluid pressure [1]. Sucker rod pumping systems are among the most popular and low cost artificial lift means with their simple hydraulic pumping designs [2,3]. As one of the earliest inventions for oil industries, sucker rod pumping systems have been studied extensively [4,5]. Over 80% of the wells world wide still depend on this oldest and most extensively used mechanism [6,7]. Yet, despite of their extended history and popularity, many engineering and operational issues still exist [8][9][10]. Despite the effectiveness and efficiency of this seemingly simple pumping mechanism amply demonstrated over decades [11,12], in practice, it has been reported that the leakage is related to the clearance C, or rather two times the gap δ, defined as the difference between the inner diameter of the barrel D b and the outer diameter of the plunger D a [13,14]. It is widely reported that large clearances often cause significant leakage which often renders the entire sucker rod pump systems inefficient.
In this paper, for the purpose of illustrating the need for perturbation studies before full-fledged computational fluid dynamics (CFD) simulations, we focus on the so-called leakage issue related to the necessary gap between the plunger and the barrel for the smooth operation of the pump mechanisms [15,16]. As illustrated in Figure 1, a typical sucker rod pump system includes a moving contact area between the traveling unit (plunger) and a fixed tube (barrel). Further more, a traveling ball valve is attached with the end of the plunger and a standing ball valve is installed with the tube near the oil and gas reservoir where the oil and gas are pushed through the fissure formation [17,18] by the so-called formation pressure. To facilitate the reciprocal motions of sucker rod pump systems, a clearance between 3 to 8 mills (one thousandth of an inch) is often reserved to avoid the direct contact (abrasion) between the plunger and the tube [19,20] and to facilitate the transport of sands and particles in oil and gas mixtures [21,22].
Although engineering designs are becoming more and more software driven, physical insights with both temporal and spatial scales are still important in particular during the preliminary design stages. We would like to demonstrate, through the analytical studies of the so-called leakage issues in sucker rod pumps prevalent in petroleum industries, the need for in-depth understanding of contributions of the geometries with extreme spatial ratios and the flow characteristics with specific relaxation time. These insights will help to understand the complexity and many design aspects and enable more efficient use of computational resources with challenging spatial aspect ratios and temporal scales.
In this study, the gap size between the outer surface of the plunger and the inner surface of the barrel is measured with a mill (one thousandth of an inch) whereas the plunger axial length L p is measured with inches or even feet. In the computational approaches, the extreme aspect ratio between the plunger length L p and the plunger outer radius R a (D a = 2R a ), or the barrel inner radius R b (D b = 2R b ), as well as the gap size δ, defined as R b − R a , or a half of the clearance C, does present significant challenges. Furthermore, in engineering practice, the plunger and the barrel are often not concentric and have various levels of eccentricity. The temporal scales, namely relaxation times, are estimated with both expansions in Bessel functions for the annulus flow region and expansions in Fourier series when such a narrow circular flow region is approximated with a rectangular one. These engineering insights derived from the perturbation approaches have been confirmed with the use of full-fledged CFD analyses with sophisticated computational tools. With these initial confirmations published in an early paper [16], new perturbation studies on the sucker rod leakage issue with eccentricities have been presented in this paper. The volume flow rate or rather leakage due to the pressure difference is calculated as a quadratic function with respect to the eccentricity, which matches with the early prediction and publication with comprehensive CFD studies [16]. In short, a healthy combination of ever more powerful modeling tools along with the physics, mathematics, and engineering insights with dimensionless numbers and classical perturbation approaches may provide a balanced and more flexible and efficient strategy in complex engineering designs with the consideration of parametric and phase spaces.

Analytical Approaches
As illustrated in Figure 1, during the operation of the sucker rod pump systems, within the gap between the plunger outer surface and the barrel inner surface, two fundamental fluid flows exist, namely, a pressure difference driven Poiseuille flow and a boundary motion induced Couette flow. Both Poiseuille and Couette flows are idealized quasi-static laminar flows and could be governed by the following equations where the plunger length is L p , the fluid density is ρ, and the pressure gradient ∂p ∂z is expressed as − ∆p L p .
where the plunger length is L p , the fluid density is ρ, and the pressure gradient ∂p ∂z is expressed 109 as − ∆p L p . Let's assume the plunger velocity as U p , namely, u(R a ) = U p . Furthermore, to be consistent with the leakage definition which is from the top to the bottom, the pressure difference ∆p is positive when the upper region pressure p h is higher than the lower region pressure p l , moreover, the plunger velocity U p is positive when it is moving downwards. Consider a plunger and a barrel pump system with an inner radius R a which is a half of the outer diameter (OD) of the plunger and an outer radius R b which is the half of the inner diameter (ID) of the barrel. The gap denoted as R b − R a = δ or the clearance C, which is equivalent to 2δ is small in comparison with both R a and L p . According to the specific geometries of typical sucker rod pump systems and their operating conditions, we could first ignore the inertial and time dependent effects and denote the pressure difference p h − p l as ∆p, the Navier-Stokes equation in the cylindrical coordinate system can be simplified as where p h and p l represent the pressure on the top of the plunger and on the bottom of the plunger, or rather within the sucker rod pump, namely, between the traveling valve and the standing valve, µ refers to the dynamic viscosity of the fluid. From Equation (2), we derive where C 1 and C 2 can be decided based on the boundary conditions. In order for us to understand the selection of these two constants C 1 and C 2 , let's continue with this steady linear partial differential Equation (3). For Newtonian viscous fluid, with the linear superposition principle, we can solve the Couette flow and the Poiseuille flow separately. For the Poiseuille flow, on the inner surface of the pump barrel and the outer surface of the plunger, we have the kinematic conditions u(R a ) = 0 and u(R b ) = 0. The velocity profile within the annulus region expressed as Equation (3) has For small clearance, with the Taylor's expansion, we derive hence, the solution of Equation (3) can be expressed as Moreover, we can easily establish the flow rate through the annulus region as Within O(δ 3 ), the flow rate due to the pressure difference, namely Poiseuille flow, Q p is established as with the perturbation ratio = δ R a .
Consequently, the viscous shear force acting on the plunger outer surface in the direction from the top to the bottom can be calculated as Likewise, for the Couette flow, with the moving outer surface of the plunger, together with the sucker rod, since the barrel is stationary, thus, for Newtonian viscous fluid, the fluid velocity at the barrel inner surface is zero. Hence, we have the kinematic boundary conditions u(R a ) = U p and u(R b ) = 0, and the flow field can be expressed as with Using the Taylor's expansion, we have the simplified expression for the flow field, Notice that the gradient of the velocity profile at the plunger surface matches with the approximation with respect to the thin gap. Moreover, the flow rate due to the shear flow, namely, Couette flow, Q c can still be written as yet the flow direction is the same as the plunger velocity U p , namely, from the bottom to the top when the upper region pressure p h is higher than the lower region pressure p l .
Again, using the Taylor's expansions in Equation (5), we have As a consequence, the viscous shear force acting on the plunger surface in the direction from the top to the bottom can be calculated as It is clear that the leading term of the shear force F c is consistent with the Newtonian fluid assumption with the dynamic viscosity µ. Finally, by combining the Couette and Poiseuille flows, we can have the steady state solution for the velocity profile φ(r) expressed as where the Taylor's expansion in Equation (5) is employed for the coefficients in Equations (4) and (10) with δ = C/2. It is the cylindrical coordinate system that renders this seemingly simple problem complicated. If however, we utilize the scaling based on the physics and mathematics, for the large aspect ratio between the plunger length L p and the gap size of the annulus region δ as well as between the plunder radius R a and the gap size, we can cut open the annulus region and simplify the flow domain as a rectangular box as shown in Figure 2 with an axial length L p (z direction), a width 2πR a (x direction), and a height h (y direction) [23,24]. Notice here that even with the eccentricity which is marked with the difference 2e between the widest gap and the narrowest gap, or rather e, the distance between the center of the outer surface of the plunger and the center of the inner surface of the barrel, there exists a mid symmetrical axis at x = πR a and the flow regions with x ∈ [0, πR a ] and [πR a , 2πR a ] are identical. Once we recognize the symmetry, we only need to consider one half of the annulus region with eccentricity and the half of the perimeter is denoted with x ∈ [−πR a /2, +πR a /2], as depicted in Figure 2. Of course, for the concentric sucker rod pump, we have a uniform gap with h = δ. However, with eccentricities, such a height will be a function of x which will be discussed separately in Section 3. For narrow annulus regions, the governing Equation (1) for the Poiseuille and Couette flows can be simplified as where w is the velocity component in the axial or z direction and the pressure gradient in z direction is still constant. Again, for the Poiseuille flow, on the inner surface of the pump barrel at y = h and the outer surface of the plunger at y = 0, we have the kinematic conditions w(0) = 0 and w(h) = 0. Hence, the velocity profile within the annulus or rather simplified rectangular region can be expressed as Furthermore, we can easily establish the flow rate Q p through the concentric annulus region with h = δ as δ 0 2πR a w(y)dy. The flow rate due to the pressure difference Q p is established as It is not difficult to confirm that the leading term in Equation (7) matches with the simplified expression in (17). Consequently, the viscous shear force acting on the plunger outer surface in the direction from the top to the bottom can be calculated as It is again confirmed that the leading term in Equation (8) matches with the simplified expression in (18). Note that the viscous shear force acting on the plunger outer surface in the direction from the top to the bottom. Moreover, due to the small gap size, it is reasonable to assume the shear force acting on the outer surface of the plunger is the same as that of the inner surface of the barrel [23,24]. Thus, these two surface shear forces will balance the total normal force due to the pressure difference over the plunger length, namely, where F p stands for the viscous shear force acting on the plunger outer surface due to Poiseuille flow. It is clear that Equation (19) is consistent with Equation (18) and the leading term in Equation (8). In fact, in engineering practice, the dominant term is often sufficient. It is obvious that with the help of the physics and mathematics insights [25,26], the simplified rectangular domain is much easier to handle than the annulus region and this advantage will be more crucial when we discuss the relaxation time and the case with eccentricities in Section 3.
Similarly, for the Couette flow, on the inner surface of the pump barrel at y = h and the outer surface of the plunger at y = 0, we have the kinematic conditions w(0) = U p and w(h) = 0. Hence, the velocity profile within the annulus or rather simplified rectangular region can be expressed as Furthermore, we can easily establish the flow rate Q c through the concentric annulus region with h = δ as The flow rate due to the shear motion at y = 0 (outer surface of the plunger) is established as which matches with the leading term in Equation (12). Consequently, the viscous shear force acting on the plunger outer surface in the direction from the top to the bottom can be calculated as where F c is the viscous shear force acting on the plunger outer surface due to Couette flow. In comparison with Equation (13), it is again confirmed that the leading term matches with the simplified expression in (22). Furthermore, in order for us to derive Equation (23) from a full-fledged Navier-Stokes equations, we must identify whether or not the fluid flow is in the turbulent region as well as the transient effects [27,28]. First of all, within the gap which is measured in mills, for typical oils, the kinematic viscosity ν at 100 o C is around 5.3 cSt or 5.3 × 10 −6 m 2 /s, about five times that of the water and the plunger velocity U p is no more than 40 in/sec, thus the so-called Reynolds number Re = U p δ/ν is much smaller than 100 let alone the turbulent flow threshold around 2000. Although the Reynolds number is a clear indication about the quasi-static nature of the Couette and Poiseuille flows within the narrow annulus region, in order to have some guidance with respect to the selection of the sampling time in the experimental measurements of the pressure and the displacement within the sucker rod pump unit, we must investigate further the inertia effects and other time dependent issues. Consider the overall governing equation for the viscous flow within the annulus region R a ≤ r ≤ R b as expressed as where the plunger length is L p and the pressure gradient ∂p ∂z is expressed as − ∆p L p .
Note that the pressure difference ∆p is positive when the upper region (top) pressure is higher than the lower region (bottom) pressure which is consistent with the leakage definition. Assuming the plunger velocity is U p , namely, w(R a ) = U p , by combining the Couette and Poiseuille flows, we can have the steady state solution for the velocity profile φ(r) expressed in Equation (14).
The final unsteady velocity profile w(r, t) can be written as where the governing equation for the transient partw(r, t) is expressed as with the boundary conditionsw(R a , t) = 0 andw(R b , t) = 0. Using the separation of variable method and common special functions [29], we introducew(r, t) = ψ(t)ϕ(r). As a consequence, we have ϕ(R a ) = ϕ(R b ) = 0 and the following governing equationsψ where ν is the kinematic viscosity and the time scale τ is also called the relaxation time.
For an exponentially decreasing function expressed as e −t/τ , the tangent line at the origin always provides an horizontal intercept τ and in general within 5 and 6 times the relaxation time τ, the function is considered as sufficiently close to the steady solutions. From Equations (26) and (27), we have ψ(t) = A o e −t/τ along with the characteristic expression It is clear that the characteristic solution expressed in Equation (28) is based on the combination of the first and the second kind of Bessel function satisfying Equation (27) and the characteristic time τ defined as the relaxation time in Equation (26) can be determined by the boundary conditions of ϕ(r), namely ϕ(R a ) = ϕ(R b ) = 0. Therefore, in order for Equations (26) and (27) to have nontrivial solutions, we must have nontrivial or nonzero solutions of A 1 and A 2 , for the following linear system of equations, It is clear that we must have the determinant of the coefficient matrix A in Equation (29) equals to zero, namely Finally, combining the steady and the transient solutions, the complete velocity profile can be expressed as A k e −t/τ k g k (r) (31) where the special function can be expressed as and the coefficient A k is calculated as In order to separate the effects of the pressure difference and the plunger motion and getting ready for quasi-static or unsteady integrations, Equation (14) is rewritten as follows.
where we have As a consequence, the velocity profile can be expressed as where the coefficient E k is calculated as It is known that for viscous fluid to squeeze through a narrow annulus channel, the required pressure gradient or pressure drop given a fixed channel length will be proportional to the gap δ to the power of three, namely, O(δ 3 ). The derivation based on Poiseuille flow does confirm such a notion. Moreover, the leakage issue is mostly prominent during the upward stroke, namely, the sucking process. In the upward stroke, the plunger is moving up and as a result will drag fluid up instead of down, which is opposite to the leakage. This effect is due to so-called Couette flow which will introduce a flow rate in a different direction of the power of one, namely, O(δ). Note that in the Couette flow, we focus only on the shear flow introduced by the inner wall (plunger) translational motion whereas no rotational motions as discussed earlier literatures will be considered.
In this paper, through the solution of the partial differential Equation (25) in cylindrical coordinate systems, it is possible for us to derive the relaxation time scales and confirm the quasi-static nature in the majority of oil production systems and the related measuring technologies in particular the sampling frequency around 30 Hz. We will also delineate the transient or rather quasi-static time dependent leakage due to Poiseuille flow from that of Couette flow.
As the plunger is lifting up, the pressure difference ∆p, namely, p h -p l or the pressure gradient ∂p ∂z , namely, − p h − p l L p with L p denotes the plunge axial length will drive the fluid from the high pressure area, namely, the top of the plunger, to the low pressure area, namely, inside the pump chamber and below the plunger. This leakage is due to the so-called Poiseuille flow within the annulus region. It is clear that this leakage is due to the pressure difference. The other flow is due to the viscous shear of the plunger motion. It is important to recognize that during the upstroke, this part of flow rate is in the same direction of the plunger, in this case, from the bottom to the top, which is opposite to the leakage direction. This flow is due to the so-called Couette flow within the annulus region. Consider a plunger and a barrel pump system with an inner radius R a which is a half of the outer diameter (OD) of the plunger and an outer radius R b which is the half of the inner diameter (ID) of the barrel, in fact, the width of the rectangular domain is simply the perimeter based on the outer radius of the plunger for this first order approximation as illustrated in Figure 2. In practice, the plunger length is measured in inches or even feet whereas the gap size is measured in mills. Thus, the gap denoted as R b − R a = δ or the clearance C which is equivalent to 2δ is very small in comparison with R a and L p . As a consequence, Equation (25) in cylindrical coordinate systems can be expressed in Cartesian coordinate systems as illustrated in Figure 2 ∂w with the boundary conditionsw(0, t) = 0 andw(h, t) = 0. Notice that to prepare the discussion in Section 3 with eccentricities, we use the gap function h(x) where the coordinate x is selected to represent the entire perimeter of the outer surface of the plunder or the inner surface of the barrel. In the concentric case, we of course have the constant h which is defined as the gap size δ.
Using the separation of variable method and common special functions [29], we introducew(y, t) = ψ(t)ϕ(y). As a consequence, we have ϕ(0) = ϕ(h) = 0 and the following governing equationsψ From Equation (40) and the corresponding boundary conditions, we can easily obtain the sinusoidal expansion of the function ϕ(y), namely, sin( nπy δ ), with the modal sequence number n. Thus, we have the evaluation of the relaxation time τ expressed as Notice that the Couette effects cannot be approximated by the special function in Equation (32) because of the imposed boundary conditions [30]. Consider a practical case with D = 1.5 in and C = 5 mill, the measured plunger velocity U p can reach around 40 in/sec, namely, R a = 0.75 in and δ = C/2 = 0.0025 in. For the fluid considered in this paper, we set the density ρ as 1.0014 × 10 3 kg/m 3 and the plunger length L p as 48 in. Assume the dynamics viscosity as 0.9867 cp, we have the Reynolds number based on the clearance ρCU p /µ as 65.5 which is much smaller than 100. Consequently, we have the following graphical representation in Figure 3 of the determinant expressed in Equation (30). The eigensolutions α k with k = 1, 2, · · · , relate to the non-dimensional characteristic time or relaxation time τ k = δ 2 π 2 ν 1 α 2 k . Using the parameters and material properties for the same test case, τ 1 = 4.146 × 10 −1 ms, τ 2 = 1.037 × 10 −1 ms, and τ 3 = 4.607 × 10 −2 ms. In fact, the non-dimensional relaxation times ξ n evaluated by the roots in Figure 3 can be simply expressed as 1/n 2 . This observation is indeed confirmed by the simplification of the Bessel functions by the sinusoidal functions as R a /δ → ∞. It is clear that the typical sampling period of 8.333 ms is sufficiently larger than the five or six times the largest intrinsic relaxation time as predicted in Figure 3. Therefore, the transient or time dependent effects can be ignored. This result is consistent with the finding using full fledged Bessel functions and also confirms the rationale for the non-dimensional quantity δ 2 n 2 π 2 ν chosen for Figure 3.

Eccentricity Issues
Finally, we will consider the eccentricity effects on leakage. We hope this comprehensive study will help the industry to better understand the nature of the leakage with respect to pressure difference, eccentricity, and motion related to the plunger of typical sucker rod pump systems [31,32]. In the near future, experimental setups will be established to validate the widely used empirical Patterson Equation and provide more evidences and confirmations to our theoretical studies.
As illustrated in Figure 1, the plunger is in general closely fit within the pump barrel. In general, there will be a clearance of 3 mills for a typical oil well. There are two considerations of this type of simple fitting. First of all, the friction force is due to viscous shear instead of Coulomb friction due to the direct contact or abrasion. This friction force will be demonstrated to be fairly insignificant in comparison with the other forces. Finally, for sandy oil fields, it is also necessary to reserve some space to avoid sticking and abrasion.
The eccentricity of the plunger's lateral position in the barrel has a great effect on liquid leakage. In practice, for completely eccentric position, leakage rates 60% greater than for concentric cases can be expected. For a thorough study of eccentricity issues, we introduce the eccentricity as the ratio between the distance e of two centers of the circular outer surface of the plunger and the circular inner surface of the barrel with the gap size δ, i.e., 2e/C. As illustrated in Figure 2, we can cut the perimeters following the line between the two centers utilizing the symmetry. Therefore, we within half of the perimeter represented with x ∈ [−πR a /2, +πR a /2], instead of [0, πR a ], we can simply interpolate the gap function h(x) as where the minimum and maximum gap sizes are δ − e and δ + e, respectively. Of course, we always have the eccentricity between 0 and 1, namely, 0 ≤ e ≤ δ. Due to symmetry, we only need to consider half of the width in this case, namely, x ∈ [−πR a /2, +πR a /2], therefore, the velocity component in the axial (z) direction can be expressed as, following Equation (16) Furthermore, we can establish the flow rate Q p through the annulus region with eccentricity e as 2 The flow rate due to the pressure different Q p is then established as with = δ R a and the eccentricity is evaluated as e δ .
Due to the complexity of the problem, we also validate our analytical approaches with a full-fledged CFD simulation with various eccentricities.

Computational Approaches
In the computational simulation with ADINA, as illustrated in Figure 4, we applied the pressure differential for three out of four strokes, namely 1295 sample points or 10.792 s, leaving only the stroke during which both the traveling valve (TV) and the standing valve (SV) are closed. Notice that in this test case, we have the stroke per minute (SPM) as 5.067, thus the period of the sucker rod pumping unit is about 11.842 s. To be consistent with the rate of 120 sampling points per second, namely, the sampling period 8.333 ms, the same time step size is also adopted. Furthermore, to keep the consistent units, in ADINA, the pressure unit is psi, the dynamic viscosity unit is psi · s as well. However, since the spatial dimension is represented with inch, therefore, the conversion factors 32.147 for gravitational constant and 12 for foot and inch conversion must be adopted in the left hand side (LHS) of the Navier-Stokes equation, which means, the density must be adjusted accordingly in addition to the standard conversion from kg/m 3 and lbm/ft 3 . Therefore, the actual input to ADINA if the dimensions are in inches and the pressure differentials are in psi, thus we have 9.358 × 10 −5 for effective density and 14.505 × 10 −8 for effective dynamic viscosity. Notice the challenges the extremely small gap size (a few mills, or a few thousandths of an inch) in comparison with the plunger length (48 in) have posted to the meshing of the annulus flow region. In fact, in Figure 4, in order to visualize the flow region, we have to zoom in further, as a result only a tiny segment of the annulus region could be depicted. Clearly, for the Poiseuille flow the velocity is zero at the outer surface of the plunger and the inner surface of the barrel as depicted in Figure 4.
For the concentric case, we utilized the measured plunger velocity and pressure differential data provided by Echometer Company in Wichita Falls, Texas. As depicted in Figures 5 and 6, the simulation results with the actual plunger velocity and pressure differential measurements as well as the inertial effects are nearly identical to those based on the quasi-static analytical results in Equations (7) and (12) which are evaluated without the inertial effects. This fact along with the relaxation time scales confirms that for the typical operation conditions illustrated in this paper and in engineering practice, the inertial effects within the annulus region is insignificant and can be ignored [6,28]. Again, according to the theoretical studies on the relaxation time scales which are confirmed with both the approximated rectangular flow region approximated with Fourier series and the cylindrical flow region calculated with Bessel functions, we can conclude that the time variations of the flow rate as represented in Figs 5 and 6 are in fact quasi-static in nature. The time scales calculated in this paper also reflect the reality that in petroleum industries, the sampling period of 8.333 ms is sufficiently large in order to avoid the transient effects of flows within annulus regions of the sucker rod pump systems. Essentially, we could consider these time variations are snapshots of steady solutions since the relaxation time in the transition solution is very small in comparison with the sampling time.
Finally, to incorporate the eccentricity effects, we have resorted to the full-fledged computational fluid dynamics simulation with ADINA [16]. As shown in Fig. 5 of the recent ASME Journal of Fluids Engineering paper [16], the computational results with eccentricities were compared with the corresponding concentric case. The new perturbations results depicted in Equation (44) have been compared with the published results in Ref. [16] and documented in Figure 7. In this paper and the early ASME paper [16], the eccentricity is measured by 2e/C, where e is the shift of the center of the inner cylinder. Different from Ref. [33], instead of studying the leakage rate vs. pump depth and upward moving process, in this paper, we focus on the degree of eccentricity measured by 2e/C, or e/δ. For the case with D = 1.5 in and C = 5 mill, and the plunger length L p = 4 ft, as illustrated in Figure 7, which has also been confirmed by the recently published computational results in Ref. [16], the misalignment can sometimes introduce up to 90% elevation of the leakage which is consistent with the analytical prediction based on the perturbation approaches in Equation (44).

Conclusions
In petroleum industry, because we must limit the contact as well as the wear and tear of the sucker rod pump plunger and the pump barrel, it is not desirable to completely eliminate the gap between the outer surface of the plunger and the inner surface of the barrel thus the leakage due to the pressure difference which could be as high as a few thousand psi. Furthermore, by reserving a slight gap between the pump plunger and the pump barrel, the sucker rod pump will also be able to operate in oil field with the pollution of sands. It is thus desirable for engineers to have a better understanding along with quantifiable relationship with respect to the leakage and many relevant operational conditions such as the pressure difference and eccentricity. Leakage issue like many engineering problems is multifaceted and complex. For complicated problems such as this, simple analytical, or computational, or empirical studies alone are not sufficient. Yet, we also would like to have more design insights through deep analytical studies before full-fledged simulation tools are employed.
Taking into consideration of the excessively small gap or clearance (a few mills) and its related radius of curvature, pump barrel or plunger radius, we compared and confirmed the matching of the transient solutions with Bessel functions for cylindrical coordinates and Fourier series for Cartesian coordinates. We also compared and confirmed the matching of the relaxation times derived with the first and the second kind Bessel functions for cylindrical coordinates and Fourier series for Cartesian coordinates. The fact that the multiple of the largest relaxation time with five or six is still smaller than the sampling time step which accounts for thirty samples a second for actual cases with physical dimensions and experimental measurements confirms that the quasi-static assumption widely employed in the field is sufficiently accurate and appropriate.
The results derived from the full-fledged three-dimensional computational fluid dynamics (CFD) simulation based on the same transient (or rather quasi-static) pressure difference and plunger velocity which by itself is very challenging due to extreme spatial dimensions does match very well with the analytical solutions based on steady Poiseuille and Couette flows for both concentric and eccentric annulus regions and their rectangular domain approximates. In short, a combination of analytical and computational approaches is presented with comparison to existing empirical measurements. This is in particular evident for the calculation of relaxation times with both cylindrical and Cartesian coordinates which is consistent with the petroleum industry practice in the field. In conclusion, a healthy combination of ever more powerful modeling tools along with the physics, mathematics, and engineering insights with dimensionless numbers and classical perturbation approaches may provide a balanced and more flexible and efficient strategy in complex engineering designs with the consideration of parametric and phase spaces.