A Computational Fluid Dynamics Study for Simulation of the Unsteady State Two-Phase Flow in Vertical Separators and Plunger Pumps Simultaneously

Developing and optimizing separators are a challenging task that demands advanced laboratory facilities or costly field assessments. Computational Fluid Dynamics (CFD) simulations can be a prompt and economical alternative. This study aims to find a robust and reliable method to simulate unsteady state two-phase flow in the sophisticated vertical separators connected to plunger pumps using CFD. The unsteady nature of reciprocating pumps and vertical two-phase gas–liquid flow forces the unsteady state boundary condition to the separators. Consequently, to simulate them, the plunger pump was considered a part of geometry, and an unsteady state turbulent two-phase flow was selected for the simulation method. Due to the geometry vastness and complexity of such a problem, the mesh generation method became vital for the convergence and reliability of simulations. Among several mesh generation methods and software, a new MATLAB meshing method developed to generate structural hexahedral mesh proved robust and time-efficient. First, a set of experimental data of a simple separator was used to select, tune and validate the simulation method. Then, a real sample sophisticated separator and a plunger pump were analyzed using the method. Results demonstrated the ability of the model to predict phase separation and motions and its capability for modeling two-phase flow inside the helical, gravitational, and hybrid separators together with a plunger pump. Performing simulation for three up and two down strokes of the plunger showed that after disappearing the initial condition effects, the model represents the same results and separation efficiency for all cycles of plunger movements.


Introduction
Separators are widely used in many industries.They aim to separate and prevent the undesired phases from entering the mainstream.One field that depends on separators is the sucker rod pump (SRP) system, the oldest and most commonly used artificial lifting system with several advantages, such as high efficiency, simplicity, and minimal maintenance [1].The main limitation of this method is gas and sand control.To solve this problem, separators must be installed before the pump inlet.These separators are divided into of separator performance and some rules of thumb that are still practicable.These principles were applied in Campbell and Meridian's studies [6], which developed practical methods for tailoring separator performance to specific well production conditions.After that, Robles and Podio [7] began experimental studies in the laboratory, which included visual observation of the separator fluid mechanics using a full-size Plexiglas well and a conventional rod pump.
Patterson and Leonard [9] also investigated different downhole gas separator designs in Wyoming.Their objective was to study the operation of separators in the dewatering of low-pressure natural gas wells and coalbed methane wells.The inlet to the gas separators was smaller in these designs and, together with some baffles, allowed gas to escape from the gas separator.They could acquire good gas separation in field installation.However, the field installations demonstrated the operation of gas separators; estimating the influence of each design parameter was extremely difficult, leading to Lisigurski's [10] laboratory study of the gas separator geometry.He built an apparatus simulating a production well, with significant modifications that could handle different volumetric flow rates of water and air to evaluate the performance of several separator designs (eleven downhole gas separator models were tested).Eventually, he developed an experimental methodology.
McCoy et al. [2] also performed laboratory testing on five different gas separator designs and presented field data on two.They examined the best position for installing a downhole separator in relation to the casing perforations and the effect of separator inlet port geometry and dip tube diameter on efficiency.Another experimental analysis was performed by Bohorquez et al. [3].They studied different models of gravity-driven separators, which differed in terms of inlet port structure, dip tube size, and the location of the separator inlet ports relative to the well perforations, using a full-scale model of the well and separators in the laboratory.Based on the results of the experiments, they proposed a new separator design that used centrifugal force to separate the two-phase air-water flow.McCoy et al. [11] developed a packer-type separator which was studied by field tests and achieved good pump filling with liquid.They also described a quantitative technique for evaluating the effectiveness of downhole gas separators based on comparing the pump's liquid filling in relation to the percentage of liquid present in the casing annulus surrounding the pump.Najafi et al. [4] studied the efficiency of an oil-air helical downhole separator experimentally.They designed three different inlet structures and found that the separator's inlet geometry significantly impacts its efficiency.
Although downhole gas separators have been used and studied for many years, the industry still needs to develop guidelines or standards to evaluate the performance of downhole gas separators [11].This is because it is difficult to assess the effect of a single parameter in the field tests without being influenced by other factors.Although laboratory testing solves this problem, experiments are only partially accurate and correspond to actual field situations due to limitations.The first limit is safety issues that force operators to use fluids such as water and air instead of hydrocarbon gas in solution with oil and differences between the separation conditions of these two groups of fluids.Another reason is ignoring the influence of the well pressure in the laboratory simulation.Since it is impossible to generate the high-pressure state of the well in the laboratory, previous studies considered the atmospheric pressure [4] or a range of 5-10 Psi [2,3].In addition, the pressure changes during downhole pump movement were not considered in experimental investigations.Finally, the modeling equipment and mechanism in the laboratory are expensive, especially the full-scale ones, which are not always affordable and timeconsuming.All the reasons prove that it is advantageous to use computational methods instead of experimental ones to simulate downhole separators.While there are few numerical studies of downhole separators [12][13][14], these are mostly case studies focusing on estimating the efficiency at different gas-liquid ratios or geometry improvements.According to the separator installation position (Fig. 1a), the separator outlet is connected to a downhole plunger pump.Movements of the pump's plunger cause pressure variations at the separators' outlet, leading to phase separation inside these devices.This aspect is identical to the real working conditions of separators, while it has not been studied before and is presented in this work.Furthermore, simulation of the upstroke and downstroke of the plunger is equally important since the two-phase flow behaves differently during these two movements, which is also considered in this paper.
The geometry of this study has a long length, several narrow holes, and high curvature helical part.For grid generation, all these factors should be considered at the same time.Previously studied geometries with vanes, spirals, or helixes used a hybrid mesh for their computations [15][16][17][18].Using hybrid mesh rises the cell number significantly (for example, 4 million cells for a 1226.5 cm 3 domain with a spiral section of 0.16 of the total pipe length [19]).To avoid the increase of cells number and make the simulation possible for a domain with approximately 285,000 cm 3 volume, a new approach for generating a completely structural grid is developed using MATLAB coding.
Finally, there had been some simplifier assumptions in previous numerical works.These simplifications were considered because of the complexities related to the turbulent multiphase flow simulation.These include performing a 2dimensional simulation [20] or using a steady-state solver [17,21] for modeling downhole separators.With the aid of the self-made grid, a method for unsteady 3-dimensional simulation of downhole gas-liquid separators in the real scale is presented considering the well casing annular section and a downhole plunger pump.Also, real properties of hydrocarbons and given well conditions are implemented.The numerical method is discussed and implemented to simulate one simple separator and one sophisticated separator connected to a plunger pump.

