Numerical simulation of flow behavior of particles in a porous media based on CFD-DEM

Journal of Petroleum Science and Engineering 171 (2018) 140–152

Contents lists available at ScienceDirect

Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol

Numerical simulation of flow behavior of particles in a porous media based on CFD-DEM

T

Shuyan Wanga,∗, Haoting Lia, Ruichen Wanga,b, Ruichao Tiana, Qiji Suna, Yimei Maa a b

School of Petroleum Engineering, Northeast Petroleum University, Daqing, 163318, China Petroleum Engineering, Missouri University of Science and Technology, Rolla, 65409, USA

A R T I C LE I N FO

A B S T R A C T

Keywords: Porous media Microchannel Computational fluid dynamic (CFD) Simulation Discrete element method (DEM) Grains of sand

The movements of particles in porous media exist in various modules of oil field development and production. Flow behavior of particles in a porous media is simulated by means of CFD-DEM. The lubrication force originated from the normal “squeezing” motion of the particles is considered. By changing the porosity of the porous media, the distributions of velocity, residence time of particles, and granular temperature are analyzed within the porous media. At a low porosity, particles will be blocked and deposited in the porous media. Increasing porosity prolongs the trajectory of particles, and particles do reciprocating motion and flow out the microchannel. Simulated results indicate that with an increase of porosity, the residence time is decreased, and then increased, while the granular temperature is in reverse. The predicted axial velocity of particles and the granular temperature decrease with an increase of liquid viscosity, while the mean residence time and contact force increase with an increase of liquid viscosity. Simulations indicate that the granular temperature and the contact force are increased, and the residence time of grains is reduced, with an increase of liquid velocity. The main contribution here is that the hydrodynamic characteristics of particles in a porous media are revealed based on the detail information of flow behaviors of particles.

1. Introduction Petroleum plays a vital role in modern industry. Today petroleum is predominantly processed to provide transportation fuels and numerous petrochemicals and materials, including pharmaceuticals, solvents, fertilizers, pesticides, and plastics (Parra, 2004). Crude oil generation is followed by migration from the source rocks into carrier beds and along faults, until an appropriate trap is reached. The movement of crude oil leads to petroleum accumulations (McAuliffe, 1979). Petroleum occurs naturally in sedimentary rocks that contain abundant clay minerals, sandstone particles, fine particles of rock and quartz crystals (Tannenbaum et al., 1986). The existence of fine solid particles influences oil velocity and free path, with the change porosity of the microchannel altered by fine particle deposition, therefore an irreversible and permanent impairment of reservoir porosity and permeability are caused (Ilyas et al., 2018). Moghadasi et al. (2004) used glass beads and sand beds to study the general behavior of fine particle movement and deposition in porous media. And a mathematical model is presented that simulates the porosity impairment by particle movement and deposition. Wu et al. (2012) reviewed the progress in the understanding of the role of clay minerals in crude oil formation, migration and

∗

accumulation. Studies showed that clay minerals influenced the migration of crude oil in underground petroleum reservoirs. Parnell et al. (1996) indicated that the petroleum fluid had a substantial aqueous component to entrain quartz during the migration of petroleum. The transport of grains of sand was one of the important flows in natural rock fractures. Crandall et al. (2010) used fracture meshes to develop relationships between the observed roughness properties of the fracture geometries and flow parameters that were of importance for modeling flow through fractures in field scale models. Rock cracks are important conduits for petroleum migration. Therefore, an understanding of petroleum migration in and along cracks is critical to successful appraisal of the petroleum systems and exploration targets within these basins, where oil and gas occurs mainly along cracks, suggesting petroleum migration and accumulation is closely related to the cracks activity (Cao et al., 2010). Mathematical modeling and computational fluid dynamic (CFD) simulations of the flow in the rock cracks give detailed information regarding the local phase velocity, spatial distributions, liquid phase flow patterns, particle residence time, and particle motion trace. This is especially important in the regions where measurements are either difficult or impossible to be obtained. Such information can be useful in

Corresponding author. E-mail address: [email protected] (S. Wang).

https://doi.org/10.1016/j.petrol.2018.07.039 Received 31 March 2018; Received in revised form 8 July 2018; Accepted 12 July 2018 0920-4105/ © 2018 Elsevier B.V. All rights reserved.

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

solid flow, and is relatively well documented in the literature (Xu and Yu, 1997; Kang et al., 2009; Washino et al., 2013); the salient features of the model equations used will be summarized below. In fact, the sand grains are the different sizes and shapes, and the size and shape of sand gains will influent the flow behaviors of grains and liquid, but the sand grain volume fraction is negligibly small. Initially, the sand grains of 500 were put into the channel, so the sand grain volume fraction is 0.025. For simplify, it is assumed that the solids phase consists of monosized particles with the same diameter and density.

