© L.T. Akanji & J. Chidamoio, published by IFP Energies nouvelles, 2020

Nomenclature

V : Volume, m³

F : Right-hand term

N : Interpolation function

P : Pressure, Pa

W : Weighting function

μ : Fluid viscosity, Pa s

ψ : Parabolic function

ρ : Density, kg/m³

ω : Domain boundary

: Unit normal to the boundary

g : Acceleration due to gravity, m/s²

U : Fluid velocity, m/s

h : Maximum width of geometry, m

L : Length of geometry, m

D : Diameter of geometry, m

ε : Flow field mismatch

Subscript

x, y, z:

Superscript

T : Matrix transpose

e : Element

1 Introduction

Fluid flow through pipes has attracted a great deal of scientific and engineering discipline interests due to their prevalence in biological systems such as blood flow in human body and multi-phase in flow in oil and gas producing wells. Flow in porous media is often modelled as periodically assembled tubes (e.g., Galdi and Robertson [1], Lahbabi and Chang [2], Bernabé and Olson [3], Payatakes et al. [4, 5]). In this case, the total pore-space can be represented by a network of tubes with each tube being assigned a hydraulic conductance C, defined by analogy with the electrical conductance [6]. In oilfield applications, several authors have studied the effect of production tubing diameter on oil production and vertical lift performance (e.g., Chidamoio et al. [7]). Hernandez [8] reported that large production tubing found in old oil fields where reservoir pressure has declined to very low values may be replaced by small diameter tubing in order to produce at a stable rate.

Several solutions exist for creeping flow of a Newtonian fluid through vertical pipes. Approximate solutions have been obtained using the Hagen–Poiseuille (i.e., lubrication approximation) for flow in pipes by integrating the equation at each axial location along the pipe to obtain the overall pressure drop (e.g., Al-Atabi et al. [9]). This has allowed for modelling and simulation of Newtonian and non-Newtonian fluid flow in networks of interconnected pipes for oil transportation, production tubing design or a simplified model of porous media (e.g., Al-Atabi et al. [10], Sochi [11], Sabooniha et al. [12]). Navier–Stokes equation has also been applied using analytical methods such as asymptotic series solution (e.g., Sisavath et al. [13]) and numerical methods such as the collocation method (e.g., Tilton and Payatakes [14]), geometric iteration method (e.g., Deiber and Schowalter [15]) and the boundary integral method (e.g., Hemmat and Borhan [16]).

Although, results obtained from analytical solutions to the Navier–Stokes equation generally offer improvement over the lubrication approximation with regard to the details of the flow field, numerical approaches are preferable in cases where detailed flow structure is desired [13]. Numerical methods are usually in good agreement with experimental data and have generally produced very good results. However, computational costs may be prohibitively high.

The Direct Numerical Simulation (DNS) methods have been used to obtain accurate turbulent flow field in a pipe (e.g. Sirisup et al. [17] Ge and Xu [18], Ould-Rouiss et al. [19], Tamano et al. [20], Zhu et al. [21], Li et al. [22, 23]). Ge and Xu [18] developed a numerical scheme capable of extending the scope of the spectral method without solving the covariant and contravariant forms of the Navier–Stokes equations in the curvilinear coordinates. They represented the primitive variables by the Fourier series and the Chebyshev polynomials in the computational space. Li et al. [22] used an explicit coupling scheme between particles and the fluid, which considers two-way coupling between the particle and the fluid for the investigation of particle-laden turbulent flows in low Reynolds number axisymmetric jet. Zhu et al. [21] concluded that the sub-grid scale model may be neglected at a spatial resolution of about 10⁶ and a duct flow at the particular friction Reynolds number of 600 will be appropriate. Laboratory experiments usually give a detailed structural behaviour of flow in production tubing, but, numerical approaches are cheaper and preferable in preliminary investigation and analysis.