Theory
An unsteady turbulent multiphase flow model should be set to simulate the performance of a downhole gas separator and a plunger pump.To model the plunger motion and perform realistic simulations, the transient solver was used for all the simulations.Given that the fluid domain was large enough and the heat transfer was rapid due to the turbulence of the flow, with a suitable approximation, the temperature variation could be assumed to be negligible, so solving the energy equation was unnecessary.Regarding the turbulence model, the realizable k-ε model [22] was used as it showed better convergence over other 2-equation turbulence models due to its more accurate formulation for turbulent viscosity.After analyzing the Eulerian-Eulerian multiphase flow approaches, including Volume of Fluid (VOF), Mixture, and Eulerian models, the VOF model was chosen.This model can accurately predict separation, tracks the interface between fluids, and tracks the movement of large bubbles in liquids that were not acquired using other multiphase models [23].

Governing Equations
In the VOF model, each phase is represented by its volume fraction, and tracking the interface between two phases is accomplished by solving a continuity equation for the volume fraction of the secondary phase.
where α L is the volume fraction of the liquid phase, ρ L is the liquid density, and − → v is the mixture velocity.The volume fraction of the primary phase, gas here, is calculated as follows.
where α G is the volume fraction of the gas phase.The momentum equation is the next governing equation and is solved by considering that the phases share the velocity field.The momentum Eq. ( 3) is dependent on the volume fractions of all phases by the properties ρ (density) and μ (viscosity), as shown in Eqs. ( 4) and (5).
In Eq. ( 3), p is the pressure, − → g is the gravity, and − → F is the surface tension force between two phases.In this study, surface tension effects were considered negligible compared to inertial forces, arising due to velocity (W e ρv 2 l σ >> 1).Regarding turbulence governing equations, k-ε realizable equations were applied.All the parameters are kept in their default mode.Details of equations are not presented here to avoid repetition and are available in [24].

Separator Efficiency
The performance of the separator can be evaluated by its efficiency.The efficiency can be expressed as follow: where α G, Inlet is the gas volume fraction at the casing inlet and α G, Pump barrel is the gas volume fraction at the separator outlet (in the downhole pump barrel), which are presented in Eqs.(7) and (8).(7) α G, Pump barrel the volume of gas inside the pump barrel the total volume of the pump barrel (8) In Eq. ( 7), Q, each phase's flow rate is given as the inlet boundary condition.Also, there is another method for evaluating separators' performance.This method describes separator outlet liquid fraction (OLF) as the fraction of the liquid flow rate leaving the separator to the total flow rate [2].

