A review of methods, applications and limitations for incorporating fluid flow in the discrete element method

Tuo Wang, Fengshou Zhang,*, Jason Furtney, Branko Damjana

a Key Laboratory of Geotechnical and Underground Engineering of Ministry of Education, Tongji University, Shanghai, 200092, China

b Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Shanghai, 200092, China

c Itasca Consulting Group, Inc., Minneapolis, MN, 55041, USA

Keywords:Hydro-mechanical process Fluid/discrete element method (DEM)coupling Geomechanics Numerical modeling

A B S T R A C T

1. Introduction

In recent years, the coupling technologies of fluid-solid have made great progress in the field of resource exploitation and geotechnical engineering due to the increase of engineering challenges and the enhancement of computer performance. The interaction between fluid and geomaterials has relevance to many geomechanics-related engineering projects, such as enhanced geothermal system (Xu et al., 2008), geological disposal of radioactive waste (Wang et al., 2006), sequestration of CO2(Rutqvist et al., 2007), sand production (Rutqvist et al., 2012), internal erosion (Bonelli and Marot, 2011; Cheng et al., 2018), slope failure(Dugan and Flemings, 2000; Zhao et al., 2017) and so on. Nevertheless, a comprehensive understanding of the hydro-mechanical coupling processes in geomechanics remains a great challenge.First, the flow of fluids through the pores and fractures in geomaterials is complex. Second, there are strong interactions between deformation of rock or soil and fluid flow, and these processes take place at very different time and length scales. Finally,because of the complexity of tectonic history,mineral composition,pore structure,and fracture networks,there is large uncertainty in the properties of geomaterials, especially in deep underground engineering with limited access to the engineered rock formation.Therefore, coupled hydro-mechanical processes in geomechanics will continue to be a subject of research within industry and academic community.

Numerical modeling has been a vital tool to investigate the coupled hydro-mechanical processes for geomaterials. The discrete element method(DEM),originally developed by Cundall and Strack(1979), has become widely accepted and used for the modeling of rock masses and granular soils,including the deformation,damage,fracturing, and stability coupled with fluid and heat transport processes. There are several advantages of using DEM to model hydromechanical processes. DEM uses individual discrete particles, governed by Newton’s law of motion and simple contact laws, to represent geomaterials.The representation can either be directly at the grain scale or use an engineering upscaling approach. This particle-based representation is similar to the basic constituents of rocks and soils. DEM exhibits emergent mechanical behaviors that arise from microscale grain(particle)packing and a set of grain and cement properties.Phenomenological damage laws are not required,since fractures or damages evolve naturally with breakages of bonds between particles.Material inhomogeneity can be naturally applied by adjusting the properties of the grains and bonds.There is no need to make any assumption about the geometry of induced fractures.Three-dimensional (3D) non-planar fractures, including their intersections and branching, are part of the model solution and are explicitly represented.Microseismicity can be predicted as a consequence of fracturing and joint deformation in the model. The integration of DEM-predicted synthetic microseismicity and actual microseismic data from the field can provide model validation.Such validation increases prediction confidence.

Fig.1. Microscale images of rock and sand with different porosities and cementations:(a)Marble with interlocking grains(Bandini et al.,2012);(b)Sandstone with broad grain gap(Haimson, 2007); and (c) Quartz sand (Brzesowsky et al., 2014).

The subject of computational fluid dynamics(CFD)has also made great progress since the 1970s(Anderson,1995).Numerical schemes have been developed to extend the Navier-Stokes and Euler equations to include turbulence (Hanjali′c et al., 2004), two- phase flow(Sussman et al.,1994), and other features. Recently, methods have been proposed to resolve meshing challenges, including meshless algorithms (Shadloo et al., 2016), dynamic and unstructured grid techniques (Liu et al., 2010), and the coupled lattice-Boltzmann method (LBM) (Chen and Doolen,1998). With the development of DEM and CFD, a number of modeling techniques are available to couple fluid-flow models with DEM for representing hydromechanical processes. However, these methods are applied to a broad spectrum of applications and no single fluid flow/DEM coupling method is suitable for all applications.A review of coupled fluid flow/DEM models was given by Furtney et al.(2013).This paper aims to provide an updated and comprehensive review of the theory and applications of DEM modeling of hydro-mechanical processes.And relevant highlights from some previous papers on every method are provided. Given their importance, the availability or unavailability of best practice guidelines is outlined. Different fluid flow/DEM coupling methods are introduced in this paper,including their implementation principles, application scenarios, and advantages and disadvantages. This paper is organized as follows.Section 2 introduces DEM and the specific fluid/DEM coupling schemes.Sections 3-7 describe the coupling schemes based on traditional continuumbased CFD,the LBM,dynamic fluid mesh(DFM)method,smoothed particle hydrodynamics (SPH), and pore network model (PNM),respectively.Section 8 describes some simple approximate methods to incorporate fluid flow effects in DEM models.Finally,a summary with speculation on future development trends is given in Section 9.

2. Overview of DEM

2.1. The principles of DEM

In traditional finite element method (FEM) or finite difference method (FDM), rock and soil are often regarded as continuous media on the macroscale. While in essence, they are discrete materials composed of a series of particles,pores,and fractures(Fig.1).Numerical schemes based on continuum mechanics have some limitations in the representation of discreteness, fracture and damage processes in rock and soil.The DEM provides an alternative for solving these problems using a micro-mechanics approach.DEM is a numerical method that simulates the movement and interaction of stressed assemblies of particles (typically polygonal or circular in two dimensions (2D), and polyhedral or spherical in 3D). The method approximates soil and rocks as particles packed and cemented together, naturally representing the microstructure and mechanical properties of rock and soil (Liu et al., 2021). The concept underlying the modern method was originally developed by Cundall (1971) as the distinct element method and was originally based on rigid polygons.The early application of the method was to analyze the stability of blocky rock slopes.The method was extended to 3D and different particle (block) shapes, including circular or spherical rigid particles. After nearly 50 years of development, DEM is widely used in the fields of geological, geotechnical, petroleum, and mining engineering.

Fig. 2. DEM calculation cycle (Itasca Consulting Group Inc, 2015).

