Advances in Numerical Analysis

Volume 2017 (2017), Article ID 2672438, 20 pages

https://doi.org/10.1155/2017/2672438

## Numerical Simulation of Barite Sag in Pipe and Annular Flow

Department of Drilling Engineering, China University of Petroleum (East China), Qingdao 266580, China

Correspondence should be addressed to Mingbo Wang

Received 15 April 2017; Accepted 14 May 2017; Published 10 August 2017

Academic Editor: Guillaume Galliero

Copyright © 2017 Patrick Kabanda and Mingbo Wang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

With the ever increasing global energy demand and diminishing petroleum reserves, current advances in drilling technology have resulted in numerous directional wells being drilled as operators strive to offset the ever-rising operating costs. In as much as deviated-well drilling allows drillers to exploit reservoir potential by penetrating the pay zone in a horizontal, rather than vertical, fashion, it also presents conditions under which the weighting agents can settle out of suspension. The present work is categorized into two parts. In the first part, governing equations were built inside a two-dimensional horizontal pipe geometry and the finite element method utilized to solve the equation-sets. In the second part, governing equations were built inside a three-dimensional horizontal annular geometry and the finite volume method utilized to solve the equation-sets. The results of the first part of the simulation are the solid concentration, mixture viscosity, and a prediction of the barite bed characteristics. For the second part, simulation results show that the highest occurrence of barite sag is at low annular velocities, nonrotating drill pipe, and eccentric drill pipe. The CFD approach in this study can be utilized as a research study tool in understanding and managing the barite sag problem.

#### 1. Introduction

With the ever increasing global energy demand and diminishing petroleum reserves, recent advances in drilling technology have resulted in numerous directional wells being drilled as operators strive to offset the ever-rising operating costs. Deviated-well drilling allows operators to exploit reservoir potential by penetrating the pay zone in a horizontal, rather than vertical, fashion. With consideration of eliminating drilling problems such as stuck pipe, torque and drag, wellbore instability, and low rates-of-penetration, these wells are being drilled increasingly with invert-emulsion drilling fluids.

In the drilling industry, the term “barite sag” refers to the settling of weighting materials in drilling fluids under flowing (pumping) or no flowing (no pumping) conditions. In as much as there exist substitutes such as iron titanium oxide, calcium carbonate, and manganese tetra oxide, barite is used extensively since it provides high density with wide accessibility and favourable economics and it is ecofriendly. Notwithstanding mentioning, the API 13D defines barite sag (in the field) as the change in drilling fluid density observed when circulating bottoms-up; it is almost always characterized by drilling fluid having a density below nominal, being followed by drilling fluid densities above nominal, and being circulated out of the annulus [1].

Modelling and simulation studies employ a combination of both mathematical formulations and numerical techniques where the governing equations for a problem definition are solved by either writing a computer program (code) following a discretization scheme such as finite difference method (FDM) or using CFD software to obtain a solution following a numerical method such as finite element method (FEM) or finite volume method (FVM). The solution obtained is referred to as a numerical solution since a particular numerical procedure is followed and the whole process generally called numerical simulation. Numerous researchers have studied the barite sag phenomenon by employing modelling and simulation studies.

Paslay et al. [2] presented a fundamental theoretical effort to describe dynamic sag in drilling fluids by considering the behaviour of a system of particles in a non-Newtonian Bingham fluid. The authors employed continuum mechanics to develop a model of dynamic sag in an inclined annulus in terms of the fluid and particle properties. They focused on the period when pumping and rotation are stopped. They indicated that the particles settle for a few minutes immediately following the cessation of rotation and pumping when the gelling properties of the stationary drilling fluid have not been fully recovered. The analysis upon which the results were deduced is laminar flow and the predictions appeared to be consistent with the operational guidelines presented by Dye et al. [3].

Nguyen et al. [4] described a fundamental mathematical approach (based mainly on the continuum methodology) to examine the sedimentation of barite particles in shear flow of Newtonian fluids; they established a numerical method to solve four partial differential equation-sets describing dynamic barite sag in pipe flow of Newtonian fluids. They calculated the concentration of solid in both the axial and radial directions as a function of time by using an explicit FDM method.

Hashemian et al. [6, 7] performed a study on barite sag by () modelling of laminar flow of non-Newtonian fluid in annulus to obtain velocity profile and () consideration of the solid particles in the fluid to predict the particle traveling path and time. The simulation was based on a proposed particle tracking method called “Particle Elimination Technique.” It is important to note that the simulation approach employed is much dependent on a parallel experimental study. The estimation of unknown parameters (, mass rate of barite bed back to suspension, , average velocity of the barite bed, and , viscosity of the barite bed) that are input parameters in the simulation is based on the experimental results. This is rather a shortcoming of this approach as it cannot be performed independently and later make a comparison with experimental data.

