A New Methodology to Analyze Fluid Pound in Sucker-Rod Pump Systems: Phenomenological Approach

— Sucker-rod pump system is the most used artificial lift method in the world. Although used for more than a century at industrial scale, its working is not completely understood. The main tool employed to analyze operation is the dynamometric charts which is obtained with equipment already producing. A major problem found in these systems is the fluid pound that indicates a partial barrel filling with the producing fluid. This specific fluid pound problem is detectable by this tool but after the event occurs. In order to offer a tool at design stage, this study proposes an alternative method based on phenomenological flow modeling and rigorous geometry details. The method is based on Fluid Structure Interaction (FSI) technique in which both Computational Fluid Dynamics (CFD) and Rigid Body Analysis are solved simultaneously. The geometry used was comprised in the standing valve part of the subsurface pump. At the present stage, the model considers isothermal, incompressible single-phase turbulent flow. The aim of the methodology is to calculate the volume of fluid admitted into barrel at upstroke phase. It is proposed a fill parameter, which represents the fraction of filled barrel volume, to indicate the existence or not of fluid pound problem. The methodology allows to investigate several operational parameters. In the present paper, the stroke length, pressure intake and cycles per minute were investigate. The tendencies previewed by the methodology were in agreement with field observations, demonstrating the potential of this tool.


INTRODUCTION
Artificial lifting methods are extremely important to optimize and maintain well production, even with declining pressure over the years. Depending on the characteristics of the wells and the fluids produced some methods of artificial lift are more indicated than others [1].
One of the most widely used lifting methods in the world is the Sucker-Rod Pump (SRP) mainly due to its simplicity of operation and widespread technical knowledge.( [2], [3], [4]) In the onshore production of Northern Espírito Santo State (Brazil) this method has a great representativeness of use due to the fact that it has, in the majority, production of heavy oils besides being low depth wells. In addition, due to the need for cyclic steam injection, this is the most appropriate artificial elevation method [5]. However, when the temperature of the reservoirs decreases and the viscosity of the oil increases, this method presents some deficiencies that in some cases makes even the operational continuity impossible. The dynamic level (fluid level above the pump suction) should be as small as possible, being standardized at 50 meters, which is approximately equivalent to a pressure of 5 kgf/cm². In relation to fluids, the viscosity is quite high (especially at low temperatures) and BSW varies greatly throughout the productive life due to the injection of steam.
One of the major problem points in the SRP, when operating with viscous oils, is the standing valve. Normally, due to the difficulty of flow, the pump inlet diameter is maximized aiming at a reduction of the load losses favoring a complete filling of the pump jacket. This maximization most often causes oversizing of the pump and sub-surface components. When high viscosities and high pump speeds are combined, one of the phenomena detrimental to the system, the "fluid pound", is identified. This phenomenon occurs due to the complete non -filling of the pump barrel, often caused by this over-sizing of the pump. Thus, in the downward pumping cycle the travelling valve remains closed until the maximum level of fluid inside the barrel is reached. At this point, abrupt opening of the valve occurs and the ball can damage the travelling valve cage or damaging the ball itself [6].
Other phenomena occur in the standing valve, such as the deposition of sediments, making it difficult to close properly, causing leaks and inefficiency of the system. This problem can be mitigated by constructively altering the standing valve components, including the ball diameter and weight. Understanding the points of pressure drop and stagnation pressure also avoids another problem caused by the "gluing" of the ball in the seat, thus hindering the ingress of the fluid by the intake of the pump [7].
One of the main tools for analyzing the wells is the surface and bottom chart which provides important information about the behavior of the loads throughout the course of the pump. Through these charts it is possible to identify if the well is producing correctly or if there is an operational or failure problem. One of the major problems in this type of analysis is that it cannot be extended to other wells since the chart is obtained from field data ( [4], [8], [9], [10]). In this way, a methodology that predicts the behavior of the pump in function of the operational data would be of great utility.
A first reported attempt to introduce fluid flow phenomenological based model to analyze SRP was made by Li et al. [11]. They propose the use to this tool to design and research purposes.
Thus, in the face of the challenges of the complete understanding of the operation of this type of device, this work presents a methodology for the analysis of wells operating with SRP through the mathematical modeling of the incompressible and isothermal fluid flow produced inside the standing valve coupled to the model of movement of the ball as a rigid body. The coupling is accomplished by integrating the forces of the flow on the surface of the sphere. The three-dimensional geometry is detailed so that the standing valve components of the SRP subsurface pump are faithfully modeled. The set of governing equations are solved using the finite volume method using the ANSYS FLUENT v14.5 computational tool [12]. The opening and closing times and speeds of the ball were determined as a function of stroke, intake pressure and pump speed (cycles per minute or CPM) in order to create an operating map summarizing the stroke length/ CPM / intake pressure to prevent "fluid pound".

