1 Introduction

Accurate and efficient acid placement is a major concern in matrix stimulation of carbonate formations. During acidizing operations, the stimulation fluids will flow preferentially through the paths of least resistance, penetrating the zones with the highest porosity and permeability and then further increasing their conductivity while bypassing the lower-quality zones where stimulation is needed most. Moreover, in the presence of complex natural fractures and vugs in carbonates, the reservoir pressure and permeability may vary tremendously among zones. This makes achieving uniform stimulation and stable production levels a challenging task. To overcome these issues, diversion techniques have been developed with the aim of achieving effective placement of the stimulation fluids, creating a deep and uniform network with highly conductive flow paths (wormholes) between the reservoir and the wellbore, without causing reservoir damage (Asiri et al. 2013; Bukovac et al. 2012; Cohen et al. 2010; Garrouch and Jennings 2017; Ghommem et al. 2015; Kam et al. 2007; Qiu et al. 2018; Siddiqui et al. 2006).

Several alternatives are possible (and commercially available) to effectively place and divert stimulation fluids, and these can be grouped into two main categories: mechanical and chemical diversion (or the combination of both) (Alsaba et al. 2017; Alshubbar et al. 2018; Calcada et al. 2015; Cohen et al. 2010; Davoodi et al. 2018; Droger et al. 2014; Garagash et al. 2019; Gomaa et al. 2011; Kam et al. 2007; Li et al. 2018; Malik et al. 2017; Sau et al. 2015; Solares et al. 2008; Thompson and Gdanski 1993; Wang et al. 2015a, b; Zhang et al. 2018). Mechanical diversion can be deployed for multistage stimulation completions. It is mainly performed by injecting ball sealers made of nylon, hard rubber or biodegradable materials such as collagen as part of a stimulation pumping sequence to plug and shut off the perforations and large fractures that are receiving most of the injected fluids (Alshubbar et al. 2018; Calcada et al. 2015; Cohen et al. 2010; Droger et al. 2014; Li et al. 2017, 2018; Sau et al. 2015; Solares et al. 2008; Wang et al. 2015a, b; Zhang et al. 2018). Alternatively, a straddle packer arrangement can be used to isolate the required treatment interval.

Due to the completion and financial constraints, chemical diverters such as self-diverting acids (SDA) and viscoelastic diverting acid (VDA) (Schlumberger 2018b) are commonly used in carbonate acid stimulation. The use of these fluids has shown great potential for improving acid placement in carbonate reservoirs due to the rheological properties of the fluids that enable self-diversion to zones of lower injectivity (Asiri et al. 2013; Bukovac et al. 2012; Cohen et al. 2010; Garrouch and Jennings 2017). The SDA is a polymer-based in situ gelled acid. It consists of HCl mixed with a gelling agent and a pH-sensitive crosslinker. Initially, the SDA has a low viscosity. As the acid reaction with the formation takes place, the polymer crosslinks at pH = 1, which dramatically increases the viscosity. The resulting highly viscous acid slows down the fluid penetration rate and diverts the following acid stages to other intervals of the reservoir. At pH = 3.5, the gelled acid breaks and reverts to low viscosity. The narrow pH window limits the controllability of the SDA and therefore restricts its efficiency and its broad use. Moreover, polymer-based diverters leave residue after the gels break, potentially compromising the stimulation result by damaging the matrix or the surface of the wormholes, especially if the volume is improperly designed (Garrouch and Jennings 2017).

To overcome the residue issue associated with polymer-based diverters, the viscoelastic surfactant-based diverting acid (VDA) was developed and introduced to the field (Garrouch and Jennings 2017). The viscosity increase in the VDA relies on the byproduct of the acid reaction with the formation rather than the fluid pH. The surfactant-based gel does not contain particles, which makes it nondamaging to the formation. However, if incorrectly applied or tested, the gel may not break in gas reservoirs or may form an emulsion with oil.

