# Numerical Evaluation of CFD-DEM Coupling Applied to Lost Circulation Control: Effects of Particle and Flow Inertia

Authors: Marcos V. Barbosa; Fernando C. De Lai; Silvio L. M. Junqueira
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

In the present study, the transport and deposition of solid particles to mitigate the loss circulation of fluid through a fracture transversely placed to a vertical channel is numerically investigated. These solid particles (commonly known in the industry as lost circulation materials—LCMs) are injected into the flow during the drilling operation in the petroleum industry, in hopes to control the fluid loss. The numerical simulation of the process follows a two-stage process: the first characterizes the lost circulation flow and the second the particle injection. The numerical model comprises an Eulerian–Lagrangian approach, in which the dense discrete phase model (DDPM) is combined with the discrete element method (DEM). A parametric analysis is done by varying the vertical channel Reynolds number, the particle-to-fluid density ratio, and the particle diameter. Results are shown in terms of the particle’s bed geometric characteristics, focusing on the location inside the fracture where the particles deposit, and the particle bed length, height, and time spent to fill the fracture. Also monitored are the fluid loss reduction over time and the fractured channel bottom pressure (which can be related to the fracture pressure). Results indicate that using a slow/intermediate flow velocity, associated with heavy particles with small diameters, provides the best combination for the efficient mitigation of the fluid loss process.


---

## Body

## 1. Introduction

The fundamentals of flow of fluid and solid particles, or fluid-solid two-phase flow, are of great interest because of several important engineering applications, such as sediment transport, hydrocyclones, fluidized beds, filtration processes, sedimentation, cuttings removal, and wellbore stability [1].Singular among the existing practical problems of two-phase fluid-solid flows is the loss of drilling fluid through the porous formation during the drilling process of oil exploration. The loss occurs by the leakage of drilling fluid through pores or discontinuities (such as fractures) [2] of the rock formation in the wellbore as drilling progresses. Lost circulation events have always been a significant cause of nonproductive drilling time, increasing well construction costs. Overall, such events could be classified as “light” loss when loss rates are below 50 bbl/h, moderate for losses that range between 50 and 150 bbl/h, and severe for losses over 150 bb/h. The industry spends millions of dollars a year to combat lost circulation and its detrimental effects, such as stuck pipe and kicks [3, 4]. Commonly related to formations with high permeability (caves and vugs), the loss of circulation is intensified in the presence of fractured zones, whose effects are notably more dramatic when the fractures are connected to a network of natural fissures [5].Historically, materials to combat the lost circulation, so-called LCM (lost circulation material), have been used by the oil industry to minimize drilling fluid losses. Such materials are typically classified as flakes, fibers, and granulars (calcium carbonate, mica, etc.). Recent strategies list a range of applications such as graphite, cross-linked pills, high fluid loss compression, or a blending of these materials in different concentrations [6, 7]. Because of the high cost of drilling fluids, oil companies resort to a “sealing flow” technique by adding particles of selected granulometry to the drilling fluid itself, in hopes that the particles will settle along the fracture and eventually seal the leaking path, stopping the loss of drilling fluid [4, 8–10].The study of this sealing flow is complex and frequently difficult to perform experimentally, particularly due to limitations of available equipment (harsh flow conditions) and difficulty in performing flow visualization (opacity of the medium). Nevertheless, the alternative of performing numerical modeling and simulations of sealing flow, issue to be addressed in the present work, is not without challenges of its own. The main limitation of a numerical approach has been the large computational power necessary to resolve the flow and interactions of all the individual particles flowing along with the fluid. Notwithstanding, numerical simulations of sealing flow have recently been brought to the forefront (i.e., they became practical) due to the recent advances in the available hardware and software.When modeling fluid-particle flow in sealing problems, the model needs to be able to address not only fluid motion itself, but also the particle-fluid and the particle-particle interactions in an efficient fashion [11]. In this respect, the combination of computational fluid dynamics (CFD) and discrete element method (DEM) is recognized as one of the most powerful combinations for the numerical study of fluid-particle interactions, especially when a four-way coupling is paramount to obtaining simulation results that are realistic and accurate.The DEM accounts for particle motion, including collisions among the particles, by solving Newton’s second law of motion for each individual particle [12]. The method considers a soft-sphere approach in which the collisions among several particles can be easily addressed in the same time step, making it a preferable choice for dense flows whether the particles are spherical or not [13].The capability and use of the CFD-DEM combination have been reported recently in the literature for a wide range of applications, including fluidized beds [14–19], filtration processes [20–24], hole cleaning and sediment transport in the oil and gas industry [25–28], hydrocyclones, vortex flow, and instabilities [29–31], and bed-load transport [32, 33]. Specifics of the CFD-DEM coupling have also been documented, including studies of different types of DEM coupling [34], numerical errors involved in the models [35], requirements for the coupling with distinct fluid mechanics numerical models [36–38], and the use of nonspherical particles [39].In the present work, a CFD-DEM combination is used to simulate the transport of particles along a simplified vertical and fractured wellbore. The objective is to learn how the particles behave along the fracture and how they deposit and eventually seal a leaking path of the drilling fluid, or at least reduce its loss. Of great interest is the determination via numerical simulations of the amount of fluid being lost through the fracture considering different flow configurations, characterized by the distinct flow Reynolds number, particle-to-fluid density ratio, and particle diameter. The outcome of such parametric exercise is a prediction of when the addition of the particles will yield the sealing of the leaking path, making the sealing flow of practical interest.


## 2. Problem Characterization