II.
METHODOLOGY The methodology proposed in this work consists of four stages: 1. transforming the CAD geometry of the standing valve into a calculation domain for the fluid flow; 2. discretization of the calculation domain taking into account the mesh around the sphere (deformable boundary); 3. resolution of the governing equations according to the operating parameters (cpm, stroke and intake pressure) in two stages (permanent and transient); 4. analysis of the volume of fluid admitted to the jacket (downstream of the cage) and the time taken to fill it (total or partial): fluid pound or not.

CO MPUTATIO NAL DO MAIN
The valves (standing and travelling) are of extreme importance in the mechanical pump, so much so that Takacs [2] considers the valves as the heart of the whole system since all the efficiency of the pump depends on the correct functioning of the same. Based on this argument for a better system efficiency, the design of the valves should be as optimized as possible, especially when it comes to pumping fluids with high viscosity where the load losses are higher.
As the problem is of an operational nature, and there are no evidences in the literature consulted of studies for him, it was decided to use a pump installed in the field, as object of study, in such a way that the results have utility in the industry. Fig. 1 shows all the valve's construction parts. The construction of the geometry for the analysis of the flow took into consideration several constructive aspects of the pump, obtained by means of measurements made in the own equipment using a caliper with precision of 0.05 mm. The valve dimensions are actual and measured in a Bolland pump, which corresponds to a 2¾"diameter pump [13]. This model was chosen due to its great use in the Northern Espirito Santo onshore wells, Brazil.
Due to the complexity of the standing valve and reduction of the computational mesh, 3 planes of symmetry were used so that the final geometry was represented by a slice of 120º, according to Fig. 2.    (left) and the constructed geometry (right). From the point of view of flow modeling, the volume that matters is the inside of the valve. Thus, the negative of the geometry is taken and the volume is obtained as shown in Fig. 4. Fig. 4A, 4B and 4C represent the constructed geometry, the constructed geometry with the part of barrel tube and the negative of the valve geometry, respectively. In addition, it is possible to observe that there is a space between the cage and the exit. This space aims to decrease the influence of the outlet on the flow in the valve.

a. CO MPUTATIO NAL MESH
Once the contours for the simulations are defined, one starts for the generation of meshes . The type of mesh used was the unstructured with tetrahedral and prismatic elements, built by the Octree method. This type of mesh has become advantageous considering the complexity of the geometry under study. Five meshes were generated for the test of independence of meshes, the number of elements is presented in Table 1.  . 5 shows the volumetric elementsthat crosses plane half height of the cage for the 5 meshes generated.

Fig. 5 -Cut Plane location (left); Comparison between the volumetric elements of the five meshes generated (right).
It is worth mentioning that the mesh deformation occurs only in the places near the sphere, that is, below the sphere there is a stretching of the mesh whereas in the upper part there is a compression of the elements of the mesh.

IV. MATHEMATICAL MODELING a. GO VERNING EQ UATIONS
The governing equations are the equation of continuity, (1), and conservation of linear momentum (Navier-Stokes equations), (2), describe the Newtonian fluid motion with specific mass ρ and viscosity μ constants : whereu is the velocity vector, p, the static pressure, t, time, g is the acceleration of gravity. Two comments are valid for the set of equations (1) and (2): first, the equations and therefore their variables are averages in time; second, as a consequence of the above statement, an effective viscosity arises from the average process, rather than the dynamic viscosity μ,and is defined as The above comments reflect the URANS (Unsteady Reynolds Averaged Navier-Stokes) approach used in this work, based on the concept of turbulent viscosity [14]. The definition of turbulent viscosity is the second term on the direct side of (3). In this definition, a new parameter (Cμ = 0.09) and two variables (k, turbulent kinetic energy and ε, k dissipation) are used to be obtained by the following two transport equations [15] where the constants of the model are: σk = 1,00; σε = 1.30; C1ε = 1.44; C2ε = 1.92. The term source Gk of equations (4) and (5) is defined as The motion of the ball (metallic sphere)is modeled through Newton's second law applied to rigid bodies. From this perspective we can write that whereae is the acceleration vector of translation, m, the mass of the sphere and fe are forces acting on the sphere. In the context of this work, restrictions were imposed on this model. 1. The velocity vector (and hence the acceleration) has only one non-zero component (axial); 2. There is no rotation of the sphere; 3. The forces acting on the sphere are gravity and drag.
Thus, we can write the above vector equation as a scalar equation wherefg is gravity and fd is the drag force. The drag force is calculated directly from the solution of the conservation equations by integrating the pressure and viscous forces on the sphere area, according to the following equation wherep is the static pressure and τw is the viscous stress on the wall obtained from the integration of the velocity field.
Once the acceleration is calculated, the translational velocity of the sphere is known to multiply by the time step. The displacement of the center of gravity of the body, in turn, is known by multiplying the velocity by the time step. This calculation is always performed at the end of the time step.