In reservoirs that contain highly permeable streaks, such as natural fractures or vugs, the use of viscous gels may be not sufficient to provide efficient diversion, and fibers have been recently used as a diversion material. Degradable polymeric fibers are mixed with VDA to form a new generation of flow diversion agents (Schlumberger 2018a). The acid rapidly viscosifies as it spends in the formation while the fibers form a porous medium that reduces or entirely blocks acid flow into thief zones (with high permeability) and redirects stimulation fluids to lower-permeability zones. Case histories of this new diversion fluid showed that the stimulation of highly heterogeneous carbonate reservoirs improved significantly with the use of this novel acid-based fluid containing degradable fibers (Asiri et al. 2013; Bukovac et al. 2012; Cohen et al. 2010). For instance, as reported in a case study (Schlumberger 2011), despite performing five conventional acid stimulation treatments of a naturally fractured southern Mexico well (with a permeability contrast of more than 300:1) over a period of 3 years, a continuing decline in the production was observed. Furthermore, a production log showed that less than one-quarter of the interval was contributing to hydrocarbon production. Applying the degradable fiber diversion system provided a uniform zonal coverage, and the oil production stabilized at 254 m3/d, that is, 700% increase in the oil production.

As presented above, there are several types of diverting agents for stimulating carbonates relying either on temporary plugging or viscosity control mechanisms (Alsaba et al. 2017; Alshubbar et al. 2018; Asiri et al. 2013; Bukovac et al. 2012; Calcada et al. 2015; Cohen et al. 2010; Davoodi et al. 2018; Detournay 2016; Droger et al. 2014; Garrouch and Jennings 2017; Gomaa et al. 2011; Kam et al. 2007; Li et al. 2018; Malik et al. 2017; Osiptsov 2017; Solares et al. 2008; Sau et al. 2015; Thompson and Gdanski 1993; Wang et al. 2015a, b; Yu 2015; Zhang et al. 2018). Although these methods are often successful when applied during stimulation treatments, the industry still needs more reliable operational guidelines and a deeper understanding of the fundamental aspects of the diversion mechanism. The current empirical design methodology is based on experimental observations and mainly corrected by best practices. For this reason, an acid stimulation job may result in formation damage due to incorrect diverter volume, improper breaking and insufficient clean up energy. In these cases, the maximized post-stimulation benefit that cannot be fulfilled is short lived or even not reached at all. Furthermore, when dealing with a fiber-based product, it is very hard for the operators to determine its merit and implement it with the highest confidence. Therefore, an important step is to develop a qualitative and quantitative understanding of the interactions between the diverting agents and suspending fluids in the formation under different operating conditions. This understanding can be achieved by integrating advanced mathematical modeling and computer simulations with experimental work. This will help to recommend reliable diverting agents and operational guidelines for specific carbonate reservoir geology (lithology and heterogeneity).

The key mechanism behind the diverting capability of fibers is their ability to flocculate and bridge when interacting with suspending fluids. Since the introduction of fiber-based techniques for well intervention and stimulation applications, several laboratory-scale experimental studies (Alsaba et al. 2017; Alshubbar et al. 2018; Calcada et al. 2015; Cohen et al. 2010; Li et al. 2018; Sau et al. 2015; Wang et al. 2015a) have been devoted to examining the fiber bridging phenomenon. The bridging experiments are conducted in pipe flow loops with slots of different shape and size mimicking fractures, as shown in Fig. 1. The bridging ability is tracked by monitoring the pressure signal and observing the fiber plug at the end of the experiment (Cohen et al. 2010; Sau et al. 2015). Li et al. (2018) proposed a novel fiber suspension system. They investigated the impact of the fiber length, fiber concentration, injection rate and permeability on the plugging efficiency. Their experimental study demonstrated the accessibility and good plugging effectiveness of the proposed fiber suspension system, and then they concluded that it can be efficiently applied in high-permeability reservoirs. Wang et al. (2015a) conducted a series of laboratory experiments using a triaxial stress system to investigate the diversion mechanisms of fiber-laden fluid systems for fracturing applications. They found that the dominant factors are mainly the fracture width, fiber concentration, heterogeneity and anisotropy of a reservoir, fluid viscosity and the injection rate. In a recent paper, Yang et al. (2019) performed experimental investigation to study the plugging mechanism of mixed fibers and particulates in a mimicked hydraulic fracture. To enable superior plugging performance, they found that fiber injection into the hydraulic fracture needs to be followed by injected particulates. Furthermore, the ratio between the fibers and particulates to be injected is suggested to be 1:1. Zhang et al. (2019) conducted a series of experiments on a 3D-printed rock model using a large-scale triaxial fracturing system to investigate the fiber-assisted diverting fracturing process. They identified the optimal fiber concentration from an economic perspective, fiber length and injection rate to enable temporary plugging and successful acid diversion.

Fig. 1
figure 1

Flow loop experiment. Left, general schematic. Right, photograph taken from Schlumberger Dhahran Carbonate Research Center (SDCR) laboratory

