Micromechanics of proppant agglomeration during settling in hydraulic fractures

Abstract

The micromechanics of proppant settling in quiescent fluid in a rough and relatively narrow rock fracture is investigated. The study focuses on particle–particle and particle–wall interactions in a dense-phase particle settling. The study used a coupled discrete element method and computational fluid dynamics (DEM–CFD) code because DEM–CFD is the most suitable computational method for modeling the frequent interactions of dense assembly of rigid particles and enables modeling of two-way solid–fluid interactions. Due to frequent particle–particle interactions of grains submerged in fluids, the particle interaction model in DEM is improved by the incorporation of the effects of lubrication due to a layer of fluid surrounding particles. Results of the numerical study are compared to previous experimental and theoretical relationships. The findings of the study highlight conditions that lead to proppant aggregation due to the fluid viscosity and fracture width in relation to particle diameter ratio. In the light of the DEM–CFD results, it was found that published relationships are inadequate in describing the settling rates for proppant in a rough and narrow hydraulic fracture and high fluid viscosity. Micromechanical particle interactions during settling and erratic upward and fluid counter-flow may cause proppant trajectories that are not always in the direction of gravity in a rough fracture resulting in clogging of the fracture or forming faster settling particle agglomerates. The maximum packing density 0.3–0.4 (3.9–5.9 lbs/gal) of proppant in a narrow and rough hydraulic fracture was obtained in this study, which is lower than the usually assumed one of 0.5 (9.8 lbs/gal) for any given fracture roughness.

Introduction

Although much attention has been devoted to modeling of proppant flow and transport in hydraulic fractures, relatively little effort has been devoted to the understanding of proppant agglomeration in narrow and rough fracture zones. This study explains the effect of fluid lubrication as a new important agglomeration mechanism that has not yet been incorporated in the body of knowledge of proppant flow and transport. This study explains conditions that lead to narrow channel clogging from the micromechanical particle–fluid interaction perspective. An improved understanding of proppant behavior in rough and narrow fracture zone will assist reservoir engineers to recognize and mediate problems associated with fracture clogging during hydraulic fracturing. During the shut-in stage, when fluid pressure deceases due to the termination of the hydraulic fracturing process, proppant settling is pronounced, as well as fracture width decreases due to formation of in situ stresses. Therefore, the findings presented in this paper can be directly applied to the shut-in stage. In this study, conditions for particle agglomerations as a consequence of settling in narrow rough fracture zone are explored. In particular, parameters such as slurry density, fluid viscosity and fracture width are evaluated toward the possibility of proppant agglomeration in 2, 4 and 8 mm-wide fracture. The objective of this study is to investigate: (1) whether or not the existing proppant settling relations can be used also for narrow fracture zones and (2) under which conditions they are valid. Curve-fitting models were previously derived for predicting proppant settling rates, including the effects of the particle concentration in the fluid and the relative size of the fracture width relative to the proppant diameter. For proppant settling in pipes and narrow slots, existing experimental correlations relate average slurry settling velocity to a single particle settling velocity based on particle concentrations assuming uniform spatial distribution in the fluid (Hannah and Harrington 1981; Harrington et al. 1979; Novotny 1977; Shah 1982, 1993; Shah and Lee 1986; Shah and Lord 1990; Shah et al. 2007).

Particle agglomerations and irregular particle motions were previously observed in several experimental studies on flow of particle–fluid mixtures in pipes and slot flow apparatus. Hammond (1995) concluded that fluid mechanical mechanisms responsible for particle migration are not fully understood. For example, significant migration of solids toward the experimental pipe center was observed in non-Newtonian cross-linked polymer solutions. Particle migration occurred in 1 cm-diameter pipe of 2 m length, but under the same conditions, tests with Newtonian-based slurry showed very little migration toward the center of the pipe within 2 m length (Hammond 1995). As a consequence, the observed settling rates were greater when proppant had migrated into a central column than when it remained uniformly distributed across the fracture cross section. Aggregation of particles falling in viscoelastic fluids was studied by Joseph et al. (1994), who have related fluid normal stresses of viscoelastic fluids with particles agglomeration and chaining (Fig. 1). Joseph et al. (1994) observed experimentally that large values of the elastic stress ratio σ _n/τ (where σ _n is the fluid normal stress and τ the shear stress) are sufficient, but not necessary for strong inter-particle interactions that lead to particle agglomeration. On the contrary, another set of experiments comparing settling in Newtonian and Boger fluids showed that fluid aggregation is not caused by fluid elasticity. Gheissary and van den Brule (1996) concluded that a change in microstructure of the fluid has an effect on the settling velocity that goes beyond its effect on the viscosity alone and plays an important role in particle–particle interactions.

Fig. 1

Experimental results of particle settling in 1.25 % aqueous polyox solution (a), SI (b) and STP (c) (Joseph et al. 1994)

In another study, the aggregation of spherical particles was observed in non-Newtonian fluids during experiments where a suspension of 60–70 μm spheres in aqueous polyacrylamide solution was placed between two glass plates. Rectilinear and rotary shear was applied by moving glass plates, which resulted in different types of aggregation (Michele et al. 1977). Shear thinning plus memory of certain fluids, which create corridors of reduced viscosity, can also contribute to particle agglomeration. Figure 1 shows experimental results of particle settling in: (a) 1.25 % aqueous polyox solution, (b) SI and (c) STP (Joseph et al. 1994). SI is a solution of 5 % by weight of polybutylene in decalin plus 50 % polybutylene oil, and STP is a solution of polybutylene in petroleum oil and weakly viscoelastic fluid with small normal stresses and basically constant viscosity for shear rates <100. Feng et al. (1994) observed how channel walls affect sedimentation of a single particle and concluded that walls increase the drag on the particles and settling particle always keeps a certain distance from the wall due to repulsion in a Newtonian fluid. They indicated that this effect may be of significance in transport of particulate mixtures in fluids (Feng et al. 1994). Clark et al. (1981) reviewed different methods for computing vertical settling velocities and hindered settling, and concluded that hindered settling equations gave somewhat different values in reviewed studies. Clark et al. (1981) also investigated wall effects on proppant settling in non-Newtonian fluids using static and dynamic tests. Wall effects diminished at slot aperture versus particle diameter ratio of 25. They also observed that particle concentration effects on particle settling agree with previously published relations, but for relatively low volumetric concentrations. At higher volumetric particle concentrations, particle segregation was observed. At particle volumetric concentration lower than 0.36 in low-viscosity fluids, experimental particle settling velocities fall into the calculated boundaries, with a major deviation at higher fluid viscosity. No satisfactory explanation for deviations from the Stoke’s law exist (Clark et al. 1981). Roodhart (1985) also concluded that Stoke’s law based on Power law viscosity is insufficient for predicting particle fall rates in both flowing and quiescent fluids. Kirkby and Rockefeller (1985) observed that sedimentation behavior in fracturing fluids was dominated by particle clustering phenomenon. He also concluded that settling velocities depend on proppant concentrations, but in a very different manner than generally assumed by proppant transport models.