These aforementioned DNS methods have focussed mainly on the modelling of flow structure and less on the actual engineering of the physical problem particularly in the numerical evaluation of the velocity fields within realistic geometries. Further, numerical computation on realistic representation of geometric entities is still a challenge. Commercial software packages (such as ANSYS Fluent) and some other open-source applications are used in numerical simulation of Stokes flow in pipes. However, there are significant limitations in such computation on representative oil and gas production tubing discretised with millions of finite-elements. Such investigations are important given that production tubing is usually completed in production casing or liner perforated for oil and gas production from hydrocarbon reservoirs where inflow performance is evaluated. Such production tubing systems may also have gas-lift mandrels which have to be incorporated in the numerical model design and computation. The workflow developed herein aims to proffer solution to this problem. The workflow allows for detailed flow structure; a caveat in the analytical solutions approach, to be captured whilst providing an opportunity to simulate flow directly on physical representation of long production tubing where computational speed is desirable. In this work, flow simulation is conducted on physical geometric models of production tubing of varying diameters using Computer-Aided Design (CAD) models. The CAD models are constructed with Non-Uniform Rational B-Spline (NURBS) curves and surfaces of order 3, capable of capturing curvatures with tolerance-based level of detail. Geometries are meshed using unstructured spatially variable adaptive grid, tracking free-form entities such as NURBS. Stokes flow simulation are then conducted on production tubing of different diameters to compute pressure and associated velocity fields respectively. Computation is implemented in complex systems modelling platform (CSMP++); an application programmer interface engineered in ANSI/ISOC++. The developed model is verified against analytical solution and then validated against experimental and numerical results.

2 Methodology

A FEM-based algorithm for Stokes flow in production tubing is formulated by solving the discretised equation using an Algebraic MultiGrid Solver (SAMG); an efficient linear solver library based on algebraic multigrid system [24]. The workflow involves: (i) building the geometric model with CAD in 2D and 3D domains, (ii) meshing the geometric samples using unstructured mesh consisting of line, triangle and tetrahedron elements, (iii) assigning boundary and initial conditions, (iv) discretising the flow equations, (v) solving equations using algebraic multigrid solver and, (vi) post-processing and model visualisation. A flowchart for the computational approach developed in this work is shown in Figure 1.

2.1 Governing equations

The flow of fluid in a pipe can be described by the general form of the momentum equation, constrained by the mass conservation equation. The momentum equation can be written as:(1)where, ρ and U are the density and velocity respectively. In equation (1), the term a represents the rate of increment in momentum per unit volume; b is the change in momentum due to convection; c is the pressure gradient; d represents the viscous and turbulent contributions; e is the gravitational forces. The mass conservation equation can be expressed as:(2)

In this investigation, the pipe fluid pressure is computed and compared with analytical solutions and data obtained from laboratory measurements. Furthermore, flow is assumed incompressible and Newtonian, thus:(3)

The fully developed flow is assumed to be achieved when the inertia term is small compared to the viscous term; low Reynolds number (e.g., Guet et al. [25], Batchelor [26], Akanji and Matthai [27]), hence, the inertia term can be ignored and equation (3) can be rewritten as:(4)where,(5)is the reduced pressure gradient. We consider a steady-state flow in the production tubing where at any fixed point in space, the velocity does not vary with time; hence ignoring transient effects and assuming that the fluid density is constant, equation (4) reduces to the Stokes equation. Thus the momentum and mass conservation equations can be written as:(6)(7)

In computational fluid mechanics, the correct choice of the shape functions for velocity and pressure can be adjusted in the case of an isothermal and incompressible flow through the inf-sup compatibility condition of the Ladyzhenskaya–Babuska–Brezzi (LBB condition) [28]. The equivalence of LBB for all shape functions involving velocity, pressure, temperature and electromagnetic fields is not known for the cases of temperature deviation and electromagnetism. The focus of this work is to adopt a robust FEM-based computational technique for practical field applications involving oil, gas and/or water flow in production tubing without exploiting the LBB condition.

2.2 Discretisation of the momentum equation