The present work focuses on the modeling aspect of fiber dynamics in fluid flow. The objective is to investigate the fundamental mechanisms of wellbore bridging by fibers of various characteristics when encountering openings (e.g., natural fractures or wormholes) with different shapes and sizes. A transport model to predict fiber–fluid interactions in carbonates is developed and used to identify the key parameters influencing the bridging phenomenon and then affecting the performance of fluid placement.

In this paper, we consider fluid diversion in a plug-and-perf completion of a horizontal well with multistage acid fracturing. The idea is to inject fibers for temporary plugging previously formed fractures during acid fracturing operations. The goal is to control the uniform formation and growth of fractures over the entire well. To this end, we present a numerical investigation of the bridging performance of fibers in a mimicked fracture. Our study aims at integrating mathematical modeling and computer simulations with laboratory experiments to develop a fundamental understanding of fiber dynamics in fluid flow. Such understanding is intended to provide guidance on controlling efficiently the process of fiber bridging required to achieve accurate fluid placement in the wellbore.

2 Mathematical modeling background

2.1 DEM-based fiber model

Due to the inherent complexity of fiber dynamics in fluid flow, particle-level simulations may be a useful complement to theoretical and experimental investigations. As such, many researchers have proposed fiber suspension models, in which the fibers are represented by multi-rigid-body systems of simple particles, for instance, interconnected spheres (as shown in Fig. 2) or ellipsoids. This approach is usually referred to as a discrete element method (DEM). One important feature of this method is the ability to model fibers with various equilibrium shapes, aspect ratio and stiffness. The geometry of the fiber can be then controlled by the size and the number of the particles (Elata and Berryman 1996; Guo et al. 2013; Switzer 2002; Switzer and Klingenberg 2003).

Fig. 2
figure 2

Discrete element-based representation of a flexible fiber

The dynamics of fibers is modeled by tracking the motion of each particle (e.g., sphere) under internal interactions with the neighboring particles through bonds and external forces, such as hydrodynamic forces, contact forces with other fibers and gravitational forces. In summary, the dynamical process includes the following steps:

  • Detect contact among particles and their surroundings (e.g., walls).

  • Apply force–displacement law to compute the contact forces based on the relative motion between two entities at the contact and the contact model (as will be discussed in Sect. 2.2).

  • Apply the law of motion (Newton’s law) to each particle to update its velocity and position based on the resultant force and moment arising from the contact and any body forces acting on the particle. The governing equations obtained by applying the conservation of linear and angular momenta can be expressed as

    $$m_{\text{p}} \frac{{{\text{d}}v_{\text{p}} }}{{{\text{d}}t}} = F_{\text{s}} + F_{\text{b}}$$
    (1)
    $$F_{\text{s}} = F_{\text{d}} + F_{\text{p}}$$
    (2)
    $$F_{\text{b}} = F_{\text{g}} + F_{\text{con}} + F_{\text{coh}}$$
    (3)
    $$I_{\text{p}} \frac{{{\text{d}}\omega_{\text{p}} }}{{{\text{d}}t}} = T_{\text{c}}$$
    (4)

    where \(m_{\text{p}}\) and \(v_{\text{p}}\) are the particle mass and velocity, \(F_{\text{s}}\) is the surface force, \(F_{\text{d}}\) is the force due to drag, \(F_{\text{p}}\) is the force due to pressure, \(F_{\text{b}}\) is the body force, \(F_{\text{g}}\) is the force due to gravity effect, \(F_{\text{con}}\) is the force due to contact developed between particles and their surroundings, \(F_{\text{coh}}\) is the attractive body force due to cohesion, \(I_{\text{p}}\) is the mass moment of inertia of the particle, \(\omega_{\text{p}}\) is the particle angular velocity and \(T_{\text{c}}\) is the contact torque.

2.2 Contact, bonding and cohesion

The mechanisms that produce fiber aggregation are highly related to fiber–fiber and fiber–wall contacts. A contact model describes the physical behavior occurring at a contact (a small area where the physical touching and eventual overlap take place). Although a particle aggregation may exhibit complex dynamical behavior, this behavior may be achieved using relatively simple contact models. The contact model consists of up to three parts: (1) a contact-stiffness model, (2) a slip and separation model and (3) a bonding model. The contact-stiffness and the slip and separation models fully describe the physical behavior at all particle–wall and unbonded particle–particle contacts.

