#### Abstract

Estimation of pressure losses and deposition velocities is vital in the hydraulic design of annular drill holes in the petroleum industry. The present study investigates the effects of fluid velocity, fluid type, particle size, particle concentration, drill string rotational speed, and eccentricity on pressure losses and settling conditions using computational fluid dynamics (CFD). Eccentricity of the drill pipe is varied in the range of 0–75%, and it rotates about its own axis at 0–150 rpm. The diameter ratio of the simulated drill hole is 0.56. Experimental data confirmed the validity of current CFD model developed using ANSYS 16.2 platform.

#### 1. Introduction

In the petroleum industry, predicting frictional pressure losses and settling conditions for the transportation of drilling fluids in the annuli are important for drilling operations. Inaccurate predictions can lead to a number of costly drilling problems. A few examples of such problems are loss of circulation, kicks, blockage, wear, abrasion, and improper rig power selection. The existing empirical models become less accurate as those involve many simplified assumptions. CFD simulations help to minimize such assumptions by using the physics-based Navier–Stokes equations to model the hydrodynamics of the flow system. Current work is focused on developing a comprehensive CFD model which is capable of considering the effects of all important drilling parameters, such as fluid velocity, fluid type, particle size, particle concentration, drill pipe rotation speed, and drill pipe eccentricity.

#### 2. Literature Review

The estimation of pressure loss in an annulus is more difficult compared with pipe flow due to the complexities in hydraulics resulting from the complex geometry [1, 2]. From an empirical perspective, the issue is usually addressed by replacing pipe diameter in the pipe flow models with an “effective diameter” of annulus. A number of definitions of “effective diameter” have been proposed till date. However, it is difficult to select a definition for a field application as those were developed and/or applied empirically. A comparison of multiple definitions in predicting pressure losses is presented by Anifowoshe and Osisanya [3]. Other issues that make the estimation of pressure losses in drilling holes difficult are the eccentricity and the rotational speed of inner drill pipe. Many studies have been done on the flow of non-Newtonian fluids in annuli to introduce empirical/analytical models which allow to take these effects into account [1, 3–10]. The results of the previous studies show that the annular pressure losses for non-Newtonian (power law) fluids flowing in a drill hole depend on drill pipe rotation speed, fluid properties, flow regimes (laminar/transitional/turbulent), diameter ratio, eccentricity, and equivalent hydrodynamic roughness.

Using commercially available CFD packages like ANSYS FLUENT to predict pressure losses for the annular transportation of drilling fluids is comparatively a new approach. Sorgun and Ozbayoglu [9] demonstrated the better performance of CFD model compared with the existing empirical models in predicting frictional pressure losses. Sorgun [8] investigated the effect of pipe eccentricity on pressure loss, tangential velocity, axial velocity, and effective viscosity by using CFD. Erge et al. [6] presented a CFD modeling approach which is applicable to estimate frictional pressure losses in an eccentric annulus with inner pipe rotation while circulating yield power law fluids. However, most of these CFD studies were limited to the laminar flow of a single phase in hydro dynamically smooth annuli.

Annular flow of drilling fluids containing cuttings, i.e., slurry, has not been studied in sufficient detail. Examples of works in the field of annular slurry flow are available in references [11–14]. The focus of these studies was to understand the hydrodynamics of the slurry flow in annulus from real-time experiments and to produce empirical models based on data analysis. Recently, different researchers [15–17] have used CFD in studying the transportation of slurry in annuli. Ofei [15] examined the effect of the rheological parameters of the carrying fluids on the velocity of solid. Sorgun and Ulker [16] compared the predictions of pressure losses obtained using artificial neural network (ANN) and CFD. Both methods produced comparable results. Sun et al. [17] studied the effects of inclination, rotational speed, and flow rate on the distribution of solid concentration and the frictional pressure loss. Similar to the single-phase annular flow works, most of these CFD studies were limited to the laminar slurry flow conditions.

#### 3. Methods