the understanding of the motion of liquid and particles in rock cracks (Slider, 1976). In the present study, the cracks in the reservoirs are represented as microchannels, and the pore structure is simulated by laying out some fixed big particles, where gaps of fixed big particles are available. Grains of particles flow in gaps of fixed big particles through the microchannels with liquid carrying. In the CFD modeling of liquid-solid two-phase flow, Eulerian-Lagrangian (E-L) models are widely applied to fluid-solid flow. The E-L approach traces individual particles described by Newton's laws of motion on an individual particle scale, while the continuum liquid phase flow is described by Navier-Stokes equation on a computational cell scale (Deen et al., 2007; Zhu et al., 2007; Tsuji, 2007). The mechanism of particle-to-particle collisions may be described by the soft-sphere model (Cundall and Strack, 1979; Hoomans et al., 1996). The soft-sphere model or discrete element method (DEM), allows for multiple particle overlaps while the net contact force is obtained from the addition of all pair-wise interactions. Grof et al. (2009) predicted the stability of a random packed layer of cohesive particles under gravity and laminar cross-flow by means of CFD-DEM. Definitely, CFD-DEM can be used to investigate the local particle flow field, highlighting the motion of particles in apparently chaotic vortices continuously forming and disappearing during the process of sand production in oil exploration. Fatahi and Hossain (2016) proposed a DEM based numerical model to analysis the fluid flow through reservoir's porous media with complex characteristics, which can potentially be used to handle any reservoir conditions and complex set of discontinuity and fluid flow regime. Qiu (2015) proposed an approach to investigate fluid flow through porous media using the DEM for modeling particles dynamics and the lattice Boltzmann method for simulating the fluid motion. Wachs (2009) investigated effects of particles of arbitrary shape and actual contacts on the flow of particles in a Newtonian fluid using DEM solver. The simulations showed that particles actually touched each other, which was in contrast with the repulsive force based on the collision model, while the 2D sedimentation of isometric polygonal particles with collisions was presented. Qiu and Wu (2014) developed a hybrid scheme coupling the discrete element method (DEM) with computational fluid dynamics to model solid-liquid flows in a closed container. The motion of the solid phase was obtained using DEM, in which the particle-particle and particle-wall interactions were modeled using the theoretical contact mechanics, and the two phases were coupled through the Newton's third law of motion. At present, the researches on the hydrodynamics behaviors concentrate on one-phase flow, not considering solids particle flow (Santanna et al., 2013; Hatzignatiou et al., 2013; Ferreira and Moreno, 2017), therefore, detail investigations of flow characteristics on liquidsolid two-phase in rock cracks are still lacking. In particular, sedimentary rocks in low-permeability oilfield are rich of clay particles, quartz crystals and sandstone particles. The existence of fine particles will affect the path and direction of the migration of underground crude oil along the pores of rock. Furthermore, the internal structure of porous media will lead to the differences in the migration of particles, which causes the deposition of solid particles in varying degrees, and then modified the porosity of rocks, thus affecting the migration of crude oil. Therefore, the detailed flow behaviors of particles are crucial to reveal oil migration law. In the present study, the flow behavior of grains of sand, the solids phase, is simulated by means of CFD-DEM in the rock cracks of varying porosity. The hydrodynamics in rock cracks, modeled as microchannels, is simulated. By changing the velocity and viscosity of liquid, porosity of microchannel, the distributions of grains are studied within the microchannel. The variations of velocity and residence time of grains, granular temperature and contact forces are analyzed at different porosities, velocities and viscosities of liquid.

2.1. Equation of motion for liquid phase Generally, in Eulerian-Lagrangian method the fluid phase flow is solved by a locally averaged approximation of the continuity and Navier-Stokes equations with an averaging scale of the order of the computational cell. The equations of conservation of mass and momentum for liquid phase are:

∂ (ρl εl ) + ∇⋅(ρl εl ul) = 0 ∂t

(1)

∂ (ρl εl ul) + ∇⋅(ρl εl ul ul) = −εl ∇P + εl ∇⋅τl + εl ρl g − Fpl ∂t

(2)

where g is the acceleration due to gravity, P is the liquid pressure, εl is the liquid volume fraction, τl is the viscous stress tensor, and ρl is the density of liquid. The coupling term Fpl between the particle and liquid phases is approximated as the sum of the drag on each particle within the corresponding fluid control volume. The stress tensor of liquid phase can be represented as

τl = μl [∇ul + (∇ul)T ] −

2 μ (∇⋅ul) I 3 l

(3)

where μl is the viscosity of liquid phase. 2.2. Equation of motion for a particle Spherical particles of uniform size are investigated in the present work. The particles are tracked individually based on the Newton's second law of motion. Each particle has two types of motion: translational and rotational. The motion of each individual particle is governed by the laws of conservation of linear momentum and angular momentum, expressed, for the i-particle, by

mi

Ip

d vi = −Vp ∇P + mi g + fd + dt

d ωi = dt

N

∑ j=1 (flj + fcj)

(4)

N

∑ j=1 Tpij

(5)

where mi and vi are the mass and velocity of a particle, and Vp is the volume of a particle. Tp is the torque arising from the tangential components of the contact force. Ip and ω are the moment of inertia and angular velocity of a particle. The terms on the right-hand side of Eq. (4) are the liquid pressure gradients, gravity, drag force exerted from the fluid, lubrication force and contact force, respectively. These interparticle forces and torques are summed over the N particles in contact with particle i. The contact force between particles is calculated based on the soft-particle method. The liquid-solid interaction force, or drag force, is determined for each particle. The drag force depends on not only the relative velocity between the solid particle and fluid but also the presence of neighboring particles, i.e., local volume fraction of the solids phase. The drag force is expressed by considering these factors as follows:

fd = 2. Eulerian-Lagrangian liquid-solid flow model

βVp 1 − εl

(ul − up)

(6)

where εl and β are the volume fraction of the fluid and an inter-phase momentum transfer coefficient. A proper drag model for the description

The CFD-DEM approach has been widely used in the field of liquid141

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

of β is vital in solid-fluid interaction problems. In present simulations, the inter-phase momentum transfer coefficient is predicted by the Huilin-Gidaspow model (Lu et al., 2004). As two particles move relative to each other with only a thin film of liquid medium separating them, hydrodynamic forces, often referred to as lubrication forces, arise due to the motion of the interstitial liquid. In the simple case where the particles are moving away from each other, the liquid has to flow into the gap between the particles; analogously, when the particles are moving toward each other, the liquid has to flow out of the gap. This motion of the liquid generates pressure gradients and viscous stresses that are cumulatively expressed by a hydrodynamic force on the particles. Because the approach of the particles is often oblique, the hydrodynamic force has a normal component (in the direction that connects the particle centers) and a tangential component which results in torque on the particles. It has been shown that when the particles are in close proximity, the prime contribution to the lubrication force is the normal component, resulting from the normal “squeezing” motion of the particles; it is only this component that is considered in this work. This force is given by the following equation (Crowe et al., 1997):