A contact-stiffness model provides an elastic relation between the contact force and the relative displacement for unbonded particles. For instance, the damped linear spring as shown in Fig. 3 may be utilized. This model is composed of an elastic spring and a dashpot. The normal/tangential force (Fn/Ft) is generated by the repulsive elastic spring, and the damping force is applied by the dashpot, which dissipates energy during the contact. A more complex inelastic collision can be represented by the Hertz–Mindlin contact model (Elata and Berryman 1996).

Fig. 3
figure 3

Schematic of a damped linear spring model

The slip model allows two entities in contact to slide relative to one another due to friction. They can also separate if a tensile force develops between them, and they are not bonded together. The slip condition occurs when the shear component of force reaches the maximum allowable shear contact force. This force is usually taken to be the minimum friction coefficient of the two contacting entities multiplied by the magnitude of the compressive normal component of the force.

The bonding model makes the fibers elastically deformable. This model introduces attractive interparticle forces to the fiber. It uses the concept of a massless bar connecting a pair of bonded particles (see Fig. 4). The bar can transmit force and torque between neighboring particles and is subject to breakage under the maximum load. The bond radius controls the overall fiber flexibility. Connecting the spherical particles through larger bond radius yields fibers that are more rigid.

Fig. 4
figure 4

DEM interactions: particle bonding

Fibers in fluid flow may interact in two different ways: collisions and cohesive interactions. Collisions occur when the fibers are in physical contact and are usually described by contact-stiffness models discussed above. Cohesion is the property of particles (of the same material) to stick to each other due to mutual attraction. The cohesive interactions may be governed for instance by the square-well potential in which the attractive forces come into play after the separation distance between two particles reaches a specific distance dc. A typical interaction sequence is shown in Fig. 5. After the particles are close enough to each other, particles move in a straight-line trajectory, which leads to an inelastic collision. Because of the repulsive forces, particles travel in a straight path until the departing-cohesive interaction, which occurs at a separation distance of dc. If the magnitude of these forces is not large enough, an internal reflection of the particles from the well occurs, which is referred to as a capture–cohesive interaction and yields to fiber aggregation. If the repulsive forces are significant enough to overcome the cohesive well, an escaping cohesive interaction occurs and then the particles separate (fiber dispersion). We note that for purpose of simplicity, the interaction in Fig. 5 is shown for only two particles and without incorporating other effects (e.g., wall, hydrodynamic). This physical aspect requires a careful investigation as being one of the key mechanisms of fiber flocculation and dispersion.

Fig. 5
figure 5

Fiber interaction sequence

2.3 Fiber–fluid coupling

The motion and interaction of fibers in a suspension/floc/cake are strongly influenced by hydrodynamic forces exerted by the suspending fluid. As such, incorporating this aspect is important to analyze the bridging performance. The motion of an incompressible fluid phase in the presence of a secondary solid phase (fibers) is governed by a modified set of Navier–Stokes equations given by

$$\frac{{\partial \alpha_{\text{f}} }}{\partial t} + \nabla \cdot \left( {\alpha_{\text{f}} \varvec{u}_{\text{f}} } \right) = 0$$
(5)
$$\frac{{\partial \alpha_{\text{f}} \varvec{u}_{\text{f}} }}{\partial t} + \nabla \cdot \left( {\alpha_{\text{f}} \varvec{u}_{\text{f}} \varvec{u}_{\text{f}} } \right) = - \alpha_{\text{f}} \nabla \frac{p}{{\rho_{\text{f}} }} + \nabla \cdot \tau - C_{{{\text{s}} - {\text{f}}}}$$
(6)

where \(\alpha_{\text{f}}\) is the volume fraction occupied by the fluid, \(\rho_{\text{f}}\) is the fluid density, \(\varvec{u}_{\text{f}}\) is the fluid velocity vector, \(\tau\) is the stress tensor of the fluid phase, \(p\) is the pressure, and \(C_{{{\text{s}} - {\text{f}}}}\) denotes the momentum exchange with the solid phase. We follow the Eulerian–Lagrangian approach to model the multiphase flow. As mentioned above, the Navier–Stokes equations are solved for the continuous phase (base fluid) while each individual particle within each fiber is tracked based on the Lagrangian approach. Such approach provides instantaneous information of the position and velocity of particles, and the forces being applied on them. It also has the capability to handle the flexibility of fibers, particle–particle collision (based on contact and cohesion models as described in the previous subsections), and the impact of walls. However, it may be computationally demanding for a large number of particles being tracked.

