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 everrising operating costs. In as much as deviatedwell 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 twodimensional horizontal pipe geometry and the finite element method utilized to solve the equationsets. In the second part, governing equations were built inside a threedimensional horizontal annular geometry and the finite volume method utilized to solve the equationsets. 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.
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 everrising operating costs. Deviatedwell 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 ratesofpenetration, these wells are being drilled increasingly with invertemulsion 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 bottomsup; 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 [
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. [
Nguyen et al. [
Hashemian et al. [
Kulkarni et al. [
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 CFDDEM approach, based on preliminary results, is proposed for future research on the study of barite sag in annular sections.
The mathematical approach is in two parts. In the first part, the governing equations are built inside a twodimensional horizontal pipe geometry and the finite element method (FEM) utilized to solve the equationsets, for studying the solid concentration distribution of a solidliquid system in pipe flow. In the second part, governing equations are built inside a threedimensional horizontal annular geometry and the finite volume method (FVM) utilized to solve the equationsets, for studying the barite sag phenomenon in an annulus under flow conditions.
Description of the EulerianEulerian approach to study the solid concentration distribution in a pipe in shear flow of a solidliquid 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.
Assuming that the mass transfer between the two phases is zero, the following continuity relations hold for the continuous and dispersed phase [
Both phases are considered incompressible, in which case (
When (
In order to control the mass balance of the two phases, the EulerEuler model interface solves (
The momentum equations for the continuous and dispersed phases, using the nonconservative forms of Ishii [
For fluidsolid mixtures, (
The fluid phases in the above equations are assumed to be Newtonian and the viscous stress tensors are defined as
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 (
Using an expression for the mixture viscosity, the default values for the dynamic viscosities of the two interpenetrating phases are taken as
On the other hand, userdefined expressions for the dispersed phase and mixture viscosity can be employed [
The drag force added to the momentum equation is defined as
For dilute flows the drag force coefficient
For fluidsolid mixtures, a model for the solid pressure,
Description of the 3D model employed to study the behaviour of barite particles in an annulus based on the CFD approach is as follows.
The local averaged NavierStokes equations describe the 3D equations of motion of viscous, unsteady, and incompressible fluid phase.
For parcels, individual particles are not tracked; instead a single parcel represents a set of identical particles, at some mean centroid
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 userdefined restitution coefficients:
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
Input data for the simulation – 2D CFD model.
Parameter  Variable  Value  Units 

Diameter of pipe 

0.0508 ( 
m (in.) 
Length of pipe 

3.65 ( 
m (ft.) 
Fluid inlet velocity 

0.1556 (30.64)  m/s (ft/min) 
Viscosity of liquid 

0.062 ( 
Pa·s (cP) 
Density of liquid 

898.78 ( 
kg/m^{3} ( 
Density of solid 

4198.92 ( 
kg/m^{3} ( 
Diameter of solid 

0.000025 ( 
m ( 
Initial concentration of the solid 

0.067  
Deviation angle 

90  ° (from vertical) 
Number of CFD elements  63,790  
CFD timestep 

0.1  s 
Physical time simulated 

36  s 
Horizontal pipe section,
Mesh for the horizontal pipe geometry.
Every single phase of the mixture behaves as if it was a lone phase, with the exception of the case when interacting with the other phase.
The equations of motion describing the mixture are similarly those for a single phase and are the consequence of the summation motion equations for the individual phases over all phases.
Flow of the liquid phase is in one direction (that is, axial), whereas that of the solid phase is in two directions (both axial and radial).
The liquid and solid densities are constant; in other words, phases are incompressible.
There is no consideration of Brownian motion.
There is no particle interaction arising from collisions and friction between the particles. Thus, solid pressure is neglected.
Wall effects are neglected.
Isothermal and laminar flow are considered.
Noslip condition at the walls.
Inlet boundary conditions, that is, velocity inlet.
Outlet boundary conditions, that is, pressure outlet.
The physical model is an annular test section as depicted in Figure
Input for the simulation – 3D CFD model; dispersed phase modelled as Lagrangian phase.
Parameter  Variable  Value  Units 

Drill string length 

1.8288 ( 
m (ft.) 
Pipe diameter 

0.0508 ( 
m (inch) 
Hole diameter 

0.1016 ( 
m (inch) 
Angle of inclination 

0–90  ° (degrees) 
Fluid inlet velocity 

0.1524 ( 
m/s (ft./min) 
Dynamic viscosity 


Density of liquid 

958.61 ( 
kg/m^{3} (ppg) 
Fluid behaviour index 

0.44  — 
Consistency index 

0.63 (1.316)  Pa·s 
Drill pipe rotation speed 

0–100  rpm 
Eccentricity ratio 

0, 1  — 
Shape of particles  spherical  
Particle diameter 

0.0000249 ( 
m (microns) 
Particle density (dry density) 

4193.92 (35.0)  kg/m^{3} (ppg) 
Coefficient of restitution 

1.0  — 
Number of CFD cells  84,102  
CFD timestep 

0.01  s 
Physical time simulated 

36  s 
Physical model for the annular section: (a) concentric geometry, (b) eccentric geometry.
Mesh for the annular section: (a) concentric section, (b) eccentric section.
Flow of the liquid phase is in one direction (that is, axial), whereas that of the solid phase is in two directions (both axial and radial).
The phases are incompressible; that is, the densities of the solid and liquid are constant.
Interaction forces such as shear lift force, drag force, and pressure gradient force exist between the liquid and solid phase.
Solid particles are spherical.
Solid particles are considered as Lagrangian phases.
The liquid phase is a generalized power law (nonNewtonian) fluid.
Isothermal and laminar flow are considered.
Noslip condition at the walls.
Inlet boundary conditions, that is, velocity inlet.
Outlet boundary conditions, that is, pressure outlet.
The 2DCFD model is compared with the mathematical model of Nguyen et al. [
Comparison between the CFD model and the mathematical model of Nguyen et al. [
Comparison between the CFD model and the mathematical model of Nguyen et al. [
Initially, when
Plot of concentration of solid profile in axial direction,
Plot of concentration of solid profile in lateral direction,
For low annular velocities (i.e., 0.1556 m/s), there is rapid formation of a barite bed at the bottom (lower) section of the pipe. The concentration of solid in this bed is not much greater than the concentration of solid in the fluid, which flows in the top (upper) side of the pipe. Put differently, the bed is actually the fluid with a greater concentration of solid and can be easily removed. This layer is what is referred to as the fluidized bed (see Figure
Solid concentration distribution along the pipe in the axial direction, prediction of the barite bed characteristics.
In brief, there exist three layers during the sedimentation of barite particles in the pipe: the clarified fluid layer (i.e., the fluid that flows upward), the fluidized bed, and the solidified bed layer. The dispersed phase volume fraction contour plots are shown in Figure
2D CFD contour plots for the dispersed phase volume fraction. (a) Fluidsolid mixture nonsegregated. (b) Fluidsolid mixture segregated; developing solidified bed and suspension layers can be visualized. The colored box represents magnitude of the volume fraction of the dispersed phase.
2D CFD contour plot for the dispersed phase volume fraction. Notice the development of the clarified fluid region, suspension layer, and the solidified bed at the bottom of the pipe. The colored box represents magnitude of the volume fraction of the dispersed phase.
The relationship between continuous phase, dispersed phase, and mixture viscosity follows the correlations presented in (
Comparison between CFD model and results of Nguyen et al. [
The minimum flow time from inlet to outlet at a velocity of 0.1524 m/s (lowest annular velocity) is 12 s and 2.4 s for a velocity of 0.762 m/s (highest annular velocity). The visualization of barite accumulation clearly illustrates the tendency for the particles to aggregate on the bottom (lower) side of the test section (particularly in the case of an eccentric annulus). This is for the case of a low annular velocity, in this case 0.1524 m/s (30 ft/min) with a stationary drill pipe. The redistribution of barite particles into the fluid stream, at a low annular velocity, is due to rotation of the drill pipe (Figure
3D visualization of the apparent barite accumulation on the lower side of a fully eccentric annulus (horizontal configuration) at
Particles (barite particles) are injected into the fluid stream, at the inlet, at a mass flow rate of 0.055 kg/s. This configuration simulates the case of high particle concentration (i.e.,
Visualization of barite particle distribution at low annular velocity and stationary drill pipe for a horizontal concentric annular section, that is,
Visualization of barite particle distribution at low annular velocity and rotating drill pipe for a horizontal concentric annular section, that is,
Particles (barite particles) are injected into the fluid stream, at the inlet, at a particle flow rate which defines both the period of injection and the injection velocity. Additionally, this configuration simulates the case of low particle concentration (i.e.,
Initiated at the start of circulation: Figure
Visualization of barite behaviour in a horizontal concentric annulus at a high fluid inlet velocity and drill pipe rotation;
Visualization of the particle tracks coloured according to the velocity of the particles in Figures
Particles (barite particles) are injected into the fluid stream, at the inlet, at a particle flow rate which defines both the period of injection and the injection velocity. Additionally, this configuration simulates the case of low particle concentration (i.e.,
Initiated at the start of circulation: Figure
Visualization of barite behaviour in a horizontal eccentric annulus at a high fluid inlet velocity and drill pipe rotation;
In as much as significant efforts have been made to address the key critical issues in numerical simulation of barite sag, some simplifications have also been made:
The 2D model configuration for the horizontal pipe does not account for pipe rotation. Additionally, the influence of mean velocity on dynamic barite sag is not accounted for.
The 3D model configuration for the horizontal annulus does not account for different shapes and sizes of particles as it assumes uniform spherical particles. Different particle shapes and sizes can be introduced by employing different particleshape models in the simulation and then observing the sag tendency.
In the 3D model, only a qualitative analysis of the results is considered. At low particle concentrations, the quantitative results are comparable to those available in published literature and at very high particle concentrations; the model suffers from a convergence problem and thus does not produce reasonable results. Therefore, there is a need to perform a convergence improvement study so as to improve on the accuracy of the results and aid in the quantitative analysis of the results.
In both models, no inclination angle other than 90° is considered. Inclination angles between 30° and 90° can be introduced in the simulation and the resulting sag tendency observed and analysed.
In the 3D model, based on the CFDDEM approach, only a theoretical background for the governing equations (see the appendix) is provided and one stage of simulation successfully performed (see Figure
Flow stream populated with barite particles: case of a concentric annular section. (a) Particle velocity scene: red band shows particles with highest velocity. (b) Residence time scene: red band shows highest travel time; particles that have reached an extreme end of a section.
Illustration of an integrated CFDDEM solver implementation scheme.
Illustration of the CFDDEM coupled implementation scheme.
Nevertheless, the 3D model developed above, to the authors’ knowledge, is the first attempt to perform an independent numerical simulation to investigate dynamic barite sag in an annular section. Furthermore, the model based on the CFDDEM approach is the first attempt to propose the use of the CFDDEM method for the investigation of barite sag behaviour, keeping in mind that the system is a liquidsolid flow. Past researchers [
A numerical simulation approach has been undertaken to investigate the different aspects of barite sag behaviour under influencing factors, such as annular velocity, drill pipe rotation speed, and eccentric drill pipe, as well as the rheology of drilling fluid, that is, Newtonian and nonNewtonian fluid. For the Newtonian fluid case, governing equations were built inside a 2D horizontal pipe geometry and the finite element method (FEM) utilized to solve the equationsets whereas for the nonNewtonian fluid case, the governing equations were built inside a 3D horizontal annular geometry and the finite volume method (FVM) utilized to solve the equationsets. Furthermore, it is worth noting that the drill pipe motion is modelled as a grid flux in the convective term, instead of a body force due to system rotation in the momentum equations.
The 2DCFD model shows that the concentration of solid increases with time at the bottom (lower) section of the pipe. For a Newtonian fluid, which has viscosity of 0.062 Pa·s, the CFD model results indicate that, at low fluid inlet velocity (i.e., 0.1556 m/s), there is rapid formation of a barite bed at the bottom (lower) section of the pipe. Additionally, the model results show that there exist three layers during the sedimentation of barite particles in the pipe: the clarified fluid layer (i.e., the fluid that flows upward), the fluidized bed, and the solidified bed layer. There exists a critical solid concentration below which mixture viscosity is independent of solid concentration and beyond which the converse is true.
The developed 3DCFD model outputs positive results in contrast to some experimentally reported observations in the study of the dynamic barite sag phenomenon:
For the case of low annular velocity (i.e., 0.1524 m/s), drill pipe rotation serves to disturb the apparent barite bed and redistribute the barite particles back into the flow stream.
The effect of drill pipe rotation has a pronounced effect for the eccentric annulus compared to the concentric scenario.
The combined effect of high annular velocity (i.e., 0.762 m/s) and drill pipe rotation (i.e., 50 rpm) results in a tremendous reduction in the barite sag occurrence. Still, the effect is more pronounced for the eccentric annulus than the concentric scenario.
Maintaining the drilling fluid circulation at a high annular velocity and with rotating drill pipe ensures that no barite sag occurs.
The DEM modelling introduces extra body force representing interparticle interaction due to particle contacts with other particles and with mesh boundaries. Thus, (
Putting together the above considerations, (
Note that the momentum transfer to the particle from the continuous phase is simply
Besides the standard Lagrangian linear momentum equation, the DEM particle equations of motion incorporate angular momentum conservation equations:
Contact forces acting on
The tangential component of the contact force,
The tangential damping force,
Accordingly, the tangential torques acting on
Note that, for particlewall collisions, the formulas stay the same, but the wall radius and mass are assumed to be
Drag torque reduces the difference in the rotational differences between a particle and the fluid in which it is immersed. The drag is a torque applied to a DEM particle [
The rotational drag coefficient,
Note that (
Computational fluid dynamics
Discrete element method
Finite difference method
Finite element method
Finite volume method
2dimensional
3dimensional
Reynolds number
Temperature shift factor
Eccentricity ratio
Drag coefficient
Length of pipe/drill string length
Diameter
Diameter of pipe
Diameter of hole
Pressure
Density
Viscosity
Volume fraction
Shear rate
Power law exponent
Consistency factor
Deviation/inclination angle
Coefficient of restitution
Angular velocity.
Normal
Tangential
Solid
Particle
Liquid
Fluid.
The authors declare that they have no conflicts of interest.
This work is financially supported by the Fundamental Research Funds for the Central Universities (Grant no. 16CX05019A) and the National Natural Science Foundation of China (Grant no. 51104171).