The displacement-related elastic contact stiffness, kn, depends on the elastic properties of the colliding particles. The spring coefficients kn and kt are calculated from the following equations (Tsuji et al., 1993). For particle-particle collisions, the spring and damping coefficients are: −1

kn =

2 − γj ⎞ 2 − γi kt = 8 ⎛⎜ + ⎟ Gj ⎠ ⎝ Gi

d i dj ⎞ 1 3πμl [nij⋅(vp, i − vp, j)] ⎜⎛ nij ⎟ 2 d i ⎝ + dj ⎠ δn, ij

ηn = ηt = α kn

fct , ij = −kt δtij3/2 − ηt vtij

(9)

δn1/2

(11)

m i mj m i + mj

δn1/4

(12)

where α relates to the restitution coefficient. The elastic contribution of the impact energy absorbed during the compression is released during the restitution phase of the impact and leads to the elastic force that separates the contact partners. The absorption of kinetic energy during the impact can be described by a restitution coefficient. The coefficient of restitution is a ratio of the square root of the elastic strain energy released during the restitution to the impact energy, i.e. the initial kinetic energy (Tsuji et al., 1993).

where nij is the unit vector connecting the centers of particles with direction from i to j, and vp,i and vp,j are the velocity vectors of the two particles. This force opposes any change of the separation between the two particles: it is repulsive when the particles are moving toward each other and attractive when they are moving away from each other. Because the magnitude of the force diverges to infinity at zero separation, i.e. when the particles are in contact, the lubrication force, if acting alone, would prevent the contact of two particles. In reality however, additional forces act on the approaching particles and the particles are rough so they come in contact before the lubrication force reaches the theoretically predicted extreme values. To simulate the effect of roughness, a cut-off distance is used to prevent the lubrication force from reaching unrealistically high values. In case of a collision between two particles, the contact force fc is divided into normal (fcnij) and tangential (fctij) components, and these forces are modeled in Fig. 1. The elastic part of the normal contact force, fcn, is represented by a nonlinear spring, where the force is proportional to the stiffness, kn, and the displacement δn3/2 in the contact, i.e. the overlap. These two forces are given by the following equations (Kim and Karrila, 1991; Crowe et al., 1997): (8)

d i dj 2(di + dj )

(10)

whereγ is the Poisson ratio, and Ep and Gp are the longitude and transverse elastic moduli, respectively. To account for visco-elastic material properties that cause energy dissipation, a damping factor η related to the normal coefficient of restitution en is included into the DEM model, as proposed by Tsuji et al. (1993):

(7)

3/2 fcn, ij = (−kn δnij − ηn vrij nij ) nij

d i dj 2(di + dj )

−1

2

fl = −

2 1 − γ j2 ⎞ 4 ⎛ 1 − γi + ⎜ 3 Ei Ej ⎟ ⎝ ⎠

ln e

α=

⎧ 2 n2 π + ln en ⎨ 1 ⎩

0 < en ≤ 1.0 en = 0

(13)

2.3. Boundary conditions and numerical methods Fig. 2 (a) shows the rock crack links and overlaps. Microfractures in the overlap region connect large fracture segments that link and amalgamate into a single continuous macrofracture (Laubacha et al., 2004). Fig. 2 (b) represents the snapshot of a microfracture. Some rock particles distribute randomly in the microfracture, where these rock particles become fixed, and the fluid then flows around these fixed particles. In order to simplify the calculation in this study, the simulated porous structure in microchannel based on the structure is shown in Fig. 2(c). The fixed particles in the macrofracture represent the sandstone particles, while the sandstone particle sizes of sandstone reservoir by 80% range from 3.9 μm to 62.5 μm in China. As for Daqing oil field in China, the sandstone reservoir with high permeability and porosity, the sandstone sizes are from 8 μm to 50 μm. Pommer and Milliken (2015) introduced the components of detrital grains in pore medium. They found that carbonate skeletal debris was dominantly coccolithic and included intact coccospheres (typically 10–20 μm), disaggregated individual plates (1–10 μm in size), spines (10–20 μm in diameter), etc. Kondla et al. (2015) summarized the results of petrographic and geochemical analysis of Middle Triassic strata from the Sverdrup Basin in the Canadian Arctic. Grain size analysis of all pellet samples in the core revealed that the average grain size is fine silt ranging from 8 to 16 μm. Here, the porosity of microchannel depends on the distribution of fixed big rock particles, where fixed big rock particle sizes range from 8 μm to 20 μm. The fixed particles are randomly scattered in the microchannel in order to reflect the real structure of the porous media. Moreover, the arrangement of fixed big particles displays heterogeneity based on the FISH language built in PFC3d program, therefore, the gaps between fixed particles are different in the computational domain, and then the non-uniform void fraction in different parts of rocks are achieved, seeing Fig. 2(c). The computational domain with different porosities is obtained through changing the number of fixed big particles, as shown

where k and η are the spring and damping coefficients, respectively. δnij and δtij are the normal and tangential displacements between particle i and particle j, respectively.

Fig. 1. Model of contact force of two particles. 142

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

Fig. 2. Simulated geometry of microchannel in crock cracks. (a. crock cracks between two macrofractures. b. microchannel in rock cracks. c. partial enlarged drawing of simulated pore structure. d. simulated microchannel model).