Kulkarni et al. [8] presented a novel method of predicting real-time sag behaviour in the wellbore by employing a comprehensive computational approach to model the sag behaviour in wellbores using fluids composition (i.e., weight-material (barite) size/concentration and oil-water ratio)/properties (i.e., rheology) and wellbore geometry (i.e., inclination and diameter)/operating conditions (i.e., temperature, pressure, and time for which the fluid is uncirculated) information. The model developed was referred to as the wellbore sag model (WSM). The authors reported that, generally, the WSM captured the appropriate characteristics of the fluids and successfully predicted their respective sag behaviour.

The objective of this paper is to present an entirely independent numerical simulation study on barite sag in pipe and annular sections by employing CFD. Additionally, a CFD-DEM approach, based on preliminary results, is proposed for future research on the study of barite sag in annular sections.

#### 2. Mathematical Approach

The mathematical approach is in two parts. In the first part, the governing equations are built inside a two-dimensional horizontal pipe geometry and the finite element method (FEM) utilized to solve the equation-sets, for studying the solid concentration distribution of a solid-liquid system in pipe flow. In the second part, governing equations are built inside a three-dimensional horizontal annular geometry and the finite volume method (FVM) utilized to solve the equation-sets, for studying the barite sag phenomenon in an annulus under flow conditions.

##### 2.1. Horizontal Pipe Section

Description of the Eulerian-Eulerian approach to study the solid concentration distribution in a pipe in shear flow of a solid-liquid system: This model assumes that solid concentration only changes in the axial and lateral directions. In addition, the model is developed under laminar flow conditions.

###### 2.1.1. Mass Balance

Assuming that the mass transfer between the two phases is zero, the following continuity relations hold for the continuous and dispersed phase [9]: Here (dimensionless) denotes the phase volume fraction, *ρ* (kg/m^{3}) is the density, and (m/s) is the velocity of each phase. The subscripts and denote quantities relating to the continuous and dispersed phase, respectively. The following relation between the volume fractions must hold

Both phases are considered incompressible, in which case (1) can be simplified as

When (3) and (4) are added together, a continuity equation for the mixture is obtained:

In order to control the mass balance of the two phases, the Euler-Euler model interface solves (4) together with (5). Equation (4) is used to compute the volume fraction of the dispersed phase, and (5) is used to compute the mixture pressure.

###### 2.1.2. Momentum Balance

The momentum equations for the continuous and dispersed phases, using the nonconservative forms of Ishii [10], areHere (Pa) is the mixture pressure, which is assumed to be equal for the two phases. The viscous stress tensor for each phase is denoted by (Pa), (m/s^{2}) is the vector of gravitational acceleration, (N/m^{3}) is the interphase momentum transfer term (that is, the volume force exerted on each phase by the other phase), and (N/m^{3}) is any other volume force term.

For fluid-solid mixtures, (7) is modified in the manner of Enwald et al. [11]

The fluid phases in the above equations are assumed to be Newtonian and the viscous stress tensors are defined aswhere *µ* (Pa·s) is the dynamic viscosity of the respective phase.

In order to avoid singular solutions when the volume fractions tend to zero, the governing equations above are divided by the corresponding volume fraction. The implemented momentum equations for the continuous and dispersed phase are as given in (10) and (11), respectively.

###### 2.1.3. Dispersed Phase Viscosity

Using an expression for the mixture viscosity, the default values for the dynamic viscosities of the two interpenetrating phases are taken as A simpler mixture viscosity covering the entire range of particle concentrations is the Krieger type [12]:where is the maximum packing limit, by default 0.62 for solid particles. Equation (13) can be applied when .

On the other hand, user-defined expressions for the dispersed phase and mixture viscosity can be employed [4]:where (SI unit: Pa·s) is viscosity and the subscripts , and represent continuous phase, dispersed phase, and mixture, respectively; is the dispersed phase concentration.

###### 2.1.4. Interphase Momentum Transfer

The drag force added to the momentum equation is defined as where is the drag force on the continuous phase, is the drag force on the dispersed phase, is a drag force coefficient, and the slip velocity is defined as

###### 2.1.5. Dilute Flows

For dilute flows the drag force coefficient can be modelled asIn this case drag coefficient can be computed from the Schiller-Naumann model in (18), where the particle Reynolds number is given as in (19):

###### 2.1.6. Solid Pressure

For fluid-solid mixtures, a model for the solid pressure, , in (11), is needed. The solid pressure models the particle interaction due to collisions and friction between the particles. The solid pressure model implemented uses a gradient diffusion based assumption:where the empirical function (21) is given by () as described in Enwald et al. [11]

##### 2.2. Horizontal Annular Section

Description of the 3D model employed to study the behaviour of barite particles in an annulus based on the CFD approach is as follows.

###### 2.2.1. Governing Equations for Fluid Flow

The local averaged Navier-Stokes equations describe the 3D equations of motion of viscous, unsteady, and incompressible fluid phase.

*Mass Conservation Equations.* The mass conservation equation is expressed aswhere is the fluid velocity, is the fluid density, is the volume fraction of the fluid phase, and is the time.