Together,Newton’s laws of motion and the relationship between force and relative displacement between particles (i.e. constitutive law)are the basic laws implemented by DEM,as shown in Fig.2.The dynamic behavior is represented numerically by a time step algorithm in which it is assumed that the velocities and accelerations are constant within each time step.The solution scheme is identical to that used by the explicit FDM for continuum analysis (Itasca Consulting Group Inc, 2015). In other words, the time step is small enough so that, within a single time step, disturbances cannot propagate further from any particle to its immediate neighbors.Thus, during the simulation, the interparticle forces exist only between the adjacent particles. The magnitude of the interparticle contact force is determined by the relationship between force and relative particle displacement,called the contact model.

Fig. 3. Clump constructed by spheres (after Zheng et al., 2019).

2.2. DEM contact models

There are different DEM contact models intended for different applications. The linear contact model is the basic model for simulating granular material,and linear elasticity(no tension)and frictional sliding are provided(Cundall and Strack,1979).The linear contact model can be conceived as two elastic springs acting in the normal and shear directions at particle-particle or particle-wall contact points. Frictional behavior is introduced via the Coulomb sliding criterion, and viscous behavior is introduced via dashpot elements. More complex contact models that include additional physical effects are built on the basic ideas of the linear contact model. A summary of DEM contact models relevant to hydromechanical modeling is given here.

The rolling resistance contact model and Hertz contact model extend the linear contact model. In the rolling resistance model, a rolling resistance mechanism is added in order to dissipate energy during the rotation of non-spherical particles in granular flow(Iwashita and Oda, 1998). In the Hertz contact model, contact stiffness varies with the relative displacement of particles, so that the normal and shear forces at the contacts can include nonlinear effects(Dong et al.,2018).The hysteretic contact model is the Hertz model with attached dashpot elements,which is used to reproduce the rheological behavior of soil(Machado et al.,2012).For modeling rocks, the linear contact bond model, linear parallel bond model,flat-joint model, and soft bond model are most commonly used(Potyondy and Cundall,2004;Potyondy,2011,2017;Ma and Huang,2018). The linear contact bond model is an extension of the basic linear model to include a finite tensile and shear strength. The contact bonds allow tensile forces to develop at a contact but do not resist relative rotation.The parallel bond contact model is regarded as two additional elastic springs with constant normal and shear stiffnesses,acting uniformly over a cross-section(rectangular in 2D,and circular in 3D) located on the contact surface and centered on the contact point. These springs act in parallel with the springs of the linear component.This allows the parallel bond contact model to transmit moment between particles (Potyondy and Cundall,2004). The flat-joint model and soft bond model are similar to the parallel bond model,but they allow partial contact breakage or softening with strain to improve the uniaxial compressive/tensile strength ratio of the model(Potyondy,2012a;Ma and Huang,2018).The Burgers contact model is able to reproduce creep mechanisms using Kelvin (Eq. (1)) and Maxwell (Eq. (2)) models connected in series in both the normal and shear directions (Itasca Consulting Group Inc, 2015). Beyond these models for continuous rock, joints and fractures in rock can be represented by the smooth- joint contact model which uses contacts oriented perpendicular to the joints (Mas Ivars et al., 2008). A summary of the commonly used contact models and their key features is shown in Table 1.

Table 1 Summary of the commonly used contact models in DEM.

where σ represents the stress, ε is the strain,tis the time,Eis the elastic modulus, and η is the viscosity coefficient.

2.3. Generation of non-spherical particle

In traditional DEM model, spheres (3D) or discs (2D) are the most common particle shapes.However,the mechanical properties obtained from spherical particulate systems are often much different from irregular particulate systems (Lu et al., 2015; Ji and Liu, 2020; Yin and Wang, 2021). Therefore, the coupling method between fluid and irregular-shaped particles is also the focus of this paper.In this section,some common construction methods of nonspherical particle in DEM are introduced.

Clump which is a rigid collection of bonded spherical particles is widely used to generate non-spherical particles(Fig.3).Clumps are able to translate and rotate (generalized velocity and angular velocity, or spin of the clump centroid)) like spheres, based on the identification criteria of particle contact and obey the equations of motion.

Another way to generate realistic-shaped particles is superquadric method, which can generate particles by the variation of aspect ratio of particles and the smoothness of the surface (Fig. 4)based on ellipsoid equation. The particle shape equations (Barr,1981) are defined as

wherea,bandcare the semi-axis lengths of the particles along the major axes;n1controls the sharpness ofx-zandy-zplanes;andn2controls the sharpness ofx-yplane.

After generation of particles, the contact recognition of nonspherical particles is mainly to locate the contact point and the nearest particle.Wang et al.(2018)developed a series of techniques to recognize the shortest distance between particles and contact point. This method can locate contact quickly and is the most widely used method.

Besides clumps and super-quadric particles,dilated polyhedron particles are also commonly used to generate realistic particles.The generation process of dilated polyhedron particles is shown in Fig. 5. This method uses a particle of regular shape (i.e. sphere or cube)sweeping the exterior surface of arbitrary polyhedron based on Minkowski sum method (Varadhan and Manocha, 2006). The advantage of this method is that the center point of the conventional particle can be used to calculate the contact,which simplifies the calculation of contact.

Fig. 4. Particles generated by super-quadric method with different semi-axis lengths and aspect ratios (Wang et al., 2018).

In addition to the three commonly methods above, there are some other ways to generate non-spherical particles,such as polysuperellipsoid particles (Zhao and Zhao, 2019), B-spline particles(Maˇsín et al.,2006)and polyhedral particles(Boon et al.,2012).To date, these methods have not been widely used in the study of fluid/DEM coupling, thus they are not described in detail in this paper.

3. DEM coupled with CFD

Parallel to the development of DEM, CFD has made great progress, based on numerical solutions of Navier-Stokes equations (Eq.(5)).CFD solvers typically discretize domains using finite volume or finite element approximations, and nonlinear flow equations are solved with an iterative method. Most CFD methods work by initially estimating the pressure field by solving a nonlinear Poisson’s equation and then iteratively applying correction terms to the pressure equation until the velocity solution converges to the continuity equation(often incompressible flow).The details of CFD solvers are beyond the scope of this review. A good overview was given by Ferziger and Peric(2012).There are broadly two methods for representing solid particles in a CFD model,i.e.Euler-Euler and Euler-Lagrange methods.In the Euler-Euler method,both the fluid and the solid phases are represented as continuous phases and are both solved in the Eulerian frame of reference. A continuum numerical method is used to predict the variation of these two phases over time.The interaction between the two phases is described by the porosity and a drag force term. The drag force term (and the equal but opposite fluid-particle interaction force)depends on the local relative velocity and local solid phase volume fraction of each phase(Deen et al.,2007;Geng and Che,2011).The Euler-Lagrange method, also known as the discrete particle model (DPM), links a continuous fluid phase solved in the Eulerian frame of reference with individual discrete particles solved in the Lagrangian frame of reference. This method is the basis of most coupled CFD-DEM simulations. DEM calculates the particle motion and collisions.Typically, the DEM particles must be small relative to the CFD element size. A flowchart of the CFD-DEM coupling method is shown in Fig. 6.