Given the solvers for fluid motion (CFD-based) and fiber dynamics (DEM-based), and a coupling interface for the interaction between the two phases (e.g., through enforcing the conservation of momentum), a scheme for calculating the dynamics of fiber suspensions can be formulated as described in Fig. 6. The geometrical and material properties of fibers along with the appropriate contact and cohesion models are first introduced. The geometry of the computational domain along with the fluid properties is specified as per the fracture and suspending fluid under investigation.

Fig. 6
figure 6

Modeling approach framework

The coupling routines are detailed in the flowchart shown in Fig. 7. These include the following sequence of steps:

Fig. 7
figure 7

Flowchart of the coupled computational model of the fiber–fluid system

  • Clumped particles (fibers) are injected into the fluid domain.

  • The position and orientation of each particle are obtained using the DEM-based approach and applying Newton’s law to each particle given the total force and torque that arise from the mechanical contact (developed between particles and their surrounding), bonding, cohesion and hydrodynamic effects.

  • For each particle, the corresponding mesh grid block is identified.

  • For each grid block, the volume fraction is determined.

  • The momentum exchange (source term in the modified set of the Navier–Stokes equations) between fluid and solid phases is calculated.

  • The flow and pressure fields are updated according to the current structure of the fiber suspension.

  • The hydrodynamic forces acting on the particles are computed and transferred to the DEM solver.

  • The same procedure is repeated. Fiber injection may be stopped when observing dominant fiber flocs and reaching the required fiber concentration. A careful selection of the timestep should be made to guarantee the convergence of the numerical solution (for both discrete fiber model and fluid flow) and avoid instability issues.

For the sake of simplicity and speeding up the numerical simulations, one can also employ the one-way coupling in which the impact of the fluid flow on fibers is accounted for but not the reverse. Other effects such as electrostatic, magnetic and heat transfer can be incorporated as well.

Having a reliable coupled numerical fiber–fluid model, one may perform a sensitivity analysis of the tendency for fibers to flocculate when varying the fiber properties and the operating conditions. In particular, the size of the fiber flocs and their role in the bridging performance can be correlated with the characteristics of the fluid (viscosity, flow rate and regime) as well as with the fracture size. Such studies are limited (if not unfeasible) when relying only on experiments.

3 Fiber dynamics in fluid flow

3.1 Fiber flocculation in sheared flow: mechanism of interest

Flocculation describes the aggregation of fibers into clumps when a highly concentrated suspension of fibers is sheared. Figure 8 presents an example of flocs obtained from industrial fibers suspended in a water-based fluid. It shows clearly the capability of fibers to entangle, flocculate and form a network structure when interacting with fluid flow. Having such feature, fibers are viewed as an attractive option for acid diversion. During the acidizing jobs, operators use mixers (as shown in Fig. 9) to prepare the fiber-based diverting fluid (obviously by applying shear) prior to its injection to the well.

Fig. 8
figure 8

Fibers flocculated in suspensions. a Dry fibers. b Fibers dispersed in water-based fluid. c Entangled fibers forming a network structure. Pictures were taken from photographs courtesy of Schlumberger

Fig. 9
figure 9

Batch mixing of the degradable diverting acid system (Asiri et al. 2013)

Currently, there is a need to gain an insight into the fundamentals of how fluid diversion works in carbonates by different diverting mechanisms and diverters. As an attempt to address this, we develop a numerical model to examine fiber dynamics in fluid flow and investigate, in particular, the impact of different parameters on the bridging performance.

3.2 Numerical model setup

As a representation of a natural fracture, we consider a geometry consisting of a cylindrical pipe with a restriction at the end. This constitutes a representative geometry of lab-scale setup and imitating a natural fracture aperture to some extent. The variation in the pipe diameter before reaching the restriction is smooth (inclined boundary). The total length of the pipe is 5.5 in., and its diameter varies from 1 to 0.25 in. The dimensions of the simulated pipe with a restriction are designed to mimic the typical openings encountered in carbonates. At the inlet face (left of the pipe), we assume a uniform and constant inlet flow rate (Dirichlet boundary condition). A fixed pressure outlet (atmospheric pressure) is considered at the end of the pipe restriction (Neumann boundary condition). At the side walls, no penetration and no slip conditions are applied. The model geometry is built and meshed with the engineering software Star-CCM+ (Siemens 2017) in which all the computation and post-processing of the subsequent numerical results are performed. This software tool has demonstrated great capability to simulate multiphase flow with good accuracy in comparison with experiments (Marafaing et al. 2018; Tang et al. 2015). The system of incompressible fluid flow governing equations is numerically solved employing the finite volume method as a space discretization technique with an algebraic segregated solver. The geometry and computational mesh are shown in Fig. 10.