Electrostatic theory of agglomeration refers to a system that consists of a large number of small particles in a suspending fluid. While particles collide with each other because of Brownian motion, they can be attracted to each other so that they stick together. When solid particles are immersed in a fluid, there is a tendency for ions of one sign to be preferentially adsorbed into the solid and for the oppositely charged ions to remain in the neighboring fluid. Inter-particle potentials are important in determining the nature of agglomerates and the kinetics of agglomeration, as well as keeping the colloid stability. Deryaguin, Landau, Verwey and Overbeek (DLVO) theory of colloid stability is based on the recognition of electrostatic repulsion which opposes aggregation and a universal attractive van der Waals force which acts to bind particles together (Allen and Smith 2001). Turbulent shear and inertial particle agglomeration models (Allen and Smith 2001; Saffmann and Turner 1988) define the agglomeration rate of aerosol particles based on their sticking efficiency, particle radius, fluid properties and turbulent energy dissipation rate. In published experiments, particle agglomeration was observed, but what causes it was not explained. For lubrication causing the agglomeration, different shear fluid behavior only promotes collision frequency, but does not change local fluid properties in normal direction between two particles that approach each other.

Theoretical and experimental models of proppant settling

Proppant settling was studied in different experimental slot apparatus in the past and curve-fitting relations were used for predicting the amount and velocity of proppant settling in hydraulic fracturing simulators. For very low proppant concentrations and when the wall effects are negligible, Clark et al. (1981) proposed the following relationship:

$$v_{\text{s}} = \frac{1}{{1 + 6.88c_{\text{v}} }}v_{0} ,$$

(1)

where v _s is the slurry settling velocity, v ₀ the single particle velocity and c _v the volumetric particle concentration in the slurry. The wall effects appear to diminish significantly at the average slot width to particle diameter ratio of 25. Daneshy (1978) proposed the following equation:

$$v_{\text{s}} = v_{0} \left[ {\frac{{\left( {1 - c_{\text{v}} } \right)}}{{10^{{1.82c_{\text{v}} }} }}} \right].$$

(2)

Gadde et al. (2004) presented an empirical correlation fitting the experimental data curve (Gadde et al. 2004):

$$v_{s} = v_{0} \left( {2.37c_{\text{v}}^{2} - 3.08c_{\text{v}} + 1} \right).$$

(3)

Figure 2 shows a comparison of predicted settling velocities from published relations. As the particle volumetric concentration increases, the ratio of settling rates of concentrated slurry to a single particle settling decreases. The maximum packing is found to be around c _v = 0.6 or less.

Fig. 2

Correlations for the effects of solid concentration on settling velocity

The empirical correlation for the effect of fracture walls on proppant settling was obtained from the large vertical slot experiments for proppant flow in hydraulic fracturing research (Liu and Sharma 2005). The experimental data appear to locate on two lines: one when the ratio of particle diameter to cell width is <0.9 and the other when the ratio exceeds 0.9. The empirical correlations are:

$$\frac{{v_{\text{w}} }}{{v_{0} }} = 1 - f(\mu )\frac{r}{B};\left( {\frac{r}{B} < 0.9} \right),$$

(4)

$$\frac{{v_{\text{w}} }}{{v_{0} }} = g(u)\left( {1 - \frac{r}{B}} \right);\left( {\frac{r}{B} \ge 0.9} \right),$$

(5)

$$f(\mu ) = 0.16\mu^{0.28} ,$$

(6)

$$g(\mu ) = 8.26e^{ - 0.0061\mu } ,$$

(7)

where v ₀ (cm/s) is the settling velocity of a single unbounded particle, v _w the settling velocity corrected for the presence of walls, μ the fluid viscosity (or apparent viscosity for non-Newtonian fluids), r the particle radius and B the cell half-width. Figure 3 shows the comparison of the above relations with experiments (Liu and Sharma 2005).

Fig. 3

Wall effect on particle settling (Liu and Sharma 2005)

Novotny (1977) proposed the prediction for wall effect for low Reynolds numbers (Re ≤ 1):

$$\frac{{v_{\text{w}} }}{{v_{0} }} = 1 - 0.6526{\kern 1pt} (d/w) + 0.147{\kern 1pt} (d/w)^{3} - 0.131{\kern 1pt} (d/w)^{4} - 0.0644{\kern 1pt} (d/w)^{5} ,$$

(8)

where v ₀ (cm/s) is settling velocity of single unbounded particle and v _w is the settling velocity corrected for the presence of walls (Novotny 1977). The effect of cell walls on particle settling is quite significant, and it is found that it is more pronounced as the fluid viscosity is increased. The effects of fracture walls for high-viscosity Newtonian fluid and for power law shear-thinning fluid are similar. Both types of fluids show a significant reduction in the settling rate as the ratio of particle diameter to cell width increases. Wall effect corrections were experimentally derived by Shah (1982) from cylindrical tube. For creeping flow or Stoke’s region:

$$\frac{{v_{w} }}{{v_{0} }} = \left( {1 - \frac{{D_{\text{p}} }}{{D_{\text{c}} }}} \right)^{2.25}$$

(9)

and for turbulent flow or Newton’s region:

$$\frac{{v_{\text{w}} }}{{v_{0} }} = \left[ {1 - \left( {\frac{{D_{\text{p}} }}{{D_{\text{c}} }}} \right)^{1.5} } \right],$$

(10)

where D _p is the particle diameter and D _c the cylinder diameter.

The focus of the study presented in this paper is on particle interactions during settling at different concentrations in narrow and rough channel in Newtonian fluid, when inter-particle collisions frequently occur. The assumption is introduced, which postulates that particle–particle and particle–wall interactions that occur frequently in narrow channels can under certain conditions cause agglomerations at constant viscosity fluid regime. This study is performed using an improved discrete element method coupled with continuum fluid dynamics (DEM–CFD) with modified particle contact model. The novel particle contact model introduces fluid lubrication force during fully submerged particle collisions. The lubrication contact model allows for better modeling of dense particle–fluid slurries where particle collisions are frequent and affect both fluid field and particle motion. Fluid lubrication is a nonlinear force that is responsible for particle kinetic energy dissipation, while they approach each other in viscous fluid (Barnocky and Davis 1988; Davis et al. 1986). It was shown in this study that the kinetic energy dissipation is directly responsible for enhanced particle agglomerations, because it decreases particle velocity and they do not rebound, but stick next to each other in viscous fluid. This study reveals that particle clumps are not stable and may disperse with time as fluid drag forces push them apart. However, rough and wavy fracture walls can locally trap particles clump and may cause further disruption of proppant settling.