Fig. 5. The generation process of dilated polyhedron particles (Ji and Liu, 2020).

In the CFD-DEM coupling,the basic Navier-Stokes equations are modified to consider solid particles by adding the porosity and drag force terms(Zhang et al., 2020a):

where ρfis the fluid density,pis the fluid pressure,μ is the dynamic viscosity of the fluid,is the fluid velocity,nis the porosity, andbis the body force per unit volume applied by the particles to the fluid in each fluid element.Arrows in Eq.(5)represent vectors(the same below).

Combined with a continuity equation that considers the porosity effect, the pressure field and velocity field can be solved:

The magnitude of the drag force applied to each individual particle can be obtained as follows (Di Felice,1994):

whereCdis the drag coefficient,ris the particle radius,andu→is the particle velocity.

The drag coefficient is defined as

whereRepis the particle Reynolds number defined as

The drag force is defined as

wheren-χis an empirical factor to account for the local porosity,and the empirical coefficient χ is defined as

The particle-fluid interaction force includes a drag force and a pressure gradient force:

whereis the acceleration of gravity.The pressure gradient term accounts for buoyancy if the hydrostatic gradient is part of the fluid pressure field. Sometimes the hydrostatic pressure gradient is omitted for numerical stability and the buoyancy forces are added explicitly. The body force per unit volume exerted on the fluid is

whereVis the volume of the fluid element,and the sum is over the particles which overlap the fluid element.

As the most common coupling method, the CFD-DEM method has been applied in many fields. In the area of geotechnical engineering,the sedimentation of granular materials(Zhao et al.,2014;Ferdowsi et al., 2017), internal erosion and related engineering problems are widely investigated by applying this method(Tao and Tao, 2017; Zhang et al., 2019). Tsuji et al. (1993) first adopted the method to investigate the behavior of particles in a fluidized bed using a 2D model.Mesh generation,visualization of fluid flow,and the position of particles are shown in Fig.7.Recently,3D CFD-DEM models have become widely used, and there have been significant advances in the technology of mesh generation and visualization.Hu et al. (2019) presented a CFD-DEM model to investigate suffusion in gap- and well-graded samples.Their study reproduced the full suffusion process and confirmed that the coupling method can reasonably describe the behavior of the system.Cheng et al.(2018)proposed a semi-resolved CFD-DEM model which is combined with non-distributed Lagrange multiplier/fictitious domain (non-DLM/FD) (Haeri and Shrimpton, 2012) to investigate the seepageinduced fine particle migration within the gap-graded soils. This modified CFD-DEM technique predicts a high-resolution flow field around coarse particles.

Fig. 6. Flowchart of the CFD-DEM coupling engine (modified from Hu et al., 2019).

Fig. 7. Mesh generation (a), visualization of fluid flow (b), and the motion of particles (c) of CFD-DEM simulation (Tsuji et al.,1993).

In petroleum engineering, the method is used to model the migration of proppant(Zhu et al.,2018),wellbore expansion under fluid injection and sand production (Climent et al., 2013). Zhang et al. (2017a) developed a CFD-DEM model to study proppant embedment and fracture conductivity after hydraulic fracturing.The study illuminated the process of mechanical interaction between a proppant pack and the rock formation during the process of fracture closing. Grof et al. (2009) investigated the evolution of sand production in the wellbore using the CFD-DEM model. The study revealed the process of separation of small clusters of cohesive particles and wellbore expansion.

The fluid mesh in conventional CFD-DEM is fixed in space,therefore, it is not possible to deal with moving boundaries. To address this limitation, the immersed boundary method (IBM)(Peskin, 1972) was created. In the IBM, a virtual body force is introduced so that a non-zero fluid velocity boundary condition can be imposed at a moving solid boundary. Kajishima et al. (2001)adapted the method to simulate the particle movement in fluid by taking the virtual body force applied to a fluid cell as simply approximated one by the solid volumetric fraction and the relative velocity of two phases in that cell.Subsequently,Yuki et al.(2007)simplified the method further.Guo et al.(2012)used CFD-DEM-IBM method to simulate slightly compressible two-phase gas-particle flows with complex boundaries and illustrated the advantages of introducing IBM into the traditional CFD approach. Delaney et al.(2012) used non-spherical super-quadric particles in the coupling scheme through equivalent spherical radius of the particle and a drag coefficient for non-spherical particle developed by Hölzer and Sommerfeld (2008).

The conventional CFD-DEM coupling requires a relatively coarse CFD grid because the particles must fit within the fluid elements.This limitation makes it difficult to study flow structures that are at or below the length-scale of the particles. One example is high Reynolds number turbulent flow that requires relatively fine meshes. For these cases, there are generally three numerical methods which are based on conventional CFD for approximating turbulent flows, i.e. direct numerical simulation (DNS) method(Moin and Mahesh, 1998), large eddy simulation (LES) method(Chin-Hoh, 1984), and Reynolds averaged Navier-Stokes equation(RANS) method (Li et al., 2009). There have been several related efforts on coupling DNS and LES solvers with DEM (Zhou et al.,2004;Berrouk and Wu,2010;Gui et al.,2010;Yoshida et al.,2013).

DNS resolves turbulent fluid flow directly using a high- resolution mesh. The Navier-Stokes equations are solved directly at the length scale necessary to capture all ranges of turbulence. The largest scale is in the order of flow width and the smallest scale is the Kolmogorov length scale (Schneiders et al.,2017):

where ηkis the Kolmogorov length scale, υ is the kinematic viscosity of the fluid,and κ is the dissipation rate of turbulent kinetic energy.

Yoshida et al. (2013) developed a coupled DNS-DEM procedure to predict the slurry viscosity of fine suspended particles and the results were in good agreement with the experiments and values obtained by empirical formulas.

Gui et al. (2010) investigated the particle collisions in swirling jets with a coupled DEM-DNS method. A mesh with resolution three time finer than the Kolmogorov length scale was used to represent the fluid. According to the work of Moin and Mahesh(1998), this setup is sufficient to resolve the relevant details of the turbulent flow.