particles and wall are set to −0.90 -0.95 and −0.95, respectively. At the inlet, the uniform distribution is given by superficial velocity of liquid phase, while the outlet is set to pressure outlet. Liquid is injected only in the axial direction. At an impenetrable wall, the tangential and normal velocities of liquid phase are set to zero (no slip condition). Initially, the sand grains of 500 randomly distribute in the inlet, where velocity of the sand grains is zero. Numerical simulations are performed by means of CFD-DEM code in which the governing equations of motion for the solid and liquid phases are simultaneously solved to obtain the velocity of the particle phase, velocity of the liquid phase, and local volume fraction. The equation of motion for the liquid phase is solved by the finite difference method using the software SIMPLE (Patankar, 1980). Newton's equations are solved algebraically. In the finite difference method, the flow domain is divided into cells, the size of which is smaller than the macroscopic motion of particles in the system but larger than the particle size. The differential three-dimensional equations of liquid motion are considered in an Eulerian framework, in which the fluid cells are fixed to a reference frame. The Lagrangian particle motion equations, taking into account the liquid-solid interaction, are simultaneously solved with the motion equations of the liquid phase to provide particle positions and velocities. In Eulerian-Lagrangian flow mapping, each fluid cell is considered to contain a group of particles interacting with the fluid. The position and velocity of individual particles are calculated using the Lagrangian equation of particle motion. The simulation results are shown in the forms of particle movement and velocity vectors. In addition with Lagrangian simulations, the tracking of tracer particles in the microchannel is used to describe the flow of particles. Based on this iteration, the CFD solver iterates over the next time step until the flow field again converges to a stable solution. Three dimensional particle flow program Particle Flow Code (PFC3D) is used to solve particle motions, which is commercial software. While the equations for the liquid are solve through adding user define functions based on FISH programming language. In this simulation, a constant time step of

Table 1 Parameters used in simulations. Particle diameter

0.004 (mm)

Particle density

Channel diameter Liquid velocity Liquid viscosity

channel length Liquid density Stiffness

Damping coefficients

0.2 (mm) 0.1/0.2/0.5 (mm/s) 0.001/0.005/0.009 (Pa s) 0.05

Time step

1 × 10−7s

Friction coefficient

2500 (kg/ m3) 1 (mm) 850 (kg/m3) 800 N/m 0.3

in Fig. 2(d), where the porous structure is no display. As for the 3D computational domain, the x-axial direction is the transportation direction, namely the axial direction of the cylindrical computational domain, and the y-axial and/or z-axial direction is radial direction. There are four kinds of three-dimensional microchannels with different porosities in rock cracks used in the present numerical simulations. The microchannel is 0.20 mm in diameter and 1.0 mm in length, where the fixed big rock particles distribute randomly across the microchannel. Generally, the fracture can be divided into large fracture (≥1 mm), middle fracture (0.5–1 mm), small fracture (0.1–0.5 mm), and microfracture (≤0.1 mm) based on the width of the fracture (Fang et al., 2010). Liu and Liu (2017) stated that the fracture dimensions were the width ranging from 0.01 to 2 mm, and the length in millimeter or centimeter level. In present study, the geometry of a small fracture is adopted, in which the typical hydrodynamic characteristics have been discussed in the process of waterflooding. As for the 3D computational domain, the x-axial direction is the transportation direction, namely the axial direction of the cylindrical computational domain, and the y-axial and/or z-axial direction is radial direction. The computational parameters are listed in Table 1. The grain of sand diameter and density are 2 μm and 2500 kg/m3. The restitution coefficients of the sand grains, fixed big particles and the wall are assumed to be 0.90, 0.95 and 0.95, respectively. The roughness coefficients of sand grains, fixed big 143

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

Fig. 3. Motions of particles in the microchannel at porosity of 0.55.

1.0 × 10−7 s is used. The CPU time required for the calculation of 1.0 s is about 8.0 h on a PC (1 TB hard disk, 8 GB RAM and of 2.4G Hz CPU × 2).

microchannel at the porosities of 0.55, 0.60, 0.65 and 0.70, respectively. The liquid velocity and viscosity are 0.2 mm/s and 0.001 Pa s. The effect of microchannel porosity on motion of representative particles is obvious. The representative particles are originally locate at the different positions due to the particles flowing randomly into the microchannel. At porosities of 0.55 and 0.60, the motions of representative particles are restricted, and some representative particles are stranded due to blocking by fixed big particles located in the microchannel. It can be apparently seen that the number of particles stranded in the microchannel at a low porosity is more than a high porosity. From these figures, it is evident that with the increase of porosity, the movements of representative particles become more free. When the porosity decreases to a certain value (ε = 0.60), the movements of some representative particles will be blocked and stranded within the microchannel. This implies that the trajectories of particles in the low porosity are fewer than those in the high porosity. It corresponds with the laws of oil migration as a lateral drainage to the nearest macrofracture zones (Lopatin et al., 2003). The reason lies in the low porosity, namely, the more fixed big particles are located in the microchannel. These barriers influence on the trajectories of particles, so the porosity of microchannels in the rock cracks is a key effect on migration of oil underground. The porosity of the microchannel also influences the time of sand grains crossing through the channel. The time consumed by moving particles from the entrance to the outlet of the model is defined as the

3. Simulation results and discussions 3.1. Flow behavior of grains of sand with different porosities The motion of the particles is fully controlled by various forces acting on it, including the fluid drag, lubrication force, particle-particle interactions and particle-wall interactions. Firstly, we focus on the flow pattern of particles in the microchannel. Fig. 3 describes the motion of particles in the microchannel at the porosity of 0.55. Particles are marked as red in color, and particles are traced in the whole transportation process. Initially, the red particles presented are introduced in the inlet. It can be observed that most of particles are pushed forward along the microchannel by means of the liquid, while some particles lag behind the main stream of particle as might be expected, and a few of particles are stranded in the pore. Observation of the trajectory of particles reveals that most of particles flow out from the bottom part of microchannel, indicating that the gravity of particles influents the motion mode of particles. In short, most particles in the microchannel move forward, and that a few of particles rest on the pore, which will modify the pore structure, and then it is adverse to the oil migration. Fig. 4 shows the trajectories of representative particles in the 144

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

Fig. 4. Trajectories of representative particles in the microchannel.