Fundamental research regarding the behavior of two particles falling in quiescent Newtonian fluid has been performed experimentally and theoretically (Fig. 4). Two spheres dropped side by side, which are initially separated by a small gap or no gap, are found to separate from each other (Joseph et al., 1994). Using direct numerical simulation results, this behavior is attributed to hydrodynamic mechanisms that cause circular particles to rotate and drift away from each other in a Newtonian fluid. Theoretical analysis on the interaction of two spheres settling in a fluid found that there were second-order attractive forces between two settling spheres (Brunn, 1977). Based on the findings for two side-by-side spheres that never attract in Newtonian fluid, the aggregation of spheres in multi-particle systems in a Newtonian fluid is usually assumed to be negligible. However, when two spheres are subjected to end-to-end settling, the well-known “drafting–kissing–tumbling” behavior is observed (Gheissary and van den Brule 1996; Liu and Joseph 1993), which indicates that the particle–particle interaction exists.

Fig. 4

Falling of two spheres in Newtonian fluid: a experimental results, b simulation results (Joseph et al. 1994)

Numerical model and validation

DEM–CFD model with improved contact lubrication model

Discrete element method coupled with computational fluid dynamics (DEM–CFD) as implemented in the commercially available Particle-Flow Code PFC^2D (Itasca 2004) is used in this study. The code is improved with a user-defined particle contact model that captures fluid lubrication effects during particles collisions. Two-dimensional geometry for the particles and flow field can be predictive for the three-dimensional case, because the DEM–CFD code has out-of-plane dimension equal to the particle diameter, in this case d = 0.6 mm. Particles are modeled as real spheres and the solid volume fraction is represented rigorously for the given model depth, i.e., one layer of particles. Vertical channels are considered, which have 2, 4 and 8 mm widths and 20 cm length. Particles are uniform in diameter, and different initial concentrations are used in the study. Initial particle concentration is achieved so that a number of particles confined in a square space are set to fall from above the x–x axis (but still submerged in fluid). An algorithm is written to add more of the same particles when the last particle crosses the x–x axis. This updated code is suitable to study micromechanical behavior of dense-phase slurry composed of proppant settling in Newtonian fluid in a narrow channel. DEM uses a finite difference scheme to compute individual particle trajectories, velocities, accelerations and forces within an explicit time-stepping procedure (Cundall and Strack 1979).

Particle motion is calculated using the law of motion and superposition of all the relevant forces that act on a particle at a given time step: contact forces, gravity, buoyancy and fluid drag. Fluid flow calculation in PFC^2D uses a numerical solution of the Navier–Stokes equations. Fluid drag, gravity and buoyancy forces act on a particle center of mass, while particle or wall contact and frictional forces, which induce particle movements and rotations, act on a particle surface. An important limitation of the DEM–CFD code is that fluid shear stress does not affect particle rotation, and non-Newtonian fluids cannot as yet be modeled. Fluid flow is modeled with CFD formulation that is based on modified Navier–Stokes equations and adds the effect of a particulate solid phase mixed into fluid, where particle radius is small compared to the length of a single fluid cell. Fluid and particle flow and transport are modeled as a two-way phase coupling, where forces are applied from the fluid to the particles and vice versa (Itasca 2004).

For modeling dense-phase flows with particle volumetric concentration higher than 0.1, a two-way phase coupling based on (Crowe et al. 2011) is recommended. The SIMPLE fluid coupling scheme (Patankar 1980) for incompressible viscous flow is applied to fluid cells within a fixed rectangular geometry. The effect of solid phase on fluid flow is modeled in a simplified manner as an average effect over many particles, instead of modeling the details of fluid flow between the particles. The effects of solid particles in the fluid flow are introduced to the numerical solution of Navier–Stokes equations in terms of porosity and coupling force (Bouillard et al. 1989). For porosities e < 0.8, the following fluid flow equation is implemented in PFC^2D (Bouillard et al. 1989):

$$\rho_{\text{f}} \frac{{\partial e\vec{v}}}{\partial t} + \rho_{\text{f}} \vec{v}\rho \nabla (e\vec{v}) = - e\nabla p + \mu \nabla^{2} (e\vec{v}) + \vec{f}_{\text{b}} ,$$

(11)

$$\frac{\partial e}{\partial t} + \nabla (e\vec{v}) = 0,$$

(12)

where ρ _f is the density of the fluid, $\vec{v}$ the fluid velocity, t the time, ρ the particle density, $\nabla$ the del operator indicating partial spatial derivative operator (with respect to the x and y directions in two-dimensional space), f _b the body force per unit volume, and p the pressure and μ the fluid dynamic viscosity. The fluid velocity can be observed in these equations as the interstitial velocity $\vec{v}$ and the macroscopic Darcy velocity defined as e $\vec{v}$. The interstitial velocity is the actual velocity that a parcel of fluid has as it moves through the pore space. The macroscopic velocity is the volumetric flow rate per unit cross-sectional area. This is a non-physical velocity, since it assumes that the flow occurs across the entire cross-sectional area where flow actually only occurs in the pore space. Equation (11) gives the fluid momentum equation with included effects of varying porosity and Eq. (12) is the conservation of mass (continuity) equation for an incompressible fluid in a porous medium. It can be seen that both time and space derivatives of porosity are present in the formulation. The porosity e is defined as:

$$e = 1 - \frac{{V_{\text{p}} }}{V},$$

(13)

where V _p is the total volume of all the particles contained within the fluid element and V the volume of the fluid element. Boundary conditions are applied as a slip or non-slip velocity function at the outside of the fixed grid geometry and initial conditions are applied as a fixed normal velocity or a fixed pressure and a fixed or a free shear velocity. The drag force applied by particles to fluid in each fluid element is defined as:

$$\vec{f}_{\text{b}} = \beta \vec{U},$$

(14)

where $\vec{f}_{b}$ is the drag force per unit volume, β is a coefficient and $\vec{U}$ is the average relative or slip velocity between the particles and the fluid, defined as:

$$\vec{U} = \bar{\vec{u}} - \vec{v},$$

(15)

where $\bar{\vec{u}}$ is the average velocity of all particles in a given fluid element and $\vec{v}$ is the fluid velocity. Different expressions for coefficient β are given for porosities with values higher and lower than 0.8 (Ergun 1952):