Fig. 10
figure 10

Geometry of the computational domain: pipe with a restriction at the end

As for the DEM-based fiber model, the particles are spheres with a diameter of 1.6 mm and a density of 1.4 g/L. More details on the particle properties and relevant model parameters are presented in Table 1.

Table 1 DEM-based fiber model properties

3.3 Results and discussion

3.3.1 Fiber calibration

We evaluate the fiber Young’s modulus using a simply supported beam under a central point loading test (see Fig. 11a) while varying the radius of the bond that connects two neighboring particles and controls the fiber flexibility. Based on the Euler–Bernoulli beam theory, the maximum bending of the fiber can be expressed in terms of the fiber length \(L\), the applied force \(F\), Young’s modulus \(E\) and the second inertia moment \(I\) as

$$w_{L/2} = \frac{{FL^{3} }}{48 EI};\quad I = \frac{\uppi }{4}r^{4}$$
(7)

where \(r\) is the particle radius. Therefore, by computing the central fiber bending, one could extract Young’s modulus. The variations of Young’s modulus with the bond radius multiplier (bond radius/particle radius) are shown in Fig. 11b. Clearly, increasing the thickness of the fibers of the bond connecting two neighboring particles within the fiber increases its stiffness. One could tune then the bond radius to obtain the fiber flexibility of interest. Note that Young’s modulus of the fibers simulated in the present study is in the order of few MPa.

Fig. 11
figure 11

Mechanical test to assess fiber rigidity. a Simply supported beam under central point loading and equivalent DEM-based discrete model. b Estimated fiber Young’s modulus as a function of bond radius multiplier

3.3.2 Fiber clustering and bridging

We perform a numerical experiment to study the fiber clustering and bridging in incompressible fluid flow. We first solve numerically for the fluid flow, check its convergence to steady state and then start injecting flexible fibers. Details on the relevant properties of the fibers and the suspending fluid are summarized in Table 2. The results are obtained while considering one-way coupling, that is, only the fluid flow affects the fiber motion and deformation. We present in Fig. 12 the set of snapshots corresponding to successive time steps during the fiber clustering and the flow of the formed fiber flocs along the pipe. As would be expected, the fluid flow is accelerated when reaching the restriction to ensure the conservation of mass. Color levels denote the size of the fiber clusters. Fibers are injected randomly from the inlet and undergo artificial attractive forces in the sphere placed near the inlet to form fiber flocs before reaching the restriction. Some of the fibers interact with the walls and slow down. The fibers of the first floc enter the restriction all together with a sliding motion between each other and lead to temporary and rapid bridging (t = 0.48 s). We observe a strong rearrangement of the fibers as they arrive at the restriction. The following fiber floc with bigger size (hitting the restriction at about t = 0.8 s) gets stuck at the restriction and leads to permanent bridging. The subsequently injected fibers continue to accumulate after bridging. In this case, the concentration of fibers is high enough to provoke many collisions between flowing fibers and prevent any rearrangement.

Table 2 Fiber–fluid system setup (numerical simulations shown in Fig. 12)
Fig. 12
figure 12

Snapshots corresponding to successive timesteps during the flow of fiber flocs. Fibers are injected randomly from the inlet and undergo artificial attractive forces in the sphere placed near the inlet to form fiber clustering before reaching the restriction. These forces produce the flocculation behavior of the fibers when they get clustered to reproduce qualitatively the experimental observation as shown in Fig. 8

3.3.3 Sensitivity analysis

To characterize the fiber floc size, we introduce the mean floc diameter as a function of the fiber length \(L\), the crowding number \(N\) and the fiber aspect ratio \(c\), as given by Switzer (2002):

$$D_{\text{floc}} = \frac{L}{2}\left( {1 + \frac{N}{c}} \right);\quad N = \frac{2}{3}\phi_{\text{f}} c^{2}$$
(8)