##### 3.1. CFD Simulation

In the current work, the CFD simulation model of the annular slurry flow is developed using ANSYS Fluent 16.2 platform. Following previous works [18–20], a multi-fluid granular model is used to describe the flow behavior of a fluid-solid mixture. The granular version of Eulerian model is selected as the multiphase model (Appendix C). This is because high solid volume fraction is expected to be used for this study and the granular version captures the hydrodynamics of high concentration slurries consisting of varying grain sizes. It allows modeling of multiple separate but interacting phases. The phases can be liquids, gases, or solids in nearly any combination. The Eulerian treatment is used for each phase, in contrast to the Eulerian–Lagrangian treatment that is used for the discrete phase model.

The description of multiphase flow as interpenetrating continua incorporates the concept of phasic volume fractions, which represent the space occupied by each phase. Each phase satisfies the laws of conservation of mass and momentum individually. The conservation equations are modified by averaging the local instantaneous balance for each of the phases [21] or by using the mixture theory approach [22]. A detailed description is available in Appendix A.

For Eulerian multiphase calculations, the phase-coupled SIMPLE (PC-SIMPLE) algorithm is used for the pressure-velocity coupling. PC-SIMPLE is an extension of the SIMPLE algorithm to multiphase flows [23, 24]. The velocities coupled by phases are solved in a segregated fashion. The block algebraic multigrid scheme used by the density-based solver is used to solve a vector equation formed by the velocity components of all phases simultaneously [25]. Then, a pressure correction equation is built based on the total volume continuity. Pressure and velocities are then corrected so as to satisfy the continuity constraints.

To ensure stability and convergence of iterative process, a second-order upwind discretization was used for momentum equation, and first upwind discretization was employed for volume fraction, turbulent kinetic energy, and its dissipation. Upwinding refers to the face value derived from quantities in the cell upstream, or “upwind,” relative to the direction of the normal velocity. When first-order accuracy is desired, quantities at cell faces are determined by assuming that the cell-center values of any field variable represent a cell-average value and hold throughout the entire cell; the face quantities are identical to the cell quantities. Thus, the face value is set equal to the cell-center value of the upstream cell when first-order upwinding is selected. In contrast, when second-order accuracy or second-order upwinding is desired, quantities at cell faces are computed using a multidimensional linear reconstruction approach [26]. In this approach, higher order accuracy is achieved at cell faces through a Taylor series expansion of the cell-centered solution about the cell centroid. The total simulation process is shown in Figure 1.

##### 3.2. Turbulence Model Selection

Turbulent quantities for fluid flow are computed using Reynolds stress model [27–29] (Appendix D). Abandoning the isotropic eddy-viscosity hypothesis, the RSM closes the Reynolds-averaged Navier–Stokes equations by solving transport equations for the Reynolds stresses, together with an equation for the dissipation rate [30]. Here, five additional transport equations are required in 2D flows and seven additional transport equations must be solved in 3D (please see Appendix B for further details). The turbulence model was selected for the current work, analyzing the relative performance of different models. A typical example of the analysis is presented in Figure 2. In this figure, % difference refers to the percentile difference between the experimental measurements of pressure loss and the corresponding CFD predictions. In most cases, the difference was less than 10% when RSM was used. A list of the important results is presented in Table 1.

##### 3.3. Length Independence Study

The length of the flow domain is considered long enough to achieve the fully developed flow. Minimum entrance length considered for the flow development is 50*D*_{h}, where hydraulic diameter *D*_{h} = OD–ID [36, 37]. An example of the length independent test is shown in Figure 3. The simulation results were not dependent on length after 3 m from the inlet.

##### 3.4. Mesh Analysis

The computational grids for an annular section are generated using ANSYS Fluent, and the meshing is finalized on the basis of proper mesh independency check. Multiple layers of inflation near wall are added from both inner and outer walls to compute the characteristics of different parameters near wall more precisely. Shear stress between wall surface and fluids is much higher, and this inflation helps to create denser meshing near wall. An example of computational grid distribution and mesh independence test is shown in Figures 4 and 5. Mesh independent results could be produced for more than 800000 number of nodes. All the results presented in the current work were obtained using around 900000 nodes.