OLF
Q L, separtor−outlet Q L, separtor−outlet + Q G, separtor−outlet (9) This fraction can also be calculated using the volume of liquid filling the pump barrel to the total volume.

O L F
the volume of liquid insid the pump barrel the total volume of the pump barrel (10) 3 Software and Tools In this study, MATLAB R2020b was used for writing grid generation codes.Ansys Meshing, Fluent Meshing, Fluent, and Results component systems of Ansys 2020 R2 were used for mesh generation, analysis, and postprocessing, respectively.In addition, pressure, volume, and temperature (PVT) calculations for defining fluid properties were performed by HYSYS Process 2.2 All simulations were performed using 20 parallel processors and considering double precision.The transient simulations were carried out using a 3.70 GHz Intel® Core™ i9-10900 K CPU with a 128 GB RAM system.

Grid Generation Methods
Grid generation for unsteady turbulent two-phase simulation of a plunger pump and a separator inside the casing simultaneously was a challenging problem.The total volume of the domain was vast (15,000 cm length, largest diameter of 18 cm, and approximately 285,000 cm 3 volume), while it had 276 narrow diameter holes (7 mm) and high curvature parts that had to be considered at the same time.In addition, in the pump section, where the geometry had a moving part, there was another challenge due to the limitation of the Fluent software in simulating the dynamic layering mesh.This type of dynamic mesh simulation is only available with a structural grid.
Therefore, a smaller and less complex geometry were initially considered for grid generation.This geometry (Fig. 3) consisted of a hollow cylinder with a helix attached to its inner tubing.Different mesh generation software was studied to generate a mesh with the highest quality and accuracy while minimizing the solution time.
The selected method was used for the grid generation of this study's simple and sophisticated separators.

Fluent Meshing
Fluent Meshing is a new, intelligent mesh generation software that can generate meshes without partition creation.This software can recognize geometry areas with smaller dimensions and create a more refined mesh without requiring user definition.The weakness of this software is that it cannot generate a structure grid in the entire geometry [25].

Ansys Meshing
Ansys Meshing is another software used to generate the mesh.This software can create structured and unstructured meshes.It is necessary to subdivide the geometry into cubic blocks, specify the number of subdivisions or the element size for each edge of each block and then specify the mesh type for each volume to create a structural grid in Ansys Meshing.This process can be challenging and time-consuming in large geometries with small parts.

Structural Mesh with MATLAB
Performing a simulation with a structural grid can reduce the solution time and required memory [26].There are many software programs for generating 3D meshes, some of which can generate structural meshes.The programs examined in this study were Flac3D, Gmsh, Salome, Autodesk, CFX, CF-MESH, Abaqus, FLUIDYN, and MATLAB.After investigating how each of them works, it was found that some Fig. 4 Steps of the structural mesh generation using MATLAB, a Grid generation for a circle and b a cylinder, c Section view of the final mesh, d 3-D view of the helix and inner tubing have limitations, like generating 3D meshes only for 3D surfaces or requiring precise geometry slicing, which, like Ansys Meshing, was time-consuming.It was found that it is possible to create a structural and three-dimensional mesh using MATLAB coding.
There were several obstacles in the process of mesh generation using MATLAB.The first issue was that MATLAB could not read geometry files.For this reason, the geometry should be simplified and generated step by step by coding in MATLAB.
To generate mesh for the geometry of Fig. 3, an attempt was made to create a structural square mesh within a circle (Fig. 4a).Then, by adding the third dimension, cubic blocks were designed to shape a cylinder (Fig. 4 (b)).In the next step, the Cartesian coordinates were converted to Cylindrical.Thus, the inner walls and helix positions were identified and removed from the initial cylinder (Fig. 4c and d).Finally, the structural grid for the fluid domain of the geometry of Fig. 3b was accomplished.
The generated mesh had to be imported into Fluent software in the next step.Because Fluent cannot read the format of the MATLAB code, a function was added to the code to convert the MATLAB format into one of the readable formats of Fluent (Tecplot).Then, the fully structured grid was transferred to Fluent software for analysis.
At that point, the imported mesh had an entire exterior wall and an interior zone.Therefore, inlet and outlet surfaces had to be specified.To this end, the walls were separated in the Fluent, and the mesh was prepared to implement the solution.
Additionally, the developed MATLAB meshing method allows the generation of different-size grids within a geometry.The MATLAB code was defined based on the cylindrical coordinate system.In this study, the number of divisions in the tangential coordinate was considered the same for the whole geometry.However, approaching from the perimeter to the center of the cylinder, the tangential size of the grids became finer.Regarding the radial size of the grids, the number of radial divisions increased in the central part and sections with complex geometry.Also, due to the long length of this study's geometries, the grids' size in the height direction was considered the same.
To evaluate the independency of results from the grid and choose the appropriate cell density, four grids with 294,660, 196,440, 98,220, and 69,540 structural hexahedral cells were generated using MATLAB coding for the geometry of Fig. 3.
The same simulations (details in Table 2) were performed for these grids.The results of the tangential velocity, axial velocity, oil volume fraction, total pressure, and run-time are shown in Fig. 5.As shown in Fig. 5, increasing the number of cells from 98,220 to 294,660 did not cause a notable effect on the simulation results; however, it decreased the computational efficiency.Therefore, the MATLAB-generated grid with 98,220 cells was selected as the most accurate and efficient one to be compared with other mesh types.