In the FEM technique, the continuous variables, velocity U and pressure P are approximated by the variables and respectively. Here, nodal fluid pressure P_i is computed through simple functions of spatial variable known as shape functions N_j and velocity fields are approximated at the barycentre of each finite element. Computational steps involve: multiplication of the residual of the Partial Differential Equations (PDE) by a weighting function; integration by parts; representation of the approximate solution as a linear combination of polynomial basis functions defined on a given mesh; substitution of the functions in the weak formulation; solving the resulting algebraic system for the nodal variables (see for instance Smith et al. [29]). In this Bubnov–Galerkin FEM, the weighting functions are taken as the same interpolation functions.

We adopt a two-step segregated solution method (or pressure correction type method) in solving the Stokes equations (6) and (7). Basically, a Poisson-like equation:(8)is solved for a function ψ and mapped on the discretised geometric model. The obtained function ψ is used in solving the pressure-correction equation:(9)to obtain the pressure field along the production tubing. The velocity field is then obtained thus:(10)

Flow in the vertical direction along the z-axis is considered and the numerical integration equivalent to the discretised form of the PDEs (6) and (7) can be written thus:(11)and,(12)

Similar numerical integration expressions written for equations (8)–(10) and adopted in this work are shown in Figure 1. Where, the index i is the coordinate in physical space and a summation convention applies and,(13)and,(14)

Equations (11) and (12) can be written in matrix form as:(15)

Further details on discretisation of the above PDEs can be found in [27]. The matrix expressed by the discretisation of the above PDEs (Eqs. (8)–(10)) is conditioned in PDE operator as described in Matthäi et al. [30] and highlighted in Figure 1.

2.3 Numerical solution and code set-up

The piecewise finite-element integrals that resulted from the governing equations combine into systems of linear algebraic equations of the form A·x = b. The discretisation of the geometries will often require millions of finite-element nodes and the matrix A will contain millions of equations. This therefore means that an efficient matrix solver will be required if details of the flow behaviour are to be adequately captured in a computationally effective manner. In order to meet this requirement, Matthai and Roberts [31], Garcia et al. [32] and Akanji and Matthai [27] have demonstrated that the FEM form of the fluid pressure equation can be solved rapidly by Algebraic Multigrid Methods (AMG) because A is sparse, symmetrical and positive definite. In this solution approach, vector x is initialised with a trial solution . A·x = b is then restricted to coarser grids and the trial solution is smoothened until A·x = b can be solved directly using LU matrix decomposition. The obtained solution is then interpolated back onto finer grids to obtain an improved and accurate trial solution through repeated smoothening, coarsening and interpolation. A key feature of AMG is the ability to reduce the error components on both the low and high part of the eigenspectrum of matrix A efficiently. Further, it leads to a convergence rate which is independent of the mesh size thereby allowing for optimal execution times for systems of linear problems with billions of unknowns. In this work, the pressure field in the Stokes equation is computed using the state-of-the-art Algebraic MultiGrid Solver (SAMG) [24].

2.4 SAMG set-up parameters

SAMG’s initial dimensioning is configured by setting some default primary parameters (see [24]). The secondary parameters are usually accessed by using the variable iswtch with the sub-parameter ndefault (i.e. default switch). Table 1 displays typical values which are assigned to the secondary parameters by default to all the numerical computation set-up. A default value of 10–13 (e.g. ndefault = 10) is used if there are no “critical” positive off-diagonal entries. Interpolation is carried-out in the simplest manner. The larger the second digit of ndefault the more aggressive coarsening becomes. Aggressive coarsening basically reduces the memory requirements but at the expense of a slower convergence. The set-up time (t_s), per-cycle time (t_c) and total time (t_t) are typically indicative of the computing time associated with the models A–D (described in Sect. 2.5 below).

2.5 Geometric model construction

Four cylindrical pipe geometries of varying diameters were constructed using a Computer-Aided-Design (CAD) package tool. The CAD tool allows for Non-Uniform Rational B-spline (NURBS) curves and surfaces of order 3 to be used in order to accurately capture the curvatures with tolerance-based level of detail, at the edges bottom and top of the production tubing [27, 33]. Absolute tolerance of 1 × 10⁻⁹ m; relative tolerance of 1 × 10⁻⁷ percent and angle tolerance of 1 × 10⁻³ degrees were applied in order to differentiate all the discernible features in the models. The dimensions of the cylindrical pipes geometries of varying diameter are presented in Table 2.

