Regular Article
Numerical study on particle transport and deposition in rough fractures
^{1}
School of Petroleum Engineering, China University of Petroleum (East China), Qingdao
266580, China
^{2}
College of New Energy, China University of Petroleum (East China), Qingdao
266580, China
^{3}
School of Civil and Environmental Engineering, University of Science and Technology Beijing, Beijing
100083, China
^{4}
Department of Oilfield Exploration & Development, SINOPEC, Beijing
100029, China
^{*} Corresponding authors: lgong@upc.edu.cn,
^{*}
liyang@sinopec.com
Received:
7
November
2019
Accepted:
3
March
2020
The transport and deposition of particulate materials through fractures is widely involved in environmental engineering and resource development engineering. A 3D Computational Fluid DynamicsDiscrete Element Method (CFDDEM) coupling method was used to investigate the particle and fluid flow. The Gauss Model was applied to construct the rough surfaces. First, the numerical results were compared with the previous results and reasonable agreements were obtained. Second, the results indicated a novel flow pattern of particles in rough fractures. Then, a comprehensive particle sedimentary analysis indicated that the deposition distance of particles was inversely proportional to the particle size and density ratio. In addition, the particle deposition rates were increased by the mean roughness and there was an uneven sediment distribution impacted by roughness. Reasons for this uneven sediment distribution were analyzed in detail. Moreover, the bridge plugs of particles considering the closure of fractures were simulated as well. A part of particulate materials would be filtered at the inlet due to size effect and the transport distance of entered particles decreased significantly when the particle was large. A critical particle radius (R < 0.27 mm) that can flow through closure fracture in this work was found. This work can provide a clear insight into the migration and deposition characteristics of particles in the rough fractures underground.
© X. Wang et al., published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
: Velocity of the particle i, m/s
: External forces for per unit fluid, N
Greek symbols
α _{f} : Fluid volume fraction
β : Drag coefficient of the particle
ρ _{f} : Fluid density of phase f, kg/m^{3}
δ _{ n } : Overlap distance between particles i and j, m
: Angular velocity of the particle i, rad/s
1 Introduction
With the development of unconventional reservoirs, the flow mechanism in pores and fractures have attached the attention of scholars (Li and Voskov, 2018; Sun et al., 2017; Yang et al., 2019c; Zhu et al., 2016). Numerous publications have been reported focusing on porescale numerical investigations of the gas flow in porous media and many of the results are helpful to study reservoir permeability (Song et al., 2016; Wang et al., 2019a; Zhang et al., 2019; Zhao et al., 2020). Such as Tianyi Zhao et al. indicated gas flow characteristics in the shale matrix by Lattice Boltzmann Method (LBM) and the effects of slippage, adsorption layer, residual water saturation and surface area all were considered (Zhao et al., 2018). Especially the dynamic sorption was included in LBM schemes which are robust and efficient by Zhang and Sun (2019). However, except for the fluid flow mechanism the permeability damage results from the plugging of particulate materials in pore or fractures are studied as well because the loss of the particulate material also is a thorny issue for underground engineering and resource development engineering (Lijun et al., 2019; Oliveira Jr. et al., 2012; Wang et al., 2019c). The improper operation or engineering accidents will cause particles to leak into fractures underground and then result in significant economic losses and damage to the reservoirs (Sun et al., 2017; Yang et al., 2019b). As the basis of this particulate pollutant control, it is necessary to investigate particle flow and deposition characteristics in fractures. In addition, particle migration and deposition in fractures are common challenges in many other environmental and industrial projects such as the transport of proppant particles in fracturing engineering, treatment and plugging of slurry losses and structure optimization of granular flowing experiments (Haddadi and Di Carlo, 2017; Patankar et al., 2000; Yan et al., 2019; Zeng et al., 2019).
There are many efforts have been made to prevent this hazard and many research groups investigated the particle flow underground by experimental and numerical methods (Haddadi and Di Carlo, 2017; Hemmati and Rafee, 2018; Li et al., 2020; Yang et al., 2019a). It was indicated that the fractures were the important transport channels of particulate materials (Wang et al., 2019b; Yan et al., 2016). The previous results also showed that the characteristics of particle flow were complicated and involved the inherent structure complexities of fracture surfaces, multiple properties of particles, and the difficulties of coupling calculation (Han and Cundall, 2017; Wang et al., 2017; Zhang et al., 2017). However, not many field measurements and visual analysis for particle transport and deposition through fractures underground have been reported, since it is infeasible to detect directly in realtime. Fortunately, the numerical method has been considered as an effective method to further study these complex processes that cannot be realized by existing experimental techniques. Generally, the simulation strategies for particle flow are divided into the EulerianEulerian coupling method and the EulerianLagrangian coupling method (Basu et al., 2015; Han et al., 2007b). A EulerianEulerian coupling method is a macroscopic approach treating fluid and particle as continuum media which is not suitable for tracking movement characteristics. In contrast, the particle is regarded as a discrete phase in the EulerianLagrangian coupling strategy. In this strategy, particle flow is computed by the Discrete Element Method (DEM), fluid flow is computed by the Lattice Boltzmann Method (LBM) or Finite Volume Method (FVM). In addition, the Immersed Boundary or Moving Immersed Boundary Method is introduced in the resolved coupling method and drag force formula is used in unresolved coupling conditions (Han et al., 2007a; Kuruneru et al., 2017; Zeng et al., 2016; Zhang et al., 2017). Some semiresolved Computational Fluid DynamicsDiscrete Element Method (CFDDEM) approaches with advanced particle volume fraction calculation methods are introduced recently as well (Cheng et al., 2018; Wang et al., 2019d; Zhang et al., 2018). The CFDDEM coupling method is adopted in this manuscript because it has both computational efficiency and particle computing accuracy.
In the field condition, there are many factors that should be taken into consideration for particle flow simulation. The collisions between particles, the nonuniform sizes and densities of mixed particles, and the different amounts of every kind of particle all would affect the distribution of deposited particles in fractures (Huang et al., 2018; Zou et al., 2015). In relevant work, DEM was proposed as a suitable method to fully consider these factors (Wu et al., 2017; Zeng et al., 2016, 2018; Zhou et al., 2018). Compared with other particle tracking methods, such as DPM (Discrete Phase Model), the DEM has more advantages for its richness of contact force model and its accuracy of computing motion of particles with multi sizes, shapes, and densities. Especially, the DEM can compute particle deposition reasonably from the perspective of energy losses rather than using rigid boundary judgment conditions.
Random and complex morphology is one inherent feature of underground fractures, especially the surface roughness is increased due to acid corrosion or treatment when the calciteenrichment area is dissolved in priority (Li et al., 2016; Ogilvie et al., 2006; Yang et al., 2018). Because of the fluctuation characteristics of fracture surface accord with normal distribution, the Gauss Model has been widely accepted for constructing surface profiles that more tally with the actual situation (Guo et al., 2015). The roughness has an important influence on the flow and deposition pattern of dispersed particles in fracture and extensive relative study has been done in experiments and numerical simulation aspects. On the one hand, the visual experiments were conducted. It was found that particle size, mass loading, and the roughness degree play roles on the particle velocity profiles (Sommerfeld and Kussin, 2004). The uneven distribution of particles in the rough fracture was studied in experiments as well (Huang et al., 2018; Tan et al., 2017; Zou et al., 2015). On the other hand, various numerical simulation works obtained more detailed information on particle deposition characteristics, such as the effects of the micro shape of the rough channel (Dorari et al., 2015; Hemmati and Rafee, 2018; Hong et al., 2018). However, the compressive numerical investigation is still worth implementing.
The aim of this work is to investigate the deposition characteristics of particulate materials in rough fractures underground. There are following three sections in this manuscript. Section 2 presents the details of the CFDDEM coupling model for the particlewater flow and Gauss model for surface roughness. Section 3 validates the model compared with previous work and discusses the results of this manuscript. The influence of multiple factors on particle transport and deposition is analyzed comprehensively such as the surface roughness, fracture closure, particle injection height, particle density, particle size, and particle number difference of mixed particles. The conclusions are summarized in Section 4 eventually.
2 Methodologies
2.1 CFDDEM coupling model
The particle and water flow in the rough fracture are investigated by the CFDDEM coupling method. The CFD equation set is computed in a commercial software ANYSY Fluent to obtain the information of fluid field containing particles. The movement of particles is computed by the Discrete Element Method in the EDEM software. The coupling between these two phases is implemented by a User Defined Functions (UDFs) code in Fluent. Through this UDF, the numerical simulation results between two software can be shared and invoked in realtime. The whole flowchart about the coupling method is presented in Figure 1.
Fig. 1 Flowchart about the CFDDEM coupling method. 
The movement of particles is computed by the Discrete Element Method and the interaction between contacted particles is characterized by the overlap as presented in Figure 2. The particle translation and rotation are computed by equations (1)–(3):(1) (2) (3)where represents the gravity; represents the translational velocity of the particle i; represents the angular velocity of the particle i; represents the total contact force on particle i; represents the normal direction between the center of particle i and j; represents the particle rotational inertia; represents the vector sum of the torques on the particle i.
Fig. 2 Softsphere overlap model. 
The governing equations of fluid flow based on NavierStokes description are:(4) (5) (6)where α _{f} represents the fluid volume fraction; ρ _{f} represents the density of fluid; represents the velocity of fluid; represents the viscous stress tensor; represents the source term of the external forces; represents the pressure gradient force; represents the drag force and represents the lift force.
The fluid forces on the particles are calculated in the coupled CFDDEM model and the detailed fluid force calculation formulas are presented in the Appendix.
2.2 Construction of fracture with rough surfaces
Roughness often was ignored by the traditional fracture model using two parallel plates due to the complexity of geometry but it influenced the micro flow characteristics in fracture a lot. Most of the morphology information of rough surfaces such as peak size and distributions of roughness peaks are consistent with normal distribution, thus rough surface can be computed and reconstructed by Gauss model statistically. In other words, these peaks map information is concluded in one mathematical model equivalently for fracture roughness description.
In the Gauss (GS) Model a Fourier analysis is used for generating the surface roughness z (x, y) in equations (7) and (8):(7)where x and y represent the number of nodes in the direction X and Y of fracture geometry region respectively; l _{ x } is the surface length in the direction X and l _{ y } is the surface length in the direction Y; cl _{ x } is the correlation length in the direction X and cl _{ x } is the correlation length in the Y direction; z represents the height of every rough peak; the ift is the function that returns the twodimensional inverse Discrete Fourier Transform (DFT) computed by the Fast Fourier Transform (FFT) algorithm and F _{g} is the Gaussian filter:(8)
The geometry information of rough fracture (Fig. 3), the properties of the fluid and particle used in this work are summarized in Table 1.
Fig. 3 The rough fracture geometry and the investigated region. 
Simulation settings and parameters.
2.3 Validation
To verify the reliability of this model to simulate the dense particle and fluid twophase flow in rough fracture the profile of packed particles in fracture simulated by the model in this work was compared with the reference result (Zeng et al., 2016). The simulation setting and parameters of cases for comparison are presented in Figure 4 and Table 2. The deposited particles are in red color and the fluid is in blue color in Figure 4. The particle dune in the reference is outlined by the dotted line in black. The height and shape of the particle dune in both results are constant with a slight difference.
Fig. 4 The particle dune profile in the Eulerian view. The dune outline of the reference result (Zeng et al., 2016) was outlined by a dotted black line. 
The properties of particle and fluid.
In addition, to validate the accuracy of particle tracking of the CFDDEM coupling model the particle flow through the channel was simulated and compared with previous work (Haddadi and Di Carlo, 2017). A group of particles is injected from the inlet and flow out from the outlet of the micro channel as shown in Figure 5. The trajectory of a representative particle is tracked and compared with that in reference. Thus, the particle tracking results computed by the CFDDEM method have a reasonable match with the reference results.
Fig. 5 The particle trajectory comparison between the simulation and reference results (Haddadi and Di Carlo, 2017). 
3 Results and discussion
3.1 Effects of the rough bottom
The influence of fracture surface features on particle flow patterns and deposition characteristics has attracted the attention of many scholars recently. In this paper, some novel migration characteristics of suspending particles (density ratio was 1) in rough fracture were investigated.(9)
The particle trajectories in the rough fracture (RP1 and RP2) are different from those in the smooth parallel plates (S) as shown in Figure 6. The vibration amplitude on x, y and z directions (AX, AY and AZ) of the particle trajectory was computed to describe the particle flowing fluctuation. It was indicated that a certain particle moves forward in the threedimensional channel and experiences up and down motion during swaying around. The overall movement is a spiral trajectory that cannot be visualized by both twodimensional simulations and threedimensional simulations without considering rough surfaces before. The velocity vector of particles shown in Figure 7 indicates that the particles would wind around the peak of the rough bottom, while the other particles with less kinetic energy deposit in the valley of the rough bottom after a period of spiraling motion. In addition to the effects on the particle trajectories, the uneven distribution of deposited particles results from the surface roughness was found as well. The roughness increased the probability of particleparticle and particlewall collisions. After a series of collisions, the kinetic energy of some particles was reduced to a level which was not enough to escape vortexes. Eventually, particles deposit (Fig. 7) and the uneven distribution of particles is formed that more particles tend to deposit on the valleys of the rough bottom.(10)
Fig. 6 (a) Comparison of particle movement characteristics in rough and smooth fracture. (b)–(d) Particle trajectory on the XZ section (view along the direction of movement), YZ section (side view), and XY section (top view). 
Fig. 7 The velocity vector of particles near the rough bottom and the decrease of kinetic energy due to collisions. 
Furthermore, the effect of the mean surface roughness on the particle deposition ratio was computed. The mean roughness of fracture surfaces in this work is 0.264 mm, 0.539 mm, 0.787 mm and 1.080 mm, respectively. The length (l) of the investigated area is 80 mm and the particle deposition rate is defined in equation (11). It was found that fluid maximum velocity and deposition rate increased as mean roughness grown as shown in Figure 8.(11)
Fig. 8 The velocity of fluid and particle deposition rates under different roughness. 
3.2 Transport and deposition of mixed particles
Considering the actual mixture condition of particulate materials, the effects of density ratio and size on particle transport and deposition were investigated.
3.2.1 Effects of density ratios
The particles with various density ratios were injected into the fracture in the same way and each kind of particle has the same size and amount. The density ratio is 1.5, 2.5, 3.5 and 4.5 respectively. The particle radius is 0.25 mm and the number of every kind of particle is 100. It can be found that particles experienced separation due to the particles with a larger density ratio had smaller horizontal deposition distance. Furthermore, the relationship between the maximum deposition distance and density ratio is consistent with the parabola law as shown in Figure 9. Additional simulations with the density ratio of 2, 3 and 4 are implemented for validating this interpolation pattern and these data (in red) is fitting well with the dotted line in Figure 9.
Fig. 9 Maximum deposition distance of particles with various density ratios. 
According to the law of deposition distance for particles with different density in Figure 9, the amount adjustment or difference of each kind of particle would change the eventual distribution of deposited particles in fracture. The effects of the socalled amount difference were investigated in this work. The detailed particle number of every kind of density is listed in Table 3 and Case 1 is regarded as the reference case that all kinds of particles have the same amount. The deposition area was divided into area A, B and C and the number of deposited particles in three areas were analyzed.
Amount distribution of four kinds of particles.
As shown in Figure 10c, there are much fewer lowdensity particles deposited at area C in Case 2 compared with reference Case 1. Due to the blue particles is the only kind in area C, the amount difference leads to a reduction in the total number in this area as well. In Figure 10b, Case 3 has fewer particles in green and Case 4 has fewer particles in red which both lead to sedimentation reduce in middle section B. Case 5 has fewer particles in the near entry area A than other four cases as shown as Figure 10d due to the dominance of the dense density particle in black.
Fig. 10 Particle amount distribution in areas A, B, and C. 
3.2.2 Effects of particle sizes
Because it is hard to keep the radius of particle the same and many kinds of sizes were always mixed during the transport of particulate materials in fractures. The uneven size becomes one significant factor that will affect the transport process and distribution results in a fracture. The effects of particle sizes and their injection heights were investigated in this paper.
The particles were injected randomly from the inlet area, whose size had a normal distribution (Fig. 11). The horizontal migration distances of every kind of particle are recorded in Figure 12 and a previous separation phenomenon is found. It was indicated that the horizontal transport distance is inversely proportional to the particle diameter but it is proportional to the injection heights.
Fig. 11 Particle radius distribution. 
Fig. 12 The influence of particle diameter and injection position on the deposition distance of particles (the three rows of particles in orange were observed). 
The row of particles with various radius injected from high position experienced sorting. Smaller particles broke away from the cluster of particles gradually during the sedimentation and had farther deposition distance. Particles from medium and low injection height had a similar law of sorting during the settlement process. However, because of the shorter settling time from lower injection height, the separation was not obvious. Besides, big particles injected from a higher position and small particles from the lower position would land on the bottom at the same time and had close horizontal migration distance. Consequently, the particles from different initial layers were mixed together after settlement.
3.3 Effect of rough sides
Particle migration and deposition characteristics in fracture with rough top and bottom walls were investigated and potential influence from uneven density, amount, size and injection height was analyzed in detail. However, tortuosity of fracture, in other words, the rough sides should be considered as well. It is a significant feature of fractures and has effects on the flowing and deposition patterns of particulate materials. The particle migration and deposition process through fracture with rough surface and sides were simulated. It was predicted that a part of particles would deposit in the corner of rough sides in priority.
On the one hand, the zigzag corner of rough sides affected the fluid flow such as inducing some local vortices. The motion of inertia particles was closely related to the fluid flow characteristics and a part of particles would go through these corners. Under the effects of fluid forces and collisions, the kinetic energy of particles experiences a fluctuation and decreases to zero eventually. The particles deposited when kinetic energy reduced to a level that was not enough to escape vortexes. To analyze the particle movement in this fracture a single particle trajectory from the top view was computed. As presented in Figure 13, in the region inside the blue dashed box corresponding to region 1, the particle dives along the channel under the drag of fluid; in the green region corresponding to region 2, particle still has enough kinetic energy to escape the sharp corner of the rough side in a parabolic path; in the region 3 corresponding to the region with orange color, particle swirls in the vortices of the fluid then it still moves forward; in region 4 that inside the red dashed box, the kinetic energy of the particle has been consumed a lot after collisions and the particle has a tendency to stop at the corner. Particles finally reach the region 5 indicated by the purple box where kinetic energy consumption further and the particle velocity gradually approached 0 after swirling and collisions. The particle deposit in region 5 eventually in Figure 13.
Fig. 13 (a) The trajectory of certain selected particles. (b) The decrease process of average kinetic energy. 
Moreover, the impacts of rough bottom and sides on the particles were investigated compressively. The roughness induced local fluid vortices, increased particle collisions, and kinetic energy dissipation. Driven by the fluid, most of the particles flowed into the valleys of the bottom and zigzag corners of sides and lost kinetic energy after collisions. Consequently, a number of particles tended to deposit in the valley and zigzag corners in priority. This uneven sedimentary distribution can be observed clearly from simulation results in Figure 14.
Fig. 14 Distribution of particles at the bottom (top view). 
3.4 Effect of closure of the rough fracture
When considered the closure of rough fracture, the openings of different parts of the inlet are different and there are some contact spots due to the initial apertures between rough surfaces are uneven. A part of particulate materials would be filtered due to size effect. For the particles that entered fracture, some flowed around the rough protruding peaks while others formed the bridge plugging. The fingerlike particle distribution was formed in rough fractures rather than the particle layers. It was found that both the number of entering particles and transport distance decrease as particle size increases in Figure 15. When the particle is rather small (R < 0.27 mm in this work), the closure just has a weak impact on the number of particles entered fractures. The impact becomes significant for cases with bigger particles (R > 0.27 mm in this work) and particle transport distance is short due to the bridge plugs.
Fig. 15 The plugged particle distribution in closed fractures. 
4 Conclusion
As the basis of particulate pollutant control, it is necessary to investigate particle flow and deposition characteristics in fractures. The CFDDEM coupling numerical method was used to compute the migration and deposition process of particles in rough fractures. In addition, the Gauss Model was applied to construct the rough walls of fracture. These numerical simulation results are summarized as followed:

The simulated trajectories showed that the particles within the rough fracture held different flow patterns from those in smooth condition. The particles experienced up and down motion, swaying around and overall trajectories were spiral.

A comprehensive sedimentary analysis indicated that the deposition distance of particles was inversely proportional to the particle size and density ratio. The separation and deposition of mixed particles were simulated efficiently as well.

It was found that both rough surfaces and sides of fracture have impacts on the deposition and distribution of particles. Bottom mean roughness increased particle deposition rates. Particles tended to settle in the valleys of the rough bottom and the zigzag corners of rough sides. Possible reasons for this uneven sediment distribution were analyzed: firstly, local fluid vortices were induced by surface roughness; then particles flow into the valleys and corners driven by fluid and lost kinetic energy gradually after collisions; consequently, most of the particles deposited in the valley and corners in priority.

The bridge plugs in closed rough fractures were simulated efficiently. A part of the particle would be filtered and stopped form entering fracture due to size effect and the transport distance decreased significantly when the particle was large. A critical particle radius (R < 0.27 mm) that can flow through closure fracture in this work was found.
Overall, the trajectory of sparse particles and the distribution of deposited particles were influenced by wall roughness while the law of particle separation was more related to the particle parameters. This work gave a clear insight into the transport and deposition process of particles in fractures which could provide guides for hazardous and pollutant particle treatments underground. To investigate the particle transport and deposition in the environment and engineering further, the improved related numerical model and experiments are still necessary. Such as the model of this work is in the labscale and more efforts are needed to achieve the simulation in the field scale of hydraulic fracturing. In addition, the fracture propagation and the leakoff are not considered in this study and more effects should be done in the future.
Acknowledgments
This work was supported by the Major Projects of the National Science and Technology (No. 2016ZX05061, No. 2016ZX05060, and No. 2017ZX05009), the National Natural Science Foundation of China (No. 51936001, No. 51674208, No. 51504276, No. 51504277, and No. 51904321), the Shandong Provincial Natural Science Foundation (ZR2019JQ21), the Fundamental Research Funds for the Central Universities (No. 17CX02008A, No. 17CX05003, No. 18CX02031A, No. 18CX07012A, and No. 19CX05002A), Key Research and Development Plan of Shandong Province (2018GSF116009).
Appendix
The drag force is given by equations (12)–(14):(12) (13) (14)
The lift force is given by equation (15):(15)
The pressure gradient force is given by equation (16):(16)where μ represents the fluid viscosity, V _{p} represents the particle volume, C _{L} represents particle the lift force coefficient, d _{p} represents the particle diameter, represents the particle velocity, represents the fluid pressure gradient, β is the drag coefficient and C _{D} is defined as:(17)
References
 Basu D., Das K., Smart K., Ofoegbu G. (2015) Comparison of EulerianGranular and discrete element models for simulation of proppant flows in fractured reservoirs, in ASME 2015 International Mechanical Engineering Congress and Exposition, American Society of Mechanical Engineers Digital Collection. [Google Scholar]
 Cheng K., Wang Y., Yang Q. (2018) A semiresolved CFDDEM model for seepageinduced fine particle migration in gapgraded soils, Comput. Geotech. 100, 30–51. [Google Scholar]
 Dorari E., SaffarAvval M., Mansoori Z. (2015) Numerical simulation of gas flow and heat transfer in a rough microchannel using the lattice Boltzmann method, Phys. Rev. E 92, 6, 063034. [Google Scholar]
 Guo L., Xu H., Gong L. (2015) Influence of wall roughness models on fluid flow and heat transfer in microchannels, Appl. Therm. Eng. 84, 399–408. [Google Scholar]
 Haddadi H., Di Carlo D. (2017) Inertial flow of a dilute suspension over cavities in a microchannel, J. Fluid Mech. 811, 436–467. [Google Scholar]
 Han K., Feng Y., Owen D. (2007a) Coupled lattice Boltzmann and discrete element modelling of fluid–particle interaction problems, Comput. Struct. 85, 11–14, 1080–1088. [Google Scholar]
 Han K., Feng Y., Owen D. (2007b) Numerical simulations of irregular particle transport in turbulent flows using coupled LBMDEM, Comput. Model. Eng. Sci. 18, 2, 87. [Google Scholar]
 Han Y., Cundall P. (2017) Verification of twodimensional LBMDEM coupling approach and its application in modeling episodic sand production in borehole, Petroleum 3, 2, 179–189. [CrossRef] [Google Scholar]
 Hemmati Y., Rafee R. (2018) Effects of the shape and height of artificial 2D roughness elements on deposition of nano and microparticles in the turbulent gas flow inside a horizontal channel, J. Aerosol Sci. 122, 45–58. [Google Scholar]
 Hong W., Wang X., Zheng J. (2018) Numerical study on particle deposition in rough channels with different structure parameters of rough elements, Adv. Powder Technol. 29, 11, 2895–2903. [Google Scholar]
 Huang H., Babadagli T., Andy Li H., Develi K. (2018) Visual analysis on the effects of fracturesurface characteristics and rock type on proppant transport in vertical fractures, in SPE Hydraulic Fracturing Technology Conference and Exhibition, Society of Petroleum Engineers. [Google Scholar]
 Kuruneru S.T., Sauret E., Saha S.C., Gu Y.T. (2017) A coupled finite volume & discrete element method to examine particulate foulant transport in metal foam heat exchangers, International Journal of Heat and Mass Transfer 115, 43–61. [Google Scholar]
 Li J., Qiu Z., Zhong H., Zhao X., Huang W. (2020) Coupled CFDDEM analysis of parameters on bridging in the fracture during lost circulation, J. Petrol. Sci. Eng. 184, 106501. [CrossRef] [Google Scholar]
 Li L., Voskov D. (2018) Multilevel discrete fracture model for carbonate reservoirs, in ECMOR XVI – 16th European Conference on the Mathematics of Oil Recovery, European Association of Geoscientists & Engineers, pp. 1–17. [Google Scholar]
 Li N., Dai J., Li J., Bai F., Liu P., Luo Z. (2016) Application status and research progress of shale reservoirs acid treatment technology, Nat. Gas Ind. B 3, 2, 165–172. [CrossRef] [Google Scholar]
 Lijun L., Jun Y., Hai S., Zhaoqin H., Xia Y., Longlong L. (2019) Compositional modeling of shale condensate gas flow with multiple transport mechanisms, J. Petrol. Sci. Eng. 172, 1186–1201. [CrossRef] [Google Scholar]
 Ogilvie S.R., Isakov E., Glover P.W. (2006) Fluid flow through rough fractures in rocks. II: A new matching model for rough rock fractures, Earth Planet. Sci. Lett. 241, 3–4, 454–465. [Google Scholar]
 Oliveira Jr J.A.A., Zago J., Fontes C., Waldmann A.T.A., Martins A.L. (2012) Modeling drilling fluid losses in fractured reservoirs, in SPE Latin America and Caribbean Petroleum Engineering Conference, Society of Petroleum Engineers. [Google Scholar]
 Patankar N.A., Singh P., Joseph D.D., Glowinski R., Pan T.W. (2000) A new formulation of the distributed Lagrange multiplier/fictitious domain method for particulate flows, Int. J. Multiph. Flow 26, 9, 1509–1524. [CrossRef] [Google Scholar]
 Sommerfeld M., Kussin J. (2004) Wall roughness effects on pneumatic conveying of spherical particles in a narrow horizontal channel, Powder Technol. 142, 2–3, 180–192. [Google Scholar]
 Song W., Yao J., Li Y., Sun H., Zhang L., Yang Y., Sui H. (2016) Apparent gas permeability in an organicrich shale reservoir, Fuel 181, 973–984. [CrossRef] [Google Scholar]
 Sun H., Yao J., Cao Y.C., Fan D.Y., Zhang L. (2017) Characterization of gas transport behaviors in shale gas and tight gas reservoirs by digital rock analysis, Int. J. Heat Mass Transfer 104, 227–239. [CrossRef] [Google Scholar]
 Tan Y., Pan Z., Liu J., Wu Y., Haque A., Connell L.D. (2017) Experimental study of permeability and its anisotropy for shale fracture supported with proppant, J. Nat. Gas Sci. Eng. 44, 250–264. [Google Scholar]
 Wang D., Yao J., Chen Z., Song W., Sun H. (2019a) Imagebased corescale real gas apparent permeability from porescale experimental data in shale reservoirs, Fuel 254, 115596. [CrossRef] [Google Scholar]
 Wang M., Feng Y., Pande G., Chan A., Zuo W. (2017) Numerical modelling of fluidinduced soil erosion in granular filters using a coupled bonded particle lattice Boltzmann method, Comput. Geotech. 82, 134–143. [Google Scholar]
 Wang X., Yao J., Gong L., Sun H., Yang Y., Zhang L., Li Y., Liu W. (2019b) Numerical simulations of proppant deposition and transport characteristics in hydraulic fractures and fracture networks, J. Petrol. Sci. Eng. 183, 106401. [CrossRef] [Google Scholar]
 Wang Y., Geng F., Yang S., Jing H., Meng B. (2019c) Numerical simulation of particle migration from crushed sandstones during groundwater inrush, J. Hazard. Mater. 362, 327–335. [Google Scholar]
 Wang Z., Teng Y., Liu M. (2019d) A semiresolved CFD–DEM approach for particulate flows with kernel based approximation and Hilbert curve based searching strategy, J. Comput. Phys. 384, 151–169. [Google Scholar]
 Wu H., Gui N., Yang X., Tu J., Jiang S. (2017) Numerical simulation of heat transfer in packed pebble beds: CFDDEM coupled with particle thermal radiation, Int. J. Heat Mass. Trans. 110, 393–405. [CrossRef] [Google Scholar]
 Yan X., Huang Z., Yao J., Song W., Li Y., Gong L. (2016) Theoretical analysis of fracture conductivity created by the channel fracturing technique, J. Nat. Gas Sci. Eng. 31, 320–330. [Google Scholar]
 Yan X., Huang Z., Yao J., Zhang Z., Liu P., Li Y., Fan D. (2019) Numerical simulation of hydromechanical coupling in fractured vuggy porous media using the equivalent continuum model and embedded discrete fracture model, Adv. Water Resour. 126, 137–154. [Google Scholar]
 Yang Y., Liu Z., Yao J., Zhang L., Ma J., Hejazi S.H., Luquot L., Ngarta T.D. (2018) Flow simulation of artificially induced microfractures using digital rock and lattice Boltzmann methods, Energies 11,8, 2145. [Google Scholar]
 Yang Y., Li Y., Yao J., Zhang K., Iglauer S., Luquot L., Wang Z. (2019a) Formation damage evaluation of a sandstone reservoir via porescale Xray computed tomography analysis, J. Petrol. Sci. Eng. 183, 106356. [CrossRef] [Google Scholar]
 Yang Y.F., Wang K., Zhang L., Sun H., Zhang K., Ma J.S. (2019b) Porescale simulation of shale oil flow based on pore network model, Fuel 251, 683–692. [CrossRef] [Google Scholar]
 Yang Y., Yang H., Tao L., Yao J., Wang W., Zhang K., Luquot L. (2019c) Microscopic determination of remaining oil distribution in sandstones with different permeability scales using computed tomography scanning, J. Energy Resour. Technol. 141, 9, 092903. [Google Scholar]
 Zeng D., Zhang E., Yao G., Zhu H., Xian Q., Shi T., Ding Y., Yi Y. (2018) Investigation of erosion behaviors of sulfurparticleladen gas flow in an elbow via a CFDDEM coupling method, Powder Technol. 329, 115–128. [Google Scholar]
 Zeng J., Li H., Zhang D. (2016) Numerical simulation of proppant transport in hydraulic fracture with the upscaling CFDDEM method, J. Nat. Gas Sci. Eng. 33, 264–277. [Google Scholar]
 Zeng Q.D., Yao J., Shao J. (2019) Study of hydraulic fracturing in an anisotropic poroelastic medium via a hybrid EDFMXFEM approach, Comput. Geotech. 105, 51–68. [Google Scholar]
 Zhang G., Gutierrez M., Li M. (2017) A coupled CFDDEM approach to model particlefluid mixture transport between two parallel plates to improve understanding of proppant micromechanics in hydraulic fractures, Powder Technol. 308, 235–248. [Google Scholar]
 Zhang J., Ma G., Dai Z., Ming R., Cui X., She R. (2018) Numerical study on pore clogging mechanism in pervious pavements, J. Hydrol. 565, 589–598. [CrossRef] [Google Scholar]
 Zhang L., Jing W., Yang Y., Yang H., Guo Y., Sun H., Zhao J., Yao J. (2019) The investigation of permeability calculation using digital core simulation technology, Energies 12, 17, 3273. [Google Scholar]
 Zhang T., Sun S. (2019) A coupled Lattice Boltzmann approach to simulate gas flow and transport in shale reservoirs with dynamic sorption, Fuel 246, 196–203. [CrossRef] [Google Scholar]
 Zhao J.L., Kang Q., Wang Y., Yao J., Zhang L., Yang Y. (2020) Viscous dissipation and apparent permeability of gas flow in nanoporous media, J. Geophys. Res. Solid Earth e2019JB018667. [PubMed] [Google Scholar]
 Zhao T., Zhao H., Li X., Ning Z., Wang Q., Zhao W., Zhang J. (2018) Pore scale characteristics of gas flow in shale matrix determined by the regularized lattice Boltzmann method, Chem. Eng. Sci. 187, 245–255. [Google Scholar]
 Zhou K., Hou J., Sun Q.C., Guo L.L., Bing S.X., Du Q.J., Yao C.J. (2018) A study on particle suspension flow and permeability impairment in porous media using LBM–DEM–IMB simulation method, Transport in Porous Media 124, 3, 681–698. [Google Scholar]
 Zhu G.P., Yao J., Sun H., Zhang M., Xie M.J., Sun Z.X., Lu T. (2016) The numerical simulation of thermal recovery based on hydraulic fracture heating technology in shale gas reservoir, J. Nat. Gas Sci. Eng. 28, 305–316. [Google Scholar]
 Zou Y., Ma X., Zhang S., Zhou T., EhligEconomides C., Li H. (2015) The origins of lowfracture conductivity in soft shale formations: an experimental study, Energy Technol. 3, 12, 1233–1242. [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Flowchart about the CFDDEM coupling method. 

In the text 
Fig. 2 Softsphere overlap model. 

In the text 
Fig. 3 The rough fracture geometry and the investigated region. 

In the text 
Fig. 4 The particle dune profile in the Eulerian view. The dune outline of the reference result (Zeng et al., 2016) was outlined by a dotted black line. 

In the text 
Fig. 5 The particle trajectory comparison between the simulation and reference results (Haddadi and Di Carlo, 2017). 

In the text 
Fig. 6 (a) Comparison of particle movement characteristics in rough and smooth fracture. (b)–(d) Particle trajectory on the XZ section (view along the direction of movement), YZ section (side view), and XY section (top view). 

In the text 
Fig. 7 The velocity vector of particles near the rough bottom and the decrease of kinetic energy due to collisions. 

In the text 
Fig. 8 The velocity of fluid and particle deposition rates under different roughness. 

In the text 
Fig. 9 Maximum deposition distance of particles with various density ratios. 

In the text 
Fig. 10 Particle amount distribution in areas A, B, and C. 

In the text 
Fig. 11 Particle radius distribution. 

In the text 
Fig. 12 The influence of particle diameter and injection position on the deposition distance of particles (the three rows of particles in orange were observed). 

In the text 
Fig. 13 (a) The trajectory of certain selected particles. (b) The decrease process of average kinetic energy. 

In the text 
Fig. 14 Distribution of particles at the bottom (top view). 

In the text 
Fig. 15 The plugged particle distribution in closed fractures. 

In the text 