Selecting the Mesh Generation Method
Simulations were carried out with different mesh types to choose the best grid.A transient turbulent two-phase simulation was performed under the same conditions with different grids created with MATLAB, Ansys Meshing, and Fluent Meshing for the geometry of Fig. 3.Note that in Fig. 3, the total helix length is the same in both (a) and (b) geometries.Details of mesh types are presented in Table 1, and these grids can be seen in Fig. 6.Also, polyhedral and polyhexacore grids were generated with Ansys Meshing, which did not show accurate results and thus is not reported here.In these simulations, the flow was chosen to enter the geometry    2. These grids were compared using two criteria: orthogonal quality and aspect ratio (Table 3).The run-time with different grids was also evaluated to achieve the convergence accuracy of 10e − 4. In addition, the velocity contours and liquid phase contours (Fig. 7) were checked to determine the accuracy of the grids in detecting separation.
As can be seen from Table 3, both quality criteria of MAT mesh type are significantly higher compared to other grids, with this grid having the lowest number of cells.The highest orthogonal quality of 0.271 was achieved by generating multiple meshes with Ansys Meshing, and the mesh generated with Fluent Meshing had the maximum orthogonal quality of 0.21.In contrast, this quality was 0.997 and almost equal to one in the structured mesh generated with MATLAB, representing the best quality.The same pattern can be seen for the aspect ratio index.Therefore, the MATLAB mesh represented the highest quality compared to other meshes.
The solution time was another critical factor in the grid selection since all solutions in this study were unsteady, and the simulation's applicability depended on the duration of each solution.Because the mesh generated with MATLAB was structural, it had a significantly shorter solution time than other meshes.Therefore, this factor also made MATLAB Mesh the best option.
It can be seen from the contours in Fig. 7 that the accuracy of AM1 and AM2 mesh types was not optimal since the circular shape of the helix was not captured due to the small number of meshes (approximately 100,000).The number of cells was increased to solve this problem.Contour AM3 in Fig. 7 shows a 4 million cell grid with the helix curvature problem solved.This grid had the same smallest and largest mesh size as the MAT mesh, and the simulation results with it were similar to MAT mesh, while the number of cells was 40 times, and the solution time took 77 times longer.
Although the FLM mesh type did not have the curvature problem in the helical section and represented it well, its results were less accurate than MAT and AM3 mesh types.Furthermore, according to Fig. 7, MAT and AM3 mesh types accurately represented the separation and fluid velocity within this geometry.Due to the similarity of the results of these two meshes, MAT mesh was chosen for the simulation with less computational effort, solution time, and higher quality.

Numerical Modeling
Details of the numerical model are represented in this section.
The pressure-based transient solver was used to perform the simulation.The gravitational force was considered because it separates the heavier phase from the lighter one.In the VOF model, the explicit formulation for the volume fraction was chosen because it does not require an iterative solution of the transport equation during each time step, as is required for the implicit scheme.In addition, the implicit body force option was enabled to make the solution more robust.When body forces like gravity are present in a multiphase flow, segregated algorithms converge poorly unless a partial balance of pressure gradient and body forces is considered.The implicit body force treatment can account for this effect and allows the flow to reach a realistic pressure field very early in the iterative process.Furthermore, the k-ε realizable model was used as the turbulence model.The standard wall functions were selected among different wall functions available with the k-ε realizable model.The standard wall functions reasonably predict most high-Reynolds numbers in wall-bounded flows.The non-equilibrium wall functions are not selected since they are applicable in flows involving separation from the wall, reattachment, and impingement, like airfoil flows.Since, in the simulations of this study, the flow conditions are not departing too much from the ideal conditions (see [25] for the definition of unusual conditions), there is no need to use enhanced wall functions.The standard wall functions in ANSYS Fluent are provided as a default option and have been most widely used in industrial flows.
In all simulations, the mass flow condition was used for the inlet boundary, and the pressure condition was used for the outlets.In addition, intensity and hydraulic diameter for the turbulence specification were defined in all boundaries.Boundary conditions and material properties are discussed in more detail in each simulation setting.
The Pressure-Implicit with Splitting of Operators (PISO) scheme was used for pressure-velocity coupling.After solving the pressure correction equation, the PISO algorithm satisfies the momentum balance by performing two additional corrections, neighbor correction, and skewness correction.Also, this scheme requires slightly more CPU time per solver iteration; it can drastically reduce the number of iterations needed for convergence, especially for transient problems, compared to other segregated methods.
The least-square cell-based scheme for the gradient, PRESTO! for the pressure and second order upwind for the momentum, the turbulent kinetic energy, and the turbulent dissipation rate were chosen for discretization.The geometric reconstruction scheme was also used to discretize volume fractions.This scheme represents the interface between fluids using a piecewise linear approach and is the most accurate scheme in ANSYS Fluent.In addition, a first-order implicit method was used for the transient formulation.

