Particle-fluid flow and transport within rough fractures

Proppant injection is an important part of a hydraulic fracturing programs in which fluid-particle slurry is injected into rock fractures. Injected particles are lodged between fracture surfaces during wall close-in thereby propping open the fracture, improving connectivity and production. This paper investigates behaviour of proppant particles within artificially generated rock fractures, providing insight into transport behavioural differences caused by realistic surface roughness. Better understanding of proppant behaviour within more realistic rough fracture conditions provides greater understanding of proppant transport as compared to past works where smooth walled fracture configurations were utilized. A clearer understanding is important in providing more accurate evaluation of realistic proppant flow and distribution and improving injection design. In this study a roughened surface, analogues to actual rock fracture surface, is artificially generated based on a rock surface’s fractal dimension and asperity height standard deviation. Computational representation of the rock surfaces and flow domain is generated. Resolved Discrete Element Method coupled with computational fluid dynamics (DEM-CFD) is implemented in this study to evaluate proppant particle transport behaviour within the fractures. This work highlights importance of considering fracture surface roughness in evaluating proppant flow and transport and more generally the impact of rough boundary conditions of particle-fluid systems.


Introduction
Proppant slurry inclusion is an integral part of many rock fracture system enhancement projects. The process consists of injecting natural or ceramic granular particles as a slurry into rock fracture systems under pressure. The proppant particles travel down the fracture systems and are lodged into place upon removal of pressure and fracture closure. To successfully design a proppant slurry treatment schedule, clear comprehension of the proppant particles' conveyance and settling behaviour is needed to accurately evaluate proppant distribution and travel. Majority of past studies of proppant behaviour have relied on simplified, smooth walled experimental and numerical configurations. Further, evaluation of settling based on simplified Stokes law behaviour has been used in fracture modelling software [1].
In-field fractures faces however contain rough surface asperities and present tortuous flow paths through them [2]. Further, accuracy of simple Stokes law-based behaviour for describing particle settling in fractures has contrasting findings. Early study of rough fracture face effects on particle settling rates was explored by Novotny [3]. Settling behaviour of particles in smooth and rough wall slot apparatuses were compared. It was concluded in that work that rough fracture walls presented no further effects beyond those otherwise observed due to typical wall effects. However, the empirical work of Liu and Sharma [4] showed greater hindered particle settling rates for rough-walled fracture configurations. Further, study by Luo and Tomac, [5] observed trapping of particles in narrow rough surface configurations. Numerical investigation utilizing computational fluid dynamics (CFD) by Huang et al. [6] concurred with the finding that settling rates decreased in rough fracture wall conditions. Huang et al. [6] and Liu and Sharma [4] also concluded that horizontal particle velocity, in the direction of flow, decreased due to rough wall conditions. Due to contrasting findings, further investigation into rough fracture wall effects is warranted. In this work, evaluation of a single particle's behaviour within smooth and rough walled fracture domains is studied. Numerical modelling via resolved, coupled Discrete Element Method with computational fluid dynamics (DEM-CFD) is employed providing clear particle and fluid behaviour data. Resolved DEM-CFD additionally allows for the ability to resolve fracture features to size smaller than the proppant particle while conforming the CFD mesh to the surface geometry, providing a truer representation of fracture surfaces relating to particles.

Particle settling
Particle terminal settling rate per Stoke law is evaluated as [7]: where is the terminal settling velocity of a single particle in an unbounded infinite fluid, ρ and ρ are the densities of the particle and the fluid, respectively, is gravity (i.e. 9.81 m/s), is the particle's diameter, and μ is the fluid's dynamic viscosity. However, Stokes law is only accurate for particles with low particle Reynolds number values, specifically below approximately 0.1, which is typically designated as the 'Stokes range'. The particle Reynolds number, , is defined as: with and equaling the fluid's and particle's velocities, respectively. Evaluation of terminal settling rates beyond the Stokes range of values, requires evaluation of drag on the particle using Newtonian drag formulation. The resulting equation for terminal settling velocity applying this formulation is [8]: with equaling the empirically determined drag coefficient. For this formulation, the drag coefficient in the Stokes range can be expressed as a linear trend by , = /24. Typical proppant values fall in the range of values of ≤ 100 [9]. Thus, evaluation outside of the Stokes range is needed. Figure 1

Wall Effects
Novotny [3] evaluated effects of smooth walls on proppant settling. Formulation to describe the hindered settling velocity is: where is the fracture opening width and is the hindered settling velocity due to fracture opening width wall effects. Liu and Sharma [4] also presented equations to describe the hindered settling velocity between smooth walls, as described by: with μ in poise. Liu and Sharma [4] additionally provided expressions for the modified horizontal particle velocity due to fracture opening width effects in smooth walled fractures: with equal to the modified horizontal particle velocity and ⃗ equal to the mean cross-sectional fluid velocity.