$$\beta = \frac{(1 - e)}{{\bar{d}^{2} e^{2} }}\left( {150(1 - e)\mu + 1.75\rho_{\text{f}} \bar{d}\left| {\vec{U}} \right|} \right);e < 0.8,$$

(16)

$$\beta = \frac{4}{3}C_{\text{d}} \frac{{\left| {\vec{U}} \right|\rho_{\text{f}} (1 - e)}}{{\bar{d}e^{1.7} }};e \ge 0.8,$$

(17)

$$C_{\text{d}} = \left\{ {\begin{array}{*{20}l} \begin{gathered} \frac{24}{{Re_{\text{p}} }};Re_{\text{p}} \le 1 \hfill \\ \frac{{24\left( {1 + 0.15Re_{\text{p}}^{0.687} } \right)}}{{Re_{\text{p}} }};\quad 1 < Re_{\text{p}} \le 1000 \hfill \\ \end{gathered} \hfill \\ {0.44;Re_{p} > 1000} \hfill \\ \end{array} } \right.,$$

(18)

$$Re_{p} = \frac{{\left| {\vec{U}} \right|\rho_{\text{f}} e\bar{d}}}{\mu },$$

(19)

where $\bar{d}$ is the average diameter of the particles occurring in the element and C _d is the turbulent drag coefficient defined in terms of particle Reynolds number Re _p. A force equal and opposite acts on the particles in each fluid element. The body force term $\vec{f}_{\text{b}}$ has a unit of force per unit volume. This force is applied to each particle proportional to the volume of each particle. The drag force applied to individual discrete particles is:

$$\vec{f}_{\text{drag}} = \frac{4}{3}\pi r^{3} \frac{{\vec{f}_{\text{b}} }}{(1 - e)},$$

(20)

where $\vec{f}_{\text{drag}}$ is the fluid drag force on the particle and r particle radius. The total force exerted by the fluid on the particle is the sum of the drag force and the buoyancy force if the gravity is turned on in the simulation. The particle motion in the fluid is then integrated over time according to the default PFC^2D scheme, where mechanical and fluid forces are added together by superposition.

To improve the submerged particle collision behavior in DEM, a lubrication force is introduced to the contact model using the elasto-hydrodynamic theory. The elasto-hydrodynamic deformation of a solid elastic sphere, which is immersed in a viscous fluid and in close motion toward another sphere or a plane solid surface, is modeled using the theoretical two-particle collision model (Davis et al. 1986). The model is based on the criteria for predicting whether two solid particles will stick together or rebound from each other subsequent to impact while immersed in the fluid. Lubrication force is derived from the Navier–Stokes equations for continuous fluid layer between particles and small elastic deformation (Davis et al. 1986):

$$F(t) = \frac{{6\pi \mu r^{2} v}}{x},$$

(21)

where F(t) is the particle contact force, μ the fluid dynamic viscosity, r the particle radius, v the relative velocity of two approaching particles and x the distance between particle surfaces. The lubrication force is introduced in PFC^2D for normal component of contact force as particles approach each other at small distance, using the apparent radius of the area around the particle for contact activation (Fig. 5).

Fig. 5

Schematic of the apparent (RI) and real (RT) radii and the approaching distance r _ij

The elastic rebound depends on the overlap of two particles if they are in contact with the real radii and the lubrication damping force acts upon contact when it is activated. The lubrication force is activated at that moment and its magnitude depends both on the approaching velocity and the distance between particles (or overlap). During the time-stepping procedure, if the particles are close enough that they overlap with their real radii, then the elastic rebound is activated (as well as the particle to particle friction). The contact force logic can be written here as:

$$r_{ij} = r_{i} + r_{j} - d_{ij}$$

(22)

if r _ij  > 2r _c = r _cr,

$$F_{\text{c}} = m\textit{\"{a}} = - 6\pi \mu a^{2} \frac{{v_{ij} }}{{x_{ij} }} = - l_{\text{c}} \frac{{v_{ij} }}{{(r_{ij} - 2r_{\text{c}} )}}$$

(23)

if r _ij  ≤ 2r _c = r _cr and

$$F_{\text{c}} = m\textit{\"{a}} = k(2r_{\text{c}} - r_{ij} ) + cv_{ij} ,$$

(24)

where F _c is the contact force, d _ij is the overlap of the particles, r _ij the distance between the particles centers, r _c the real particle radius (r _cr = 2r _c), k the spring stiffness, c the dashpot constant and l _c the lubrication constant. The particle rebound can be elastic (i.e., only the elastic contact stiffness k is active) or it can have some damping included (both k and the damping coefficient c are active), but this is only the behavior that occurs upon the collision of the particles. The expression for fluid lubrication in Eq. (23) has a singular point when the particle surface distance x _ij approaches zero (Balakin et al. 2013; Crowe et al. 2011). However, in the DEM code the singularity is avoided by activation of spring–dashpot rebound behavior just before the real particle surfaces come in contact. Since Eq. (21) is discretized over small time steps, the particle velocity gradually decreases while particles approach each other.

The tangential component of the immersed particle contact force is simplified and modeled such as for dry particle collisions. A friction coefficient of 0.2 typical for dry sand is used (Itasca 2004). Joseph and Hunt (2004) concluded that tangential component of particles submerged in fluid at the contact point is similar to a dry system, but with a lower friction coefficient due to the lubrication effects. However, for collisions of rough particles at increasing tangential Stokes number, the friction coefficient rises to a value that approximates the dry friction coefficient (Joseph and Hunt 2004). However, because of the lubrication force that increases as particles approach each other, complete particle surface contact is almost impossible to occur for particles submerged in fluid with higher viscosity when tangential friction is activated. This principle, for example, governs the lubrication of ball bearings (Davis et al. 1986). Turbulent flow is accounted for in the DEM–CFD code and depending on the CFD mesh size particles are caught by different velocity streamlines, which can further promote collision frequency and particle agglomeration. The particle radius is crucial for the relationship between drag force and normal particle contact force with lubrication.

Particle sphericity (or roundness) parameter is not introduced in the current formulation for fluid drag and lubrication contact force. Therefore, the code is valid only for particles that have round or well-rounded edges. For the fluid drag force, reference is made on previous investigations that have been conducted experimentally on settling quartz sand particles. Surface roundness plays a significant effect in particle settling velocity, decreasing it for about 50 % when the Corey shape factor (CFS) decreases to 0.3. Therefore, for typical quartz sand with CFS = 0.7, the settling velocity decreases by 15 % (Dietrich 1982). The conventional DEM contact behavior as well as the lubrication force contact formulation proposed here can be used only for spheres. The roundness parameter would be hard to introduce into the soft sphere DEM contact model. The particle roughness problem can be generally solved using particle clumps composed of several spheres, for which the lubrication contact force formulation can be adopted or using a non-spherical DEM code. The fines that can be found in natural proppants during hydraulic fracturing are also not considered in the DEM–CFD model, since particle size is assumed uniform. Percentage of fines, if present in the fracture, would further increase the fluid turbulence and particle collision frequency.

Model verification

The user-defined DEM contact model behavior with lubrication force is tested using the coefficient of restitution of particle–wall collisions and the results are compared with published data. Both particle–particle and particle–wall normal collisions are represented with user-defined contact model described above. Yang and Hunt (2006) investigated experimentally particle–particle immersed collision behavior, and they concluded that the same correlation between coefficient of restitution and Stokes Number can be applied for both particle–wall and particle–particle collisions. The singe-particle and wall collision immersed Stokes Number is defined as:

$$St = \frac{{mv_{1} }}{{6\pi \mu r^{2} }}$$

(25)

where m is the particle mass, μ the fluid dynamic viscosity, v ₁ the particle approaching velocity and r the particle radius. Stokes number represents a relation between particle inertia forces and fluid forces and, in case of single particle Stokes number, the numerator provides a measure of available momentum in the solid phase that sustains particle motion through the liquid. To measure coefficient of restitution, a particle is dropped from an initial height (h ₁) and the rebound high (h ₂) that particle reach after impact is measured. The coefficient of restitution e _r is defined as:

$$e_{\text{r}} = \frac{{v_{2} }}{{v_{1} }} = \frac{{\sqrt {2gh_{2} } }}{{\sqrt {2gh_{1} } }} = \sqrt {\frac{{h_{2} }}{{h_{1} }}} ,$$

(26)

where v ₂ is the velocity of the particle after the impact with the wall and v ₁ the velocity of the particle before the impact with the horizontal wall, while g is gravity.

Experimental verification of the immersed particle collision model developed by Davis et al. (1986) was performed on steel and glass spheres (Barnocky and Davis 1988), who found that the theoretical model was able to capture the main features of the experimental work, The critical Stokes number is a critical value for which the particle does not rebound back to the air from the horizontal surface and as a result the coefficient of restitution is zero. Although the Davis et al. (1986) equation was directly implemented PFC^2D, the predicted critical Stokes number is smaller than experimentally obtained (Fig. 6). A discrepancy related to the critical Stokes number value between theoretical and experimental values was found due to the elastic properties and asperity interactions of the spheres (Joseph and Hunt 2004; Joseph et al. 2001; Zhang et al. 2005). Motion of particles after collision is divided into normal and tangential components. The user-defined contact model in PFC^2D directly employs the theoretical model of Davis’ et al. for normal component of collision. In PFC^2D, the same model is used for particle–particle and particle–wall collision. Figure 6 shows a non-linear curve of coefficient of restitution versus Stokes number for a small sand particle and a range of fluid types considered in this study. Figure 6 also compares the results obtained in PFC^2D with semi-empirical model developed by Zhang et al. (2005). The comparison shows that Davis’s model has sufficient accuracy without implementing the surface roughness of a single sand particle and the types of fluid in this study. The applicability of this model to the real particle settling situation has several drawbacks. The lubrication theory was originally tested on steel and glass spheres, with several studies that take into account particle surface roughness (Zhang et al. 2005). In terms of elastic response, it can be assumed that sand particles behave similar to steel and a significant plastic deformation is not likely to occur during proppant placement. However, the angularity of natural sand is not taken into consideration.

Fig. 6

PFC results using the user-defined immersed contact model of wall and sand particle (r = 0.3 mm) after theory developed by Davis (1986) and Zhang et al. (2005) model a) h _min/h ₀ = 1/5, b h _min/h ₀ = 1/10 and c h _min/h ₀ = 1/20, where h _min is the particle surface roughness and h ₀ is the particle radius

Settling of a single particle in a quiescent Newtonian fluid was modeled using the DEM–CFD model to verify model behavior as function of Stokes’ law. Table 1 shows the results of one-particle drop tests with the particle lubrication user-defined contact model with three different fluids. The numerical results are compared with theoretical predictions commonly used for proppant flow and transport (Novotny 1977). The results indicate very good suitability of the DEM–CFD model to simulate particle settling in quiescent fluid. The discrepancy between Novotny (1977) solution and PFC^2D for intermediate particle Re numbers (1 < Re < 1,000) is caused by different equations that those two models use in this region.

Table 1 Settling velocities of one particle in quiescent fluid

The user-defined lubrication model employed in PFC^2D simulation is tested against theoretical and experimental results on the example of two medium sand particles settling side by side in infinite Newtonian fluid. It is confirmed that two particles separate initially and then settle side by side (Fig. 7). Three different fluid dynamic viscosity values were used for the simulations (μ = 0.05 0.01 and μ = 0.05 Pa s) and particle separation was observed in all cases. Figure 8 shows a simplified multi-particle settling case, where particle separation was observed as well. It is important to notice that Figs. 7 and 8 represent simulations where particle–particle contacts did not take place and particles were always separated. Therefore, fluid drag forces, gravity and buoyancy govern particle behavior.

Fig. 7

Two particles settling in Newtonian fluid (μ = 0.1 Pa s) side by side

Fig. 8

A set of particles settling in Newtonian fluid (μ = 0.01 Pa s)

Results

The fractures modeled are generic and can be representative of rough fractures in different rock types including sedimentary rocks. Fracture roughness causes more disturbance in proppant settling, because it promotes contact frequency between wall and proppant. For the purposes of the study presented in this paper, rock fracture roughness was as determined from published data. Barton and Choubey (1977) provided joint roughness coefficient (JRC) of rough fracture asperities in the form of an irregular zigzag line in real coordinates for visual determination of the JRC for different rock types (Barton and Choubey (1977); Barton and de Quadros 1997). The fracture profiles were digitized and used for modeling rough fractures as series of segments of straight lines forming the DEM walls. It is assumed in this study that there is no offset or lateral displacement between two surfaces of the fracture. The effect of particle collisions on settling parameters in a narrow channel with dense-phase flow is studied using the modified PFC^2D code with the model for fluid lubrication. Frequent particle collisions and presence of wall affect the settling rates, and fluid lubrication in particular adds to the increased particle agglomeration probability. Clark et al. (1981), Daneshy (1978) and Gadde et al. (2004) proposed relationships for the prediction of average proppant settling velocities based on settling velocity of a single particle in a quiescent infinitive fluid media. The results are presented using relative particle velocity parameter, which is the ratio of average proppant settling velocity and a single particle terminal settling velocity. The velocity of a single particle was calculated from Stokes’ law and modifications were made for higher Reynolds numbers using the relationship proposed in Novotny (1977).

For particles settling in a narrow channel, where the channel width is 3.3 times larger than particle diameter, micromechanical particle distribution of sand proppant particles was observed during proppant settling for initial particle volumetric concentrations of c _v = 0.1 (2.0 lbs/gal), 0.2 (3.9 lbs/gal) and 0.3 (5.9 lbs/gal). Table 2 gives the proppant concentration conversion from volumetric units to field units, to give the readers a better idea of studied proppant concentrations. Figures 9, 10 and 11 are oriented by initial particle volumetric concentration and sub-figures are oriented by fluid dynamic viscosities. Figure 9 gives an insight into the micromechanics of proppant settling in a very narrow fracture of 2 mm width for low initial proppant concentrations where c _v = 0.1 (2.0 lbs/gal). Figure 9a shows proppant settling in μ = 0.005 Pa s fluid after t = 0.5 s, Fig. 9b in μ = 0.01 Pa s fluid after t = 2.05 s, Fig. 9c in μ = 0.02 Pa s fluid after t = 0.8 s and Fig. 9d in μ = 0.1 Pa s fluid after t = 2.05 s. It can be seen that particles agglomerate after a while and continue to settle in clumps. Figure 10 shows that for proppant with higher concentration where c _v = 0.2 (3.9 lbs/gal). Figure 10a shows particles settling in μ = 0.005 Pa s fluid after t = 2.0 s and Fig. 10b shows particles settling in μ = 0.01 Pa s fluid after t = 3.0 s. Higher viscosity fluids were also investigated but are not shown in Fig. 10, because the fracture clogged after a very short time. It can be seen by comparing Figs. 9 and 10 that particles form even larger clumps at higher initial particle volumetric concentrations in Fig. 10 than in Fig. 9. If the proppant concentrations increase even more to c _v = 0.3 (5.9 lbs/gal), particles form larger agglomerates and tend to clog the fracture, while fluid has little or no space to flow upward around the clump (Fig. 11). Figure 11 shows an example of proppant bridging in a fracture for the case of quiescent fluid (μ = 0.005 Pa s fluid after t = 0.5 s), which may occur during shutdown in hydraulic fracturing. When the fluid is also moving in a horizontal direction, the fluid drag is larger than in quiescent fluid, which can lead to a different proppant configuration. In narrow fracture, it was not possible to successfully perform setting simulations with proppant volumetric concentrations above c _v = 0.3 (5.9 lbs/gal). It is interesting to note that for a narrow fracture, increase in fluid viscosity does not prevent particles from clumping as it is usually assumed in hydraulic fracturing models. The main mechanism that forces particles toward each other in a narrow channel is the upward flow of fluids, which pushes particles toward the center of the fracture.

Table 2 Conversion of proppant concentration units

Fig. 9

Micromechanics of proppant settling in w = 2 mm rough fracture for c _v = 0.1 (2.0 lbs/gal), where a shows proppant settling in μ = 0.005 Pa s fluid after t = 0.8 s, b in μ = 0.01 Pa s fluid after t = 2.05 s, c in μ = 0.02 Pa s fluid after t = 0.8 s and d in μ = 0.1 Pa s fluid after t = 2.05 s

Fig. 10

Micromechanics of proppant settling in w = 2 mm rough fracture for c _v = 0.2 (3.9 lbs/gal), where a shows particles settling in μ = 0.005 Pa s fluid after t = 2.0 s and b shows particles settling in μ = 0.01 Pa s fluid after t = 3.0 s

Fig. 11

Micromechanics of proppant settling in w = 2 mm rough fracture for c _v = 0.3 (5.9 lbs/gal), where a shows particles settling in μ = 0.005 Pa s fluid after t = 0.5 s

The blue arrows in Fig. 12 represent vectors of average fluid velocity in each fluid cell at a given time step. From several time snapshots shown, it can be seen that particle trajectories in a narrow channel did not follow the direction of gravity and that there is a certain competition for space between particles settling and fluid flowing upward. When particles are submerged in the unrestricted volume of quiescent Newtonian fluid, during downward motion due to gravity, particles settle and move away from each other, as shown experimentally and numerically in Figs. 4, 7 and 8. However, in a narrow fracture, the amount of fluid is restricted and the relative volume of fluid is the same order of magnitude as volumetric particle content in the slurry. Due to the conservation of mass, as gravity drives individual particles which are denser than the surrounding fluid downward, under the assumption that fracture walls are not moving apart, fluid inevitably flows upward and sideways at the same time. Therefore, from rationale and observations a conclusion is drawn that particles settling squishes fluids away and upward along the narrow fracture. The resulting fluid flow field in a narrow fracture significantly affects particle motion. Particle motion is governed by gravity, buoyancy and local fluid drag in different directions and results in erratic paths. Second, because of the lubrication effect, particle kinetic energy decreases when they approach each other or a wall. As a result, while particles slow down going toward each other in a narrow channel, they form clusters. It is observed that those clusters are temporary and may break apart or settle in the same formation for a while. There is no real attraction force present between particles and their subsequent movement is always governed by fluid drag or gravity. From the results, a strong tendency toward random channel clogging caused by particle agglomeration is related to increase in particle concentrations and fluid dynamic viscosity as a consequence of lubrication effect on particles contacts. Naturally, increased particle concentration increases the probability that particles can come in close contact, and increased dynamic viscosity adds to larger lubrication effect and particle kinetic energy damping.

Fig. 12

Micromechanics of proppant settling in w = 4 mm rough fracture for c _v = 0.1 (2.0 lbs/gal), where a shows particles settling in μ = 0.005 Pa s fluid after t = 1.5 s, b shows particles settling in μ = 0.01 Pa s fluid after t = 0.8 s and c shows particles settling in μ = 0.02 Pa s fluid after t = 0.2 s

Figures 12, 13, 14 shows model results of proppant settling in a 4 mm-wide fracture. Figure 12 shows micromechanics of proppant settling in w = 4 mm rough fracture for c _v = 0.1 (2.0 lbs/gal), where (a) shows particles settling in μ = 0.005 Pa s fluid after t = 1.5 s, (b) shows particles settling in μ = 0.01 Pa s fluid after t = 0.8 s and (c) shows particles settling in μ = 0.02 Pa s fluid after t = 0.2 s. Figure 13 shows micromechanics of proppant settling in w = 4 mm rough fracture for c _v = 0.2 (3.9 lbs/gal), where (a) shows particles settling in μ = 0.005 Pa s fluid after t = 1.2 s and (b) shows particles settling in μ = 0.01 Pa s fluid after t = 1.4 s. Figure 14 shows micromechanics of proppant settling in w = 4 mm rough fracture for c _v = 0.3 (5.9 lbs/gal), where (a) shows particles settling in μ = 0.005 Pa s fluid after t = 1.5 s, (b) shows particles settling in μ = 0.01 Pa s fluid after t = 1.75 s, (c) shows particles settling in μ = 0.02 Pa s fluid after t = 0.7 s and (d) shows particles settling in μ = 0.1 Pa s fluid after t = 1.0 s.

Fig. 13

Micromechanics of proppant settling in w = 4 mm rough fracture for c _v = 0.2 (3.9 lbs/gal), where a shows particles settling in μ = 0.005 Pa s fluid after t = 1.2 s and b shows particles settling in μ = 0.01 Pa s fluid after t = 1.4 s

Fig. 14

Micromechanics of proppant settling in w = 4 mm rough fracture for c _v = 0.3 (5.9 lbs/gal), where a shows particles settling in μ = 0.005 Pa s fluid after t = 1.5 s, b shows particles settling in μ = 0.01 Pa s fluid after t = 1.75 s, c shows particles settling in μ = 0.02 Pa s fluid after t = 0.7 s and d shows particles settling in μ = 0.1 Pa s fluid after t = 1.0 s

Proppant forms agglomerate randomly while settling, but the aggregation is more visible for higher particle concentrations and in higher fluid viscosities. Trends of particles to agglomerate are similar for w = 4 mm and w = 8 mm fractures (Figs. 15, 16, 17). Figure 15 shows micromechanics of proppant settling in w = 8 mm rough fracture for c _v = 0.1 (2.0 lbs/gal), where (a) shows particles settling in μ = 0.005 Pa s fluid after t = 1.2 s, (b) shows particles settling in μ = 0.01 Pa s fluid after t = 0.9 s and (c) shows particles settling in μ = 0.02 Pa s fluid after t = 0.9 s. Figure 16 shows micromechanics of proppant settling in w = 8 mm rough fracture for c _v = 0.2 (3.9 lbs/gal), after t = 0.6 s, where (a) shows particles settling in μ = 0.005 Pa s fluid, (b) shows particles settling in μ = 0.01 Pa s fluid and (c) shows particles settling in μ = 0.02 Pa s fluid. Figure 17 shows micromechanics of proppant settling in w = 8 mm rough fracture for c _v = 0.3 (5.9 lbs/gal), where (a) shows particles settling in μ = 0.005 Pa s fluid after t = 1.2 s and (b) shows particles settling in μ = 0.01 Pa s fluid after t = 1.6 s. As the particle concentration increases, proppant particles come in contact with each other and form clumps that settle together. Fluid lubrication at particle contacts affects particle interactions and bigger clumps are observed at higher fluid dynamic viscosities.

Fig. 15

Micromechanics of proppant settling in w = 8 mm rough fracture for c _v = 0.1 (2.0 lbs/gal), where a shows particles settling in μ = 0.005 Pa s fluid after t = 1.2 s, b shows particles settling in μ = 0.01 Pa s fluid after t = 0.9 s and c shows particles settling in μ = 0.02 Pa s fluid after t = 0.9 s

Fig. 16

Micromechanics of proppant settling in w = 8 mm rough fracture for c _v = 0.2 (3.9 lbs/gal), after t = 0.6 s, where a shows particles settling in μ = 0.005 Pa s fluid, b shows particles settling in μ = 0.01 Pa s fluid and c shows particles settling in μ = 0.02 Pa s fluid

Fig. 17

Micromechanics of proppant settling in w = 8 mm rough fracture for c _v = 0.3 (5.9 lbs/gal), where a shows particles settling in μ = 0.005 Pa s fluid after t = 1.2 s and b shows particles settling in μ = 0.01 Pa s fluid after t = 1.6 s

Figure 18 shows a comparison between average particle settling velocities in w = 2, 4 and 8 mm-wide channels and previous theoretical and experimental relationships. The ratio of average proppant settling velocity and the settling velocity for single particle in infinite media was calculated to follow the previously proposed theory (Eq. 8). Uniform proppant size is used with sand mesh 30 particles, with particle diameter d = 0.6 mm. The ratio of fracture width w and particle diameter size d is w/d = 3.3–13.3. The effect of fluid viscosity is examined and plots show proppant settling results in different fluids. The average velocity of proppant settling is a function of the proppant concentration. The average proppant settling velocity generally (without proppant agglomeration) decreases with increasing concentration. However, the extensive particle agglomeration that occurs at higher fluid viscosities and proppant concentration increases the proppant tendency to settle in clumps. The non-monotonic behavior of proppant settling velocity can be seen in Fig. 18 and shows that proppant concentration increases are caused by agglomeration. Larger agglomerates settle faster than a single particle, resulting in the average proppant settling velocity increase. It can also be seen that the increase in fluid viscosity contributes and accentuates non-monotonic behavior as particle concentration increases. For wide fracture opening, generally higher particle concentrations are possible, but due to the interactions with walls which direct particles toward each other causing agglomerations, average settling velocities are high compared to general experimental predictions (Fig. 18).

Fig. 18

Ratio of average particles settling velocities and one-particle velocity versus the channel width, w is the fracture width and d is the particle diameter

Figure 19 shows settling of high proppant concentration in a w = 8 mm-wide fracture. In this case, the proppant migrates very quickly toward the center of the fracture. Particle motion toward the fracture center can be explained by observing the fluid flow field. When fracture is wide enough relatively to the fracture surface asperity size, particle–wall interactions are less frequent at the beginning while the slurry is uniform. When particles start moving downward due to gravity, the fluid flow field adjacent to the fracture wall is stagnant due to friction and particles move away from the wall falling toward the center of the fracture. After the proppant has migrated toward the center of the fracture, average settling velocities increased. Increased velocity was observed in wider fractures for all fluid viscosities. Figure 20 shows results of relative particle velocities with different volumetric concentrations and fluid dynamic viscosities in narrow and wide channels. Results are compared to the velocity of single settling particle in a quiescent fluid. The relative settling velocities decrease with increasing particle concentration, which is in agreement with physical experiments, but in some cases higher concentrations cause proppant agglomerations and relatively faster settling of newly formed particle clumps. Figure 20 compares the effect of fluid viscosity on relative settling velocity for different widths of rough fracture with the JRC = 7 (Barton and de Quadros 1997). Particle size is the same in all cases with uniform proppant diameter of d = 0.6 mm (mesh 30). Comparing plots in Fig. 20a–c, it can be seen that as the fracture opening increases, the effect of fluid viscosity is more visible and lines that connect relative settling velocity for different fluid viscosities fall more apart. Narrow fractures inhibit agglomeration and, therefore, result in slower settling velocities.

Fig. 19

Migration of particles toward center in w = 8 mm-wide fracture with c _v = 0.4 (7.9 lbs/gal) in μ = 0.005 Pa s fluid after t = 0.6 s (a), t = 1.2 s (b) and t = 1.6 s (c)

Fig. 20

Relative particle velocities versus proppant volumetric concentration in different channel widths, relative particle velocity is the ratio of average particle velocities and single particle settling velocity

Figure 21 shows the effect of fracture width on proppant settling for different fluid dynamic viscosities. Figure 21a–c shows the wall effect for different proppant concentrations. The numerical results are compared with the previously proposed theory in Eq. (8) by Novotny (1977). The aim of comparing the settling velocity with a single particle settling is to remove the effect of fluid dynamic viscosity and particle size from the relations. However, model results show that proppant settling velocity depends on fluid viscosity (Fig. 21). In spite of dividing the average settling velocity with the single particle settling velocity value, the lines connecting model results differ for the same concentrations and d/w ratios. For example, for high fluid dynamic viscosities, the relationship proposed by Novotny (1977) underpredicts the settling velocity of the model, while for smaller fluid viscosities the values from the same relationship are larger than the model. The results in Fig. 12 indicate that particle interactions and particle–wall interaction affect the settling velocities.

Fig. 21

Effect of fracture width on proppant settling

Conclusions

Conditions for particle agglomerations as a consequence of settling in narrow rough fracture zone were investigated in this study. Relatively little effort has been devoted in the past to the understanding of proppant agglomeration in narrow and rough fracture zones. Therefore, an improved understanding of proppant behavior in rough and narrow fracture zone provided in this paper can assist reservoir engineers to recognize and mediate problems associated with fracture clogging during hydraulic fracturing. In particular, parameters such as slurry density, fluid viscosity and fracture width were evaluated toward the possibility of proppant agglomeration in 2, 4 and 8 mm-wide fractures. The objective of this study was to investigate: (1) whether or not the existing proppant settling relations can be used also for narrow fracture zones and (2) under which conditions they are valid. It was found that proppant settling relations cannot be used for narrow fracture zones when particle volumetric concentrations are higher than c _v = 0.2 and in fluid viscosities higher than μ = 0.005 Pa s.

The means to meet the stated objectives are via micromechanical observations of particle settling in a narrow channel using a modified DEM–CFD approach as implemented in PFC^2D with an improved particle contact model that accounts for lubrication effects on particle collisions. The combined fluid flow algorithm and improved particle collision contact model showed individual particles trajectories and fluid flow streamlines around individual particles and groups of particles. The micromechanical behavior of multi-particles dense systems settling in quiescent Newtonian fluids was studied from the perspective of non-linear particle collision behavior caused by a lubrication of thin fluid layer between colliding particles in DEM code. The significance of fluid lubrication force on particle motion approaching each other depends on lubrication parameters and the fluid drag, gravity and buoyancy, which also act on particles. Lubrication effect dramatically impacts particle motions when a combination of parameters is unfavorable. When low particle approaching kinetic energy is combined with high fluid dynamic viscosity and larger particle size, the lubrication force increases. It is possible, and as shown in this manuscript, that in certain conditions such as rough and narrow fracture, fluid lubrication slows down the velocities of approaching particles such that they start forming agglomerates and settle faster in clumps. As a result, in cases where lubrication dominates over other forces, agglomeration is more likely to occur. Migration of particles toward the center of the channel is more visible in wider channel with higher proppant concentrations. Fluid viscosity in a relatively narrow channel promotes particle agglomeration due to inter-particle and particle–wall collisions. Temporary clumps formed of several particles significantly increase the average proppant settling rates. Settling of proppant in higher fluid viscosities and lower fluid viscosities do not follow completely previously published relationships. Generally, it is observed that average proppant settling velocities decrease with increasing particle volumetric concentration in the slurry, which is consistent with previously published experimental results for wider slots with smooth walls. However, as particle–particle and particle–wall collisions get more frequent due to denser slurry, agglomerations are observed and more intense at higher fluid viscosities. One of the main contributions in this study is the observation of how the average proppant settling velocity increases instead of decreases at relatively high particle volumetric concentrations due to agglomeration processes.

In low fluid dynamic viscosity μ = 0.005 Pa s proppant settling velocities are still similar to published relations for 2, 4 and 8 mm-wide rough fractures. As the fluid dynamic viscosity increases to μ = 0.01–0.02 Pa s, the fracture width to particle diameter ratio plays more significant role than at lower fluid dynamic viscosity, and it can be concluded that as fracture width to particle diameter ratio increases, proppant settling velocities approach the values given by previously published relationships. At highest simulated fluid dynamic viscosity μ = 0.1 Pa s, the average proppant settling velocities are significantly higher in 2, 4 and 8 mm-wide rough fractures than published relationships, except for the case of 8 mm fracture and very low proppant volumetric concentrations (c _v = 0.1–0.2 or 2–3.9 lbs/gal). Particle volumetric concentration increase in combination with higher fluid viscosity, shown in some cases, trend toward inconsistent proppant settling velocity increase due to formation of irregular agglomerates. In such cases, proppant–fluid mixture cannot be considered as uniform slurry.

Fluid viscosity was found not to be a normalizing parameter and simple comparison with settling of single particle in quiescent fluid is not sufficient to predict the behavior of settling of concentrated slurries. At higher fluid viscosities and with increase in particle concentration, the proppant shows tendency to agglomerate and cause clumps. The micromechanics of particles settling reveals that during the sedimentation process, a combination of fluid Stokes drag upward flow and particles settling results in erratic particles paths and enhances collisions. Particles and fluid constantly interchange their positions in space in a vertical channel.

The elasto-hydrodynamic model for lubrication as implemented in PFC^2D is a step forward from the conventional DEM default, which mostly uses the Kelvin–Voigt particle contact model (spring–dashpot) because the lubrication model introduces the contact criteria for stick or elastic rebound of two spheres at their contact based on the surrounding fluid properties rather than only damped elastic rebound. In multi-particles system, if the balance of the lubrication force and the fluid approaching velocities is such that the particle slows down enough to near zero, the particles may stick next to each other and get trapped with the fluid and start agglomerating to form clumps.