Although DNS can directly resolve turbulent flow, it is computationally expensive and often impractical for engineering problems. To address this limitation, LES was introduced by Smagorinsky (1963). A coarse mesh is used in LES and the influence of sub-grid scale turbulent motion is approximated.Typically,LES models are calibrated to experimental observations.Compared with DNS, LES uses less memory and is faster. Schmeeckle (2014)used LES to simulate turbulent flow during transport of medium sand.In this work,lift and added mass are not included and smallscale acceleration and turbulence due to flow around individual particles are neglected.The author demonstrated that the model is capable of predicting bed load sediment transport rates, but the method needs further verification.

As the most widely used fluid-particle coupling method, the conventional and modified CFD-DEM methods have been adapted to many application areas. However, the traditional coupled CFDDEM method has some fundamental limitations related to the requirement of a coarse mesh and the computation cost. Details relevant to geotechnical applications like pore-scale porosity and permeability are not easily studied with this method. The convergence of CFD solvers in some CFD-DEM methods is greatly reduced due to strong body force terms and rapid porosity change.Also,the calculation of interfacial interactions among the multiphase flow with suspended particles is a challenge.

4. DEM coupled with the LBM

The LBM is a discrete mesoscopic approach to solve incompressible viscous flow. The method involves representing fluid as discrete parcels moving in a fixed lattice and subjected to simple collision rules. These discrete systems reproduce incompressible viscous flow without solving Navier-Stokes equations.Compared to traditional CFD that needs to solve the continuous partial differential equations (PDEs), a LBM solver is simpler. In addition, the LBM-DEM coupling method can naturally represent fluid flow in the pores formed by large particles and can accommodate multiphase flow.

4.1. The basic principles of LBM

LBM methods are classified according to the spatial dimension and lattice structure (discrete velocity number), and the computational accuracy is different for different forms of the method.An overview of the basic lattice-Boltzmann method is given first, followed by a description of LBM/DEM coupling. In this work, lattice Boltzmann D3Q19 (i.e. spatial dimension is 3 and discrete velocity number is 19) scheme is used to introduce the basic principles. As shown in Fig. 8, the space is divided into a series of cubic cells

Fig. 8. The schematic of LBM cell of the D3Q19 (the directions of the 19 discrete velocities are shown in the figure, and the center point is the direction of No. 0 discrete velocity) (modified form Brumby et al., 2015).

where each cell has a set of probability distribution functionsfirepresenting the density of fluid particles going through each one of the 19 discrete velocities.Traditional macroscopic properties of a flow including pressure and velocity can be recovered from these probability distribution functions.

In the lattice,the density of fluid particles at each cell position is

whereis the position of the center of the cell.

The velocity at each cell position is expressed as

whererepresents theith discrete velocity.

The D3Q19 discrete velocities can be written in a matrix form as

The velocity in each direction is given by

wherec= δx/δtis the characteristic lattice velocity defined by the grid spacing δxand the time step δt.

During each time step, the particles at all nodes move to the neighboring nodes along the lattice directions, colliding and changing the particle distribution. The change of the density distributions before and after collisions is typically computed through a single relaxation time approximation. There are two computational operations at each time step, streaming and collision operations. The streaming operation is expressed by

The collision operation is expressed by

where τ is the dimensionless relaxation time andfeqiis the equilibrium distribution function that depends on the macroscopic variables (fluid density ρfand velocityvat positionxand timet).

The two operations can be combined as

For the D3Q19 model, the equilibrium function is

where ωiis a weighting factor given by

The fluid pressurepis related to the fluid density byp=,wherecsis the speed of sound given in terms of the lattice speedcascs=. The kinematic viscosity of the fluid is given by

4.2. LBM-DEM coupling

DEM spherical particles can be well coupled with the LBM. In the scheme,LBM lattice is used to simulate the fluid flow.Inside the DEM particles,the lattice cells are fully occupied by the solid phase.Outside the sphere, the volume fraction of solid phase is zero. For cells intersecting the boundary of the sphere, the volume fraction occupied by the solid phase is between 0 and 1. A coupling law must ensure a smooth transition as the sphere enters and exits the LBM cells during its movement. The collision equation is modified as (Kruggel-Emden et al., 2016; Tran et al., 2017):

whereBnis a weight coverage function,and Ω is a collision operator representing the change of momentum due to collision of the DEM particle with the LBM cells. The weight coverage functionBnis expressed as (Brumby et al.,2015):

where εnis the volume occupation fraction for thenth cell. The collision operator for LBM cells and DEM particle is

The forces applied on the particle by the fluid are calculated based on the momentum exchange between the fluid and the particle(Han and Cundall,2013).The total forceFand torqueTover the DEM particle are given by

The Verlet integration (Verlet, 1976) used in DEM provides a clear limit for the integration step:

whereMminis the minimum value of the masses of all particles,andknis the stiffness of the contact. However, the above limit of the DEM time step is not always enough,especially for particles with a shape significantly different from sphere, such as flat particles.Experience has shown that the time step required for an accurate DEM simulation is one or two orders of magnitude below this critical limit.

The time step of LBM is given by

The time step of the coupled method should be smaller than the time steps of the individual DEM and LBM solvers.

Han and Cundall (2013) described and verified the LBM-DEM coupling system and demonstrated that it is a useful tool to simulate fluid-solid interaction in porous media.The authors used the model to simulate sand production and upward seepage flow in a dam.

Fan et al. (2019) developed a coupled LBM-DEM model to investigate the evolution of fracture conductivity and non-Darcy flow in proppant-supported hydraulic fractures.

Ding and Xu (2018) proposed coupling the LBM and DEM methods for solving problems of multiphase fluid in granular materials.In their study,a coupled model for a fluid-solid-gas system is described and the scheme is validated via classical fluid-solid coupling examples. According to their studies, LBM is an alternative, flexible and efficient approach to solve Navier- Stokes equations.

Kruggel-Emden et al. (2016) developed a fluid-particle- thermal lattice-Boltzmann scheme. In their research, fluid flow is modeled by a multiple-relaxation-time (MRT) scheme and the thermal flow by a Bhatnagar-Gross-Krook(BGK)collision model.In the research, the coupled problems of fluid flow and heat transfer were addressed by combining distribution functions for the density and internal energy.The work demonstrated that this approach can solve multi-physics problems.