The drilling operation consists, in a very simplified view, of the partial removal of a substrate material in order to reach an oil and gas reservoir. The substrate is removed by the action of a drilling column which carries a drill bit at the end. Due to the friction effect against the substrate, the drill bit tends to overheat when in operation; to keep it within an acceptable working temperature range, the bit is cooled by the action of a drilling fluid, injected at the top of the drilling column. This fluid flows down the drilling column, and it returns to the surface through the annular space comprised between the column and the formation.During drilling, the operator must observe, among other factors, the so-called operational window and the equivalent circulating density. The operational window is a combination of two pressures inside the wellbore: the pore pressure and the formation pressure. The pore pressure is defined as the pressure of fluids within the pores of a reservoir; the formation pressure is the pressure at which the formation is mechanically fractured. The equivalent circulating density (ECD) is the effective density of the circulating (drilling) fluid, taking into account the pressure drop in the annulus [40].Observing the operational window and the ECD, the drilling process may be classified as underbalanced or overbalanced (see Figure1). Underbalanced drilling occurs when ECD is lower than the pore pressure. Conversely, an operation is overbalanced when the ECD is higher than the pore pressure. Therefore, the lost circulation tends to become more prevalent during overbalanced drilling because of the tendency of the drilling process in this case to cause fractures in the formation. It is important to notice that the lost circulation phenomenon may be intensified in the presence of fractures already present in the formation (natural fractures) in addition to those caused by the drilling process (when ECD is higher than the formation pressure) [41, 42].Figure 1
Representation of the operational window, its upper and lower limits, and ECD.When facing lost circulation in overbalanced drilling, the operator may attempt to mitigate the losses by adding particles to the circulating fluid, as mentioned previously. Remarkably, such process is feasible during operation, especially in the cases of seepage or partial losses. This is so because the particles are transported to the fluid loss region by the drilling fluid itself, depositing, accumulating, and clogging the flow path.In the present study, the fluid loss region is represented by a discrete fracture, which mimics a naturally occurring fracture or even induced by the drilling process [43, 44]. Figure 2 shows a 2D representation of the wellbore formation with a horizontal fracture of a typical vertical drilling operation. The figure shows a longitudinal (top) and a transversal cutting plane (bottom), the fracture aperture (eFR) and its depth (zFR), and also the fluid flow direction in each region of the wellbore-fracture-formation.Figure 2
Representation of a wellbore with discrete fracture cut planes: longitudinal (top) and transversal (bottom).The geometric representation of the idealized fracture is presented in Figure3(a). Notice that zFR is significantly larger than eFR, providing a symmetric inlet at the fracture region. In Figure 3(b), the red region marked as (1) is the inlet section of the fluid; (2) represents the end of the fracture; (3) symbolizes the wellbore top region; and (4) indicates the injection surface from which the particles are released. The annulus is represented by half of its hydraulic diameter (hCH) which is the space between the drilling column and the wellbore wall. In relation to the fracture, whose aperture is eFR, the figure also shows the upstream and the downstream lengths, lUP and lDW, respectively.Figure 3
3D domain representation and scale of the fracture plane related to the wellbore (a); fractured wellbore simplified representation (b).
(a)
(b)Another simplification adopted in this study is to consider the formation (including the fracture surfaces) impermeable. Therefore, the lost circulation analysis is limited to the fracture flow region. Consequently, the fracture filling process is analyzed without the percolation effect that typically occurs in porous reservoirs. As the porosity can be very small in these regions, the simplification effect is considered minor in respect to the flow and deposition of the particles.To characterize the fluid loss, a two-step procedure is proposed by changing the boundary conditions applied in the problem whose goal is to observe a predetermined flow condition inside the wellbore. The flow is represented by a Reynolds number (equation (1)) defined in terms of the velocity at the wellbore inlet, Uβ,CH,i, set at surface (1) in Figure 3(b).(1)Uβ,CH,i=ReμβρβhCH,where the subscripts β, CH, and i refer to the fluid phase, the wellbore, and the inlet region, respectively. The circulating fluid density is indicated by ρ and the viscosity by μ.A special procedure is necessary for determining the flow boundary conditions in the domain because the fluid loss in a real wellbore is due to the pressure difference between the annulus and the formation, as explained earlier by the operational window, and these pressures are usually unknown. Hence, a specific amount of fluid loss is first set at the fracture outlet,qloss, at surface (2), while the remainder of the circulating fluid returns to the surface, with qreturn set at surface (3) of Figure 3(b). With all boundary conditions set, the problem is simulated until a fully developed flow condition is obtained at the annulus.When this state is reached, the pressure at the wellbore inletpβinlet and the respective pressures at surfaces (2) and (3) are collected, respectively, ploss and preturn. The second step is implemented by maintaining the velocity inlet boundary condition at surface (1) and changing the boundary conditions of surfaces (2) and (3) to the related pressures collected in the first step. Once the flow, considering the fluid loss, reaches a steady regime, the particles are released from the injection surface (4) at the same flow velocity as that of the circulating fluid to minimize any effect they might otherwise have over the fluid loss velocity field.


## 3. Numerical Model