Figure 2 shows the CAD models constructed for the purpose of this investigation. The height of the models is 3.2 m and inner diameters are (a) 0.04 m (b) 0.06 m (c) 0.075 m and (d) 0.09 m. S1, S2, S3 and S4 are the axial positions of the reference slits corresponding to the location of the actual pressure sensors placed on the physical experimental rig developed as part of a complementary pilot scale research project (see Chidamoio [34]). It is important to note that one of the advantages of this developed technique is in the ability to adequately resolve geometric entities with adequate degree of realism.

2.6 FEM discretisation of geometric models – meshing

The geometric models are discretised using unstructured mesh consisting of lines, triangles or quadrilaterals for the surface section and tetrahedral or hexahedra for the volumetric inner space. The unstructured grids can fit free-form geometrical entities, such as NURBS, with spatially variable refinement and they can also be generated automatically. For realistic meshes of free-form geometry, the quality of the resulting mesh can be evaluated by using the element-to-node ratio. A value close to 2 can be obtained for hybrid meshes when compared with 5–6 for pure tetrahedral meshes, [33]. The number of nodes and elements in each part of the meshed geometric models are shown in Table 3 and a zoom into the lower part of the geometries is shown in Figure 3.

The mesh quality is characterised by the orthogonal quality and corresponding aspect ratios. Mesh qualities are improved by using high-level diagnostic smoothening and modification algorithm in Integrated Computer Engineering and Manufacturing (ICEM) meshing tool. Typical problems that may be associated with low mesh quality include single, multiples edges, triangle boxes, overlapping elements, non-manifold and unconnected vertices. The sample mesh quality indicators for the models shown in Figures 2 and 3 are presented in Table 4 and the corresponding histogram of sample mesh quality is shown in Figures 4 and 5. Values closer to 1 are of the highest mesh quality.

3 Numerical simulation

The numerical simulation workflow developed in this work is shown in Figure 1. Samples of 2D and 3D geometric models were constructed and used as part of the verification process in both 2D channels and 3D cylindrical geometries. In all cases geometric models are built with CAD, meshed and fluid flow computation carried out based on imposed boundary and initial conditions assignment. The flow computation is carried-out based on imposed boundary and initial conditions assignment in CSMP++ platform where material properties and initial conditions are assigned with Interrelations subclass and used in computing variable coefficients. Obtained results are then post-processed and evaluated. Computed fluid pressure and velocity fields are displayed and output to VTK for visualisation.

3.1 Verification of the single-phase flow model in 2D vertical pipes

In order to verify the degree of accuracy of the solution approach developed in this work, numerical computations are carried out on increasingly refined grid cells in vertical pipes. Numerical computations involving discretisation of geometric and mathematical models are approximations with errors which generally decrease as the grid is refined. The order of the approximation is a measure of the degree of accuracy of the numerical computation; therefore, evaluation of the mesh quality is an important aspect of this process [35]. In order to verify the code developed for this single-phase flow model, verification exercise is carried out using 2D cylindrical geometries shown in Figure 6. This verification is constrained by the imposed boundary and initial conditions.

3.1.1 Boundary and initial conditions in 2D geometries

A typical 2D geometric model sample constructed using CAD models with Dirichlet boundary conditions applied to the bottom and top boundaries is presented in Figure 6a. In this case, the bottom boundary is the inlet and the top boundary is the outlet and no flow boundary conditions are assigned to the left and right sides of the geometry. Several geometries are constructed and meshed with, for instance, 8380 elements and 4189 nodes (see Fig. 6b). The fluid used for this test is water with physical properties of (ρ = 998 kg/m³; μ = 1.0e⁻³ Pa s) [36]. The computed numerical solution of the pressure and velocity fields are presented in Figures 6c and 6d respectively.

3.1.2 Comparison with analytical solution