**(a)**

**(b)**

**(c)**

**(d)**

The values of dimensionless wall distance (*y*^{+}) were checked during near wall cells generation in consideration of the convergence requirement of *y*^{+} for cells adjacent to the wall. The value of *y*^{+} depends on the wall shear stress, fluid density, hydraulic diameter, and molecular viscosity:where is the distance from the wall to the cell center; , the molecular viscosity; , the fluid density; and , the wall shear stress. Eventually, *y*^{+} depends on the mesh resolution and the flow Reynolds number. Default standard wall functions are generally applicable if the first cell center adjacent to the wall has a *y*^{+} value larger than 30 [38]. In view of the minimum requirement (*y*^{+} > 30), the value of *y*^{+} was maintained above 45 in our study.

###### 3.4.1. Convergence Rate Analysis

An optimum convergence rate of 10^{−5} was selected for the termination of iteration. Figure 6 shows an example of the analysis used to find out the optimum convergence rate within 10^{−6}–10^{−4}. The simulation results varied when the convergence value was greater than 10^{−5}. However, the results did not change at all for the values less than 10^{−5}.

##### 3.5. CFD Model Validation

As the preliminary step, the CFD modeling approach is validated with respect to the data available in open literature. Few examples of the validation are presented as follows.

##### 3.6. Single-Phase Flow through Annuli

Sets of experimental data from Kelessidis et al. [33] and Camçi [32] are compared in Figure 7 with the proposed CFD model. The geometry is taken from Kelessidis et al. (2011), where the inner diameter (ID) is 0.04 m, outer diameter (OD) is 0.07 m, and length is 5 m (horizontal concentric annuli). Fluid is taken as water and wall material is Plexiglas (hydro dynamically smooth wall, *ε*_{a} = 0). In Camçi (2003), the inner diameter is 0.0432 m, outer diameter is 0.123 m, and length is 5 m. Wall material is aluminum (smooth wall).

In Figure 7, the graph represents log-log scale. The pressure gradient increasing rate in log-log scale is almost linear for both cases. The average percentage errors of simulation results from Kelessidis et al. (2011) and Camçi (2003) are 9.88% and 8.46%, respectively, which indicates a very good agreement (it was estimated that the error of the experimental data is approximately ±10%).

##### 3.7. Two-Phase (Solid-Liquid) Flow through Annuli

Pressure gradient (Pa/m) profile of water-sand slurry flow through vertical concentric annuli is compared with Ozbelge and Beyaz [35] experimental data in Figure 8. For the CFD simulation, the liquid phase is considered as water (density 9982 kg/m^{3} and viscosity 0.001003 kg/m-s) and the solid phase is taken as feldspar (mean particle diameter 0.23 mm and mean density 2500 kg/m^{3}). Length 5 m, outer diameter 0.125 m, inner diameter 0.025 m, inlet velocity range of 0.0738–0.197 m/s, and overall slurry volumetric concentration range of 1.0%–1.8% with 0.23 mm grain size (*d*_{p}) are considered as boundary conditions. A smooth pipe of stainless steel (density 8030 kg/m^{3}) is used for the simulation. The pipe is assumed to be vertical, i.e., gravity effect is included, and gravity acceleration is directed opposite to the outlet. No slip condition for liquid and solid phases is used at the walls. Figure 9 shows the comparison of simulated and experimental two-phase frictional pressure drop through vertical annuli at different mixture velocity and at different volume concentration of slurry for 0.23 mm mean sand particle diameter (*d*_{p}). Simulated results are in good agreement with experimental values with 2.62% average error.

#### 4. Results and Discussion