Zhou et al. (2020) used a gap-graded granular material to investigate the initiation of internal erosion under different stress levels and hydraulic gradients. Brumby et al. (2015) investigated the erosion process of cohesive particles, in the context of marine sediments.Their studies showed that a coupled DEM-LBM model is suitable for representing materials other than uniformly sized loose particles.

For non-spherical particle shape, Galindo-Torres (2013) proposed a method to couple fluid and realistic-shaped particles based on LBM and the dilated polyhedral particles developed by Pournin and Liebling (2005). The basic principle is to search the minimum distance between geometric feature of DEM sphero- polyhedron and cell center.After that,a sphere of radius equal to the particle’s sphero-radius is defined and the algorithm that follows the coupling law for spheres is used. Chen and Wang (2017) used the same coupling scheme to simulate hydrofracturing process with triangular DEM particles and immersed moving boundary (IMB).The study revealed that fracture geometry is much more complex in heterogeneous samples.

Overall, the coupled LBM-DEM model captures the physics of fluid-particle interaction very well and even enables multiphase flow analysis at the pore scale. However, it is computationally inefficient compared with the conventional CFD-DEM scheme,especially in 3D cases. Further, the method is limited to low Reynolds number flows.

5. DEM coupled with a DFM

DEM coupled with a remeshing strategy, i.e. the DFM, was introduced by Zhang et al.(2020b)to study large-displacement soil suffusion problems in gap-graded soils. The numerical scheme is computationally efficient and resolves the fluid flow at the scale of the pores of the large particles.The dynamic remeshing allows the flow solution to remain accurate for large-strain movements of the soil.A flowchart describing the DFM-DEM model is shown in Fig.9.A DEM material sample is prepared first,and particles are created in a gap-graded (bimodal) size distribution. An initial fluid mesh is generated based on the configuration of coarse particles in the sample.The initial mesh is a tetrahedralization of the points given by the large particle centroids and the contact points between the large particles and the walls.The porosity and permeability for each fluid element are calculated according to the position and volume of both the coarse and fine particles. The fluid in the model is assumed to be incompressible. The fluid pressure is estimated by solving a Laplace’s equation,and the local flow rate in each pore is found by Darcy’s law. For the fine particles, the interaction force applied by the fluid is calculated according to the relative velocity of the particle and fluid as in the CFD-DEM coupling scheme described in Section 3. The fluid mesh, permeability, and fluid particle interaction force are updated at predetermined intervals. Detailed description of each step is given in the following sub-sections.

Fig. 9. Workflow of a DFM-DEM coupling model.

As the particles in the DEM model, the centroids of coarse particles and the contact points between walls and coarse particles are used to generate a new tetrahedral mesh (Fig. 10), this method provides a good representation of flow in gap-graded soils. The mesh is updated dynamically at predetermined time intervals.

After the fluid mesh isgenerated,a spatial searchalgorithm isused to find the closest tetrahedron centroid to each fine particle. Every fine particle is initially taken as completely within the tetrahedron with the nearest centroid. The volume of the sub-tetrahedra composed of the fine particle and each combination of three tetrahedranodes are obtained.If the sum of thesefourvolumes isthe same as the volume of the closest tetrahedron,this fine particle belongs to the closest tetrahedron.If not,the same process is repeated on adjacent tetrahedra until the nearest tetrahedron is found.

To obtain the porosity of each tetrahedron, the volumes of the coarse and fine particles in each tetrahedron are summed. For the coarse particles, the volume of a tetrahedron embedded in the particle (Vin Fig.11) is calculated first. Then a volume correction factor is applied to the embedded tetrahedron in order to add the small additional volume(Fig.11)near the spherical surface.For the fine particles, the total volume is recorded in the tetrahedron in which the fine particle centroid is located.Then the porosity,n, is

Fig. 10. Fluid mesh generated based on the skeleton of coarse particles in the DEM model. The yellow and blue balls represent the coarse and fine particles, respectively(Zhang et al., 2020b).

Fig. 11. Fluid meshes made up of coarse particles. V represents the volume of tetrahedron embedded in the particle, and S represents the area of a particle that receives the fluid pressure from the adjacent fluid mesh (Zhang et al., 2020b).

whereVtetis the volume of a tetrahedral mesh element;andVcoarseandVfinerepresent the volumes of the tetrahedrons embedded in the coarse and fine particles, respectively. The permeability is obtained by using the Kozeny-Carman equation (Bear,1972):

wheredmis the average diameter of fine particles in the mesh.

Low Reynolds number fluid flow in porous medium is described by Darcy’s law:

whereKis the matrix permeability.

The compressibility of the fluid is small enough to be neglected,leading to the incompressibility condition:

Combining Eqs. (35) and (36) results in a Laplace’s equation:

This equation is solved for pressure in each fluid element by using a standard cell-centered finite volume(FV)technique(Silpa-Anan and Hartley, 2008). Constant pressure boundary conditions are applied to the top and bottom of the sample. Fluid velocity at every element face is calculated via Darcy’s law. A simple interpolation scheme is used to estimate the fluid velocity at the center of each element:

The hydro-mechanical forces on the coarse and fine particles are calculated separately.For the coarse particles,the force that a single fluid element applies to the particle is calculated by multiplying the fluid pressure with the surface area of the coarse particle included in the fluid element(Fig.11). The resultant force that is exerted by the fluid on the coarse particle is

wherePiis the pressure,Siis the surface area of the particle in the fluid mesh, and the sum is over the tetrahedra attached to each coarse particle.

For the fine particles, the total force exerted by the fluid on a particle is the sum of the drag force and the buoyancy force. The calculation equations and related parameters are shown in Eqs.(7)-(13).

Zhang et al.(2020b)proposed the DFM method and combined it with DEM to investigate suffusion in gap-graded soil.The research revealed that the coupled DFM-DEM model can reproduce the expected variation of fluid flow and the migration of fine particles within pore network of coarse particles.

In coupled DFM-DEM method,a simpler linear equation is solved instead of the more complex Navier-Stokes equation.The calculation is fast and free of convergence problems that can occur in CFD-DEM.The coarse tetrahedral mesh can represent any geometry of interest.Significantly,local variations in porosity due to soil skeleton deformation and fine particle migration are described well.However,the method can only be used to solve incompressible porous media flows.Moreover,the current scheme has only been applied to specific gap-graded particle assemblies,where pores are formed from coarse particles. It has the potential to be applied to more general cases of particle assembly(e.g.uniform size distribution).

6. DEM coupled with SPH