The number of elements required to obtain realistic parabolic velocity profile is tested. In this case, five numerical experiments are carried out using a 2D geometry of 0.06 m × 1 m flow model presented in Figure 6 for steady state single phase flow. The numerical solution of equations (6) and (7) with density assumed equal to 1 and fluid viscosity kept equal to 0.001 Pa s is compared with the analytical solution of the velocity represented by equation:(16)where, y denotes the width of the 2D geometry, h is the maximum width of the geometry, dP/dx is the pressure drop along the pipe length. For the purpose of this analysis, five (5) sensitivities on mesh refinement along a reference slit placed at a designated point on a geometric pipe of L/D = 16.7 (see Fig. 7) were carried out. The number of elements at the reference slit varies from 5, 10, 16, 20 and 32 with an inlet and outlet pressures of 2.5 and 1.01 Pa respectively. Numerical computation of axial velocities on each linear element along the slit is then evaluated and compared with equivalent analytically computed velocities.

The computed velocities using linear finite elements are piece-wisely constant within each element and they show a stair-step parabolic profiles in contrast to the smooth analytical solution. An adequate resolution can be achieved based on how well the numerical flow velocity profile matches the analytical solution. Radial profile of the axial velocity and the corresponding pressure fields for each mesh sensitivity are presented in Figure 8.

The flow field mismatch between the analytical and numerical methods is estimated thus:(17)where, U_a is the flow field from the analytical solution to equation (16) and U_n is the flow field from the numerical solution to equations (6) and (7) in 2-dimensions.

As observed from Figure 9, the velocity field mismatch decreases with increasing mesh refinement with an exponential decay rate which can be expressed as:(18)where, x is the number of elements along the reference slit, a and b are constants; determined for this analysis as a = 27.638 and b = −0.071. Figure 9 shows the representation of this fitting with a correlation coefficient of R² = 0.9613. This trend is in agreement with the study of [27] where a maximum and minimum velocity mismatch of ≈22% and ≈1% were observed with 7 and 21 hybrid finite-elements (mixed elements consisting of triangles and quadrilaterals) across the reference slits respectively for the 2D channel geometry of 30 × 10 μm. In this study, a maximum and minimum velocity mismatch of 17.4% and 2.9% were observed with 5 and 32 non-hybrid finite-elements (only triangles) along the reference slit respectively in a 2D geometry of 0.06 × 1 m. The benefits of the use of hybrid finite-elements become more pronounced at higher refinements where a slight error reduction may be observable. Similar observation has been reported in [37] where the use of hybrid grid provided a slightly lesser computational time compared to non-hybrid grid.

3.2 Verification of the single-phase flow model in 3D vertical pipes

The single-phase algorithm that solves the Stokes equation is tested on 3D geometric pipes with dimensions of 0.06 m diameter and 1.0 m length. Figures 10a and 10b illustrates the CAD geometries and the corresponding FEM meshes consisting of 310 228 elements and 56 128 nodes. The fluid used for this test is water with physical properties of (ρ = 998 [kg/m³]; μ = 1.0 × 10⁻³ [Pa s]) [36]. Dirichlet boundary conditions were applied to the bottom and top sides of the cylinder, while no slip boundary conditions were assigned on the perimeter of the pipe. An inlet pressure of 2.5 Pa and an outlet pressure of 1.01 Pa were assigned to the model as shown in Figures 10c and 10d. The numerically computed simulation output of pressure field and the corresponding velocity fields along the pipe are presented in Figures 10c and 10d respectively. The numerical simulation results show that the fluid pressure and velocity fields are consistent with the boundary and initial conditions assigned and correspond to direction of fluid flow along the pipe.

4 Model validation and benchmarking with commercial software package

In order to validate the developed model, numerical computation was carried out on 2D geometries using the numerical solution approach described in Section 2 and results compared with corresponding computations obtained from commercial software package (ANSYS Fluent). Two comparisons with commercial simulator involving (i) computation of axial velocity field profiles along the pipe length and (ii) computation of parabolic velocity field profiles at predefined reference slits along the pipe length, were carried out. Further, validation was conducted by comparing computed results with experimentally measured data. The laboratory measured data were obtained as part of a complementary investigation on flow in pipe research project (see Chidamoio [34]).