After achieving good validation of proposed model with experimental data, a parametric analysis is done to observe the effect of the changing flow rate, drill pipe rotation, eccentricity, and particle size (detailed data tables are presented in Appendix (A and B)). The parameters used for the analysis are presented in Table 2.

##### 4.1. Effect of Flow Rate

The effect of fluid flow on maximum bed concentration is analyzed in Figures 10–13. Water with sand particle mixture (slurry) is used as operating fluid. Four different conditions are taken into consideration with fix sand inlet concentration (20%) and sand particle size (0.1 mm). The conditions are as follows:(i)Concentric annuli with stationary inner pipe(ii)Concentric annuli with 150 rpm rotating inner pipe(iii)50% eccentric annuli with stationary inner pipe(iv)50% eccentric annuli with 150 rpm inner pipe

From each case of the analysis, it is observed that bed concentration near the bottom wall decreases when the flow rate increases. Because of gravitational force and horizontal orientation, sand particles have the tendency to gather near bottom wall and create a flow blockage. Decreasing flow enhances the process. From the analysis, it is found that below Reynolds number bed concentration near bottom wall is more than 25% in all cases and when it goes down , the percentage is above 50%. In Figure 14, sample contour distribution of sand particles is shown at different flow rates. Operating conditions are taken from Figure 10. Contour distribution of annuli cross section at 3 m distance from inlet displays the effect of fluid flow and gravitational force on concentration distribution clearly.

**(a)**

**(b)**

**(c)**

**(d)**

##### 4.2. Effect of Drill Pipe Rotation and Eccentricity

Inner pipe of annuli can affect pressure loss by changing its eccentricity and rotation. The effect of inner pipe rotation and eccentricity is presented in Figures 9 and 15. Single-phase water is used as fluid in this analysis. From Figure 15, it can be visualized that pressure loss increases with the increase of rotational speed, and this trend is same at different flow rates. Stationary, 50 rpm, 100 rpm, and 150 rpm rotation of inner pipe are taken into consideration. Increasing trend is higher at high flow rates. Because of rotational speed particle–particle and particle–wall collision and rebound increase, which results in higher pressure loss.

Figure 9 shows the effect of inner pipe eccentricity on pressure loss. Concentric, 25% eccentricity, 50% eccentricity, and 75% eccentricity are analyzed in this figure. At a fixed flow rate (Re = 50000), with the increase of inner pipe eccentricity, pressure loss decreases. This trend is applicable during inner pipe stationary condition. With inner pipe rotation, the trend is opposite. Because of extra collision added by pipe rotation, this change occurs.

##### 4.3. Effect of Particle Size

In slurry flow, sand particle size has an effective role on particle blockage near bottom wall of horizontal annular pipe. Figure 16 analyzes the effect of sand particle size at concentric and stationary inner pipe with 5% inlet sand concentration. With the increase of sand particle size, particle deposition near bottom wall increases because individual particle weight increases, which eventually results in this blockage. From previous analysis, we have found maximum bed concentration decreases with the increase in the flow rate (Figures 8, 10–13), but for particle size 0.005 mm, maximum bed concentration near bottom wall is almost constant at different flow rates. That means, for smaller particle size (<0.01 mm), the effect of the flow rate on particle deposition in negligible.

To be noted, these mean particle sizes are selected from particle size distribution (PSD) charts considered during CFD analyzes. One of the PSD charts is shown in Figure 17, with mean particle size 0.1 mm, where 0.1 mm particle size is selected from eight different sizes based on cumulative weight (%) of each particles. To simplify the parametric analysis process, only mean particle sizes are shown.

#### 5. Conclusions