In contrast to traditional CFD methods that use the Eulerian frame of reference, SPH uses a Lagrangian representation of the fluid.SPH does not require a traditional computational mesh,thus it can naturally handle complex geometries and moving boundaries. The method was developed by Lucy (1977) and Gingold and Monaghan (1977), and was originally used in studies in astrophysics, but is now applied in many fields including electromagnetism, heat transfer, and fluid mechanics (Libersky et al., 1993).The basic SPH method is introduced first, and then DEM/SPH coupling is described.

6.1. The theory of SPH

Kernel interpolation theory is the foundation of SPH. This transforms the fluid flow PDEs into integral equations where a kernel estimate of the field variables at an arbitrary point is given.

Information related to the flow field is known only at the position of discrete points,and the integrals are calculated as a weighted sum of neighboring points.In this way,spatial derivatives operating on the physical quantities are transformed into operations on the interpolation kernel. The values required for the calculation are provided by the discrete points and an interpolation kernel,and no fluid mesh is needed. The conventional kernel estimate is

whereWrepresents the kernel function,his a smoothing length,Dis the support of the kernel function,is a position vector,and′is a neighbor vector.

A spline is often used as the kernel function (Monaghan and Lattanzio,1985):

Eq.(40)is a continuous function and can be approximated by a summation over a discrete set ofNpoints:

wheremjis the mass of smooth particlej; ρjis the summation density; andf() is the value associated with particlej, which is located at′.

The equations of mass and momentum conservation in particle form are

whererepresents the velocity of the particlei,is the velocity difference between particlesiandj,pjis the pressure at the particle,andis the sum of the body forces (including gravity).

In order to improve the boundary deficiency and the tensile instability, Chen et al. (1999) proposed a corrective smoothed particle method by applying a Taylor series expansion to the kernel estimate:

whereURi,ρi,andpiare on the left side of the contact point;URj,pj,and ρjare on the right side of the contact point;andciandcjare the speeds of sound at the respective densities.

At fixed boundaries, the fluid velocity is fixed to zero and the fluid pressure of particles is given by the solution. For moving boundaries, the velocities and positions of the moving solid boundary are updated by independently solving the equations of motion of the rigid bodies. For modeling storm surge waves, a special moving boundary condition is used to create waves of a given size:

where η0is the target wave surface,ηsis the local wave surface,X0is the wave amplitude, ω is the frequency of the waves,His the wave height,Ais a transfer function,kwis the wave number, anddis the depth.

The repulsive force between the boundary and the fluid particle pairs is

whereiandkdenote the fluid and boundary particles,respectively;andpkis the pressure on the boundary particle. This form is the basis for SPH-DEM coupling.

6.2. SPH-DEM coupling

In a coupled SPH-DEM model, a layer of boundary particles between the SPH and DEM particles is used to calculate repulsive forces.The schematic diagram is shown in Fig.12 and the principle is as follows.

Fig. 12. Schematic diagram of coupled SPH-DEM method (modified from Wu et al.,2016).

The interaction force between fluid and particle is as follows:

where the subscript f represents the fluid, and b represents the solid.

The total external force and torque on a particle are given as

Ren et al. (2014) presented a 2D coupling SPH-DEM model to study the interactions between waves and rubble-mound breakwaters. Their model was validated by comparison with physical experiments. Fakhimi and Lanari (2014) proposed a coupled SPHDEM model for the simulation of rock blasting. In the simulation,the gas particles can flow through the microcracks and contact points between the grains,but are not allowed to flow into the rock grains.

Liu et al. (2020) employed repulsive force method to simplify the interaction between SPH particles and dilated polyhedra for modeling rock dumping process because that support domain can be truncated by the surface of the dilated polyhedron(Markauskas et al., 2017). Cleary et al. (2020) used super-quadric particles in a coupled SPH-DEM method to predict the behavior of slurry grinding due to medium and coarse rock interactions in a 3D pilot semi-autogenous (SAG) mill.

The coupled SPH-DEM method works well for modeling freesurface flows and dynamic boundary problems,such as storm surge waves,wave flume tests,and rock blasting where clear boundaries exist between the fluid and solid phases.

Fig.13. Schematics of the pore network model: (a) Pore domains formed from close chains of bonded particles; and (b) The drag force on a particle, Ffluid, as a resultant from the pore pressures of surrounding pore domains (Huang et al., 2019).

7. DEM coupled with a PNM

As early as the 1950s, Fatt (1956) represented flow in porous media as a discrete network of flow paths connecting pore spaces.Many researchers have used this approach to study the mechanics and fluid properties at the pore scale and the approach has been widely used in coupling with DEM. The pore network method is a natural counterpart to DEM and gives good insights into hydromechanical effects since fluid pressures can be applied directly to discrete grains. The flow network is made up of pores which are formed by the discrete particles(the flow network is the dual of the particle contact network). The pressure is constant in each pore(Fig.13a). In 2D sphere packings,the pores can be defined exactly,but in 3D,they must be approximated.The flow network is updated as the solid particles deform under the influence of fluid flow or external force. The effect of the fluid on a given particle is simply the resultant force of the pressures of all the pores enclosing that particle(Fig.13b).The flow pathways used to connect the pores are idealized as parallel-plate flow pipes at each DEM contact.Conceptually, these flow pipes are pore throats (Fig. 13). The aperture of the pipe can be determined by the location of two contact particles.The fluid is assumed to flow only in these nominal pipes and the flow rate,q(volume per unit time) is given by the Hagen-Poiseuille equation:

whereais the aperture of the flow pipe,p1andp2are the fluid pressures in the two neighboring pores,andLPis the length of the flow pipe.

The increment of pore pressure due to flow and volumetric strain is

whereKfis the fluid bulk modulus,Vpis the current pore volume,Nfpis the number of flow pipes around the pore,and Δtis the fluid time step. This system is sensitive to large changes in aperture as pressure drop is proportional to aperture cubed.To avoid numerical problem, the aperture is often taken as a small residual value plus the contact separation,and the aperture is typically capped at some upper value.

At first, the method was used to simulate the variation of permeability under stress in coupled fluid/DEM.Bryant et al.(1993)presented a network model of a disordered porous medium without assumptions regarding pore structure. Bruno (1994)developed a coupled fluid/DEM model with a deformable fluid flow network to calculate stress-induced permeability chage and damage in sedimentary rock.Compared to experimental results, it has been demonstrated that the method provides good insight into the physical mechanisms that control the stress-induced permeability change.

