CFD Analysis of Turbulent Flow of Power-Law Fluid in a Partially Blocked Eccentric Annulus

: This study focuses on analyzing the turbulent flow of drilling fluid in inclined wells using the Computational Fluid Dynamics (CFD) technique. The analysis is performed considering an annulus with a fixed eccentricity of 90% and varying fluid properties, diameter ratio, and bed thickness to examine velocity profile, pressure loss, and overall wall and average bed shear stresses. CFD simulation results are compared with existing data for validation. The pressure loss predicted with CFD agrees with the data. After verification, predictions are used to establish a correlation that can be applied to compute bed shear stress. The established correlation mostly displays a discrepancy of up to 10% when compared with simulation data. The correlation can be used to optimize hole cleaning and manage downhole pressure.


Introduction
Hydrocarbon production from horizontal and inclined wells is now a widespread practice. One of the major drawbacks of drilling horizontal and inclined wells is that the inclination causes wellbore eccentricity. This eccentric wellbore geometry causes a nonuniform fluid velocity profile in the annulus creating low-velocity stagnant regions to develop. The formation of the stagnant zones makes the cuttings removal process inefficient, leading to the development of stable solids beds. The formation of the bed restricts fluid flow and increases the pressure loss and bottom hole pressure, which needs to be maintained within the operating pressure window to prevent drilling problems.
Failure in ensuring a good hole cleaning leads to several drilling issues such as formation damage, fluid loss, lost circulation, high torque and drag, or stuck pipe. Therefore, a wellbore hydraulic study is vital in establishing a precise model to compute the pressure loss in the annulus and optimize hole cleaning. Several studies [1][2][3][4][5][6][7] have been performed to examine the flow in an annular region and develop effective hole cleaning solutions and manage wellbore pressure within the operating window. For an effective hole cleanup, the enhancement of bed shear stress and local fluid velocity on the bed surface is vital. If the flow rate is increased, the solids bed thickness would decrease. However, an increase in pressure due to the flow rate increase can offset the decline in frictional resistance achieved because of bed thickness reduction. Therefore, an accurate bed shear stress model is required for hydraulic optimization. The model should account for the bed shear stress variation occurring because of bed thickness variation.
Common flow parameters in the wellbore such as bed shear stress and pressure loss are affected by mudflow rate, fluid properties, and wellbore geometry (eccentricity, bed height, and diameter ratio). The development of solids beds makes the wellbore flow geometry complicated; consequently, analytical solutions are not possible to establish when a stationary solids bed is formed. A CFD investigation is frequently preferred approach to study the connection among relevant flow parameters (pressure loss, bed shear stress, and flow rate) in eccentric annuli. Previously, CFD models [7][8][9][10][11][12] for partially blocked annuli were established considering laminar flow conditions. Similar studies [13][14][15][16][17] were performed to analyze laminar flows in concentric and eccentric annuli.
Based on the CFD analysis of eccentric annular flow, Rojas et al. [7] showed the impact of the level of annular blockage on the bed shear stress under laminar flow conditions. The bed shear stress exhibited significant reduction with a decrease in the bed thickness. Results also demonstrated substantial variation in the bed shear stress in the lateral direction, which causes non-uniform bed erosion. Similar results were obtained for a concentric annulus [8]. The results indicated the reduction of local fluid velocity and bed shear stress with an increase in the blockage level. These studies demonstrated that the impact of blockage on relevant flow parameters such as pressure loss and bed shear stress can be predicted with reasonable accuracy with the application of CFD techniques.
The flow behavior of most drilling fluids used in the industry is best described using the power-law rheology model. The flow regime in inclined wells can be laminar or turbulent depending on fluid viscosity, borehole geometry, and fluid velocity. Even the effects of bed formation on the hydraulics of annular flows have been studied recently [7][8][9][10][11][12], the impact of bed thickness on the hydraulic parameters (pressure loss and bed shear stress) of turbulent flow in a partially blocked eccentric annulus has not been previously investigated. The goal of this investigation is to study the turbulent flow of power-law fluid in wellbores considering eccentric annuli with blockage and develop a relationship between bed shear stress and the level of blockage.