Modeling of a Simple Separator for Validation
A simple gravity separator was simulated to evaluate the accuracy of the presented methodology.was modeled.They examined separators in different installation positions and for various air and water flow rates.Since installing the separators so that their entry ports are above the casing perforations is the most common position in field installations [2], this experiment was chosen and simulated.
The field performance of this separator was also assessed in McCoy et al. [2] study, and the field data showed good agreement with the laboratory results.

Geometry
A separator geometry based on the Patterson-type separator in the study by McCoy et al. [2] was designed for the simulation, shown in Fig. 8a.As shown in Fig. 8b, in the simulation, the separator was considered in a 6-inch casing with a height of 54 ft, just like the lab tests.

Boundary Conditions
As seen in Fig. 8, the bottom surface of the casing was determined as the inlet and its top as the outlet.The inlet boundary condition at the casing inlet was set to mass flow inlet, and this amount for each phase was entered into the Fluent software.As given in the experimental data, the pressure outlet boundary condition was determined for the casing outlet and separator outlet.The amounts of boundary conditions are presented in Table 4.

Fluid Properties
Air and water were used as the two-phase flow in the laboratory testing.Therefore, these materials were used to simulate the simple separator, and their properties are represented in Table 5.

Grid Generation
The new MATLAB coding method generated a hexahedron structural mesh for the entire domain.The main dimensions of the geometry were changed negligibly (0.5%), making it possible to develop hexahedral meshes.A grid with 2,646,232 cells for the simple separator inside the casing was generated, representing the separator performance in good agreement with the experimental data.This grid is shown in Fig. 9.

Solver Setting and Modeling Summary
The unsteady simulation of the simple separator was performed using the proposed model and settings presented in Sect. 5.The model components for the simulation of the simple separator are demonstrated in Table 6.

Results and Validation
In the experiments of McCoy et al. [2], all results were measured after the flow had stabilized means the liquid height in the casing became almost constant.As seen in Fig. 10, this point was considered in the simulation, and the solution continued until the liquid level inside the casing stabilized.

The comparison between experiment 4 of McCoy et al. [2]
and the simulation results is shown in Table 7.
The most critical parameter showing the presented numerical model's accuracy is the correct separator efficiency prediction.In the study by McCoy et al. [2], the separator efficiency was defined as the volume fraction of the liquid at the separator outlet (Eq.( 9)).The results showed good agreement with the experimental values with a 2.4% deviation which is reasonable and acceptable.The robustness of the structural MATLAB meshing method was also proved.
In addition to the results presented in Table 7, the simulation was continued for an additional 4.5 s to ensure that the model could show the simultaneous upward movement of the separated gas and the downward movement of the liquid phase after the flow had stabilized.The upward path of the separated gas phase within the separator and its exit from the top holes in these 4.5 s are shown in Fig. 11 by a circle.This figure does not show the casing annulus to see the separator better.Also, the volume fraction contours at the entrance to the separator's dip tube are compared to the photo presented in the McCoy et al. [2] study in Fig. 12.It can be seen that in the simulation, as in the experiment, there is a significant amount of gas at the inlet of the separator's inner tube.Based on the results in Figs.11 and 12, it can be concluded that using the steady-state model for simulating the separation is inaccurate.While the flow becomes stable a few seconds after the start of the solution, the contour of the phases changes at any moment and can only be captured using the transient solution.
Fig. 10 The rise of fluid level inside the casing and its stabilization after 35 s of solution Fig. 11 The upward path of the separated gas phase within the separator and its exit from the top holes in 4.5 s after the flow stabilization