Resolved DEM-CFD
DEM can be originally attributed to the work of Cundall and Strack [11]. The model tracks and updates individual Lagrangian particles' motion in space and time. Particle collisions and particle assembly deformation are represented by overlapping of the individual particles or by particle overlap with walls in space. Motion and contact physics are resolved through Newton's second law for motion resulting from force application and transfer of energy between kinetic and potential during contact via force displacement law with spring interactions between interacting bodies. In this work, Hertzian non-linear contact model is used in collision interactions. Dissipation of energy from contact is provided by a dashpot interaction during contact, with damping dictated by the prescribed coefficient of restitution. DEM equations of motion for each particle, i, can be expressed as: where is the particle mass, is the particle velocity, is the particle rotational velocity, are the forces acting on the particle, is the torques acting on the particle, and is time. LIGGGHTS software package [12] is utilized for the DEM implementation in this work.
Computational fluid dynamics provides a numerical evaluation of fluid behavior by discretization of the Navier-Stokes equations over a solution domain. Equations describing the behavior of an incompressible Newtonian fluid are: where is the fluid pressure. The finite volume CFD software, OpenFOAM [13], is used in this study.
Coupling of the DEM and CFD occurs through the CFDEM®coupling toolkit, developed by Goniva et al. [14]. More specifically, the immersed boundary method (i.e. resolved DEM-CFD) is utilized [15,16]. In this method, particles are of size larger than CFD mesh cells. Fluid-particle interaction behavior is then resolved about each of the particles' surfaces. The resolved method can be contrasted with unresolved DEM-CFD where particles are of size smaller than the CFD cell mesh. Fluid-particle interactions are in that case are provided by empirical and theoretical coupling equations.
In DEM-CFD coupling, the fluid and solid phase simulations pass pertinent information calculated in each respective components' iteration about a shared time step. The exchanged information is used to determine the fluid influence on the particle and vice versa. Contribution of the particles' motion on the fluid phase in resolved DEM-CFD is accomplished by 1) solving the fluid behavior initially without considering the influence of the particles, giving a velocity field, , 2) imposing the motion of the particles at their respective locations in the CFD mesh, correcting velocity for those cells covered by a particle, giving a velocity field, , and 3) correcting the fluid velocity field to provide a divergence free solution [15,16]. The final step results in a velocity field described by: where is the divergence free, corrected velocity field, and is the correction factor.
Evaluating the continuity equation, Eq. 9, with the corrected divergence free velocity field and substituting in Eq. 10, results in: which is used in correcting the velocity field. The corrected pressure field is then evaluated as: Drag force from the fluid on particles is evaluated as: where is the fluid drag contribution on the particle, is designation of a respective CFD cell, represents the set of cells with solids (i.e. particle in cell region), and is the cell volume. Further details and derivations of the fluid-particle interactions can be found in the work of Hager [15] and Hager et al. [16]. Effects of buoyancy force, . , from hydrostatic pressure is also added as a particle force by: It is noted that the resolved DEM-CFD method has the limitation of high computational costs due to the large number of CFD cells required in resolving behavior around each particle. Hager and Hager et al. [15,16] note that at least 8 CFD cells across a particle's diameter are needed to achieve mesh independent results. Due to the high computational cost, simulations are limited in domain size and number of particles to allow for tractable numerical evaluation. Simulation domain size used in this work is discussed in subsection 3.3.

Synthetic Fracture Generation
Real rock surfaces can be well described with synthetic representation, using the surface's fractal dimension and root-mean-square roughness for a referenced length scale [2]. Fractal dimensions are used to describe a shape that is topologically defined between Euclidian spatial definitions, i.e. between 1-D, 2-D, and 3-D [17]. Further information regarding fractal dimensions can be found in the work of Mandelbrot [18]. Typical rock surface fractal dimension values are in the range of 2.0 to 2.5 per Brown [2]. Smaller values correspond to flatter surfaces where large correspond to rougher surfaces. An initial estimate of standard deviation based on fractal dimension can be determined as specified by Vuopio and Polla [19]: where σ is the standard deviation, is the sample's length, is the th bisection, taken as = 1, and is the equivalent profile fractal dimension which is evaluated as = − 1 [20].
Brown [2] proposed a method for producing computational representations of rock surfaces. The surfaces are constructed by considering a probability density function for heights of surface features, in this case Gaussian distribution, and an autocorrelation function or power spectrum describing spatial correlation of the surface heights. The slope of the power spectrum in a log-log plot is related to the fractal dimension. Using the fractal dimension along with the standard deviation of the surface roughness heights and sample size, a synthetic representations of rock surfaces can be generated. Implementation of this method is provided in SynFrac software [21]. For further details of this method the reader is directed to the works of Brown [2] and Ogilvie et al. [21].
In this work, surfaces are assumed to be matched (i.e. no offset in faces). To further confirm surface generation matching desired fractal dimension, surfaces were checked using the method described by Qiu et al. [22]. Standard deviation was adjusted if needed until surfaces measured to within 0.05 % of the desired fractal dimension. After surface generation, STL files were generated using a Python script and CFD mesh generated with Helyx-OS, a GUI for the SnappyHexMesh OpenFOAM software.