Theoretical Background
For non-Newtonian fluid, the Reynolds number in a concentric annulus is defined as a function of apparent viscosity, which is determined by equating the non-Newtonian pressure loss equation with that of a Newtonian fluid. Thus: * = (1) where is the effective diameter, which is determined from the hydraulic diameter of the annulus as: = 0.81 . is fluid consistency index, is average fluid velocity and is the apparent viscosity of the fluid. The hydraulic diameter for any non-circular duct is generally expressed as: = 4 ⁄ , where is the wetted perimeter of the annulus, and represents the cross-sectional area of the flow. The apparent viscosity of power-law fluid ( = ) in a concentric annulus is computed as: Equation (2) is adopted for partially blocked annulus using the hydraulic diameter concept. The Fanning friction factor ( ) is determined from the mean fluid velocity ( ) and the overall wall shear stress ( ̅ ) as: where is the fluid density. Applying the momentum balance, the average/overall wall shear stress is calculated from the pressure gradient as: According to Kozicki et al. [18], dimensionless geometric parameters can be utilized to determine the Reynolds number to estimate the pressure loss in a non-circular duct of any arbitrary cross-section. For a partially blocked eccentric annulus, Rojas et al. [7] presented correlations for determining the geometric parameters. Kozicki et al. [18] proposed a generalized Reynolds number for non-circular duct flows. Thus: where and are the dimensionless parameters and is the fluid behavior index. In inclined wellbores, the development of solids beds influences the wall shear stress and velocity profiles [9]. Thus, the bed shear stress (i.e., the average shear stress acting on the bed) could significantly vary from the overall wall shear stress. The bed shear stress ( ̅ ) governs the movement of bed particles by affecting the hydrodynamic forces (lift and drag) acting on the bed particles. By determining the pressure loss, the overall wall shear stress ( ̅ ) can be predicted and hydraulic optimization can be performed to ensure efficient hole cleaning in highly deviated wells [7,8,19,20]. Previous studies [10,12] were conducted to investigate the impact of blockage on annular flows. The bed shear stress [7,21] is a relevant flow parameter for accurate hydraulic modeling of the wellbore. The dimensionless bed shear stress ( ) is defined as [7]: The bed shear stress is defined in this form to reduce the number of flow parameters that affect its magnitude.

Modeling Flow in Partially Block Eccentric Annulus
The research methodology utilized to achieve the principal objective of this study includes CFD analysis of a partially blocked eccentric annular flow using commercial software (ANSYS Fluent) to establish a relationship between the bed shear stress and the blockage thickness. To conduct the analysis, several 3D models of annular geometries (varying geometric dimensions and blockage level) are created to simulate and examine the bed shear stress and frictional pressure loss. The simulation is performed considering turbulent flow conditions and varying flow behavior of the fluid and the annular geometry. Using results obtained from CFD simulations, a simplified correlation is developed to calculate the bed shear stress in eccentric annuli as a function of bed height. The correlation is valid for power-law fluids flowing in highly eccentric annuli (approximately 90%) when the bed thickness is between zero and 100 percent of the inner pipe diameter.
To develop a theoretical model, the following assumptions are made: (i) the flow is considered as a fully developed turbulent flow under steady-state condition; (ii) the fluid is considered incompressible and flowing under isothermal conditions, (iii) the flow behavior of the fluid is described using the power-law fluid rheology model; iv) non-rotating inner pipe; (v) no-slip condition at the surface of the bed and other wetted surfaces; (vi) stable and uniform solids bed; (vii) no impact of suspended solids on the flow; and (viii) all wetted surfaces are smooth.

Governing Equation
A variety of turbulent modeling approaches have been established to determine the approximate solutions for the equations of motion. One of these models is the k-ε model that is often applied to model turbulent flows in uniform ducts such as eccentric annulus [22] and open channels. In the standard k-ε model, the turbulent kinetic energy (k) distribution is described as: Similarly, the turbulent energy dissipation rate (ε) is described using the following differential equation.
where is the turbulent kinetic energy generated because of the mean velocity gradient.
is fluctuating dilatation in incompressible turbulence to the overall dissipation rate.
is the turbulent kinetic energy developed due to the buoyancy. and are the turbulent Prandtl numbers for and , respectively. In the standard k-ε model, the Reynolds stresses are obtained using the eddy viscosity ( ) that is given by: For power-law fluid, the viscous shear stress is computed from the shear rate ( ) as: = .
The apparent viscosity ( ( ) = / ) of non-Newtonian fluids is expressed in terms of the shear rate, which is defined in a general form as: This relationship helps evaluate the stress terms in the equation of motion. , , and are time-averaged local fluid velocities in the x, y, and z directions. Using these velocities, the Reynolds Averaged Navier-Stokes (RANS) equation for incompressible fluid can be described in a generalized form as: where Fi is the body force. ui is the time-averaged local fluid velocity in the direction xi. The continuity equation for incompressible fluid is described as: The model equations (Equations (7)- (12)) are solved numerically in the CFD software using the finite volume method implemented on a computational grid system with reasonable grid size.

