Regular Article
A finiteelement algorithm for Stokes flow through oil and gas production tubing of uniform diameter
^{1}
School of Engineering, University of Aberdeen, AB24 3FX Aberdeen, UK
^{2}
Department of Chemical Engineering, Faculty of Engineering, Eduardo Mondlane University, 257 Maputo, Mozambique
^{*} Corresponding author: l.akanji@abdn.ac.uk
Received:
16
February
2020
Accepted:
28
August
2020
Stokes flow of a Newtonian fluid through oil and gas production tubing of uniform diameter is studied. Using a direct simulation on computeraided design of discretised conduits, velocity profiles with gravitational effect and pressure fields are obtained for production tubing of different inner but uniform diameter. The results obtained with this new technique are compared with the integrated form of the Hagen–Poiseuille equation (i.e., lubrication approximation) and data obtained from experimental and numerical studies for flow in vertical pipes. Good agreement is found in the creeping flow regime between the computed and measured pressure fields with a coefficient of correlation of 0.97. Further, computed velocity field was benchmarked against ANSYS Fluent; a finite element commercial software package, in a singlephase flow simulation using the axial velocity profile computed at predefined locations along the geometric domains. This method offers an improved solution approach over other existing methods both in terms of computational speed and accuracy.
© L.T. Akanji & J. Chidamoio, published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
g : Acceleration due to gravity, m/s^{2}
h : Maximum width of geometry, m
Subscript
Superscript
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 multiphase 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 porespace 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., AlAtabi et al. [9]). This has allowed for modelling and simulation of Newtonian and nonNewtonian fluid flow in networks of interconnected pipes for oil transportation, production tubing design or a simplified model of porous media (e.g., AlAtabi 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], OuldRouiss 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 twoway coupling between the particle and the fluid for the investigation of particleladen turbulent flows in low Reynolds number axisymmetric jet. Zhu et al. [21] concluded that the subgrid scale model may be neglected at a spatial resolution of about 10^{6} 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 opensource 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 finiteelements. 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 gaslift 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 ComputerAided Design (CAD) models. The CAD models are constructed with NonUniform Rational BSpline (NURBS) curves and surfaces of order 3, capable of capturing curvatures with tolerancebased level of detail. Geometries are meshed using unstructured spatially variable adaptive grid, tracking freeform 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 FEMbased 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) postprocessing and model visualisation. A flowchart for the computational approach developed in this work is shown in Figure 1.
Fig. 1 Flowchart for model construction, FEM meshing and numerical computation of fluid pressure and velocity fields in geometric conduits. Geometric models are constructed in Computer Aided Design (CAD) package and meshed using unstructured finite elements. Meshed samples are input into CSMP++ where SuperGroup objects representing aggregates of physical entities in the geometric models are formed. Material properties and initial conditions are then assigned with Interrelations subclass being used in computing variable coefficients. Integral forms of the Partial Differential Equations (PDEs) arising from the FEM are assembled termbyterm using numerical PDE operator classes (e.g., NumIntegral dNT op dN dV). Computed fluid pressure and Stokes velocity are displayed and output to VTK for visualisation. 
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 steadystate 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 infsup 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 FEMbased 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 twostep segregated solution method (or pressure correction type method) in solving the Stokes equations (6) and (7). Basically, a Poissonlike equation:(8)is solved for a function ψ and mapped on the discretised geometric model. The obtained function ψ is used in solving the pressurecorrection 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 zaxis 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 setup
The piecewise finiteelement 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 finiteelement 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 stateoftheart Algebraic MultiGrid Solver (SAMG) [24].
2.4 SAMG setup 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 subparameter ndefault (i.e. default switch). Table 1 displays typical values which are assigned to the secondary parameters by default to all the numerical computation setup. A default value of 10–13 (e.g. ndefault = 10) is used if there are no “critical” positive offdiagonal entries. Interpolation is carriedout 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 setup time (t_{s}), percycle 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).
Typical SAMG computational data for solving the Stokes equation using the default switch ndefault setting. The system solved corresponds to one particular timestep taken from a normal production run to solve the equations on geometries A–D described in Section 2.5 below.
2.5 Geometric model construction
Four cylindrical pipe geometries of varying diameters were constructed using a ComputerAidedDesign (CAD) package tool. The CAD tool allows for NonUniform Rational Bspline (NURBS) curves and surfaces of order 3 to be used in order to accurately capture the curvatures with tolerancebased level of detail, at the edges bottom and top of the production tubing [27, 33]. Absolute tolerance of 1 × 10^{−9} m; relative tolerance of 1 × 10^{−7} percent and angle tolerance of 1 × 10^{−3} 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.
The models dimensions.
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.
Fig. 2 A typical CAD construction of the proposed pipe models of height 3.2 m and inner diameters (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 location of the actual pressure sensors placed on the physical experimental rig. 
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 freeform geometrical entities, such as NURBS, with spatially variable refinement and they can also be generated automatically. For realistic meshes of freeform geometry, the quality of the resulting mesh can be evaluated by using the elementtonode 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.
Fig. 3 A zoom into the lower section of the FEM meshes of the geometric samples of diameter (A) d = 0.04 m, (B) d = 0.06 m, (C) d = 0.075 m, and (D) d = 0.09 m. 
Number of elements in each part of the discretised CAD model.
The mesh quality is characterised by the orthogonal quality and corresponding aspect ratios. Mesh qualities are improved by using highlevel 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, nonmanifold 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.
Fig. 4 Typical histograms of mesh quality for pipe geometries of diameter (A) d = 0.04 m and (B) d = 0.06 m. 
Fig. 5 Typical histograms of mesh quality for pipe geometries of diameter (C) d = 0.075 m and (D) d = 0.09 m. 
Mesh quality indicators.
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 carriedout 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 postprocessed and evaluated. Computed fluid pressure and velocity fields are displayed and output to VTK for visualisation.
3.1 Verification of the singlephase 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 singlephase 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.
Fig. 6 A typical geometric model constructed in this work showing (a) 2D CAD configuration of 0.06 m × 1 m and (b) corresponding FEM mesh with imposed boundary conditions. The mesh is composed of triangular elements and nodes. Flow is specified from inlet to outlet directions with noflow on either sides boundaries. (c) numerical simulation results of pressure field and (d) velocity field along the 2D geometry. 
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^{3}; μ = 1.0e^{−3} 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.
Fig. 7 2D geometric mesh of L/D = 16.7 with 5 elements along the reference slit with flow directions and boundary conditions indicated. Note that the orientation of the pipe is in the vertical direction. 
The computed velocities using linear finite elements are piecewisely constant within each element and they show a stairstep 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.
Fig. 8 Velocity fields for placement of 5, 10, 16, 20 and 32 unstructured finite elements across a reference slit shown in Figure 7 Computed piecewise velocity fields (dash lines) versus continuous profiles (dotted lines) from analytical solutions. 
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 2dimensions.
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^{2} = 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 finiteelements (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 nonhybrid finiteelements (only triangles) along the reference slit respectively in a 2D geometry of 0.06 × 1 m. The benefits of the use of hybrid finiteelements 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 nonhybrid grid.
Fig. 9 Mismatch between numerical and analytical solution of the axial velocity profile along the reference slit for the geometries analysed in Figure 8. The trend line represents the fitting as described by equation (18) with correlation coefficient R^{2} = 0.9613. 
3.2 Verification of the singlephase flow model in 3D vertical pipes
The singlephase 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^{3}]; μ = 1.0 × 10^{−3} [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.
Fig. 10 (a) Cylindrical geometry model and (b) corresponding volume mesh consisting of 310 228 elements and 56 128 nodes. (c) Pressure profile along the 3D vertical pipe with the maximum test inlet pressure of 2.5 Pa and outlet pressure of 1.01 Pa. Inbetween the inlet and outlet pressure is a gradient that is consistent with imposed values and the physical geometry investigated. (d) Velocity profile along the 3D vertical pipe with the imposed no slip boundary conditions on the surface of the pipe and maximum velocity profile at the middle of the pipe corresponding to the parabolic flow profile. 
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%.
Fig. 11 Velocity profile for the numerical simulation obtained from ANSYS Fluent package is shown in (a), while (b) shows the result obtained from this numerical computation approach (CSMP++). 
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.
Fig. 12 Axial velocity distribution along the geometric domain for (a) ANSYS Fluent simulation and (b) CSMP++ simulation. In both cases, (c1)–(c3) represent the axial velocity profile at positions 0.25 m, 0.5 m and 0.75 m along the production tubing. 
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 bottomhole 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 underpredicts the experimentally measured pressure by 0.9% at axial position of L/D = 0.0, overpredicts the measured pressure by 2.73% at axial position of L/D = 18.3, overpredicts the measured pressure by 0.83% at axial position of L/D = 36.7 and underpredicts 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.
Fig. 13 FEM predicted and experimentally measured pressure at 4 measurement points S1, S2, S3 and S4 along the production tubing. 
Fig. 14 Experimentally measured versus FEM predicted (M/P) pressure at 4 points S1, S2, S3 and S4 along the production tubing. 
Fig. 15 Pressure field along vertical production tubing. Model height is fixed at 3.2 m but inner diameters D are: (a) 0.04 m, (b) 0.06 m, (c) 0.075 m and (d) 0.09 m respectively. In each of (a–d) the referenced positions along the pipe correspond to the location of the actual pressure sensors (S1, S2, S3, S4) placed on the physical experimental rig. The lengthtodiameter L/D ratios along the pipes are as indicated on the geometric models. 
Fig. 16 Velocity field along vertical production tubing. Model height is fixed at 3.2 m but inner diameters D are: (a) 0.04 m, (b) 0.06 m, (c) 0.075 m and (d) 0.09 m respectively. In each of (a–d) the referenced positions along the pipe correspond to the location of the actual pressure sensors (S1, S2, S3, S4) placed on the physical experimental rig. The lengthtodiameter L/D ratios along the pipes are as indicated on the geometric models. 
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 twodimensional 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 NonITerAtive (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 PCSIMPLE 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, PCSIMPLE 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, multicore and multinode computations on single personal computer, workstations or compute nodes. Further, ANSYS Fluent uses pressurevelocity coupling algorithms whereas in this work, the computation of pressure on each node followed by velocity field computation in a postprocess 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.
Fig. 17 Computational time (s) versus number of finiteelements for CSMP++ and other solvers in ANSYS software package: FluentSimple, CFX15, CFX4.3 and FluentNITA simulations. Numerical simulations were carried out using different computing system specifications. IBM RS/6000 workstation with Power 3 II processor was used for CFX4.3 simulation; 3.2 GHZ Intel i7 processor with 32 GB RAM and 4 cores was used for FluentSimple, CFX15 and FluentNITA; Intel core v5.1, CPU 3.20GHz, and 8 GB memory, with 4 cores was used for CSMP++ simulations. Computation becomes prohibitively expensive as the number of elements increases in ANSYS Fluent simulation. 
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 noncommercial 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 gaslift mandrels and reservoir inflow performance applications will be conducted as part of future research.
References
 Galdi G.P., Robertson A.M. (2005) On flow of a navier–stokes fluid in curved pipes. Part I: Steady flow, Appl. Math. Lett. 18, 10, 1116–1124. doi: 10.1016/j.aml.2004.11.004. [Google Scholar]
 Lahbabi A., Chang H.C. (1986) Flow in periodically constricted tubes: Transition to inertial and nonsteady flows, Chem. Eng. Sci. 41, 2487–2505. [Google Scholar]
 Bernabé Y., Olson J.F. (2000) The hydraulic conductance of a capillary with a sinusoidally varying crosssection, Geophys. Res. Lett. 27, 2, 245–248. doi: 10.1029/1999GL010842. [Google Scholar]
 Payatakes A.C., Tien C., Turian R.M. (1973) A new model for granular porous media: Part i. model formulation, AIChE J. 19, 1, 58–67. doi: 10.1002/aic.690190110. [Google Scholar]
 Payatakes A.C., Tien C., Turian R.M. (1973) A new model for granular porous media: Part ii. numerical solution of steady state incompressible newtonian flow through periodically constricted tubes, AIChE J. 19, 1, 67–76. doi: 10.1002/aic.690190111. [Google Scholar]
 Koplik J. (1982) Creeping flow in twodimensional networks, J. Fluid Mech. 119, 219–247. doi: 10.1017/S0022112082001323. [Google Scholar]
 Chidamoio J., Akanji L., Rafati R. (2017) Prediction of optimum length to diameter ratio for twophase fluid flow development in vertical pipes, Adv. Pet. Explor. Develop. 14, 1, 1–17. [Google Scholar]
 Hernandez A. (2016) Fundamentals of gas lift engineering: Well design and troubleshooting, Elsevier. [Google Scholar]
 AlAtabi M., AlZuhair S., Chin S.B., Luo X.Y. (2006) Pressure drop in laminar and turbulent flows in circular pipe with baffles – an experimental and analytical study, Int. J. Fluid Mech. Res. 33, 4, 303–319. doi: 10.1615/InterJFluidMechRes.v33.i4.10. [CrossRef] [Google Scholar]
 AlAtabi M., AlZuhair S., Chin S.B., Luo X.Y. (2010) Flow of nonnewtonian fluids in porous media, J. Polym. Sci. 48, 23, 2437–2467. [CrossRef] [Google Scholar]
 Sochi T. (2015) Flow of navierstokes fluids in cylindrical elastic tubes, J. Appl. Fluid Mech. 8, 2, 181–188. [CrossRef] [Google Scholar]
 Sabooniha E., Rokhforouz M.R., Ayatollahi S. (2019) Porescale investigation of selective plugging mechanism in immiscible twophase flow using phasefield method, Oil Gas Sci. Technol.  Rev. IFP Energies nouvelles 74, 78. doi: 10.2516/ogst/2019050. [CrossRef] [Google Scholar]
 Sisavath S., Jing X., Zimmerman R.W. (2001) Creeping flow through a pipe of varying radius, Phys. Fluids 13, 10, 2762–2772. doi: 10.1063/1.1399289. [CrossRef] [Google Scholar]
 Tilton J.N., Payatakes A.C. (1984) Collocation solution of creeping newtonian flow through sinusoidal tubes: A correction, AIChE J. 30, 6, 1016–1021. doi: 10.1002/aic.690300628. [Google Scholar]
 Deiber J.A., Schowalter W.R. (1979) Flow through tubes with sinusoidal axial variations in diameter, AIChE J. 25, 4, 638–645. doi: 10.1002/aic.690250410. [Google Scholar]
 Hemmat M., Borhan A. (1995) Creeping flow through sinusoidally constricted capillaries, Phys. Fluids 7, 9, 2111–2121. doi: 10.1063/1.868462. [CrossRef] [Google Scholar]
 Sirisup S., Karniadakis G.E., Saelim N., Rockwell D. (2004) Dns and experiments of flow past a wired cylinder at low reynolds number, Eur. J. Mech. – B/Fluids 23, 1, 181–188. doi: 10.1016/j.euromechflu.2003.04.003. [CrossRef] [Google Scholar]
 Ge M., Xu C. (2010) Direct numerical simulation of flow in channel with timedependent wall geometry, Appl. Math. Mech. 31, 1, 97–108. doi: 10.1007/s104830100110x. [Google Scholar]
 OuldRouiss M., RedjemSaad L., Lauriat G. (2009) Direct numerical simulation of turbulent heat transfer in annuli: Effect of heat flux ratio, Int. J. Heat Fluid Flow 30, 4, 579–589. doi: 10.1016/j.ijheatfluidflow.2009.02.018. [Google Scholar]
 Tamano S., Itoh M., Hoshizaki K., Yokota K. (2007) Direct numerical simulation of the dragreducing turbulent boundary layer of viscoelastic fluid, Phys. Fluids 19, 7, 75–106. doi:10.1063/1.2749816. [Google Scholar]
 Zhu Z., Yang H., Chen T. (2009) Direct numerical simulation of turbulent flow in a straight square duct at reynolds number 600, J. Hydrodyn. Ser. B 21, 5, 600–607. doi: 10.1016/S10016058(08)601900. [CrossRef] [Google Scholar]
 Li D., Fan J., Luo K., Cen K. (2011) Direct numerical simulation of a particleladen low reynolds number turbulent round jet, Int. J. Multiph. Flow 37, 6, 539–554. doi: 10.1016/j.ijmultiphaseflow.2011.03.013. [CrossRef] [Google Scholar]
 Li B., Liu N., Lu X. (2006) Direct numerical simulation of wallnormal rotating turbulent channel flow with heat transfer, Int. J. Heat Mass Trans. 49, 5, 1162–1175. doi: 10.1016/j.ijheatmasstransfer.2005.08.030. [CrossRef] [Google Scholar]
 Stüben K. (2001) A review of algebraic multigrid, J. Comput. Appl. Math. 128, 1, 281–309. doi: 10.1016/S03770427(00)005161. ID: 271610. [Google Scholar]
 Guet S., Ooms G., Oliemans R.V.A., Mudde R.F. (2003) Bubble injector effect on the gaslift efficiency, AIChE J. 49, 9, 2242–2252. [Google Scholar]
 Batchelor G.K. (1967) An introduction to fluid dynamics, Cambridge University Press. [Google Scholar]
 Akanji L.T., Matthai S.K. (2010) Finite elementbased characterization of porescale geometry and its impact on fluid flow, Transp. Porous Media 81, 2, 241–259. doi: 10.1007/s1124200994007. [Google Scholar]
 Babuvška I., Rheinboldt W.C. (1978) Error estimates for adaptive finite element computations, SIAM J. Numer. Anal. 15, 4, 736–754. doi: 10.1137/071504. [Google Scholar]
 Smith I.M., Griffiths D.V., Margetts L. (2013) Programming the finite element method, John Wiley & Sons. [Google Scholar]
 Matthäi S.K., Geiger S., Roberts S.G. (2001) Complex Systems Platform: CSP3D3.0 user’s guide, ETH Zurich Research Collection. [Google Scholar]
 Matthai S.K., Roberts S.G. (1996) The influence of fault permeability on singlephase fluid flow near faultsand intersections: Results from steadystate highresolution models of pressuredriven fluid flow, Assoc. Pet. Geol. Bull. 80, 11, 1763–1779. [Google Scholar]
 Garcia X., Akanji L.T., Blunt M.J., Matthai S.K., Latham J.P. (2009) Numerical study of the effects of particle shape and polydispersity on permeability, Phys. Rev. E 80, 021304. doi: 10.1103/PhysRevE.80.021304. [Google Scholar]
 Akanji L.T., Nasr G., Matthai S.K. (2013) Estimation of hydraulic anisotropy of unconsolidated granular packs using finite element methods, Int. J. Multiphys. 7, 2, 153–166. [Google Scholar]
 Chidamoio J.F. (2018) Experimental and numerical modelling of gaslift cavitation and instabilities in oil producing wells, PhD Thesis, Petroleum Engineering Division, University of Aberdeen. [Google Scholar]
 Anderson A.E., Ellis B.J., Weiss J.A. (2007) Verification, validation and sensitivity studies in computational biomechanics, Comput. Meth. Biomech. Biomed. Eng. 10, 3, 171–184. [CrossRef] [Google Scholar]
 Shao N., Salman W., Gavriilidis A., Angeli P. (2008) Cfd simulations of the effect of inlet conditions on taylor flow formation, Int. J. Heat Fluid Flow 296, 1603–1611. [Google Scholar]
 Kroll N., Gerhold Th., Melber S., Heinrich R., Schwarz Th., Schöning B. (2002) Parallel large scale computations for aerodynamic aircraft design with the German CFD system megaflow, in Wilders P., Ecer A., Satofuka N., Periaux J., Fox P. (eds), Parallel Computational Fluid Dynamics 2001, NorthHolland, Amsterdam, pp. 227–236. doi: 10.1016/B9780444506726/500803. [CrossRef] [Google Scholar]
 Inc ANSYS (2013) ANSYS FLUENT 12.0 user’s guide, ANSYS, Canonsburg, Pennsylvania, United States. [Google Scholar]
 Fletcher D.F., McClure D.D., Kavanagh J.M., Barton G.W. (2017) Cfd simulation of industrial bubble columns: Numerical challenges and model validation successes, Appl. Math. Model. 44, 25–42. doi: 10.1016/j.apm.2016.08.033. [Google Scholar]
 Frey P.J., Alauzet F. (2005) Anisotropic mesh adaptation for cfd computations, Comput. Meth. Appl. Mech. Eng. 194, 48–49, 5068–5082. [CrossRef] [Google Scholar]
 Alauzet F., Loseille A. (2016) A decade of progress on anisotropic mesh adaptation for computational fluid dynamics, Comput.Aid. Design 72, 13–39. doi: 10.1016/j.cad.2015.09.005. [CrossRef] [Google Scholar]
 Papoutsakis A., Sazhin S.S., Begg S., Danaila I., Luddens F. (2018) An efficient adaptive mesh refinement (amr) algorithm for the discontinuous galerkin method: Applications for the computation of compressible twophase flows, J. Comput. Phys. 363, 399–427. doi: 10.1016/j.jcp.2018.02.048. [Google Scholar]
 Alauzet F., Loseille A., Olivier G. (2018) Timeaccurate multiscale anisotropic mesh adaptation for unsteady flows in CFD, J. Comput. Phys. 373, 28–63. doi: 10.1016/j.jcp.2018.06.043. [Google Scholar]
All Tables
Typical SAMG computational data for solving the Stokes equation using the default switch ndefault setting. The system solved corresponds to one particular timestep taken from a normal production run to solve the equations on geometries A–D described in Section 2.5 below.
All Figures
Fig. 1 Flowchart for model construction, FEM meshing and numerical computation of fluid pressure and velocity fields in geometric conduits. Geometric models are constructed in Computer Aided Design (CAD) package and meshed using unstructured finite elements. Meshed samples are input into CSMP++ where SuperGroup objects representing aggregates of physical entities in the geometric models are formed. Material properties and initial conditions are then assigned with Interrelations subclass being used in computing variable coefficients. Integral forms of the Partial Differential Equations (PDEs) arising from the FEM are assembled termbyterm using numerical PDE operator classes (e.g., NumIntegral dNT op dN dV). Computed fluid pressure and Stokes velocity are displayed and output to VTK for visualisation. 

In the text 
Fig. 2 A typical CAD construction of the proposed pipe models of height 3.2 m and inner diameters (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 location of the actual pressure sensors placed on the physical experimental rig. 

In the text 
Fig. 3 A zoom into the lower section of the FEM meshes of the geometric samples of diameter (A) d = 0.04 m, (B) d = 0.06 m, (C) d = 0.075 m, and (D) d = 0.09 m. 

In the text 
Fig. 4 Typical histograms of mesh quality for pipe geometries of diameter (A) d = 0.04 m and (B) d = 0.06 m. 

In the text 
Fig. 5 Typical histograms of mesh quality for pipe geometries of diameter (C) d = 0.075 m and (D) d = 0.09 m. 

In the text 
Fig. 6 A typical geometric model constructed in this work showing (a) 2D CAD configuration of 0.06 m × 1 m and (b) corresponding FEM mesh with imposed boundary conditions. The mesh is composed of triangular elements and nodes. Flow is specified from inlet to outlet directions with noflow on either sides boundaries. (c) numerical simulation results of pressure field and (d) velocity field along the 2D geometry. 

In the text 
Fig. 7 2D geometric mesh of L/D = 16.7 with 5 elements along the reference slit with flow directions and boundary conditions indicated. Note that the orientation of the pipe is in the vertical direction. 

In the text 
Fig. 8 Velocity fields for placement of 5, 10, 16, 20 and 32 unstructured finite elements across a reference slit shown in Figure 7 Computed piecewise velocity fields (dash lines) versus continuous profiles (dotted lines) from analytical solutions. 

In the text 
Fig. 9 Mismatch between numerical and analytical solution of the axial velocity profile along the reference slit for the geometries analysed in Figure 8. The trend line represents the fitting as described by equation (18) with correlation coefficient R^{2} = 0.9613. 

In the text 
Fig. 10 (a) Cylindrical geometry model and (b) corresponding volume mesh consisting of 310 228 elements and 56 128 nodes. (c) Pressure profile along the 3D vertical pipe with the maximum test inlet pressure of 2.5 Pa and outlet pressure of 1.01 Pa. Inbetween the inlet and outlet pressure is a gradient that is consistent with imposed values and the physical geometry investigated. (d) Velocity profile along the 3D vertical pipe with the imposed no slip boundary conditions on the surface of the pipe and maximum velocity profile at the middle of the pipe corresponding to the parabolic flow profile. 

In the text 
Fig. 11 Velocity profile for the numerical simulation obtained from ANSYS Fluent package is shown in (a), while (b) shows the result obtained from this numerical computation approach (CSMP++). 

In the text 
Fig. 12 Axial velocity distribution along the geometric domain for (a) ANSYS Fluent simulation and (b) CSMP++ simulation. In both cases, (c1)–(c3) represent the axial velocity profile at positions 0.25 m, 0.5 m and 0.75 m along the production tubing. 

In the text 
Fig. 13 FEM predicted and experimentally measured pressure at 4 measurement points S1, S2, S3 and S4 along the production tubing. 

In the text 
Fig. 14 Experimentally measured versus FEM predicted (M/P) pressure at 4 points S1, S2, S3 and S4 along the production tubing. 

In the text 
Fig. 15 Pressure field along vertical production tubing. Model height is fixed at 3.2 m but inner diameters D are: (a) 0.04 m, (b) 0.06 m, (c) 0.075 m and (d) 0.09 m respectively. In each of (a–d) the referenced positions along the pipe correspond to the location of the actual pressure sensors (S1, S2, S3, S4) placed on the physical experimental rig. The lengthtodiameter L/D ratios along the pipes are as indicated on the geometric models. 

In the text 
Fig. 16 Velocity field along vertical production tubing. Model height is fixed at 3.2 m but inner diameters D are: (a) 0.04 m, (b) 0.06 m, (c) 0.075 m and (d) 0.09 m respectively. In each of (a–d) the referenced positions along the pipe correspond to the location of the actual pressure sensors (S1, S2, S3, S4) placed on the physical experimental rig. The lengthtodiameter L/D ratios along the pipes are as indicated on the geometric models. 

In the text 
Fig. 17 Computational time (s) versus number of finiteelements for CSMP++ and other solvers in ANSYS software package: FluentSimple, CFX15, CFX4.3 and FluentNITA simulations. Numerical simulations were carried out using different computing system specifications. IBM RS/6000 workstation with Power 3 II processor was used for CFX4.3 simulation; 3.2 GHZ Intel i7 processor with 32 GB RAM and 4 cores was used for FluentSimple, CFX15 and FluentNITA; Intel core v5.1, CPU 3.20GHz, and 8 GB memory, with 4 cores was used for CSMP++ simulations. Computation becomes prohibitively expensive as the number of elements increases in ANSYS Fluent simulation. 

In the text 