while 2.5% of particles begin to flow out of the channel during the residence time of 0.2–0.3s at the porosity of 0.7. During the residence time of 0.3–0.4s, the probability of particle number at the porosity of 0.6 is the largest, as more particles flow out of the channel. The next

resident time of the particles. In this study, eighty representative particles are traced to study the residence time. Fig. 5 shows probability distribution of particle number as a function of residence times. It is clearly seen that all particles stay in the microchannel in the first 0.2s, 145

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

fixed big particles. The velocity of sand grains will be reduced due to the energy dissipation by collisions of grains and fixed big particles during the residence times in the microchannel. The velocity of sand grains increases with an increasing of porosity, however, the residence time is increased at a porosity larger than 0.60. The residence time is the lowest at the porosity of 0.60 in present study. The increment in porosity means the reduction of fixed big particles, which leads to more space for free movement of grains, so the lateral movement of particles becomes large. It indicates that sand grains spend more time in reciprocating motion, therefore, the residence time of particles increase. That is to say, the residence time is increased due to long distance travelled under larger porosity, although the velocity of sand grains is increased. It also reveals that if an excess of sand grains stay for a long time in the microchannel, the porosity of microchannel will be decreased, resulting in an adverse effect on petroleum migration. Furthermore, sensitivity of computed residence time of particles to the spatial grid size is tested. The qualitative impact of the grid on the distribution of residence time of particles is shown in Fig. 6. Although the distributions of residence time of particles display the same trends, the difference is, however, obvious. The most significant discrepancy is embodied in underprediction from the coarse grids. In other words, an agreement between results from the medium grid sizes and finer grid sizes is found. It can be explained as the spatial grid size is so big that small pore structure is ignored. Therefore, the medium grid sizes are used throughout computations to save the computation expense. Both particle resident time and the spatial position of stranded particles are key factor to oil migration. Fig. 7 shows distributions of stranded particle fraction along the microchannel at different porosities. The high stranded particle fraction means more sand grains stay in the microchannel, and can not pass through the microchannel. These stranded sand grains will deposit in the channel, thus reducing the channel porosity. With an increase of porosity, the number of stranded sand gains is decreased, indicating that a large porosity contributes to the motion of sand grains and migration of oil. Furthermore, the low porosity leads to impedes oil migration in the loose sandstone oil reservoir. With the increase of porosity, more particles will be carried out of the channel by the carrying fluid, and the number of stranded grains will be reduced. Observation, the locations of stranded particles vary greatly with different porosities. In order to describe the spatial distribution of stranded particles, the length of microchannel is divided into three segments, i.e. 0.0–0.3, 0.3–0.7 and 0.7–1.0 mm, respectively. The distributed spatial scale of stranded particles is extended with a decrease of porosity. At the porosity of 0.55, the fraction of stranded

Fig. 5. Distribution of probability of particle residence time as a function of porosity.

0.2s, i.e. the residence time of 0.4–0.6s, the probability of particle number is in the range of 0.15–0.265. Within the residence time of 0.6–0.8 s, the probability of particle number at porosity of 0.55 is larger than others, while that at the porosity of 0.65 is the smallest. At the next 0.2 s, i.e. at the residence time of 0.8–1.0s, the numbers of particles flowing out of the channel at porosities of 0.65 and 0.70 are nearly equal, and their numbers are larger than others. In last 0.4 s, almost one third of particles at the porosity of 0.65 flows out. This indicates that the large porosity reduces flow resistance of grains of sand, while providing a free flow path for the sand grains. Fig. 4(d) clearly shows that the distance of motion for the particles is longer than other porosities. Roughly, both low porosity and high porosity can prolong the residence time sand grains in the microchannel. In addition, there are almost particles flowing out in every time period at the porosity of 0.70, while its histograms shows more homogenous, indicating that larger porosity in pore medium provides relative uniform flow characteristic. Fig. 6 shows the averaged residence time of particles as a function of porosity. The translational motion of particles controls the residence time in the channel. The residence time depends on two factors: sand grain velocity and distance travelled. With an increase of porosity, the residence time is decreased, then increased. The results show that the relationship between the porosity and residence time is nonlinear. At a low porosity of the microchannel, the collisional frequency of sand grains and fixed big particles is increased due to the large number of

Fig. 7. Fraction of stranded particle at different porosities along the microchannel.

Fig. 6. Particle residence time at three grid sizes as a function of porosity. 146

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

granular temperature is calculated as (Gidaspow, 1994):

θ=

1⎛1 3 ⎜N ⎝

N

∑ Cxj2 + j

1 N

N

∑ Cyj2 + j

1 N

N

⎞

∑ Czj2 ⎟ j

⎠

(14)

where Cj is the random fluctuation velocity of the j sand grain defined as the difference between the instantaneous velocity and the local mean velocity, and N is the total number of sand grains in the computational cell. Fig. 11 shows the distribution of granular temperature as a function of porosity. The granular temperature of grains decreases, reaches minimum, and then increases with increasing porosity. At a low porosity, the particle collisions play an important role in the channel. The granular temperature decreases due to the decrease of the collisions of grains and the larger fixed particles. The granular temperature increases due to the increase of the mean free path at high porosity. When the porosity is equal to 0.7, the granular temperature is almost independent upon porosity due to the collision and free path tending to steady. This indicates that the fluctuation of grains is controlled by microchannel geometry. At a low porosity, the fluctuations in grain movement are dependent on the collisional probability between grains, or between grains and fixed particles. At a high porosity, the activity of grains is a key factor on granular temperature. Therefore, both high frequency of collisions and more activity of grains can improve granular temperature.

Fig. 8. Instantaneous axial velocity of a representative particle as a function of porosity.