The analysis of the particle transport in a fluid can be performed according to two approaches. The first one, the so-called Eulerian approach, consists of treating the particles as a continuum solid phase. This method is usually preferred when the tracking of the particles can be neglected, meaning that the particle concentration field is enough to describe the particle motion. Alternatively, in the Lagrangian approach, the particles are tracked individually throughout the domain, with their trajectories computed using Newton’s second law of motion [45, 46].In the study presented hereby, the fluid and the particles are described according to, respectively, the Eulerian and Lagrangian approaches via the dense discrete phase model (DDPM) of Ansys Fluent® software. The basics of the numerical model are presented here, with a more thorough and detailed analysis found in [47]. For the continuous phase (fluid), the mass and momentum conservation equations are described by equations (2) and (3):(2)∂ρβρβ∂t+∇⋅εβρβuβ=0,(3)∂εβρβuβ∂t+∇⋅εβpβuβuβ=−εβ∇ρβ+∇⋅εβμβ∇⋅uβ+εβρβg+FDPM+SDPM,where t is the time, εβ represents the continuous phase volume fraction, uβ is the velocity vector, pβ is the pressure gradient, g refers to the gravity acceleration vector, FDPM refers to the coupling force between the phases, and SDPM is the source term due to the displacement of fluid caused by the particle entering a control volume.Particle movement is determined using both the definition of velocity and Newton’s second law of motion, given by equations (4) and (5), respectively,(4)dxpdt=up,(5)mpdupdt=Fd+Fgb+Fpg+Fvm+Fls+FDEM,where xp and up represent the particle position and velocity, respectively, and the particle mass is given by mp. By applying a Lagrangian approach to the particles, several forces can be included in the analysis of particle motion. A summary of the expressions for determining the aforementioned forces is presented in Table 1, where drag (Fd), gravitational and buoyancy (Fgb), pressure gradient (Fpg), virtual (added) mass (Fvm), and Saffman’s lift (Fls) forces are considered.Table 1
Forces considered in the numerical model.
Force
Equation
Gravitational and buoyancy
F
gb
=
m
p
ρ
p
−
ρ
β
/
ρ
p
g
Drag
F
d
=
3
/
4
m
p
μ
β
/
ρ
p
d
p
2
C
D
Re
p
u
β
−
u
p
Saffman lift
F
ls
=
C
ls
m
p
ρ
β
/
ρ
p
∇
×
u
β
×
u
β
−
u
p
Virtual mass
F
vm
=
C
vm
m
p
ρ
β
/
ρ
p
D
/
D
t
u
β
−
u
p
Pressure gradient
F
pg
=
m
p
ρ
β
/
ρ
p
u
β
∇
⋅
u
β
Collision
Normal
F
n
=
k
δ
+
γ
u
12
⋅
λ
12
λ
12
Tangential
F
t
=
−
μ
a
F
n
ζ
12Regarding equation (5), a four-way coupling accounts for the particle collisions [48]. As such, FDEM = Fn + Ft, where the last two terms represent, respectively, the normal and the tangential components [49]. In this sense, the particles are able to collide with each other as well as with any solid surface on the numerical domain, be it a fracture or a channel surface. This means that once a particle reaches the boundary of the numerical domain, two situations may occur: if the boundary is a solid surface, the collision force is calculated using FDEM; in case of an open flow region, such as the fracture or channel outlet (surfaces 2 and 3 in Figure 1), the particle will simply exit (be eliminated from) the domain, as it happens in practical situations.While the drag coefficientCD is obtained from Morsi and Alexander’s [50] correlation, Saffman’s constant Cls is calculated using Li and Ahmadi’s [51] model, being the virtual mass coefficient obtained from [52].The tangential component of the collision force is usually obtained from Coulomb’s law in whichμa is the friction coefficient and ζ12 represents the tangential direction. For the proposed problem, the tangential force is neglected.The normal force is evaluated by a spring-dashpot model [53] by means of several numerical parameters. In the spring part of the model, the overlap δ is calculated with two collision partners according to their position x related to their radius δ = |x2 – x1| – (r1 – r2). The particle stiffness constant is a numerical parameter calculated by k=πu122dpρp/3σd2, where σd represents the fraction of the particle diameter that is allowed to overlap with the collision partner and u12 is the relative velocity between the collision partners; λ12 indicates the normal direction of collision.Similarly, the dashpot part modeling employs the damping coefficientγ [54], calculated from the reduced mass m12 = 2m1m2/(m1 + m2), the coefficient of restitution (η), and from the collision time (tcol) as follows: γ = −2m12ln (η)/tcol.For the simulations presented here, an extensive sensitivity analysis was performed to ensure the correct characterization of fluid motion as well as particle transport and collision. The main numerical parameters applied on the model are summarized as follows: fluid time step Δtβ = 2 · 10−2 s, particle time step Δtp = 2 · 10−4 s, and particle’s stiffness constant k = 2.0 N/m. Also, the coefficient of restitution η = 0.9 is applied for both particle-particle and particle-wall collisions. Other numerical details of the DDPM model, such as pressure-velocity coupling, pressure correction, spatial discretization, and transient formulation, are detailed in [47].


## 4. Results and Discussion

### 4.1. Verification of DDPM-DEM Capability

In order to evaluate the capability of the model to accurately predict the interaction between the solid and fluid phases, as well as to correctly determine the collision forces amongst the particles, two verification tests were carried out: the settling velocity of a single particle released in a fluid and the bouncing motion of a single particle against a static solid surface.The settling velocity problem consists of abandoning a single sphere from a resting position inside a container of fluid and measuring its velocity up to a constant value named terminal or settling velocity. The collision problem consists, at first, of the same settling velocity problem, except that this time, the particle will collide with a static wall and bounce on it until the particle achieves a zero-velocity condition. Interacting problems of one, two, or four particles are often applied in the literature to evaluate the capacity of numerical models [55–60].The settling velocity of a single spherical particle, presented in Figure4(a), was simulated and compared with the numerical and experimental results of [61]. The results correspond to a particle density and diameter of ρp = 7710 kg/m³ and Dp = 0.8 mm, respectively, settling in water (ρβ = 998.2 kg/m³ and μβ = 1.003 10−3 Pa·s).Figure 4
Comparison between the data from [61] (MP) (a) and [62] (GLP) (b).
(a)
(b)The results presented in Figure4(a) show a good agreement with both the numerical model and the experimental results presented by Mordant and Pinton [61], even at the final stages of the settling process. For the settling velocity, the relative difference is less than 1%, 0.316 m/s for the experimental results and 0.313 m/s for the numerical simulation. Thus, one can conclude that the DDPM model is indeed able to calculate the interaction forces between the fluid and a particle with the inclusion of the selected forces shown at Table 1.As the particle deposition process occurs within the fracture, particle collisions become more important. Consequently, in order to ascertain the capability of the DDPM-DEM models, outcomes from the bouncing motion of a spherical particle are compared with the experimental results from [62], as depicted in Figure 4(b). In this case, a particle with Dp = 0.8 mm and ρp = 7800 kg/m³ is moving in a fluid composed of a water-glycerin mixture (73.7% water in volume) with properties ρβ = 935 kg/m³ and μβ = 10−2 Pa·s. Prior to the first collision, the particle attains a settling velocity of roughly 0.6 m/s.Observe the DEM model is capable to represent the wall-particle interactions, since after the first collision only a slightly lower value of the maximum velocity is observed. When the second collision happens, a time shift in the numerical simulation of roughly 0.1 s is noticed, which persists for the remaining collisions, causing a delay in the time demanded to cease the particle motion. However, it is important to notice that the speeds at each collision as well as the number of collisions were predicted correctly by the simulations.As the particles are mostly carried in the horizontal direction of the fluid motion that occurs inside the fracture instead of being exclusively moved in the vertical direction, less interest is dedicated to the repetitive collision of a particle against a wall. Therefore, the first collision that occurs when the particle reaches the surface of the fracture demands special attention. From the analysis presented in this section, it is possible to assure that the DDPM-DEM model is indeed a good approach for the problem investigated here.