Here \(\phi_{\text{f}}\) is the fiber volume fraction. We define the floc–frac size ratio as \(\beta = \frac{{D_{\text{floc}} }}{D}\), where \(D\) is the restriction diameter. We note that the crowding parameter is considered in the present study to provide a better measure of fiber impact in comparison with the volume concentration of pure fibers (Lundell et al. 2011). Suspensions of flexible fibers in pipe flow are simulated for a variety of values of the system parameters. These include fiber Young’s modulus, mass flow rate and fluid viscosity. The effect of these parameters along with the floc–frac size ratio (proportional to fiber concentration) on the tendency for fibers to bridge is illustrated in Fig. 13. Flocs formed of rigid fibers tend to flocculate and bridge, whereas higher fiber concentration is required to achieve bridging when they are flexible. Being more rigid, they can resist more the pressure and drag forces at the pipe restriction, and then they can bridge more easily. This mechanism was observed experimentally; Soszynski and Kerekes (1988) reported that nylon fiber flocs dispersed when the fiber stiffness was reduced by heating above the glass transition temperature of nylon. Increasing the mass flow rate or fluid viscosity leads to higher drag and pressure forces acting on fibers and then this would provoke their rearrangements when reaching the restriction and prevent the bridging to take place. Thus, fiber flocs of larger size and higher concentration are needed to achieve bridging when dealing with more accelerated flow and/or viscous suspending fluids.

Fig. 13
figure 13

Effect of a fiber flexibility, b flow rate and c fluid viscosity on the fiber tendency to bridge

Table 3 summarizes the effects of different parameters on the tendency to bridge. The trends in the dynamical behavior of fibers under varying operating conditions are concluded based on numerical experiments. In addition to the effects discussed above, we examine the influence of the interfiber cohesion (represented by a linear cohesion model in our numerical simulations, i.e., the application attractive forces proportional to the eventual particles overlap when colliding) on the fiber dynamics when interacting with the fluid flow. Activating the cohesion yields more compact flocs that promotes the tendency to bridge.

Table 3 Key parameters affecting fiber tendency to bridge

Several models for bridging of particulates have been proposed in the literature (Garagash et al. 2019; Osiptsov 2017). In the pressure driven flow of a suspension through a fracture, the bridging is observed to occur when the ratio of the fracture width to the particle diameter \(w/d\) is within the range of 2.5–3.0. A more recent study reported a modified version of the bridging criterion in a fracture which accounts for the particle volume fraction in the suspension flow (Osiptsov 2017). The critical fracture width is given by:

$$w^{*} = \hbox{min} \left[ {b,1 + \frac{{C_{\text{p}} }}{0.17}\left( {b - 1} \right)} \right]d$$

where the bridging factor \(b\) is commonly taken equal to 2.5 and \(C_{\text{p}}\) is the particle volume fraction. It was shown that the value of the bridging factor can vary significantly based on the roughness of the walls (Garagash et al. 2019). The aforementioned bridging criteria currently used by several commercial simulation tools depend only on the particles size and concentration (Garagash et al. 2019; Osiptsov 2017). In the present study, we study the bridging mechanism of fibers in suspension flow. We found evidence of noticeable impact of the fiber stiffness, the base fluid viscosity and the inherent cohesion among fibers on the formation and stability of fiber flocs for temporary plugging.

As a first step to verify the numerical predictions against experimental observations, based on flow loop tests conducted at Schlumberger Dhahran Carbonate Research Center, and establish a level of confidence in the numerical transport model, we perform a set of simulations while reproducing most of the operating conditions (fluid and fiber properties as reported in Table 4), and Fig. 14 shows the tendency for fibers to bridge at the restriction while varying the fiber loading and fluid flow rate. The shaded zone points out the range of values for the fiber concentration and flow rate that leads to fiber bridging. Clearly, increasing the flow rate lowers the tendency for fibers to bridge. The numerical results match closely with the experimentally observed fiber loading (15 g/L) needed to achieve bridging (represented by the dashed line in Fig. 14).

Table 4 Fiber–fluid system setup (numerical simulations shown in Fig. 14)
Fig. 14
figure 14

Bridging as a function of flow rate: experimental verification

Based on the simulation results, recommendations are made for the hydraulic/acid fracturing operations in the field: (1) Stiffer fibers enable superior plugging performance at reduced fiber concentration/loading. (2) Adjusting the injection rate can be an alternative to optimize the fiber loading needed to achieve effective plugging from an economic perspective. (3) The viscosity of the base fluid needs to be accounted for when designing fiber-assisted diverting fracturing jobs given its significant impact on the plugging mechanism.

3.3.4 Fully coupled fiber–fluid system