4.1 Benchmarking of CSMP++ simulation with ANSYS Fluent simulations – axial velocity profiles

In this case, fluid pressure has been assigned as Dirichlet boundaries at both the inlet and outlet sections of the pipe. The 2D model shown in Figure 6 is used for the purpose of this analysis. The model is a 0.06 m width and 1.0 m length pipe where the inlet is located at the bottom and the outlet at the top. In ANSYS Fluent, a Dirichlet boundary condition is applied by assigning at the inlet while the outlet boundary is placed at the top of the tubing where a Neumann pressure boundary condition is applied on the outlet. The outlet pressure is set equal to one bar, meaning that the outlet is open to the environment and the only pressure acting at the outlet is atmospheric pressure [38]. The outer surfaces of the tubing are treated as wall with no slip and no penetration , where and denote the unit tangent on the boundary and unit normal to the boundary respectively. Figures 11a and 11b show the results of the velocity field computed from ANSYS Fluent simulations and from this work respectively. The two results are in very good agreement and within an average relative error of 0.9%.

4.2 Benchmarking of CSMP++ simulation with ANSYS Fluent simulations – parabolic velocity profiles

Computation of parabolic velocity field profiles at predefined reference slits along the pipe length is carried out in this section using similar geometric configuration as described in Section 4.1. Figures 12a–12c show the velocity fields and parabolic velocity profiles along predefined reference slits located at positions 0.25 m, 0.50 m and 0.75 m on the geometric domain for the ANSYS Fluent simulation and this work (CSMP++) respectively. Both codes have a similar trend on axial velocity distribution with a slight difference in variable magnitude.

4.3 Model validation – comparison with experimental measurements

To validate this work, measured pressure data along fixed positions on the vertical pipe (Chidamoio [34]) are compared with computed pressure data obtained in this work. The experimental test section of the rig consisting of 0.06 m diameter and 3.2 m length, was replicated using numerical computation approach developed in this work (see Figs. 2b and 3b). The numerical tests were carried out along the axial distance of the test sections of L/D of 0; 18.3; 36.7 and 53.3 respectively; representing the location of the pressure transducers along the physical experimental production tubing. Dirichlet inflow and out flow boundaries conditions for pressure were assigned based on the experimental measurements of bottom-hole pressure of 133 251.64 Pa (absolute pressure) using pressure transducer located at L/D = 0.0 and tubing head pressure of 102 303.90 Pa (absolute pressure) using the pressure transducer located at L/D = 53.3. Number of elements and mesh quality of the CAD are presented in Tables 3 and 4 respectively.

The values presented in Figure 12Yb were obtained at the reference slits S1, S2, S3 and S4 in the VTK pressure field outputs corresponding to the axial positions of the pressure transducers in the experimental production tubing rig. FEM under-predicts the experimentally measured pressure by 0.9% at axial position of L/D = 0.0, over-predicts the measured pressure by 2.73% at axial position of L/D = 18.3, over-predicts the measured pressure by 0.83% at axial position of L/D = 36.7 and under-predicts the measured pressure by 0.97% at axial position of L/D = 53.3. An overall pressure drop of 575.98 Pa/m was obtained for this particular geometry. FEM solution versus experimentally measured pressure values were plotted and the results presented in Figure 13. The trend of the results obtained is observed to be within reasonable degree of accuracy in comparison with the experimental results. This can be attributed to the very high mesh refinement used in the numerical simulation. As can be seen in Figure 13, a strong correlation between FEM predicted pressure profile and experimentally measured pressure data is observed with a coefficient of correlation of 0.97. Thus, it can be concluded that, the present simulation approach can accurately predict the pressure along the production tubing. Further investigation involving the impact of production tubing diameter on liquid production rate was assessed using the developed FEM numerical approach (Fig. 14). Tubing length is fixed at 3.2 m and internal diameters are (a) 0.04 m, (b) 0.06 m, (c) 0.075 m and (d) 0.09 m and the corresponding to L/D ratios are indicated on the pipes. The pressure fields for each of the pipes are shown in Figure 15 and the velocity fields shown in Figure 16.