### 4.2. Filling a Horizontal Fracture

In the simulations to be presented, the fluid loss occurring at the end of the fracture was set asqloss = 10% of the total mass flow rate. As the lost circulation phenomenon is established, the particles are released in the numerical domain from an injection surface as shown in Figure 5, with the maximum velocity attained from the fluid flow velocity profile. As the momentum exchange from particle to fluid along time occurs, the particles tend to assume the velocity of the fluid nearby. As a consequence, at a steady state, a parabolic velocity profile is verified (the color of the particle is associated to its velocity—red means faster particles). At each time step, a total of 30 particles are injected. For instance, considering 175 seconds of injection, for a time step of 2 · 10−2 seconds, 262,500 particles will be injected into the domain.Figure 5
Particle injection process for several time steps.Before choosing the particle time step size, one must be certain that no particle will be released inside each other so that the overlap parameterδ is not incorrectly calculated. This is guaranteed by associating the particle velocity with the fluid time step size in a manner that the distance covered by the particle is sufficiently larger than the particle diameter. By observing this simple procedure, an incorrect particle injection can be avoided.It is also important to understand that not all particles will enter the fracture and contribute to the filling process. In fact, depending on the particles’ properties, fluid, and flow conditions, only 2 to 5% of all injected particles will ever get into the fracture. Also, as the fluid loss is reduced throughout the process, fewer particles will enter the fracture. This situation is presented in Figure6.Figure 6
Fracture inlet region for several instants showing the variation on the amount of particle entering the fracture.The channel geometric characteristic presents an upstream length oflUP = 1.8 m, a downstream length of lDW = 0.225 m, and a width of hCH = 0.045 m. The fracture was considered as impermeable with hFR = 0.720 m in length and width eFR = 0.01 m. The fluid considered in the simulations was a water-glycerin mixture with ρβ = 1187.6 kg/m³ and μβ = 27.973 10−3 Pa·s.The flow Reynolds number, written in terms of fluid’s mean velocity (Uβ,CH,i), equation (1), set at the channel inlet, controls the flow rate. The particle-to-fluid density ratio (ρp/β) and the particle diameter (Dp) are also the main parameters of the present study.A numerical grid test, where the fluid loss stabilization was taken as the main parameter to determine mesh independence, showed that for a relative difference of 1%, no significant changes on the flow conditions and on the particle transport/deposition rate were observed. Additionally, a grid refinement, to ensure a suitable wall effect evaluation, was applied specifically for each configuration. Thus, a channel with 18000, 20000, and 22000 control volumes and respective Reynolds numbers of 250, 500, and 750 were considered. In Figure7, a mesh detail at the fracture inlet region is shown.Figure 7
Mesh display for the fracture inlet region for Re = 500.Initially, the influence of the Reynolds number over the particle bed formation as well as the fluid loss through the fracture is investigated. The simulations were accomplished for three configurations, namely, Re = 250, 500, and 750, corresponding to a channel inlet velocity ofUβ,CH,i = 0.131, 0.262, and 0.393 m/s, respectively. For each flow time step, 1500 particles (ρp/β = 2.25 and Dp = 0.5 mm) are injected into the domain from the particle injection surface. Figure 8 shows the particle bed formed inside the fracture at the moment the particles cease to enter the fracture, Δtfill.Figure 8
Particle bed formed inside the fracture at Δtfill for Reynolds number (Re) variations (ρp/β = 2.50; Dp = 0.5 mm).One can notice a strong influence of the Reynolds number especially in the particle bed position (hbed,i), measured as the distance from the fracture inlet to the beginning of the particle bed. Clearly, the bed position is influenced by the flow quantity of movement, which in turn also impacts the particle bed extension (hbed). As expected, for low Re, the particles start to deposit closer to the fracture entrance, forming a more compact bed, which also leads to an increase in the particle bed height, represented here in terms of the fracture vertical filling percentage (e%,FR), as presented in Table 2.Table 2
Particle bed geometric characteristics in terms of Reynolds number.
Re (–)
h
bed
,
i
mm
h
bed
mm
e
%
,
FR
%
Δ
t
fill
s
250
15
108
64
76
500
64
131
54
82
750
102
196
45
136As related to Re, the flow velocity will be different for all three conditions, even if the initial percentage of fluid loss (qloss = 10%) is the same for all the cases simulated. Thus, a low Re indicates a low fluid velocity in the fracture, and the reverse is also true. In this way, the flow presents more or less momentum to carry particles along the fracture as the Reynolds number is, respectively, higher or lower. Consequently, the smaller the Re is, the closer to the fracture entrance (reduction on hbed,i) and also more compact is the particle bed, implying in a shorter length (hbed) and in a larger height (e%,FR).Given the flow inability to maintain particles suspended for a long period of time, it is possible to see in Table2 that the time required to the particle deposition process (Δtfill) is smaller for lower Reynolds number. Hence, to reduce the fluid loss in an impermeable fracture, the flow rate has to be reduced to provide a faster particle deposition in the fracture.Assuming that the understanding on how the particle injection process acts on the fluid loss through the fracture is paramount for the present study, the influence ofQloss = (qloss,t/qloss) × 100%, defined as the ratio between the fluid lost along the particle injection (qloss,t) and the initial fluid loss (qloss), is investigated for different Re. Therefore, as shown in Figure 9(a), for a fluid loss of 100%, as the particle deposition process takes place, a reduction in the fluid escaping through the fracture is expected.Figure 9
Monitoring for the fluid loss (Qloss) by the fracture (a) and dimensionless pressure (Pinlet) at the bottom of the fractured channel (b) for Reynolds number variations (ρp/β = 2.50; Dp = 0.5 mm).
(a)
(b)Another monitor, defined asPinlet = (pinlet,t/pβ,inlet) × 100%, has been configured at the inlet of the main channel, with the aim to mimic the downhole pressure. In such a monitor, presented in Figure 9(b), pinlet,t and pβ,inlet are, respectively, the pressure at the fracture entrance along the particle injection and the initial pressure.Figure9(a) shows that the fluid loss is reduced to roughly 42% of the original loss for both Re = 500 and 750. However, the time spent in the particle deposition process (which is indicated when the fluid loss reaches a plateau) is much higher for Re = 750, as described previously. By reducing the Reynolds number to 250 and consequently the flow strength, the fluid loss reaches 38%, indicating a reduction of 62% of the original fluid loss.The downside of choosing a lower velocity for the fluid is a dramatic increase in the pressure measured at the channel inlet, as can be observed in Figure9(b). The particles injected at the channel contribute to the increase in pressure due to differences in the particle-to-fluid density ratio. For lower velocities, more particles are expected to accumulate on the main (vertical) channel, which causes a prompt increase of Pinlet for Re = 250. As the particles move more rapidly throughout the domain, the increase in the pressure is less prominent.In the second set of simulations, the influence of the particle-to-fluid density ratio (ρp/β) is investigated. For this, the case Re = 500 is set as the reference flow configuration and all the particles have diameter Dp = 0.5 mm. The main objective is to observe the effect of changes in the particle density over the fracture filling process. Further on, the effect of changes in the particle diameter will be presented as well. The final condition for the particles to be packed inside the fracture is presented in Figure 10. Information on size, extension, and position of each bed is presented in Table 3.Figure 10
Particle bed formed inside the fracture at Δtfill for particle-to-fluid density ratio (ρp/β) variations (Re = 500; Dp = 0.5 mm).Table 3
Particle bed geometric characteristics in terms of the particle-to-fluid density ratio.
ρ
p
/
β
h
bed
,
i
mm
h
bed
mm
e
%
,
FR
%
Δ
t
fill
s
1.50
150
284
43
147
1.75
101
183
50
111
2.00
79
161
52
94
2.25
64
131
54
82
2.50
53
131
57
77Whilst changes in Reynolds number influence the flow strength, modifications in the particle density alter the particle quantity of movement. Consequently, the flow has more difficulty to carry heavier particles. As shown in Figure10, the lighter the particles, the greater the hbed values (see the ρp/β = 1.50 case in Table 3). On the contrary, as the particles get heavier (case ρp/β = 2.50), the fluid does not have enough momentum to carry them away. With that, a tendency of the particles to accumulate near the entrance is observed, and a more compact bed is formed.The observations concerning the influence of the particle mass shown in Figure10 are also related to the particle bed length, which gets larger for lighter particles. For the ρp/β values analyzed, it seems there is a stabilization process regarding the particle bed length because no significant difference was observed between cases ρp/β = 2.25 and ρp/β = 2.50, except for a shift on the particle bed initial position hbed,i, which is related to particle quantity of movement.Regarding the fluid loss reduction, shown in Figure11(a), which is the main objective of this study, a slight difference when using heavier particles is noticed, since the final reduction of loss circulation ranges from 39% (ρp/β = 1.50) to 47% (ρp/β = 2.50). The use of heavier particles provides a significant reduction of the time required to fill the fracture. Nevertheless, they are not as effective as using the lighter particles to reduce the loss of fluid. Notice that the difference on the fluid loss reduction for ρp/β = 1.50 and ρp/β = 2.50 is only of 4%, while the time for the filling process is cut by 50%.Figure 11
Monitoring for the fluid loss (Qloss) by the fracture (a) and dimensionless pressure Pinlet at the fractured channel inlet (b) for particle-to-fluid density ratio variations (Re = 500; Dp = 0.5 mm).
(a)
(b)Concerning the pressure monitor at the fractured channel inlet, presented in Figure11(b), the use of heavier particles has the drawback of increasing the pressure in this region by roughly 60%. Different from the change in fluid velocity, here the increase in pressure between the limit cases, i.e., ρp/β = 1.50 and ρp/β = 2.50, is only about 31%. Such behavior indicates that a technician, for instance, is better off choosing heavier particles for the filling process, since in this case, the pressure does not represent any significant restriction for the operation with respect to the particle density.Finally, an alternative analysis of the particle inertia by considering different particle diameters (Dp = 0.4, 0.5, and 0.6 mm) is presented in Figure 12. As pointed out earlier, the Reynolds number is fixed as 500 while the particle-to-fluid density ratio is kept as ρp/β = 2.50. Again, due to the same flow conditions, it is possible to see that the increase of particle mass due to the increase in its diameter is responsible for changing the fixed bed position along the fracture: as particle gets heavier the fluid flow gets less capable of transporting them.Figure 12
Particle bed formed inside the fracture at Δtfill for particle diameter (Dp) variations (Re = 500; ρp/β = 2.50).Results from Figure12 show that the particle bed length increases with the particle diameter. From the geometric characteristics of each bed, presented in Table 4, one can notice that there is an apparent stability on the fracture filling time when Dp = 0.5 and 0.6 mm.Table 4
Particle bed geometric characteristics as a function of the particle diameterDp.
D
p
mm
h
bed
,
i
mm
h
bed
mm
e
%
,
FR
%
Δ
t
fill
s
0.4
77
111
59
111
0.5
53
131
57
77
0.6
44
169
48
77As occurred in the particle-to-fluid density ratio analysis, it is possible to see from Figure13(a) that particles with less quantity of movement (lighter due to a smaller diameter) have the ability to more effectively reduce the fluid loss. In the case presented here, this is also related to the porosity of the bed. As a matter of fact, a more cohesive and less porous bed is formed when the particles are smaller, which in turn helps to reduce the fluid loss.Figure 13
Monitoring for the fluid loss (Qloss) by the fracture (a) and dimensionless pressure (Pinlet) at the fractured channel inlet (b) for particle diameter variations (Re = 500; ρp/β = 2.50).
(a)
(b)Still, the increase inPinlet in Figure 13(b) is directly linked to the particle mass variation, as the flow characteristics remain the same for the three simulations. Therefore, the higher the particle mass, the higher Pinlet. The relative difference of 32% in Pinlet for the limit cases (Dp = 0.4 and 0.6 mm) may not be as significant for the drilling operation. Otherwise, when a comparison between Dp = 0.5 and 0.6 mm cases is drawn, the former presents not only a better fluid loss reduction (roughly 54%) but also a lower increase in Pinlet and a similar Δtfill, suggesting that it would be wiser to employ particles with Dp = 0.5 mm.