Geometry Creation and Mesh Generation
A precise representation of the flow geometry is essential in performing extensive CFD simulations. The first step in a CFD modeling process is the flow geometry creation, which is followed by mesh generation (Figure 1). Numerous flow geometries are created varying bed height (0 to 100%) and diameter ratio ( = / ). To cover a wide range of wellbore geometries, three diameter ratios (25, 50, and 75%) are considered. The relative eccentricity ( = ( − ) ⁄ ) considered for all the simulation cases is 90%, where and are the diameters of the inner and outer pipes, respectively; and d is the offset distance between the centers of the pipes. The second step after the creation of the flow geometry is the appropriate mesh generation. The grids (mesh) are generated to represent the entire flow domain and ensure the required numerical precision, which is determined by the number of grid points. The numerical accuracy and grid resolution increase with the number of grids. However, the number of grids needs to be optimized to avoid very long computational time. The unstructured grid system ( Figure 1) is utilized in this study to ensure accuracy with a reasonable computational time.

Boundary Conditions
Implementing the correct boundary conditions is paramount before beginning a numerical simulation. Boundary conditions are essential to represent real scenarios of flow conditions at the flow boundaries. Five named surfaces ( Figure 2) are considered as flow boundaries, which are grouped into three subcategories as follows: Inlet (Fluid Inflow): The inlet face is named as inlet boundary (velocity inlet) located at the upstream end of the simulated wellbore at which the velocity profile is assumed uniform.
Outlet (Fluid Outflow): A pressure outlet boundary condition is defined at the 3D geometries outflow, which is located downstream of the simulated wellbore.
Wall Boundaries: The no-slip boundary condition is implemented in all walled surfaces, namely Bed Surface, Inner Surface, and Outer Surface in the CFD solver. The volume enclosed by these three wall boundaries (named surfaces) forms the flow domain (simulated wellbore).

Numerical Setup in Fluent Solver
The CFD final result needs to be independent of the initial conditions. Thus, a steadystate pressure-based technique is chosen to model turbulent flow in annuli with stationary cuttings bed. The technique is one of the most reliable and commonly used methods to model and analyze turbulent fluid flow. Apart from this, the pressure gradient effect and enhanced wall treatment options are selected to maintain the required numerical accuracy at the flow boundaries. Furthermore, the Semi-Implicit Pressure Linked Equation (SIM-PLE) algorithm is implemented in combination with the second-order discretization method for momentum, turbulent dissipation rate, turbulent kinetic energy, and pressure. The convergence limit for the continuity equation is maintained at 1 × 10 −4 and the convergence requirements for momentum, kinetic energy, dissipation rate equations are set at 1 × 10 −7 [26].

Grid Sensitivity Analysis
Mesh size sensitivity analysis is performed to ensure the accuracy of CFD simulations. Input parameters for the grid-independent study are presented in Table 1. The input parameters are the specified boundary conditions. The grid system considered for this analysis is based on different grid refinement levels (0.5, 1.25, 2.5, and 5 million grids/elements). All simulations are conducted maintaining a high Reynolds number (Re* > 14,000). Figure 3 demonstrates an increase in accuracy with the number of grids that improves numerical accuracy. After maintaining a mesh of size of more than 2.2 million grids, the pressure gradient stabilizes [25].

Validation of CFD Simulation Results
To validate the accuracy of the simulation approach implemented, the CFD results are compared with available pressure gradient measurements obtained from eccentric annulus (8 in × 4.5 in) with a recorded average bed thickness [27,28]. The measurements are obtained for diverse bed heights, fluid properties (Table 2), and flow rates. Four waterbased muds are considered. During the experiment, the flow rate was between 150 and 400 gal/min. The Reynolds numbers (Re*) of the flows are determined to distinguish the flow regime. Figure 4 compares experimental measurements with CFD simulation results. The results are predominately in agreement with measurements. Under laminar flow conditions (Figure 4c), the discrepancies are considerable (around 30%). The predictions are reasonably accurate with a maximum discrepancy of about 17% under turbulent flow conditions. The Reynolds number (Re*) obtained from Equation (1) is valid for a concentric annulus. Hence, it is an approximate indicator of the flow regime when used in partially blocked eccentric annuli. Hence, some of the low Reynolds number flow data points may not represent turbulent flows even though they are more than 2100.

Parametric Study
A parametric study is conducted to examine how the annular pressure loss is affected by the diameter ratio and bed height. Simulations are conducted maintaining high Reynolds numbers and varying the relevant flow parameters including diameter ratio, bed height, and flow behavior index. All the CFD simulation cases considered in this study were performed at a constant flow rate of 5 × 10 −5 m 3 /s. The fluid density and hole diameter are fixed at 1000 kg/m 3 and 0.05 m, respectively.

Annular Frictional Pressure Loss
CFD simulations are performed to calculate the annular pressure loss for dimensionless bed height ( = ℎ/ ) ranging from 0 to 100 percent. h is the bed height measured from the bottom of the outer pipe/cylinder. A constant flow rate is maintained for all cases and the change in annular pressure gradient ( / ) is observed ( Figure 5). With increasing bed height and diameter ratio the pressure loss sharply increased. This phenomenon is primarily attributed to the reduction in flow area caused by the increase in bed height and diameter ratio. The reduction in flow area considerably increases the average velocity, and subsequently, augments the pressure loss. Furthermore, the hydraulic resistance experienced by the fluid flowing through the annulus diminishes with a decrease in fluid behavior index due to the improvement in the shear-thinning tendency of the fluid.