4.4 Discussions

Shao et al. [36] applied Volume Of Fluid (VOF) method in ANSYS Fluent to simulate the effect of inlet conditions on Taylor flow formation in 1 mm ID capillary with two gas nozzle sizes of 0.11 mm and 0.34 mm ID. The solution domain was constructed as two-dimensional axisymmetric geometry with length of 3 mm ID and the total number of elements was 20 000. A IBM RS/6000 workstation with Power 3 II processor was used to run the simulation and the average running time for each case was 5 days.

Fletcher et al. [39] carried out numerical simulation with emphasis on the various solution methods involving the use of a new Non-ITerAtive (NITA) Eulerian multiphase solver available in ANSYS Fluent 16. The model is a bubble column of 2 m height and 0.39 m diameter. The same computational mesh comprising 36 000 hexahedral cells, is used in the ANSYS CFX simulations and PC-SIMPLE solver in ANSYS Fluent. The computations were evaluated in a 3.2 GHZ intel i7 processor with 32 GB RAM and 4 cores. For this geometry, computation time of 21.4, 92.6, and 484 s were recorded for NITA solver, PC-SIMPLE solver and ANSYS CFX respectively. As part of the strategy for achieving high accuracy, convergence order and improved CPU time, mesh adaptation has been proposed by a group of authors (e.g. Frey and Alauzet [40], Alauzet and Loseille [41], Papoutsakis et al. [42], Alauzet et al. [43]).

In order to evaluate the computational resources used in this numerical simulation approach, the computing capabilities were evaluated in a single processor, Intel core v5.1, CPU 3.20GHz, and 8 GB memory, with 4 cores. The computation time for each case studied was recorded and plotted. Figure 17 shows the computation time versus number of elements for each model presented in Table 3. Maximum of 17 s of computation time was achieved for 1 536 648 elements. From this, it can be concluded that the simulation approach proposed in this work is highly efficient in comparison with ANSYS Fluent which is computationally demanding. This efficiency stems from the fact that the numerical computation carried out in this work uses SAMG solver which is a very efficient linear solver library based on Algebraic MultiGrid (AMG) system. It supports both serial, multi-core and multi-node computations on single personal computer, workstations or compute nodes. Further, ANSYS Fluent uses pressure-velocity coupling algorithms whereas in this work, the computation of pressure on each node followed by velocity field computation in a post-process operation significantly enhances the computational efficiency of this method. Computational efficiency is crucial when conducting flow simulation in long pipes discretised with millions of elements as the use of ANSYS Fluent and CFX in such cases become computationally prohibitive.

5 Conclusion

A rigorous numerical workflow for the direct simulation of fluid flow in oil and gas production tubing was developed. This method is embedded into a novel simulation approach that begins with CAD modelling and discretisation of the geometry and ends with an efficient numeric solution of momentum and conservation equations in a non-commercial package. Contrary to existing commercial simulation packages, this technology allows for computation in physical geometric production tubing models discretised with millions of element. The following conclusions can be drawn from this work:

• The 2D and 3D single phase flow code based on Stokes equation was formulated and tested. The results from the numerical simulation test cases show that the fluid pressure field is consistent with the boundary conditions applied by the user and the velocity profile results showed that the Stokes equation parabolic profile is maintained.

• The computed pressure and velocity fields were verified against analytical solution and benchmarked with ANSYS Fluent using the same geometries and assigned variable properties. A satisfactory agreement between the codes was observed.

• The 3D flow model was validated by quantitative comparison with the experimentally measured pressure. FEM pressure field output was plotted against experimentally measured pressure along predefined reference points. Under the prevailing conditions, a strong correlation between the predicted and experimentally measured pressure was observed with a coefficient of correlation of 0.97. The impact of pipe diameter impact on liquid production rate was assessed using the formulated FEM numerical approach and for varying tubing diameters.

• Extension of the formulations to production tubing systems with gas-lift mandrels and reservoir inflow performance applications will be conducted as part of future research.

References

All Tables

All Figures