## 4.1. Verification of DDPM-DEM Capability

In order to evaluate the capability of the model to accurately predict the interaction between the solid and fluid phases, as well as to correctly determine the collision forces amongst the particles, two verification tests were carried out: the settling velocity of a single particle released in a fluid and the bouncing motion of a single particle against a static solid surface.The settling velocity problem consists of abandoning a single sphere from a resting position inside a container of fluid and measuring its velocity up to a constant value named terminal or settling velocity. The collision problem consists, at first, of the same settling velocity problem, except that this time, the particle will collide with a static wall and bounce on it until the particle achieves a zero-velocity condition. Interacting problems of one, two, or four particles are often applied in the literature to evaluate the capacity of numerical models [55–60].The settling velocity of a single spherical particle, presented in Figure4(a), was simulated and compared with the numerical and experimental results of [61]. The results correspond to a particle density and diameter of ρp = 7710 kg/m³ and Dp = 0.8 mm, respectively, settling in water (ρβ = 998.2 kg/m³ and μβ = 1.003 10−3 Pa·s).Figure 4
Comparison between the data from [61] (MP) (a) and [62] (GLP) (b).
(a)
(b)The results presented in Figure4(a) show a good agreement with both the numerical model and the experimental results presented by Mordant and Pinton [61], even at the final stages of the settling process. For the settling velocity, the relative difference is less than 1%, 0.316 m/s for the experimental results and 0.313 m/s for the numerical simulation. Thus, one can conclude that the DDPM model is indeed able to calculate the interaction forces between the fluid and a particle with the inclusion of the selected forces shown at Table 1.As the particle deposition process occurs within the fracture, particle collisions become more important. Consequently, in order to ascertain the capability of the DDPM-DEM models, outcomes from the bouncing motion of a spherical particle are compared with the experimental results from [62], as depicted in Figure 4(b). In this case, a particle with Dp = 0.8 mm and ρp = 7800 kg/m³ is moving in a fluid composed of a water-glycerin mixture (73.7% water in volume) with properties ρβ = 935 kg/m³ and μβ = 10−2 Pa·s. Prior to the first collision, the particle attains a settling velocity of roughly 0.6 m/s.Observe the DEM model is capable to represent the wall-particle interactions, since after the first collision only a slightly lower value of the maximum velocity is observed. When the second collision happens, a time shift in the numerical simulation of roughly 0.1 s is noticed, which persists for the remaining collisions, causing a delay in the time demanded to cease the particle motion. However, it is important to notice that the speeds at each collision as well as the number of collisions were predicted correctly by the simulations.As the particles are mostly carried in the horizontal direction of the fluid motion that occurs inside the fracture instead of being exclusively moved in the vertical direction, less interest is dedicated to the repetitive collision of a particle against a wall. Therefore, the first collision that occurs when the particle reaches the surface of the fracture demands special attention. From the analysis presented in this section, it is possible to assure that the DDPM-DEM model is indeed a good approach for the problem investigated here.