Modeling of a Real Sophisticated Separator and a Plunger Pump
This work is the first step of a major project to design and optimize more sophisticated separators connected to downhole plunger pumps.Due to the variety of downhole separators and scarcity of field or experimental data, we decided to validate our method with a simple separator described above and then use it for simulating the designed separators.This section describes the validated method's application for simulating a real sample sophisticated separator/plunger pump and the simulation results.

Geometry
The geometry of the sophisticated separator is shown in Fig. 13.The two-phase gas-liquid flow enters the sophisticated separator via its portals (Fig. 13a).Then, as can be seen from Fig. 13b, the flow passes through the helical section where the two-phase fluid is thrown toward the outer wall of the separator due to the centrifugal force; thereby, phase separation starts.After that, the fluids move downwards, and the separation process continues due to gravity and the difference in phases' density.The fluids are then divided more by passing through the baffles, where they become more turbulent.The separated gas sticks to the inner wall of the separator and, due to its lower density, moves upwards and leaves the separator through upper portals.In contrast, the liquid phase Since the separator portals are for the entry of two-phase fluid and exit of separated gas, the casing must also be simulated with the separator.In this study, a 7-inch casing with a length of 33 ft (10 m) was modeled.The casing was modeled for five more meters after the separator's portals to observe the liquid level inside the casing annulus during the up and down courses of the plunger.In the previous laboratory studies [2,3,11], separators were also assessed inside the casing.
In addition, the sub-plunger part of the pump was modeled to capture its movement effect on the separator performance.The geometry of the simulated part of the downhole pump is shown in Fig. 14.The entire fluid domain of this study, including the separator, casing, and sub-plunger part of the pump, along with inlet and outlet boundaries and the flow path, can be seen in Fig. 15.

Boundary Conditions
As seen in Fig. 15, the bottom surface of the casing was determined as the inlet and its top as the outlet.The inlet boundary condition at the casing inlet was set to mass flow inlet, and this amount for each phase was entered into the Fluent software.The pressure outlet boundary condition was determined for the casing outlet.In addition, the surface named the moving plunger in Fig. 15 was defined as the wall, and the dynamic layering mesh with constant speed was determined for it.The length of the plunger course is considered 7.5 m, and its speed is set to 3 strokes per minute (SPM).These parameters agree with the field operational conditions of the long-stroke SRP units.The amounts of boundary conditions are presented in Table 8.

Fluid Properties
In the simulation of the sophisticated separator performance, a given real well operating conditions and fluid properties   9 demonstrates some characteristics of a typical well, A. As seen from Table 9, fluid properties were available at bottom hole conditions, while for simulations, they were required at separator installation depth.Therefore, HYSYS Process software was used to calculate PVT and find parameters at the desired depth.The calculated viscosities and densities of the gas and liquid phases at a depth of 1400 m are presented in Table 10.

Grid Generation
Different software was used to generate mesh for the sophisticated separator.As generating mesh using Fluent Meshing was the most straightforward method, several tetrahedron grids with different dimensions and qualities were produced for the sophisticated separator without considering the pump, as shown in Table 11.
After running the solution with the introduced system and grids FLM1-5 (Table 11), it was observed that the run was very time-consuming as each 10-s simulation took about three months.Increasing the mesh dimensions increased the running speed, but the solution diverged after a few seconds of simulation.More grids were made, none of which helped solve the problem.
After that, attempts were made to generate a structural grid using Ansys Meshing.Therefore, the sophisticated separator was redesigned, and its portals and inner tubing holes were changed from cylindrical to cubic with the same crosssection, allowing partitioning into cubic blocks.Note that this initial model did not include the helix and pump sections.The number of subdivisions required for the geometry to achieve cubic blocks equaled 380 subdivisions, as shown in Fig. 16.Implementing size and mesh type for each block were timeconsuming, and there was no way of coding or algorithm for automatic grid generation; consequently, using this software did not seem viable.Finally, the grid was generated using MATLAB coding.The process of grid generation was discussed in Sects.3 and 4. The steps of the mesh generation and grids in different parts of the sophisticated separator are shown in Figs. 17 and 18, respectively.The unsteady turbulent two-phase simulation of the complete domain of the sophisticated separator was a timeconsuming task as it should be repeated several times.Therefore, to study the grid independency for this geometry, the helical section was selected as it is the most complicated part compared to the whole domain.Also, a great part of phase separation occurs in the helix.The generated grids for this part were discussed in detail in Sect. 4.
To generate cells for simulation, 6 degrees divisions in the tangential direction, variable radial divisions with a maximum of 2 mm, and constant 8 mm divisions in height were performed.Coarser grids were generated in the casing, while finer ones were produced inside the separator.Due to the long length of the geometry, constant divisions were performed in its height direction.
Therefore, 1,907,820 structural cells were created for the sophisticated separator inside the casing (Fig. 15) and increased to 3,383,640 cells after activating the dynamic mesh for the plunger part in Fluent.In reshaping for structural grid generation of the sophisticated separator, the total crosssections and lengths remained identical to the geometry of Fig. 13.