*Momentum Conservation Equations.* The momentum conservation equation is expressed aswhere is the fluid pressure, is the viscous stress tensor, and is the volume-averaged (on a cell) interaction forces (interphase momentum transfer source term between the particles and fluid). The stress tensor is defined as [13] where is the solvent viscosity which is dependent on the shear rate, is the shear rate, and is the rate of deformation tensor (25a). In (25b), is an identity matrix and the equation is valid for laminar flow. The shear rate is calculated from the second invariant of the rate of deformation tensor this is defined in (27).And , where is the velocity field.

*Rheology Model.* The relationship between the shear stress and shear rate in the fluid phase is described by the non-Newtonian generalized power law model. The mathematical representation is shown in [13]where is the temperature shift factor (for an isothermal flow, = 1), is the yield stress threshold, is the yielding viscosity, is the minimum viscosity limit, is the maximum viscosity limit, is the power law exponent, and is the consistency factor. The value of determines the class of the fluid; that is, for , the fluid is Newtonian, , the fluid is shear-thickening (dilatant), and < 1, the fluid is shear-thinning (pseudoplastic) or viscoelastic.

###### 2.2.2. Governing Equations for Particle Flow

*(I) Dispersed Phase Modelled as Lagrangian Phase*

*Basic Equations of Motion.* The most basic particle description involves only its position and velocity . These two quantities relate through the equation of motion [13]: The grid velocity is evaluated at the particle position ; its appearance in (29) indicates that the convection is that is the absolute velocity of the particle, whereas is the position of the particle with respect to the frame of reference.

For parcels, individual particles are not tracked; instead a single parcel represents a set of identical particles, at some mean centroid . The velocity of the parcel is assumed to be the same as its constituent particles; hence its equation of motion is

*Mass Balance for a Material Particle.* The equation of conservation of mass of a material particle iswhere is the mass of the particle and the rate of mass transfer to the particle.

*Mass Transfer.* The rate of mass transfer to a single particle from the continuous phase is .

*Momentum Balance for a Material Particle*. The generic form of the equation of conservation of momentum of a material particle iswhere represents the forces acting on the surface of the particle, and the body forces. These forces are in turn decomposed intowhere is the drag force, is the pressure gradient force, is the gravity force, and is the user-defined body force.

*(a) Drag Force*. The equation for drag force is [14]where is the drag coefficient, is the density of the fluid (continuous phase), is the projected area of the particle, and is the particle slip velocity (and is given as .

*(b) Drag Coefficient*. The Schiller-Naumann correlation [14] is suitable for spherical solid particles (and liquid droplets and small-diameter bubbles). For a viscous continuous phase, the correlation is as defined in (18).

*(c) Pressure Gradient Force*. The equation for the pressure gradient force is [14]where is the volume of the particle, and is the gradient of the static pressure in the continuous phase.

*(d) User-Defined Body Force. *The equation for the user-defined body force iswhere is the user body force (per unit volume).

*(e) Particle Shear Lift Force. *This applies to a particle moving relative to a fluid where there is a velocity gradient in the fluid orthogonal to the relative motion. Saffman [12] gives the lift force aswhere the direction is the direction of the velocity gradient. The 3D version of (38) iswhere is the curl of the fluid velocity and is the shear lift coefficient. can be defined using either of two published definitions. Saffman [12] provides a definition that recovers the original Saffman asymptotic solution for low Reynolds numbers (39) and on the other hand, Sommerfeld [14] provides a definition for a broader range of Reynolds numbers (41).in which and Reynolds number for shear flow is .

*Momentum Transfer. *The rate of momentum transfer to a single particle from the continuous phase is , where is as defined in (33a).

*Boundary Interface Mode. *It is important to formulate the Lagrangian phase boundary interaction mode [13].

Rebounding particles remain active in the simulation; the mode is distinguished by its treatment of the particle velocity. The rebound velocity relative to the wall velocity is determined by the impingement velocity and user-defined restitution coefficients:The superscripts and denote rebound and impingement, respectively; the subscripts and denote wall-normal and tangential, respectively. Since the left hand side of (42) can be split into orthogonal and components, it can be split into two equationswhich serve to emphasize the definition of the restitution coefficients as the constants of proportionality between impingement and rebound velocities. Both coefficients may range from 0 to 1; the latter is “perfect” elastic rebound. The tangential velocity of a wall boundary is zero unless a value is explicitly prescribed through a wall sliding option. In other words, it is nonzero only at no-slip walls.

*(II) Dispersed Phase Modelled as DEM Particles. *The detailed description is shown in the appendix.

#### 3. Model Configuration

##### 3.1. 2D Model: Horizontal Pipe Section

For the physical model of a horizontal pipe section with ID 0.0508 m (2 in.) and length 3.65 m (12 ft) see Figure 1. For the discretized representation of the computational domain which the physics solvers use to provide a numerical solution, the mesh is shown in Figure 2. Table 1 lists the inputs used in the 2D configuration. Additionally, the assumptions and boundary conditions are listed below.