## 4.2. Filling a Horizontal Fracture

In the simulations to be presented, the fluid loss occurring at the end of the fracture was set asqloss = 10% of the total mass flow rate. As the lost circulation phenomenon is established, the particles are released in the numerical domain from an injection surface as shown in Figure 5, with the maximum velocity attained from the fluid flow velocity profile. As the momentum exchange from particle to fluid along time occurs, the particles tend to assume the velocity of the fluid nearby. As a consequence, at a steady state, a parabolic velocity profile is verified (the color of the particle is associated to its velocity—red means faster particles). At each time step, a total of 30 particles are injected. For instance, considering 175 seconds of injection, for a time step of 2 · 10−2 seconds, 262,500 particles will be injected into the domain.Figure 5
Particle injection process for several time steps.Before choosing the particle time step size, one must be certain that no particle will be released inside each other so that the overlap parameterδ is not incorrectly calculated. This is guaranteed by associating the particle velocity with the fluid time step size in a manner that the distance covered by the particle is sufficiently larger than the particle diameter. By observing this simple procedure, an incorrect particle injection can be avoided.It is also important to understand that not all particles will enter the fracture and contribute to the filling process. In fact, depending on the particles’ properties, fluid, and flow conditions, only 2 to 5% of all injected particles will ever get into the fracture. Also, as the fluid loss is reduced throughout the process, fewer particles will enter the fracture. This situation is presented in Figure6.Figure 6
Fracture inlet region for several instants showing the variation on the amount of particle entering the fracture.The channel geometric characteristic presents an upstream length oflUP = 1.8 m, a downstream length of lDW = 0.225 m, and a width of hCH = 0.045 m. The fracture was considered as impermeable with hFR = 0.720 m in length and width eFR = 0.01 m. The fluid considered in the simulations was a water-glycerin mixture with ρβ = 1187.6 kg/m³ and μβ = 27.973 10−3 Pa·s.The flow Reynolds number, written in terms of fluid’s mean velocity (Uβ,CH,i), equation (1), set at the channel inlet, controls the flow rate. The particle-to-fluid density ratio (ρp/β) and the particle diameter (Dp) are also the main parameters of the present study.A numerical grid test, where the fluid loss stabilization was taken as the main parameter to determine mesh independence, showed that for a relative difference of 1%, no significant changes on the flow conditions and on the particle transport/deposition rate were observed. Additionally, a grid refinement, to ensure a suitable wall effect evaluation, was applied specifically for each configuration. Thus, a channel with 18000, 20000, and 22000 control volumes and respective Reynolds numbers of 250, 500, and 750 were considered. In Figure7, a mesh detail at the fracture inlet region is shown.Figure 7
Mesh display for the fracture inlet region for Re = 500.Initially, the influence of the Reynolds number over the particle bed formation as well as the fluid loss through the fracture is investigated. The simulations were accomplished for three configurations, namely, Re = 250, 500, and 750, corresponding to a channel inlet velocity ofUβ,CH,i = 0.131, 0.262, and 0.393 m/s, respectively. For each flow time step, 1500 particles (ρp/β = 2.25 and Dp = 0.5 mm) are injected into the domain from the particle injection surface. Figure 8 shows the particle bed formed inside the fracture at the moment the particles cease to enter the fracture, Δtfill.Figure 8
Particle bed formed inside the fracture at Δtfill for Reynolds number (Re) variations (ρp/β = 2.50; Dp = 0.5 mm).One can notice a strong influence of the Reynolds number especially in the particle bed position (hbed,i), measured as the distance from the fracture inlet to the beginning of the particle bed. Clearly, the bed position is influenced by the flow quantity of movement, which in turn also impacts the particle bed extension (hbed). As expected, for low Re, the particles start to deposit closer to the fracture entrance, forming a more compact bed, which also leads to an increase in the particle bed height, represented here in terms of the fracture vertical filling percentage (e%,FR), as presented in Table 2.Table 2
Particle bed geometric characteristics in terms of Reynolds number.
Re (–)
h
bed
,
i
mm
h
bed
mm
e
%
,
FR
%
Δ
t
fill
s
250
15
108
64
76
500
64
131
54
82
750
102
196
45
136As related to Re, the flow velocity will be different for all three conditions, even if the initial percentage of fluid loss (qloss = 10%) is the same for all the cases simulated. Thus, a low Re indicates a low fluid velocity in the fracture, and the reverse is also true. In this way, the flow presents more or less momentum to carry particles along the fracture as the Reynolds number is, respectively, higher or lower. Consequently, the smaller the Re is, the closer to the fracture entrance (reduction on hbed,i) and also more compact is the particle bed, implying in a shorter length (hbed) and in a larger height (e%,FR).Given the flow inability to maintain particles suspended for a long period of time, it is possible to see in Table2 that the time required to the particle deposition process (Δtfill) is smaller for lower Reynolds number. Hence, to reduce the fluid loss in an impermeable fracture, the flow rate has to be reduced to provide a faster particle deposition in the fracture.Assuming that the understanding on how the particle injection process acts on the fluid loss through the fracture is paramount for the present study, the influence ofQloss = (qloss,t/qloss) × 100%, defined as the ratio between the fluid lost along the particle injection (qloss,t) and the initial fluid loss (qloss), is investigated for different Re. Therefore, as shown in Figure 9(a), for a fluid loss of 100%, as the particle deposition process takes place, a reduction in the fluid escaping through the fracture is expected.Figure 9
Monitoring for the fluid loss (Qloss) by the fracture (a) and dimensionless pressure (Pinlet) at the bottom of the fractured channel (b) for Reynolds number variations (ρp/β = 2.50; Dp = 0.5 mm).
(a)
(b)Another monitor, defined asPinlet = (pinlet,t/pβ,inlet) × 100%, has been configured at the inlet of the main channel, with the aim to mimic the downhole pressure. In such a monitor, presented in Figure 9(b), pinlet,t and pβ,inlet are, respectively, the pressure at the fracture entrance along the particle injection and the initial pressure.Figure9(a) shows that the fluid loss is reduced to roughly 42% of the original loss for both Re = 500 and 750. However, the time spent in the particle deposition process (which is indicated when the fluid loss reaches a plateau) is much higher for Re = 750, as described previously. By reducing the Reynolds number to 250 and consequently the flow strength, the fluid loss reaches 38%, indicating a reduction of 62% of the original fluid loss.The downside of choosing a lower velocity for the fluid is a dramatic increase in the pressure measured at the channel inlet, as can be observed in Figure9(b). The particles injected at the channel contribute to the increase in pressure due to differences in the particle-to-fluid density ratio. For lower velocities, more particles are expected to accumulate on the main (vertical) channel, which causes a prompt increase of Pinlet for Re = 250. As the particles move more rapidly throughout the domain, the increase in the pressure is less prominent.In the second set of simulations, the influence of the particle-to-fluid density ratio (ρp/β) is investigated. For this, the case Re = 500 is set as the reference flow configuration and all the particles have diameter Dp = 0.5 mm. The main objective is to observe the effect of changes in the particle density over the fracture filling process. Further on, the effect of changes in the particle diameter will be presented as well. The final condition for the particles to be packed inside the fracture is presented in Figure 10. Information on size, extension, and position of each bed is presented in Table 3.Figure 10
Particle bed formed inside the fracture at Δtfill for particle-to-fluid density ratio (ρp/β) variations (Re = 500; Dp = 0.5 mm).Table 3
Particle bed geometric characteristics in terms of the particle-to-fluid density ratio.
ρ
p
/
β
h
bed
,
i
mm
h
bed
mm
e
%
,
FR
%
Δ
t
fill
s
1.50
150
284
43
147
1.75
101
183
50
111
2.00
79
161
52
94
2.25
64
131
54
82
2.50
53
131
57
77Whilst changes in Reynolds number influence the flow strength, modifications in the particle density alter the particle quantity of movement. Consequently, the flow has more difficulty to carry heavier particles. As shown in Figure10, the lighter the particles, the greater the hbed values (see the ρp/β = 1.50 case in Table 3). On the contrary, as the particles get heavier (case ρp/β = 2.50), the fluid does not have enough momentum to carry them away. With that, a tendency of the particles to accumulate near the entrance is observed, and a more compact bed is formed.The observations concerning the influence of the particle mass shown in Figure10 are also related to the particle bed length, which gets larger for lighter particles. For the ρp/β values analyzed, it seems there is a stabilization process regarding the particle bed length because no significant difference was observed between cases ρp/β = 2.25 and ρp/β = 2.50, except for a shift on the particle bed initial position hbed,i, which is related to particle quantity of movement.Regarding the fluid loss reduction, shown in Figure11(a), which is the main objective of this study, a slight difference when using heavier particles is noticed, since the final reduction of loss circulation ranges from 39% (ρp/β = 1.50) to 47% (ρp/β = 2.50). The use of heavier particles provides a significant reduction of the time required to fill the fracture. Nevertheless, they are not as effective as using the lighter particles to reduce the loss of fluid. Notice that the difference on the fluid loss reduction for ρp/β = 1.50 and ρp/β = 2.50 is only of 4%, while the time for the filling process is cut by 50%.Figure 11
Monitoring for the fluid loss (Qloss) by the fracture (a) and dimensionless pressure Pinlet at the fractured channel inlet (b) for particle-to-fluid density ratio variations (Re = 500; Dp = 0.5 mm).
(a)
(b)Concerning the pressure monitor at the fractured channel inlet, presented in Figure11(b), the use of heavier particles has the drawback of increasing the pressure in this region by roughly 60%. Different from the change in fluid velocity, here the increase in pressure between the limit cases, i.e., ρp/β = 1.50 and ρp/β = 2.50, is only about 31%. Such behavior indicates that a technician, for instance, is better off choosing heavier particles for the filling process, since in this case, the pressure does not represent any significant restriction for the operation with respect to the particle density.Finally, an alternative analysis of the particle inertia by considering different particle diameters (Dp = 0.4, 0.5, and 0.6 mm) is presented in Figure 12. As pointed out earlier, the Reynolds number is fixed as 500 while the particle-to-fluid density ratio is kept as ρp/β = 2.50. Again, due to the same flow conditions, it is possible to see that the increase of particle mass due to the increase in its diameter is responsible for changing the fixed bed position along the fracture: as particle gets heavier the fluid flow gets less capable of transporting them.Figure 12
Particle bed formed inside the fracture at Δtfill for particle diameter (Dp) variations (Re = 500; ρp/β = 2.50).Results from Figure12 show that the particle bed length increases with the particle diameter. From the geometric characteristics of each bed, presented in Table 4, one can notice that there is an apparent stability on the fracture filling time when Dp = 0.5 and 0.6 mm.Table 4
Particle bed geometric characteristics as a function of the particle diameterDp.
D
p
mm
h
bed
,
i
mm
h
bed
mm
e
%
,
FR
%
Δ
t
fill
s
0.4
77
111
59
111
0.5
53
131
57
77
0.6
44
169
48
77As occurred in the particle-to-fluid density ratio analysis, it is possible to see from Figure13(a) that particles with less quantity of movement (lighter due to a smaller diameter) have the ability to more effectively reduce the fluid loss. In the case presented here, this is also related to the porosity of the bed. As a matter of fact, a more cohesive and less porous bed is formed when the particles are smaller, which in turn helps to reduce the fluid loss.Figure 13
Monitoring for the fluid loss (Qloss) by the fracture (a) and dimensionless pressure (Pinlet) at the fractured channel inlet (b) for particle diameter variations (Re = 500; ρp/β = 2.50).
(a)
(b)Still, the increase inPinlet in Figure 13(b) is directly linked to the particle mass variation, as the flow characteristics remain the same for the three simulations. Therefore, the higher the particle mass, the higher Pinlet. The relative difference of 32% in Pinlet for the limit cases (Dp = 0.4 and 0.6 mm) may not be as significant for the drilling operation. Otherwise, when a comparison between Dp = 0.5 and 0.6 mm cases is drawn, the former presents not only a better fluid loss reduction (roughly 54%) but also a lower increase in Pinlet and a similar Δtfill, suggesting that it would be wiser to employ particles with Dp = 0.5 mm.