Solver Setting and Modeling Summary
The unsteady simulation of the sophisticated separator was performed using the proposed model.Separate modeling   of the plunger pump and separator leads to dismissing the pump's effect on the separator's performance and efficiency.Also, the separator outlet and the pump intake share the same variable boundary condition, which can only be captured by simultaneously modeling the pump and separator.This makes the simulations more realistic compared to using a constant flow rate or pressure for the separator outlet/pump inlet boundary conditions.The simulation was performed for 50 s, including three upward and two downward plunger courses.The plunger's upward movement was modeled using a dynamic mesh with the layering method, and its downward motion was simulated by deleting meshes by time.
In the downhole pumps, during the downstroke, fluid discharges from the region between two valves (Fig. 14).Since determining an outlet boundary to move is not available in Ansys Fluent, it was not possible to simulate the downstroke automatically.Dividing the downstroke into 15 steps (0.5 m), deleting the meshes of the pump barrel in each step, and continuing the run was a new approach to solving this issue.
The model components for the simulation of the sophisticated separator are demonstrated in Table 12.

Results
This section outlines the numerical results of the simultaneous simulation of the sophisticated separator and the plunger pump.As discussed, the simulation was run for three up and two down courses.Since the set initial condition influenced the results of the first up and down strokes, these were not taken into account, and the simulation continued for the second upward course.
Additionally, to gain confidence in the accuracy of the second upstroke results, the second downstroke and third upstroke were simulated.It was observed that the results of the second and third upstrokes were identical.Since each stroke lasts 10 s, the results are presented for 20-50 s of the simulation.The animation of the phase separation during 50 s of the simulation is available in the online version of this article.The 50 s simulation took 705.6 h, including 163.2 h for each upstroke and 108 h for every downstroke solution.

Velocity and Pressure Results
Figure 19 shows the velocity contour at the end of the upstroke in all parts of the studied domain.The velocity contours in all time steps during the upstroke were mostly the same.However, the velocity contour was entirely different on the downward stroke, as shown in Fig. 20.This is because there was no suction at the separator outlet during the downward stroke (Fig. 21), and the flow rate through the separator therefore dropped.
The pressure change during the upstrokes (20-30 and 40-50 s) and a downstroke (30-40 s) are shown in Fig. 22.As can be seen from this diagram, the total pressure at the casing outlet remained constant due to the consistent pressure boundary condition.However, there was a minor change at the casing inlet due to the liquid level changes and a significant one at the pump inlet (separator outlet).
The sudden drop in pressure at the start and during the upstroke is undesirable as it causes more gas to be separated from the liquid at the pump inlet, where there is no way to prevent gas from entering the pump, leading to inefficiency in the pumping system.On the other hand, the pressure recovered as the plunger moved down and remained largely stable during the downstroke.This point demonstrates the importance of considering the plunger pump in the simulation of separators simultaneously.Modeling the plunger pump leads to consideration of the variable boundary condition at the separator outlet.This makes the simulations more realistic compared to using a constant flow rate or pressure for the separator outlet boundary condition.

Phase Separation Results
The results of gas-liquid separation in different sections of the sophisticated separator are shown in Figs.23 and 24.Note that the casing is not shown in these figures to avoid confusion.As seen in these Figs, the gas-liquid flow entered the separator during the upstroke, and the helix and baffles forced it to break.However, some gas bubbles separated and escaped from the separator; the gas layer formed near the inner tubing wall of the separator and moved downward.When the plunger motion reversed, the liquid level rose inside the casing (Fig. 25).Less fluid would likely to flow to the separator, allowing bigger bubbles to form.These bubbles moved upward toward the separator portals to escape to the casing, which shows the importance of modeling the downstroke.At the end of the downstroke, there was mostly no gas near the inner tubing holes of the separator.The flow streamlines at the end of the upstroke are also shown in Fig. 24c, which shows the gas phase's simultaneous upward and downward motions.
When the plunger moved upward, the liquid was suctioned to the pump inlet to fill the barrel of the downhole pump.As a result, the liquid level came down inside the casing, and a reverse flow of gas formed at the casing outlet.This amount of gas exited the casing during the downstroke and did not change the volume fraction of the two-phase flow.