In an effort to develop a widely accepted CFD model of multiphase flow through drilling annuli, the current work relates the overall plan and the initial progress of the project. In summary, this study can be recounted as follows:(i)A CFD modeling methodology to predict frictional pressure loss and settling conditions is validated. The validation is presented with examples which demonstrate its applicability in complex drilling conditions.(ii)The effects of following important drilling parameters on pressure loss and settling conditions are tested: fluid flow rate, rotational speed, and eccentricity of drill pipe and solid particle size.(iii)Fluid flow rate has the maximum impact on particle deposition. In all analyzed conditions, with the decrease of fluid flow, particle deposition near the bottom wall increases. Specific deposition velocity depends on specific requirement of maximum deposition near wall which can be calculated using our approach.(iv)With the increase of rotational speed and eccentricity of inner pipe, the energy loss (pressure loss) of fluid flow increases. However, at stationary condition, pressure loss decreases with the increase of eccentricity.(v)Preliminary results of the current research program are presented with Newtonian fluid flow. The project is expected to produce a comprehensive CFD model capable of considering all important drilling parameters with Newtonian and non-Newtonian fluid flow. Analysis with non-Newtonian fluid flow is ongoing.

#### Appendix

#### A. Parametric Study Chart with Single-Phase Fluid

The results of the parametric study chart with the single-phase fluid are given in Table 3.

#### B. Parametric Study Chart with Two-Phase Fluid

The results of the parametric study chart with the two-phase fluid are given in Table 4.

#### C. Description of Multiphase Model

The volume of phase, , is defined bywhere

The effective density of phase iswhere is the physical density of phase.

The equations for fluid-fluid and granular multiphase flows are presented here for the general case of an *n* -phase flow. The volume fraction of each phase is calculated from a continuity equation as below.where is the phase reference density or the volume averaged density of the phase in the solution domain, characterizes the mass transfer from the to phase, and characterizes the mass transfer from the to phase.

The conservation of momentum for a fluid phase iswhere is the acceleration due to gravity, is the phase stress-strain tensor, is an external body force, is a lift force, and is a virtual mass force.

The conservation of momentum for the solid phase iswhere is the solids pressure, is the momentum exchange coefficient between fluid or solid phase and solid phase , and is the total number of phases.

#### D. Description of Turbulence Model

The exact transport equation of Reynolds stress () is as follows:where is the summation of the changing rate of and transport of by convection. is the production rate of . is the diffusion transport of . is the rate of dissipation. is the pressure strain correlation term. is the rotation term.

The diffusion term used in simulation is as follows:where , .

Production rate of can be expressed as

The pressure-strain term has or slow pressure-strain term also known as the return-to-isotropy term, or rapid pressure-strain term, and as wall-reflection term. It can be presented aswherewhere and .

The wall-reflection term is responsible for the normal stresses distribution near the wall. This term is modeled aswhere is the horizontal component of the component normal to the wall, is the shortest distance to the wall, and, where and is the von Kármán constant (=0.4187).

The dissipation rate or the destruction rate of is modeled aswhereand = fluctuating deformation rate.

Rotation term is expressed aswhere = rotation vector. = alternating symbol, +1, −1, or 0 depending on *i*, j, and *k*.

#### Nomenclature

CFD: | Computational fluid dynamics |

RSM: | Reynolds stress model |

SST: | Shear stress transport |

RANS: | Reynolds averaged Navier–Stokes |

SIMPLE: | Semi-implicit method for pressure linked equations |

d_{m}: | Sand particle diameter |

C_{v}: | Sand volumetric concentration |

ID: | Pipe inner diameter |

OD: | Pipe outer diameter |

D_{H}: | Hydraulic diameter |

v: | Fluid velocity |

: | Water velocity |

: | Sand velocity |

: | Wall roughness |

: | Reynolds number |

: | Friction factor. |

#### Data Availability

The “ANSYS Simulation” data used to support the findings of this study are available from the corresponding author upon request.

#### Disclosure

A conference paper [39] by the same authors was presented with the primary datasets which are updated and modified in this article.

#### Conflicts of Interest

The authors also declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

The authors are thankful to the Faculty of Engineering and Applied Science of Memorial University of Newfoundland, Qatar National Research Fund (the grant NPRP10-0101-170091), Texas A&M University at Qatar, and ANSYS Inc., for helping and funding this research work.