In addition to permeability variation,the coupled PNM-DEM is a useful tool to simulate the evolution of hydraulic fracturing in naturally fractured reservoirs. Thallak et al. (1991) pioneered the application of this method to the simulation of hydraulic fractures.Later, Zhao and Young (2012) developed a coupled PNM-DEM model to study fluid stimulation in a reservoir containing natural fractures. The model predictions qualitatively agree with the laboratory and field observations. These studies have shown that coupled PNM-DEM models can reproduce the evolution of complex networks of natural and induced fractures. Recently, Zhang et al.(2017b) used a 2D fully coupled hybrid discrete-continuum method to study the interaction of hydraulic and natural fractures(Fig.14a). Li et al. (2020) investigated supercritical carbon dioxide(SC-CO2)fracturing in both intact and fractured rocks using a PNM.Their study validated numerical predictions of the propagation of SC-CO2-induced fractures via a comparison with the Khristianovic-Geertsma-de Klerk (KGD) analytical solution.

Beyond application to cemented rocks, the PNM-DEM model has been applied to fluid injection into granular media represented by unbonded DEM particles. Zhang et al. (2013) used a 2D PNMDEM model to simulate the process of fluid injection into a dense granular medium. The numerical results revealed the evolution process of fracture in loose granular material and indicated that the methodology is capable of reproducing phenomena consistent with those observed in laboratory injection experiments (Huang et al.,2012) (Fig. 14b). Further, the numerical results provided insight into the effects of fluid viscosity and injection rate on the evolution of fracturing. Progress has also been made on practical 3D fieldscale numerical modeling of hydraulic fracturing (Damjanac et al.,2010; Shen et al., 2020).

Catalano et al. (2014) performed a pore-scale finite volume(PFV)method to couple DEM.Their study applied PNM-DEM in 3D model. The method is consistent with the theory of poroelasticity and is able to solve transient problems in the quasi-static regime which is verified by an oedometer test simulation and in good agreement with Terzaghi’s analytical solution.

Fig.14. The simulation of fracture evolution by PNM-DEM:(a)The interaction of a hydraulic fracture and a natural fracture(Zhang et al.,2017b);and(b)The evolution of fracture in granular material (Huang et al., 2012).

In their study, regular triangulation and its dual Voronoi graph are used to discretize the void space (Fig.15). The pressure field is defined by the values of pressure in each tetrahedral element located at the vertices of the Voronoi’s graph. The mass balance equation integrated on the poreiand recast into a surface integral by using the divergence theorem gives

whereVfiis the total pore volume, andu-vis the velocity of the fluid relative to that of the solid phase.The integral on facetSijgives a flux exchanged between adjacent tetrahedra, noted asqij.

The total forceFfiexerted by the fluid on a particleiresults from the pressure and viscous stress acting at the surface:

where ∂Γidenotes the solid surface of the particlei,pais the absolute pressure, and τ is the viscous shear stress tensor.

Although the coupled PNM-DEM method represents the hydromechanical process more directly than other coupled methods, it has only been used in the simulation of small-scale problems. For the scale of practical field engineering, the simulation of fracture network evolution is still a challenging task, especially for 3D simulations (Huang et al., 2018).

8. Other methods of hydro-mechanical coupling in DEM

In addition to the methods described above, other methods of hydro-mechanical coupling in DEM have also been developed.These methods capture some key elements of fluid-solid interaction, but do not require intensive computational efforts. The fluid flow in some methods is maybe typical and not explicitly modeled,but the interaction between fluid and particles is realized in an approximate manner.

8.1. Ice load on floating structure in broken ice field

Liu and Ji(2018)used DEM with dilated polyhedral elements to simulate the interaction between ice floes and the floating structure (Fig.16). A 2D-Voronoi tessellation algorithm is employed to generate polyhedral ice floes in broken ice field. In their study, a one-way hydrodynamic force is considered. The buoyancy is calculated by meshing tetrahedrons onto the envelope polyhedron,and the drag force and moment are determined by functions of relative motion between the objective body and fluid.

Fig.15. (a) Regular triangulation and (b) its dual Voronoi diagram (Catalano et al., 2014).

8.2. Shining-lamp algorithm and pressure chain algorithm

The shining-lamp algorithm, described by Potyondy (2012b), is used for applying a confining fluid pressure to the potentially complex and moving boundary of an impermeable body.The basic principle is to suppose that many light sources shine on a body from infinitely far away. Space is divided into cells and each light source illuminates a small area (Fig.17a). The outermost particles for each area are identified in the positive and negative of each coordinate direction and forces are applied to these particles(Fig.17b).The magnitude of the force is specified by using pressure multiplied by the area of intersection.

In a 2D model, a similar effect can be achieved by identifying particle connections called pressure chain algorithm. Emam and Potyondy (2010) developed a 2D DEM model to simulate cutting rock confined by heavy drilling mud.As shown in Fig.18,the shape of the rock surface changes constantly and bears the pressure of mud during cutting process.

Wang et al. (2020) developed a DEM-continuum model to investigate the slot-shaped breakout in high porosity sandstone.Fluid pressure is applied to the DEM particles via shinning-lamp algorithm to represent the uniform confinement effect. Compared to applying loads with a DEM wall, the shinning-lamp algorithm allows local boundary deformation and applies pressure to a changing cross-sectional area. The strength of the method is the adaptability to complex evolving surface shapes.

Fig.16. Snapshots of different ice concentrations (40% and 80%) by DEM simulations (Liu and Ji, 2018).

Fig.17. The schematic of the shining-lamp algorithm: (a) An imaginary plane light source; and (b) The schematic of illumination area (Potyondy, 2012b).

8.3. Oscillatory flow

Bed load transport under the action of oscillatory flow is a common form of motion of marine material in the offshore environment(Dumas et al.,2005).Drake and Calantoni(2001)used the equation of motion for sediment particles described by Madsen(1991) to model particle transport in oscillatory sheet flow:

Fig.18. The model of cutting rock in drilling mud. Applied external forces by mud are shown in upper-left image (Emam and Potyondy, 2010).

where ω = 2π/T，φ is the waveform parameter,andiis the number of sinusoidal superpositions.

8.4. Proppant flowback

In order to better understand the proppant flowback in hydraulic fractures,Asgian et al.(1995)developed a simplified fluidparticle coupling method to reveal the evolution of the failure of proppant in hydraulic fractures(Fig.20).The efficacy of the model is shown by comparing with the experimental results (Milton-Tayler et al.,1992).