Separator Efficiency
The comparison between gas and liquid flow rates through different sections of the studied domain can be seen in the diagrams of Fig. 26.The solid black lines in Fig. 26a and b represent the gas and liquid flow rates exiting the separator and entering the pump barrel, respectively.Also, the inlet flow rates of liquid and gas are blue lines of Fig. 26.Therefore, the sophisticated separator efficiency in this specific operating condition can be calculated as follow: This method of calculating the efficiency of the downhole separator, considering the total amount of gas entering the pump, can be more accurate than using the average gas flow rate to leave the separator.Because as can be seen from the solid line in Fig. 26a, the gas flow entering the pump has some fluctuations, and the average value would not represent it accurately.Therefore, the total volume of gas exited from Also, comparing the diagrams of Figs.22 and 26 reveals the important impact of the plunger motion and, therefore, pressure variations on the separator efficiency.In lower pressure magnitudes, more separation at the pump intake happened.The gas flow rate entering the pump barrel increased, so the separator's efficiency decreased.This phenomenon is identical to the separator's real operating conditions.Considering constant outlet boundary conditions for the separator ignores these variations and would not acquire reliable results.
As discussed in Sect.

Conclusion
This study demonstrated a method for numerical modeling of the unsteady state turbulent two-phase flow in complex geometries, which was also implemented for modeling one simple and one sophisticated vertical separator and a plunger pump.The simple separator simulation was used to validate the method and showed good agreement with experimental data.A real sample sophisticated separator and a plunger pump simultaneously were modeled to examine this model's ability to simulate more complex and large geometries.Considering the plunger movements together with separators led to considering variable boundary conditions for separator simulation, identical to the real operating conditions.The later simulation, which is part of a major project, showed reasonable results, and the applied method promised to be an appropriate tool for designing and optimizing such equipment.This method can also be used to simulate equipment with large and complex geometries, such as gas lifts, jet pumps, ejectors, multiphase flowmeters, and cyclones, under near-to-real conditions.
Furthermore, a new MATLAB programming method for generating structural grids was introduced and evaluated.By having a value of almost one for the orthogonal quality and

Fig. 3
Fig. 3 Helical simulated geometries a Geometry used for AM1-3 and FLM mesh type simulations.b Geometry for MAT mesh simulation

Fig. 5
Fig. 5 Comparison of the results of MATLAB grids in a cross-section of the helix 800 mm away from the inlet

Fig. 6
Fig. 6 Generated mesh types, shown in the middle section of the geometry

Fig. 7
Fig. 7 Results of grid study after 9 s of simulation and flow stabilization, a Contour of oil volume fraction b Velocity contour (m/s)

Fig. 8 aTable 4 Table 5 × 10 − 2 Fig. 9
Fig. 8 a Designed separator based on Patterson type separator in McCoy et al. [2] study b Separator inside the casing and boundaries

Fig. 12
Fig. 12 Results of volume fraction contour at the entrance of the separator's dip tube, a photo reported in McCoy et al. [2], b simulation results in the same section

Fig. 13
Fig. 13 Sophisticated separator a Outside view b Transparent view

Fig. 14
Fig. 14 Geometry of the simulated part of the plunger pump

Fig. 15 a 8 Fig. 16
Fig. 15 a The entire fluid domain of this study, including the separator, casing, and sub-plunger part of the downhole pump.b The path of the liquid and not separated gas inside the domain

Fig. 17 Fig. 18
Fig. 17 Steps of structural mesh generation for the sophisticated separator using MATLAB.Generating a a circle and b a cylinder.c Deleting the wall of the separator, inlet, and outlet holes, d baffles e the helix from the initial cylinder

Fig. 19 Fig. 20
Fig.19 Velocity contour when the plunger reached the end of the upward stroke (30th second of the simulation)

Fig. 23
Fig.23 Results of the gas-liquid separation in the portals and the helix section of the sophisticated separator during second up and down strokes 2.1.1,the separator effectiveness can also be represented in the fraction of the pump barrel filled with liquid.The percentage of the pump barrel filed with each phase is also shown in Fig. 27.

Fig. 24 Fig. 25
Fig. 24 Results of the gas-liquid volume fraction in the baffles and inner tubing of the sophisticated separator.a During the second downstroke, b Third upstroke, c Streamlines colored by Z velocity contour in the 50th second

Fig. 26 a
Fig. 26 a Gas and b liquid flow rates through different sections of the studied domain

Table 1
Details of generated grids for mesh study

Table 2
Simulation conditions for comparing the mesh types

Table 3
Comparison of the details of generated grids and their simulation results

Table 6
Summary of the model settings used to simulate the simple separator

Table 7
Comparison of laboratory data and simulation results

Table 9
Well a properties

Table 11
Meshes generated with fluent meshing

Table 12
Summary of the sophisticated separator simulation model