c. BO UNDARY CONDITIONS
The boundary conditions are shown inFig.6. The fluid inlet (INLET) is located to the left representing the fluid produced being admitted to the pump. The ball (SPHERE) is represented as a non-slipping wall condition, as well as WALL representing the outside of the valve. The condition SYMMETRY are planes of symmetry where normal flow is not allowed. Condition OUTLET is the outlet of the fluid towards the pump jacket. The numerical values adopted in the boundary conditions in this study are summarized in Table 2. As mentioned in the methodology, the simulations are divided into two stages: permanent and transient. The permanent stage has its initial condition of fluid at rest inside the pump and with the ball 1 mm from the seat fixed in space. The transient stage (where the sphere can travel) has as its initial condition the velocity, pressure and turbulence field resulting from the permanent stage.

e. PHYSICAL PROPERTIES
The fluid used was from a real oilfield well in the northernEspirito Santo, Brazil, and its main characteristics are in Table 3. The turbulent isothermal flow model considers the fluid to be incompressible and Newtonian.The ball has a mass of 400 g and is considered perfectly smooth. If filling parameter is greater than 1, the barrel is completely full and there is no fluid pound. If filling parameter is less than 1, the barrel is no completely full and fluid pound is predicted.

V.
NUMERICAL APPROACH The numerical approach used here is presented into two parts: first part is the discretization of the conservation equations (mass, momentum and turbulence model); second part concerns the mesh deformation due to ball movement.

a. GO VERNING EQ UATIONS DISCRETIZATIO N METHO D
The governing equations are discretizedemploying the Finite Volume Method [16]using the environment of the ANSYS Fluent v14.5 commercial package [12].
The velocity coupling algorithm used is SIMPLE (Semi IMPlicit Linked Equations) [17]. The discretization schemes of the advective terms were: Standard, for pressure, second order upwind for momentum, upwind first order for turbulence equations and Least Squares Cell-Based method for the calculation of gradients .
The convergence criteria for the simulations residues was 0.001 for all equations.The time step was 10 -5 s with a simulated real time of about 20 ms depending on the operating conditions on a computer with a processor: Intel Core i5-2.6GHz; Memory RAM: 6.0 GB. For this computers configuration, the dynamic mesh simulations had a duration of approximately 12 hours in each of the 24 cases so that the sphere came out of the seat and reached the top of the cage and after that another 12 hours of simulation were made with the ball at the top in order to obtain the rate of stabilization of the fluid.

DISCRETIZATIO N
The mesh deformation algorithm is split into two steps: the equations describing the ball movement (subsection IV.b) provides the ball new position; new ball positions is diffuse through the mesh and, if the case, a new mesh is generated (Fig. 7).

c. DYNAMIC MESH ALGORITHM
With the movement of the ball, the relative position of one of the boundaries is changed and thus, the mesh must also be altered. The solver can take two actions when there is movement of contours: 1. propagate the movement of the boundary in the near inner mesh through a diffusion equation; 2. remesh. Action (i)is always performed when there is boundary movement. The equation that governs the motion of each inner node is the diffusion equation ⋅ ( ⃗⃗⃗⃗ ) = 0, (11) whereVm is the velocity of the mesh and gamma is the diffusion coefficient.
The diffusion coefficient formulation employs the concept that the further the inner node of the moving contour, the smaller the displacement of this node. Thus, gamma can be expressed as whered is the distance from the inner node to the boundary, a is a diffusivity tuning parameter. If a = 0, it means that the movement of the contour will propagate evenly throughout the mesh. On the other hand, a very large parameter (something like 10) will dampen the movement of the inner nodes close to the moving contour.
The diffusion equation is solved with the finite element method. Using this method rather than the finite volumes (used in the conservation equations) is more convenient since it directly provides the value of the displacement of the nodes. The displacement is calculated by the following equation:

International Journal of Advanced Engineering Research and Science (IJAERS)
[Vol -6, Issue-6, June-2019] https: //dx.doi.org/10.22161/ijaers.6.6.86  ISSN: 2349-6495(P) | 2456-1908(O) However, the algorithm can lead to distorted elements of poor quality. In this situation, a new generation of mesh must proceed. The solver performs the procedure automatically when low quality elements are detected (whose quality value is defined by the user). The procedure is as follows 1. the solver identifies elements that do not meet the following criteria a. skewness greater than the tolerance, b. element size greater than the specified maximu m, c. element size less than the specified minimum, 2. the solver agglomerates the elements identified in 1; 3. the solver performs its internal mesh generation algorithm only on the agglomerated region; 4. If the solver does not reach a higher quality of the new mesh, the previous mesh continues in the calculations . Due to the movement described above, the conservation equations must fit this movement. Thus, whether a volume of control v any, and a conserved scalar or vector , the transport equation must be rewritten as : where ∂v represents the faces of the control volume, v. Therefore, there is a modification in the advective term to take into account the effect of the movement of the borders on the faces of the control volumes. Another point to note is that the term transient involves the volume of the control volume and this volume may be variable from onetime step to the next. This is the second effect that must be taken into account.
Using the first order discretization for the transient term, one can write that: note that v n+1 may be different from v n . Thus, an extra calculation for volume evolution is required. This can be done by Taylor's expansion on the volume around the volume in the previous iteration, ie: In order to satisfy the conservation law, the time derivative of the volume must be evaluated as wherenf is the number of faces of the control volume and Ajis the area of face j.
In turn, the scalar product within the summation should be evaluated as the volume swept by face j in the Δt interval. , ⃗⃗⃗⃗⃗⃗⃗ ⋅ ⃗⃗⃗ = . (18)

d. REMESHING AND SMOOTHING ALGORITHMS
The smoothing method is usually used in triangular or tetrahedral meshes, but when the limit displacement is large (compared to the size of the local elements), the quality of the mesh may deteriorate. This will invalidate the mesh (such as resulting in negative cell volumes) and consequently, will lead to convergence problems when the solution is upgraded to the next time step. In order to bypass this problem and prevent the elements from violating the convergence criteria, the remeshinggenerates a new mesh locally. Thus, if the new elements satisfy the asymmetry criteria, the mesh is updated with new local elements. Otherwise, the new elements are discarded and the old mesh is retained.
For these reasons, the methods chosen for this problem were remeshing and smoothing, as they were the ones that best adapted to this particular case [12].The procedure is as follows: 1. The boundary of the cell zone moves predominantly in one direction (ie without excessive anisotropic stretching or compression of the cell zone) -Condition for using the smoothing method; 2. Movement is predominantly normal to the boundary zone -Condition for using the smoothing method; 3. Triangular or tetrahedral mesh -Condition for using the remeshing method.

VI. RESULTS AND DISCUSSION a. CASES STUDIED
The SRP used in this study is located at the Northern Espirito Santo field, Brazil. Due to the operating characteristics of the SRP selected, the available pumping stroke lengths are 86 in (218.44 cm), 106 in (269.24 cm), 124 in (314.96 cm) and 144 in (365.76 cm). Thus, for each case, a stroke lengthwas chosen and the CPM was changed to values of 5, 7 and 10. In addition to these data, the inlet pressures were changed to 1 and 5 kgf/cm². Thus totaling 24 analyzed cases.
With the increase of the stroke length, the volume that must be filled increases, so a greater amount of fluid is needed. This change also influences the CPM because the distance to be covered by the plunger increases. However, the change in CPM influences the cycle time, the up-stroke time and also the piston speed. Thus , the filling time of the barrel is also changed.
The total pump flow is influenced by the stroke length and CPM combination, and actual pump flow by efficiency.

b. VELO CITY FIELD PATTERN
The contour plot of a typical velocity field at the start and at the end of ball movement is shown in Fig. 8.