Fig. 20. The process of proppant flowback in simulation (Asgian et al.,1995).

In this method,the drag force tending to push a proppant grain in the direction of fluid flow has two components: pressure and shear.The pressure component varies with the cubic power of grain size and varies linearly with drawdown given by the analytic solution:

whereris the radius of proppant grain,andxis the distance in the direction of the fracture.

For the shear component,the analytic solution for laminar flow around an isolated grain is used:

Fig.19. Profile variations of oscillatory sheet flow (Dumas et al., 2005).

8.5. Capillary force

Particle assemblies with a small fraction of fluid in the pore space can form capillary bridges(Fig. 21) between particles. These bridges of fluid provide a cohesive force between particles. The capillary force between two particles,which are of the same size,is proportional to the radius(R) and the surface tension(γ):

Fig. 21. Capillary bridge between two moist particles.

wheredis the dimensionless force;L*is the dimensionless distance;Lis the separation distance of particles;φ is the contact angle of wetting fluid at the surface of particles;V*is the volume of the capillary bridge;riandrjare the radii of particlesiandj, respectively; anddijis the distance between particle centers.

Dimensionless force can be recovered by solving the Laplace-Young equation. Willett et al. (2000) gave a curve-fitting to the numerical solution. Grof et al. (2009) improved the approximation approach and modeled the process of particle entrainment by fluid flow. The method is only valid for small saturation values. Gras et al. (2011) investigated the link between suction and water content in granular media via DEM coupled with the Laplace-Young equation, and made comparisons with experimental results.

9. Concluding remarks

Numerical models incorporating fluid flow effects with DEM are used in a wide variety of engineering applications to predict the evolution of different processes in rocks,soils,and structures.It is a complex topic which involves many physical phenomena and a diversity of numerical calculation methods. As computer performance has improved and numerical analysis theory has advanced,it has become possible to apply these coupled fluid/DEM methods to practical problems. On one hand, DEM has developed and matured and, with the advent of more sophisticated contact models, is able to model a variety of processes in rocks and soils taking place at different scales. On the other hand, CFD has also made great progress. Multiple new calculation models and methods are proposed. More and more CFD methods are being coupled with DEM. The accuracy of these coupled methods has increased with the advances in computing power in past decade.This paper has summarized a variety of coupling methods and described the calculation principles, application scenarios, and limitations. The summary of common coupling methods is shown in Table 2. Practical experience shows that researchers and engineers must choose the appropriate approach to solve a problem according to the advantages and disadvantages of different coupling methods.

Table 2 The summary of common coupling methods.

Although coupled fluid/DEM technology has made great progress,there are still some problems to be solved in the future,which are listed below:

(1) No existing method solves efficiently and accurately the fluid/DEM coupling problem involving highly compressible fluid.

(2) With the development of computer technology, many grid generation technologies have emerged,such as the dynamic grid method, and the overlapping grid method (Desquesnes et al., 2006). However, it is still a time-consuming and challenging work to generate a suitable grid, especially in problems with complex shapes and moving boundaries.

(3) Although computing power has increased dramatically in recent years, calculation time is still a practical limitation when the grids are fine,such as those required for calculation of turbulent flow.

(4) There are few studies on multiphase flow effects on particles.

(5) The simulation of highly viscous non-Newtonian fluids(such as that used in grouting) is still difficult.

(6) Combinations of these methods in different spatial regions and multi-scale methods are challenges to be solved in the future (Furtney et al., 2013).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant Nos. 41772286 and 42077247) and the Fundamental Research Funds for the Central Universities, China.

List of symbols

Σ Stress

EStrain

TTime

EElastic modulus

HViscosity coefficient

a,b,cSemi-axis lengths of the particles along the major axes

n1Parameter that controls the sharpness of x-z and y-z plane

n2Parameter that controls the sharpness of x-y plane

ρffluid density

pfluid pressure

μdynamicviscosity ofthe fluid

nporosity

CdDrag coefficient

rParticle radius

Particle velocity

RepParticle Reynolds number

n-χ Empirical factor to account for the local porosity

χ Empirical coefficient

Acceleration of gravity

VVolume of the fluid element

ηkKolmogorov length scale

υ Kinematic viscosity of the fluid

κ Dissipation rate of turbulent kinetic energy

Position of the center of the cell

ith discrete velocity

cCharacteristic lattice velocity, andc=δx/δt

τ dimensionless relaxation time

feq

iEquilibrium distribution function that depends on the macroscopic variables (fluid density ρfand velocityvat positionxand timet)

ωiWeighting factor

csSpeed of sound, andcs=c/

BnWeight coverage function

Ω Collision operator representing the change of momentum due to collision of the DEM particle with the LBM cells

εnVolume occupation fraction for thenthcell

pPosition of the particle’s center of mass

MminMinimum value of the masses of all particles

knStiffness of the contact

δLBM

tTime step of LBM

δDEM

tTime step of DEM

nPorosity

VtetVolume of a tetrahedral mesh element

Vcoarse，VfineVolumes of the tetrahedron embedded in the coarse and fine particles, respectively

kPermeability

dmAverage diameter of fine particles in the mesh

KMatrix permeability

SjArea of each element face

WKernel function

hSmoothing length

D Support of the kernel functionPosition vector′Neighbor vector

ρjSummation density

f() Value associated with particlej

mjMass of smooth particlej

Velocity of the particle of particlei

Velocity difference between particlesiandSum of the body forces(including gravity)

η0Target wave surface

ηsLocal wave surface

X0Wave amplitude

ω Frequency of the waves

HWave height

ATransfer function

kwWave number

dDepth

NNumber of spheres on one particle

NfNumber of fluid particles interacting with the solid sphere

aAperture of the flow pipe

p1，p2Fluid pressures in the two neighboring pores

LPLength of the flow pipe

qFlow rate

KfFluid bulk modulus

VpCurrent pore volume

NfpNumber of flow pipes around the pore

ΔtFluid time step

paAbsolute pressure

τ Viscous shear stress tensor

Vsg,AVolume and projected area of the sediment grain,respectively

Cd,CmDrag and added mass coefficients, respectively

rRadius of proppant grain

L*Dimensionless distance

LSeparation distance of particles

φContact angle of the wetting fluid at the surface of particles

V*Volume of the capillary bridge

ri,rjRadii of particlesiandj, respectively

dijDistance between particle centers