particles is 9.0%, which distribute the whole length, otherwise, at the porosity of 0.70, the stranded particles are available in the first two segments, and the fraction is only 1.4%. Furthermore, the stranded particles distribute more uniform as the porosity decreases along the microchannel length. Fig. 8 shows the instantaneous axial velocity of a representative particle as a function of time at four porosities. The axial velocity of representative particles oscillates between positive and negative. The positive value means the representative particles proceed through the channel by liquid phase carrying, while the negative value means the representative particles move against the flow. This indicates that the fixed big particles in the microchannel block the small particles moving freely in the channel. With the increase of porosity, the simulated fluctuation of axial velocity of the representative particles is reduced. It means that the collisions between sand grains and fixed particles influence the velocities and path of sand grains. In order to explore the particle flow characteristic further, Fig. 9 presents instantaneous particle velocity vector at the porosity of 0.70. It can be seen that particles mainly flow forward along the axial direction of the channel, and particles pass through the small opening in pore structure. It is evident that the axial velocity of particles is far greater than lateral velocity, so the translational motion of particles dominates over the entire movement process in the channel. The contact forces are caused by collisions between two sand grains, and sand grains and fixed particles. The contact forces depend on the relative velocity between two particles and the parameters of stiffness, dissipation, and friction coefficients. It includes normal (fcn) and tangential (fct) forces. Fig. 10 shows the normal and tangential forces for a sampled sand grain (No. 19) as a function of porosities. The contact forces fluctuate violently at the low porosity, because the collision probability increases with an increase of porosity. The lack of free path and long residence time to small particles brings out an increase of collision probability. The mean value (MV) and standard deviation (SD) of contact forces are calculated. The normal component of contact force decreases from 1.42 × 10−4 to 1.33 × 10−6 N, and the tangential component of contact force is from 3.99 × 10−7 to 3.44 × 10−8 N at the porosity of 0.55–0.70. This implies that the contact force decreases with increasing porosity, the normal component of contact force is significantly larger than tangential component. Granular temperature, in analogy with the temperature in the kinetic theory of gases, is a measure of the intensity of particle collisions and depends on the fluctuating velocity of particles. From simulated instantaneous velocity of sand grains, the granular temperature is predicted in the channel. In a three-dimensional system, the mean

3.2. Flow behavior of particles with different liquid viscosities The liquid viscosity is changed with the depth of rocks because the property of fluid depends upon pressure and temperature. The high liquid viscosity results in more energy loss through the microchannel, and the velocity of fluid will be changed. The variation of the liquid viscosity has the potential to change the drag force on the particles. Fig. 12 shows the distributions of axial velocity of representative particles at the porosity of 0.55 and liquid velocity of 0.2 mm/s with different viscosities of 0.001, 0.005 and 0.009 Pa s, respectively. It can be seen that the axial velocity of sand grains is negative or positive. The axial velocity of sand grains decreases with an increase of liquid viscosity. This implies that the high liquid viscosity slow the motion of sand grains, furthermore, the loss in energy is increased with an increase of liquid viscosity. Fig. 13 shows the granular temperature as a function of the liquid viscosity. The granular temperature is increased with decreasing liquid viscosity. With the increase of liquid viscosity, more energy of liquid phase is dissipated. Thus, the fluctuation velocity of particles is reduced, and granular temperature is decreased. Fig. 14 shows the averaged residence time of sand grains as a function of liquid viscosity. It is clearly seen that the averaged resident time rises with an increase of liquid viscosity. With the increase of liquid viscosity, more energy is lost due to viscous dissipation. This will slow particles down, and increase residence time of particles in the microchannel. It is the evidence that the liquid viscosity is characterized by influencing on flow of sand grains, namely prolonging the residence time. The motions of particles in the liquid, mainly experience two types of forces, i.e. drag force exerted by the fluid along the fluid streamlines, and collisions with other grains and walls. The grain collisional effects are dominant over the fluid viscous effects. The instantaneous normal and tangential contact forces as a function of liquid viscosities are shown in Fig. 15. Both normal contact forces and tangential forces are increased with an increase of liquid viscosity. The increase of liquid viscosity leads to an increase of grain interaction, and then the contact forces are strengthened. As the liquid viscosity changes from 0.001 to 0.009Pa s, the normal contact force increases from 2.22 × 10−5 to1.42 × 10−4 N, and the tangential force increases from 5.44 × 10−8 to 3.99 × 10−7 N. The normal contact force is far greater than the tangential force. It indicates that the motion of particles is mainly 147

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

Fig. 9. Instantaneous particle velocity vector at the porosity of 0.70.

Fig. 10. Profiles of contact forces of a representative particle as a function of porosity. 148

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

Fig. 11. Distribution of granular temperature as a function of porosity.

Fig. 14. Particle residence time as a function of liquid viscosity.

field, which drives san grains acceleration in the microchannel. Fig. 16 shows the instantaneous axial velocity of representative particles as a function of time at three liquid velocities of 0.1, 0.2 and 0.3 mm/s, respectively. Simulations indicate that the axial velocities of sand grains are varied in positive and negative. This means that the particles undergo oscillatory motion in the channel due to the collision of sand grains and fixed big particles with carrying of the liquid phase. It is clearly observed that the axial velocity of sand grains increases with increasing liquid velocity. The magnitude of negative velocity is lower than that of positive velocity in three cases of different liquid velocities, therefore the particles can pass through the channel with the carrying liquid. It is important to reduce stay for small particles in the channel, because excessive particles in the channel leads to flowing space reduction and thus reduces the migration path of oil. Fig. 17 shows the distribution of averaged granular temperature as a function of liquid velocity. The granular temperature increased with an increase of liquid velocity. It indicates that an increase in liquid velocity corresponds to faster motion of the particles, while the frequency of collisions between sand grains and fixed big particles is increased, and then the oscillations of particles are increased. The granular temperature is sensitive to the liquid velocity, where the granular temperatures are 0.0147, 0.0498 and 0.0625(cm/s)2 at the liquid velocity of 0.1, 0.2 and 0.5 mm/s. The distribution of residence time of sand grains is shown in Fig. 18 as a function of liquid velocity. The residence time of sand grains decreases with an increase of liquid velocity. From Fig. 16, we can see that the axial velocity varies between positive and negative, meaning that the grains move back and forth. Therefore the travel distance of sand grains becomes long. The high residence time of sand grains attributes to collision of particles, and reduces the migration of oil. Fig. 19 shows the instantaneous normal and tangential contact forces as a function of liquid velocity. The relative velocity between the liquid phase and sand grains increases with an increase of liquid velocity, resulting in the grains suffering from large drag forces at high liquid velocity. The motion of particles is dominated by the drag force due to the liquid carrying. As the velocity of particles increases, so do the contact forces.