International Journal of Advanced Engineering Research and Science (IJAERS)
[  At the beginning of the flow there is flux through the side passages of the cage and also through the outlet at the top where the valve puller is connected. It is the moment when the lowest speeds are observed because the valve is almost completely seated in the seat. This is due to incompressibility hypothesis used in the model. In case of gas is released from solution, adjustments to the model should be made. This can be seen in Fig. 8A which shows the fluid velocity contour inside the pump.As the ball approaches the cage (Fig. 8B), high velocities are observed at the valve seat which is the smallest cross -sectional area. The velocity field forms a jet-like pattern. Depending on Reynolds number, the jet could oscillate which leads to a ball movement not aligned with the axis. In its turn, the oscillating ball movement leads to an excessive weariness.
It is noted from Fig. 9 that just above the seat a flow stagnation zone is formed. It is observed in the field that sand deposition and "oily sludge" are formed at this location. This deposition causes the phenomenon known as "gluing" of the ball. Due to this fact, it is common to operate the bottom checkingvalve, i.e. it is necessary for the plunger to touch the standing valve, causing a stirring in the ball and taking it off.This maneuver, although necessary, is extremely damaging to the pump as it damages the puller and can also damage the cage.
Another location that also occurs stagnation is at the top of the cage at the puller connecting point. Obstruction of this point makes it difficult to connect the puller to the standing valve when it is removed with the Terrestrial Production Probe. Fig. 9shows the velocity vector of the fluid within the cage at the moment the ball reaches the top of the cage. The graphic of Fig.11 shows the average fluid velocity as a function of the flow time at the pump inlet. As it is possible to see, the speed increases rapidly and stabilizes at around 737 cm/s. This value was considered for the calculation of the volumetric flow rate.

c. PRESSURE FIELD PATTERN
The pressure field obtained is shown in Fig.12. At the moment the ball touches the top of the cage (Fig.12B) it is possible to observe that a low-pressure zone is formed in the region of the puller connection. Due to this effect, a zone of stagnation of flow is formed, as it is observed in the images of the velocity profiles.A high-pressure zone is also observed below the valve. This pressure is responsible for supporting the ball at that point during the rise of the piston.  While in the graph of Fig.14, the velocity of the ballis observed as a function of its displacement time to the top of the cage. At time t = 0, the velocity of the ball is also zero, but in the following time step (0.00001 s) the velocity of the ball is approximately 143 cm/s after the velocity decreases and increases again over time.This happens because at the moment the simulation starts the ball is practically in the seat, generating a zone of high pressure just below the sphere, as it can be seen in Fig.12, after that first moment the ball moves upwards doing with which the pressure just below the ball decreases and consequently the speed of the ball also decreases. Then, the velocity of the sphere increases again until it reaches the moment it stops at the top of the cage, just before that moment the simulation with the dynamic mesh is finished and the simulation is started to stabilize the fluid flow.This time with the ball stopped at the top of the cage, as has been said previously.

e. FLUID PO UND ANALYSIS
It is determined, for the cases studied, whether or not the phenomenon of fluid pound occurs. In the following, will be presented and discussed the effects of each operational variable (CPM, stroke length and intake pressure) on fluid pound.

Effect of cpm on fluid pound
This is a paramount parameter since it is easily manipulated during production. Its effect is illustrated in Fig.15. In these results, the inlet pressure was maintained at 5 kgf/cm² and the stroke length was maintained at 124 inches. It is possible to state qualitatively that the higher the CPM the worse it is for the operation. The region above the black line indicates no fluid pound, while the region below is likely to fluid pound occur. For the operating points studied, only one (the lowest CPM) is recommended to avoid fluid bumping.

Effect of stroke length on fluid pound
The stroke of the plunger is another parameter that can be manipulated in the field. Its effect is shown in Fig.16. In these results, the intake pressure was kept constant at 5 kgf/cm² and the CPM at 7. It is possible to perceive the same trend of the previous parameter: the larger the stroke length, the greater the possibility of a fluid pound. This behavior is consistent with the intuition that the larger the barrel's volume, the more time it takes to fill it. This is a parameter that can be controlled indirectly through the fluid level in the annulus of the well. Its effect on the filling parameter is illustrated in Fig.17. The results presented have fixed values of 7 CPM and 86-inch stroke length. Unlike the other parameters, increasing the intake pressure has a positive effect on the fill. Low pressures indicate low speeds, no time to fill the barrel.