To obtain a more realistic representation of the fiber–fluid system, we perform a fully coupled simulation to examine the simultaneous effect of the fiber on the fluid flow and vice versa. To guarantee stability and convergence of the numerical solution, we reduce the time step to 10−4 s.

We present in Fig. 15a, c the transient variations of the differential pressure (between the inlet and the outlet) and the corresponding snapshots of the flow of fiber flocs obtained from the numerical model. The first snapshot (see Fig. 15c) shows the flocculation of the fibers and their flow toward the restriction. At this level, the differential pressure is relatively low. As the fibers entangle and form a plug at the pipe restriction, the pressure rises. The pressure keeps increasing as the fibers accumulate and block the fluid to flow through the restriction until reaching a peak. The pressure buildup is the manifestation of the resistance of the bridge. Then, the pressure falls sharply to stabilize at a low value when the fibers exit entirely the pipe. The numerical results are in qualitative agreement with the experimental observations of the pressure trend when injecting fibers (Cohen et al. 2010) as shown in Fig. 15b, d. The flocculation and bridging process takes a few seconds. The experiments were conducted using the laboratory-scale equipment for testing leakoff behavior. As shown in Fig. 15d, this equipment includes a tube fitted with a piston, a pump and an orifice. The acid and fibers are pushed through the piston while recording the transient pressure as the fiber-laden acid forms plugs and passes through the orifice.

Fig. 15
figure 15

Transient variations of differential pressure while injecting fibers (upper frames) and flow of fiber suspensions (lower frames). Results are obtained from the numerical model and laboratory experiments for testing leakoff behavior. a Differential pressure versus time (numerical simulations). b Pressure signal during fiber-laden acid bridging experiment (Cohen et al. 2010). c Successive snapshots of the flow of a fiber suspension. Initially, the fibers are oriented randomly. They flip and rearrange as they flow toward the restriction and exit the pipe. d Experimental setup and example of fiber plugs formed in different geometry channels (Cohen et al. 2010)

The computational time to simulate the dynamics of about 400 flexible fibers in fluid flow for 1-s physical time is approximately 12 h on 3.3 GHz processors. The simulation time is reduced to only 2 h when applying the one-way coupling and deactivating the impact of fibers on the fluid. As such, the fully coupled model is not suitable to conduct the parametric study which requires many simulations (hundreds if not thousands) to assess the bridging mechanism for different operating conditions. We have only shown few cases to demonstrate the possible tuning of the developed model to account for the two-way coupling. However, its use remains impractical for sensitivity analysis and optimization purposes.

Next, we decrease the mass flow rate and investigate the bridging mechanism. Figure 16 shows the transient variations of the differential pressure resulting from fiber injection. The results are obtained from the fully coupled simulations for a physical time of one second. The sharp increase in the differential pressure is the manifestation of the formation of fiber bridge. The stabilization of the pressure level at higher value is indicative of the occurrence of permanent fiber bridge.

Fig. 16
figure 16

Transient variations of differential pressure obtained from the fully coupled model (permanent fiber bridge)

Figure 17a shows the simulation time requirements of the case presented in Fig. 16. Figure 17b shows the corresponding simulation time when running the one-way coupled model. The results are presented in terms of elapsed time per time step. Clearly, the one-way coupled model is much less computationally demanding.

Fig. 17
figure 17

Simulation time requirements. a Two-way coupling. b One-way coupling

4 Summary and conclusions

In this paper, we developed a numerical model to simulate fiber flocculation and bridging in fluid flow. A discrete model of flexible fibers is coupled with a fluid flow solver to account for the inherent simultaneous interactions. This study aims at predicting the fiber dynamics when the fibers interact with suspending fluids and encounter restrictions that can be representative of fractures or wormholes in carbonates. Based on the present numerical analysis, we observed that the bridging capability of the fiber–fluid system could be enhanced by

  • Increasing the fiber loading

  • Increasing the fiber stiffness

  • Reducing the flow rate

  • Reducing the suspending fluid viscosity

  • Increasing the attractive cohesive forces among fibers (sticky fibers flocculate and bridge more easily)

These trends in the dynamical fiber behavior (tendency for bridging under varying operating conditions) obtained from the numerical model show consistency with experimental observations. Furthermore, the fully coupled numerical simulation gives a pressure response (usually used to track the bridging) with features similar to that measured experimentally. This modeling capability can provide a baseline and guidance to improve the efficiency of existing fiber-based systems and to design new ones to fit the needs of new potential well intervention applications.