Overall Wall Shear Stress
The average/overall wall shear stress ( ̅ ) is a hydraulic parameter used to evaluate hole cleaning performance. In Figure 6, the overall wall shear stress ( ̅ ) is presented versus dimensionless bed height for various values of fluid behavior index. The wall shear stress in the wellbore increased with bed height. This is anticipated since the wall shear stress and the pressure loss are directly related according to Equation (4). At a constant flow rate, the mean fluid velocity increases with blockage height resulting in a rise in the wall shear stress. Moreover, the wall shear stress is significantly influenced by the fluid behavior index. It reduces considerably with the improvement of the shear-thinning behavior of the fluid.  Figure 7 displays the average bed shear stress versus the dimensionless bed height for different fluid behavior indexes and diameter ratios. When the flow rate is constant, the bed shear stress increases significantly with the diameter ratio and bed height. This is predominantly due to the reduction in the flow area, which increases the mean fluid velocity. Besides this, the increase in bed shear stress is party due to the bed approaching the high-velocity core zone of the flow (Figure 8). Consistent with this observation, a recent experimental study [29] demonstrated that the cleaning of eccentric annuli becomes more challenging as the bed gets smaller and the bed shear stress diminishes. Consequently, the Reynolds stresses which help the lifting of bed particles significantly reduces as the bed thickness decreases [30]. The bed shear stress also slightly increases with the fluid behavior index.   Figure 9 presents the dimensionless bed shear stress as a function of dimensionless bed height (Hbed), fluid behavior index, and diameter ratio. The bed shear stress increases with bed height and approaches the overall bed shear stress when the Hbed reaches 100% (i.e., fully buried inner pipe case). This demonstrates the weakening of the flow stagnation effect which causes the bed shear stress to diminish once the bed surface is below the inner pipe. In addition, all the bed shear stress curves converge as the Hbed approaches 100%, showing the reduction in the influence of fluid behavior index on dimensionless bed shear stress.

Modeling Bed Shear Stress
Average bed shear stress prediction is essential to optimize hole cleaning process and manage wellbore pressure. The dimensionless bed shear stress ( ) is useful to predict the bed shear stress from the overall wall shear stress, which is calculated from the pressure gradient as shown in Equation (4). The impacts of flow and fluid parameters on the dimensionless bed shear stress are shown in Figure 9. Other fluid and flow parameters that are not shown in the figure such as fluid consistency index, flow rate, Reynolds number can be relevant; and therefore, their impacts on are investigated. These parameters are varied to perform a sensitivity analysis under turbulent flow conditions. The effect of flow rate, Reynolds number, and fluid consistency index on the dimensionless bed shear stress is minor. Applying a non-linear regression analysis, a correlation is developed for the dimensionless bed shear stress. The correlation relates the bed shear stress to the bed height, diameter ratio, and flow behavior index. The simulation data generated during this investigation is used in the development of the correlation. Thus, the dimensionless bed shear stress is expressed as:

Conclusions
A CFD study is conducted under turbulent flow conditions to investigate wellbore hydraulics for a power-law fluid in inclined and horizontal wells. Extensive simulations are conducted varying the cuttings bed height, fluid behavior index, and wellbore diameter ratio. The main contribution of this study is the development of the bed shear stress model by analyzing the simulation results. Besides this, the following conclusions can be made based on the findings of this investigation.

•
The dimensionless bed shear stress is a crucial flow parameter in eccentric annuli, which is primarily affected by the bed height, diameter ratio, and fluid behavior index.

•
The bed shear stress is a function of dimensionless bed shear stress and overall shear stress, which is affected by mean fluid velocity, wellbore size, diameter ratio, the thickness of the cuttings bed, and properties of the fluid.

•
The thickness of the cuttings bed strongly affects the annular pressure gradient. Its impact is more pronounced when the dimensionless bed height is greater than 50%. This is primarily because of the decrease in the flow cross-section area. • CFD simulation of turbulent flow provides a reasonable prediction of relevant flow parameters including pressure gradient, bed shear stress, and velocity profile.

Recommendations for Future Research
The following topics are recommended for future studies in the area of wellbore hydraulics and hole cleaning in inclined wells and better understand non-Newtonian fluid flow in partially blocked eccentric annuli: Acknowledgments: This publication was made possible by an NPRP award [NPRP11S-1228-170140] from the Qatar National Research Fund (a member of The Qatar Foundation). We would like to express our gratitude and appreciation to the University of Oklahoma and Qatar University for supporting the project.

Conflicts of Interest:
The authors declare no conflict of interest.