Model properties and conditions
Tables 1 and 2 show material properties and DEM-CFD boundary conditions. In addition to the flowing and settling between fracture walls, a particle settling in an unbounded infinite fluid case is also run. Boundary conditions for that case are specified in Table 2.  Two fracture domains are generated, one with fractal dimension of 2.0 (i.e. flat surface) and one with fractal dimension of 2.25 (i.e. rough surface). Domain for the simulation is generated by taking a subsection of the generated fracture surface equal to 10 mm (X-direction) x 3 mm (Z-direction). Aperture for both models is set to 4 mm (Y-direction). Fracture surface features are resolved to 100 micrometers. CFD cells are generated as 50 micrometer sided cubic cells. In order to allow for generated particles to achieve terminal settling or sustained flow velocities prior to entering the fracture configurations, a 5 mm long 'initiation' region is extended from the -X face of the fracture domain. Illustration of this is shown in Figure 2. Particle evaluated behavior is considered after this region, once the particle enters the fracture section of the domain. Further, DEM domain is shortened by 1mm relative to the +X face of the CFD fracture domain (i.e. DEM domain equaled 9 mm in the X direction + initiation region of 5 mm). The domain shortening serves for avoiding boundary effects at the end of the simulation domain.
Single particles are initiated within the initiation region at 1 mm from the -X face. Y and Z positions are initiated as the center of area the Y/Z face. To form the prescribed inlet velocity for the neutrally buoyant cases, an initial CFD only case is run with a small, uniform flow in the -X direction, generated from the domains' +X faces.
Resulting developed velocity at the -X face is then mapped, scaled, and adjusted to point in the +X direction for use in the simulations. An inlet mean fluid flow rate of 0.33 m/s is used in the neutrally buoyant cases. Additionally, an initially uniform velocity of 0.33 m/s within the CFD mesh is prescribed to allow for a developed flow within the fracture prior to the particle entering the evaluation region.

Results
Conveyance results for the neutrally buoyant cases compared with the evaluations of Liu and Sharma [4] are shown in Figure 3. Single particle conveyance in the neutrally buoyant, smooth fracture wall model configuration is similar to that expected by the formulation developed by Liu and Sharma [4]. However, the rough surfaced configuration simulation results are in stark contrast to the smooth wall results. From what can be seen in Figure 3, the single particle within the rough fracture in fact has higher velocity in the direction of flow. This is in contrast with the findings of Liu and Sharma [4] and Huang et al. [6] discussed in the Introduction, where it was concluded that rough surfaces resulted reduced horizontal velocity for particles.  What can be observed in the neutrally buoyant, flowing cases is that preferred regions of flow develop within the rough fracture leading to higher velocity regions compared to the smooth fracture case. Figure 4 shows a cross section through the fracture models with fluid velocity magnitudes for both smooth and rough configurations. Regions with high velocity develop within the rough fracture with magnitudes up to approximately 4 times the maximum velocities in the smooth model. The higher velocity can be attributed to the variance in cross sectional area normal to the direction of flow in the rough fracture. The smaller area leads to the higher velocity. The particle enveloped in these regions experiences an increased conveyance rate, resulting in faster travel compared to the smooth model.
Settling results compared to the evaluations of Liu and Sharma [4], Novotny [3], and settling per Eq. 3 are shown in Figure 5. The smooth walled case has fairly similar results to the evaluations of Liu and Sharma [4] and Novotny [3]. Further since wall effects are small for this case, the settling is alike to the unbounded particle in an infinite fluid modelled case and the Newtonian drag formulation from Eq. 3.  Settling within the rough fracture however is initially slightly slower due to slightly higher wall effects in the initiation region's shape. During portions of the particle's settling within the rough fracture, settling velocity approaches that of the smooth wall fracture but decreases when the particle encounters a wall asperity and collides. Eventually, the particle comes to rest at a position approximately half-way through the simulation domain where it becomes arrested on a rough asperity feature. The findings agree well with previous observations presented in the Introduction, where settling velocity was more greatly hindered in rough fractures.

Conclusions
Results from this study indicate fracture wall roughness can have a significant effect on particle settling and conveyance. The findings include: 1. Rough fracture walls do not always attenuate particle conveyance. Higher conveyance rate in a rough fracture compared to a smooth-walled fracture was observed in this study. 2. Preferential regions of high fluid velocity developed in the rough fracture configuration for the case evaluated. This resulted in the higher particle conveyance rate. 3. Settling rate of proppant in rough fracture conditions was found to be attenuated, and in this case eventually arrested, due to interference with the rough fracture walls.