Fig. 12. Instantaneous axial velocity of a representative particle with different viscosities.

Fig. 13. Distribution of granular temperature as a function of liquid viscosity.

4. Conclusions

dominated by collision of sand grains in the microchannel.

Flow behavior of grains of sand in microchannel rock cracks is simulated using CFD-DEM. The motions of grains of sand are predicted. Flow behavior of sand grains shows that the motions of sand grains changing continuously their positions and residence times within the microchannel due to collisions of sand grains and fixed big particles. Simulated results show that particle movement increases with

3.3. Flow behavior of particles with different liquid velocities The increase in liquid velocity enhances the energy input to the flow 149

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

Fig. 15. Profiles of contact forces of a representative particle as a function of liquid viscosity.

Fig. 18. Particle residence time as a function of liquid velocity.

Fig. 16. Instantaneous axial velocity of particles with different liquid velocities.

increasing microchannel porosity. While the porosity decreases to a certain value, the motion of representative particles will be blocked. With an increase in porosity of the microchannel, the residence time of sand grains is decreased, and then increased, while the mean granular temperature is in reverse. The simulated fluctuation of axial velocity and the fluctuation of liquid velocity are reduced with an increase of porosity, while the mean value (MV) and standard deviation (SD) of the contact forces are increased with an increase of porosity. The axial velocity and granular temperature of grains decrease with increasing liquid viscosity, while the mean residence time and contact force increase. The granular temperature and the contact force are increased with increasing liquid velocity, while the residence time of sand grains decreases. Note that present simulations are limited by the number and shape of grains of sand in the system and geometry of the microchannel in the rock cracks. This makes it difficult to make quantitative comparisons with experiments. Detailed simulations should provide better quantitative answers to the question of how grains of sand affect the migration of oil underground. In addition, the geometry of the microchannel should be included different shapes, sizes and arrange mode of fixed

Fig. 17. Distribution of granular temperature as a function of liquid velocity.

150

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

Fig. 19. Profiles of contact forces of a representative particle as a function of liquid velocity.

rock particles, which will be closer to the actual situation. Therefore, a thorough investigation is needed by numerical simulations and experiments, which is left as future work. Acknowledgments This work was supported by the National Science Foundation in China through grant (No.21676051) and State Key Laboratory of Heavy Oil Processing of China (2017). We also thank Dr. Matthew R. Girardi at Princeton University for polishing this manuscript. Notation Cd drag coefficient E modulus of longitudinal elasticity, N/m2 fc contact force, N fd drag force, N fl lubrication force, N fm virtual mass force, N mi mass of particle, kg g gravitational acceleration, m/s2 G modulus of transverse elasticity, N/m2 I moment of inertia, kg/m2 kn spring constant for normal direction, N/(m3/2) kt spring constant for tangential direction, N/m P liquid pressure, N/m2 T torque, N/m ul liquid velocity, m/s up solid velocity, m/s Vp volume, m3 Greek letters ρl εl τl μ β ω εp κ

γ δ η μf θ Subscripts

Poisson's ratio displacement, m damping coefficient, kg/s friction coefficient granular temperature, m2/s2

l p w n t x y z

liquid phase particle phase wall normal direction tangential direction x direction y direction z direction

References Cao, J., Jin, Z.J., Hu, W.X., Zhang, Y.J., Yao, S.P., Wang, X.L., Zhang, Y.Q., Tang, Y., 2010. Improved understanding of petroleum migration history in the Hongche fault zone, northwestern Junggar Basin (northwest China): constrained by vein-calcite fluid inclusions and trace elements. Mar. Petrol. Geol. 27, 61–68. Crandall, D., Bromhal, G., Zuleima, T.K., 2010. Numerical simulations examining the relationship between wall-roughness and fluid flow in rock fractures. Int. J. Rock Mech. Min. Sci. 47, 784–796. Crowe, C.T., Sommerfeld, M., Tsuji, Y., 1997. Multiphase Flows with Droplets and Particles. CRC Press. Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Geotechnique 29, 47–65. Deen, N.G., Van Sint Annaland, M., Van der Hoef, M.A., Kuipers, J.A.M., 2007. Review of discrete particle modeling of fluidized beds. Chem. Eng. Sci. 62, 28–44. Fang, J., Jia, C.Z., Ran, Q.Q., Zhou, X.M., Hui, G., 2010. The Study of Fracture in Volcanic Reservoir, vol. 3. Society of Petroleum Engineers – IOGCEC, pp. 2415–2420. Fatahi, H., Hossain, M., 2016. Fluid flow through porous media using distinct element based numerical method. J. Petrol. Explor. Prod. Technol. 6, 217–242. Ferreira, V.H.S., Moreno, R.B.Z.L., 2017. Impact of flow rate variation in dynamic properties of a terpolymer in sandstone. J. Petrol. Sci. Eng. 157, 737–746. Gidaspow, D., 1994. Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions. Academic Press, London. Grof, Z., Cook, J., Lawrence, C.J., Štěpánek, F., 2009. The interaction between small clusters of cohesive particles and laminar flow: coupled DEM/CFD approach. J. Petrol. Sci. Eng. 66, 24–32. Hatzignatiou, D.G., Norris, U.L., Stavland, A., 2013. Core-scale simulation of polymer flow through porous media. J. Petrol. Sci. Eng. 108, 137–150. Hoomans, B.P.B., Kuipers, J.A.M., Briels, W.J., Van Swaaij, W.P.M., 1996. Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidized bed: a hard-sphere approach. Chem. Eng. Sci. 51, 99–118.