## 5. Final Remarks

In this work, the reduction of the lost circulation in a horizontal fracture branching from a vertical channel emulating a fractured oil well was investigated. A numerical modeling accounting for the transport and deposition of solid particles inside the fracture was proposed. Numerical results obtained using the coupling of DDPM-DEM models proved to be a good alternative, not only to represent the two-phase flow but also to characterize the particle bed length, height, and position. This study provided a better understanding of using particles to seal a fracture, especially concerning to the oil and gas industry.Several observations can be pointed out from the results. When injecting a fixed type of particles, an increase in the flow velocity usually increases the time spent on the filling process because the particles tend to be carried by the flow in the fracture region for longer distances. This is reflected in the particle bed formation, whose length is inversely proportional to its height. Increasing the particle-to-fluid density ratio, the particles get heavier; therefore, they tend to accumulate closer to the fracture inlet. The same behavior is verified when the particle diameter is increased: the particles also accumulate closer to the inlet region of the fracture. In both cases, the time spent to fill the fracture is reduced. Clearly, this is related to the weight of the particle since density and diameter are directly linked to particle mass.On the contrary, particle density and diameter affect bed length and height in an opposite manner: while the particle bed length is reduced as the particle density is increased, the reverse is true for the particle diameter. This happens due to the increase in the particle mass at the annulus, providing a higher pressure and resulting in a longer particle bed.Remarkably, the effects presented here are competitive and one can choose to increase particle inertia by increasing its density or diameter. Also, the same effect may be achieved by decreasing flow velocity, for instance. Bearing in mind the operational window, the operator may alter one or more parameters to achieve a specific effect which can be represented by reducing the lost circulation to acceptable levels in the least amount of time or, if possible, curtail all losses without taking time into account.The present study represents a development of a robust and complete methodology that can help Petrobras to optimize its drilling operations in severe loss of circulation scenarios.



---

*Source: 6742371-2019-10-20.xml*