density of liquid, kg/m3 porosity viscous stress tensor viscosity of liquid, kg/ms interface momentum transfer coefficient, kg/m2 s2 angular velocity, 1/s concentration of particles surface tension 151

Journal of Petroleum Science and Engineering 171 (2018) 140–152

S. Wang et al.

Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. McGraw-Hill, New York. Pommer, M., Milliken, K., 2015. Pore types and pore-size distributions across thermal maturity. Eagle Ford Formation, southern Texas. AAPG (Am. Assoc. Pet. Geol.) Bull. 99, 1713–1744. Qiu, L.C., 2015. A coupling model of DEM and LBM for fluid flow through porous media. In: New Paradigm of Particle Science and Technology-proceedings of the 7th World Congress on Particle Technology. 102. pp. 1520–1525. Qiu, L.C., Wu, C.Y., 2014. A hybrid DEM/CFD approach for solid-liquid flows. J. Hydrodyn. 26, 19–25. Santanna, V.C., Silva, A.C.M., Lopes, H.M., Sampaio Neto, F.A., 2013. Microemulsion flow in porous medium for enhanced oil recovery. J. Petrol. Sci. Eng. 105, 116–120. Slider, H., 1976. Practical Petroleum Reservoir Engineering. (Tulsa). Tannenbaum, E., Huiainga, B., Kaplan, I., 1986. Role of minerals in thermal alteration of organic matter. II—a material balance. AAPG (Am. Assoc. Pet. Geol.) Bull. 70, 1156–1165. Tsuji, Y., 2007. Multi-scale modeling of dense phase gas-particle flow. Chem. Eng. Sci. 62, 3410–3418. Tsuji, Y., Kawaguchi, T., Tanaka, T., 1993. Discrete particle simulation of two dimensional fluidized bed. Powder Technol. 77, 79–87. Wachs, A., 2009. A DEM-DLM/FD method for direct numerical simulation of particulate flows: sedimentation of polygonal isometric particles in a Newtonian fluid with collisions. Comput. Fluids 38, 1608–1628. Washino, K., Tan, H.S., Hounslow, M.J., Salman, A.D., 2013. A new capillary force model implemented in micro-scale CFD–DEM coupling for wet granulation. Chem. Eng. Sci. 93, 197–205. Wu, L.M., Zhou, C.H., Keeling, J., Tong, D.S., Yu, W.H., 2012. Towards an understanding of the role of clay minerals in crude oil formation, migration and accumulation. Earth Sci. Rev. 115, 373–386. Xu, B.H., Yu, A.B., 1997. Numerical simulation of the gas solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics. Chem. Eng. Sci. 52, 2785–2809. Zhu, H.P., Zhou, Z.Y., Yang, R.Y., Yu, A.B., 2007. Discrete particle simulation of particulate systems: theoretical developments. Chem. Eng. Sci. 62, 3378–3396.

Ilyas, K., Hazim, A., Abdulrahman, A., 2018. Modeling cementation in porous media during waterflooding: asphaltene deposition, formation dissolution and their cementation. J. Petrol. Sci. Eng. 161, 359–367. Kang, Y.F., Yu, M.J., Miska, S., Takach, N.E., 2009. Wellbore stability: a critical review and introduction to DEM. In: Proceedings - SPE Annual Technical Conference and Exhibition. 4. pp. 2689–2712. Kim, S., Karrila, S.J., 1991. Microhydrodynamics: Principles and Selected Applications. Butterworth-Heinman Series in Chemical Engineering, Stoneham, MA, USA. Kondla, D., Sanei, H., Embry, A., Ardakani, O.H., Clarkson, C.R., 2015. Depositional environment and hydrocarbon potential of the middle triassic strata of the Sverdrup basin, Canada. Int. J. Coal Geol. 147–148, 71–84. Laubacha, S.E., Reeda, R.M., Olsonb, J.E., Lander, R.H., Bonnella, L.M., 2004. Coevolution of crack-seal texture and fracture porosity in sedimentary rocks: cathodoluminescence observations of regional fractures. J. Struct. Geol. 26, 967–982. Liu, J.D., Liu, G.X., 2017. Types, evolution of fractures, and their relationship with oil and gas migration of permian changxing formation in yuanba gas field of Sichuan basin, China. Energy Fuels 31, 9345–9355. Lopatin, N.V., Emets, T.P., Romanovb, E.A., Mal'chikhina, O.V., 2003. Petroleum migration and trapping mechanism in fine-grained source rock formation. J. Geochem. Explor. 78–79, 395–397. Lu, H.L., He, Y.Y., Liu, W.T., Ding, J.M., Gidaspow, D., Bouillard, J., 2004. Computer simulations of gas–solid flow in spouted beds using kinetic-frictional stress model of granular flow. Chem. Eng. Sci. 59, 865–878. McAuliffe, C.D., 1979. Oil and gas migration-chemical and physical constraints. AAPG (Am. Assoc. Pet. Geol.) Bull. 63, 761–781. Moghadasi, J., Müller-Steinhagen, H., Jamialahmadi, M., Sharif, A., 2004. Theoretical and experimental study of particle movement and deposition in porous media during water injection. J. Petrol. Sci. Eng. 43, 163–181. Parnell, J., Carey, P.F., Monson, B., 1996. Fluid inclusion constraints on temperatures of petroleum migration from authigenic quartz in bitumen veins. Chem. Geol. 129, 217–226. Parra, F.R., 2004. Oil Politics: a Modern History of Petroleum. L.B. Tauris & Co., Ltd. Press, New York, America.

152