Conductivity analysis of hydraulic fractures filled with nonspherical proppants in tight oil reservoir

Hydraulic fracturing is an effective way to exploit hydrocarbons in tight oil reservoir. The success of this technique lies on the creation of high conductivity channels for the oil and gas. Quantitative analysis of factors influencing the conductivity of propped fractures is significant for the design of the operation schedule. According to the definition of sphericity, cylindrical and planar proppants of different sphericity and sizes were presented. And disordered packs of nonspherical proppants in fractures of different widths were established by the discrete element method. The change of the proppant embedment and fracture width under different closure pressures were analyzed by the solid mechanics. During the fracture deformation, the pressure and velocity of the fluid flowing in fractures were simulated by Lattice‐Boltzmann method. Based on that, the effect of the proppant size, sphericity, and the fracture width on the permeability and conductivity of hydraulic fractures were investigated. The accuracy of the model was verified by the experimental data. The simulation results show that the fracture permeability and conductivity decrease with the decrease in the rock's Young's modulus, proppant size, and sphericity, and their stress sensitivity increases with the decrease in the rock's Young's modulus and the increase in the proppant size. Increasing fracture width can improve fracture conductivity more significantly than increasing fracture permeability. The permeability and conductivity of the fractures filled cylindrical proppants are higher than that of the fractures filled with planar proppants.

propped fractures can optimize the fracturing process. On the other hand, it is also the basis of oil well productivity simulation. The previous research also revealed that the stimulated volume, the conductivity of the hydraulic fractures, and the permeability of the natural fractures are the main factors influencing the development of unconventional reservoirs. [1][2][3][4][5][6] However, the transport distance of proppants is far more less than we expected, which is demonstrated by the core analysis obtained from the fractured horizontal wells. [7][8][9][10] This phenomenon lays more emphasis on the fracture conductivity.
The essential reasons for the reduction of fracture conductivity are fracture deformation, proppant embedment, and their shape variation. 11,12 However, the mechanism becomes more complex under different situations. A number of experimental and theoretical studies have been conducted to determine the relationship between the fracture flow capacity and multiple factors, such as proppant concentration and diameter, closure pressure, and temperature. Cooke established a mathematical model to estimate the variation of fracture conductivity, which proved that the proppant concentration, the fluid residue, and fracture porosity had great contribution to the fracture conductivity reduction. 13 Lacy discovered that closure pressure, proppant size, and concentration were the main factors determining the proppant embedment by a series of experiments. 14, 15 Fredd also investigated the influence of proppant concentration and strength, and properties of formation rocks on the conductivity of aligned and displaced fractures by experiments. 16 With the emerge of unconventional formations, the sparse distribution of proppants in fracture networks draws much attention from researchers. Khanna et al 17 used the elastic mechanics to calculate the narrow fracture conductivity, which obtained the relationship between the optimal concentration of proppant and the confining stress. Guo et al 18 studied the embedment and conductivity of proppant packs in sparse distribution with the contact mechanics and Carman-Kozeny formula, while Liu et al 19,20 made modification to Carman-Kozeny equation due to its deviation in the conductivity calculation of the monolayer proppant. Another factor contributing to the conductivity decline is the proppant crushing. Zheng et al [21][22][23] performed the investigation on the breakage characteristics of different proppants, which could provide more accurate prediction of fracture conductivity under high confining stress.
Based on the existing knowledge, many scholars are committed to deceasing the proppant embedment and crush, and proppants treated by SMA (surface modification agent) are consequently used extensively. An analytical model with adjustment coefficients added to the previous one was derived to compute embedment and conductivity of this kind of proppants, which also treated the proppant and formation rock as the isotropic elastomer. 24, 25 Shamsi et al simulated the dynamic conductivity of fracture filled with proppants graded differently using the Lattice-Boltzmann method. However, he did not consider the deformation of the fracture surface and proppants under the confining stress. 26,27 Recently, rod-shaped proppant has been increasingly used in the hydraulic fracturing with the alternate-slug pumping schedule. Osiptsov established a random packing model of rod-shaped proppant using the sequential deposition method, which demonstrated that this shape of proppants provided the highest permeability with the Lattice-Boltzmann method. Also, the fluid flowing in the fracture was effected by the distribution of proppants. 28 Tan designed six cases in experiments to estimate the influence of the proppant layer and type on permeability of fracture with and without proppants. 29 Wang also studied the fracture conductivity under different proppant distributions. Fractures in this research were partially and fully filled with proppants. Fracture width and proppant embedment were investigated by elastic mechanics. 30 Most of these studies, especially the theoretical investigations, treat proppants as standard spheres. And proppant particles are orderly packed in the fracture. Those assumptions are inconsistent with the real condition, which might lead to large errors in the prediction.
In this paper, nonspherical proppant particles were constructed by several spheres of equal diameters. Then, disordered packs of proppants were built with discrete element method. A coupling model of elastic mechanics and the Lattice-Boltzmann method was used to simulate the deformation of fracture surface and the fluid flowing in the fracture. Ultimately, effects of proppant size and sphericity, fracture width, and formation properties on the fracture width, permeability, and conductivity were analyzed.

| Packing model of nonspherical proppants
Proppants used in the hydraulic fracturing are mainly quartz sand and artificial ceramsite. These proppants are naturally nonspherical particles. Many parameters are adopted to quantitatively describe the shape of nonspherical particles, for example, aspect ratio or elongation, 31-33 eccentricity ratio, 34 angularity, 35 and roughness. 36 Generally, sphericity is used to describe the similarity between proppant and sphere in the petroleum industry. This parameter is one of the bases for correct selection of proppants. There are many definitions of sphericity, such as Wadell sphericity, 37,38 which is defined as follows where S is the sphericity; A e is the area of sphere equal to the volume of proppant, mm 2 ; and A p is the area of proppant, mm 2 .
Krumbein 39 also proposed a method to determine the particle sphericity with the long, intermediate, and short axes. Two parameters were introduced. They were the ratio of the intermediate to the long diameter and the ratio of the short to the intermediate diameter. The sphericity could be obtained by the chart with these two ratios as the axes.
There is also a definition in the standard for petroleum and natural gas industry of China (SY/T 5108-1997). The proppant sphericity is defined as: where d e is diameter of isovolumetric sphere of proppant, mm; d p is the circumsphere diameter of proppant, mm. This definition is adopted in this study as it is easier for us to build the nonspherical proppants.
Several spheres of equal diameters are used to construct proppants of different sphericity, as shown in Figure 1. The three-dimensional graphs on the left of the figure consist of two parts: (a) some equal-diameter spheres in the interior represent the proppant; (b) the exterior sphere is the circumscribed sphere of the proppant. Whether the proppant can pass through a mesh depends on its maximum length, that is, the diameter of the external sphere. The graphs on the right are the cross-section of the proppant. By adjusting the spherical center distance (L AB ) and radius of the sphere (r), the proppant of different sphericity can be constructed. The proppant consisting of two equal-diameter spheres is called cylindrical proppant. The proppant consisting of four equal-diameter spheres is called planar proppant. The specific parameters of different types of proppants are shown in Table 1.
To build a reasonable proppant packing model, the discrete element method is a solid choice due to its good applicability for the motion and accumulation simulation of proppant particle system. [40][41][42][43] This method disperses the whole research object into discrete, separable, and independent particles. By analyzing the kinematic and dynamic equations of particles in the system, the position and velocity information of each particle are determined. The evolution process of the whole studying object is described by the behavior of the particle system. The specific research process is as follows: firstly, the resultant force and resultant moment (including the collision between particles and the field force between particles and physical field) of particles are analyzed. The acceleration of particles is obtained by Newton's second law. The velocity and displacement of particles are obtained by time integration of particle acceleration, and then, the motion state of particles at any time is obtained; the movement of each particle is calculated cyclically until equilibrium is reached. The relation of force and displacement is applied at each contact point of particles, which is defined by a unit normal vector, n par . When particles contact with each other, there is the contact force. This force is calculated by the relative displacement and specified stiffness of contacted particles. Normally, the contact force is derived into normal and shear components with respect to the contact plane, as shown in Equation 3. The normal contact force and shear contact force could be computed by Equations 4 and 5, respectively.
where F par is the contact force of each particle, N; F n par is the normal contact force of each particle, N; F s par is the shear contact force of each particle, N; K n is the normal stiffness, N/m; k s is the shear stiffness, N/m; U n is the contact displacement in the normal direction, mm; U s par is the contact displacement of the shear direction, mm; and n par is the normal direction of each particle.
Three different proppant packs and their corresponding distribution of particle size are illustrated in Figure 2. In this study, there are nine kinds of proppant packs built to analyze the variation of fracture width, permeability, and conductivity, as shown in Table 2. The initial porosity and specific surface area are also presented in Table 2. The fracture porosity was determined by the ratio of the pore volume in fracture to the total volume of fracture. The specific surface area was obtained by the ratio of pore surface area to the pore volume in fracture. The volume and surface area of pore and the fracture volume can be measured in the simulation.

| Rock deformation model
Under the closure pressure, the embedment of proppants into the fracture surface leads to the decrease in fracture width and porosity in the fracture, which further causes the decline of fracture conductivity. In order to study the proppant embedment, the following assumptions are made: 1. There is no initial stress and strain in proppants and rock matrix. 2. The rock matrix is isotropic linear elastomer. 3. Since the ceramic proppant has higher Young's modulus than the rock, the deformation of proppant was smaller than the rock. The movement and deformation of proppants are neglected. 4. The effect of fracturing fluid filtration on rock mechanical properties is ignored. 5. The whole research process is isothermal.
Based on hypothesis (2), the stress-displacement relationship of rock matrix with homogeneous linear elastomer under external force is as follows: where ρ s is the rock density, kg/m 3 ; u is the rock displacement tensor, mm; σ is stress tensor; F v is the external force of per volume, N/m 3 ; and t is the time, seconds. The stress-strain relationship and the strain-displacement relationship are shown in Equations 7 and 8: where, ε is the strain tensor; C is the stiffness matrix. For the three-dimensional stress-strain relationship, the stiffness matrix is a tensor of 6 orders, as follows: where, E is the Young modulus, GPa; υ is Poisson's ratio.
It should be pointed out that in the process of fluid-solid coupling, the closing pressure on the proppant is the net pressure of in situ stress and fluid pressure.
The displacement of the fracture surface under the closure pressure is uneven due to the uneven distribution of proppants in the fracture. The area-weighted average displacement is used to express the change of fracture width.
where w a is the average fracture width under closure pressure, mm; w 0 is the initial fracture width, mm; u f is the displacement of fracture surface, mm; and A f is the area of fracture surface, mm 2 .

| Lattice-Boltzmann method
The traditional CFD method used to solve Navier-Stokes equation was based on continuum assumption, which is difficult and complex. The Lattice-Boltzmann method discretizes fluid into particles and establishes motion equation. [44][45][46] The physical quantities described by Lattice-Boltzmann method satisfy = C ⋅ , Navier-Stokes equation from macroscopic point of view, and it has certain advantages in modeling complex geometries and boundary conditions due to its discrete nature. 45,47,48 It has been successfully applied to fluid flow simulation in porous media. The fracture permeability was calculated by Darcy's law, which requires that fluid flowing must be laminar in the fracture. And the fluid flowing in the conductivity experiment is also laminar required by the experimental standard. Therefore, the fluid flowing in the simulation is laminar. Single relaxation model (LBGK model) and multiple relaxation model (MRT model) are commonly used LBM models. 49 In this paper, LBGK model is adopted and its evolution equation can be expressed as follows: where f i (x,t) is particle distribution function; Δt is time step, seconds; e i is discrete velocity, m/s; τ is dimensionless relaxation time; f eq i (x,t) is equilibrium state distribution function of particles; and i = 0, 1, …, m is discrete velocity direction with different lattice points.
Particle equilibrium distribution function is defined as: where ω i is weight coefficients in direction i; ρ l is density of fluid lattice points, kg/m 3 ; e s is lattice sound velocity, m/s; and v l is velocity of fluid lattice, m/s. Because this study is a three-dimensional flow model, the D3Q19 model is chosen to simulate the fluid flowing in the pore area between proppants. Its discrete velocity is as follows: where e = Δx/Δt is lattice velocity, e s 2 = e 2 /3.
For 19 discrete directions, the weight coefficients can be expressed as: The macroscopic density and velocity of fluids are: During the flow of crude oil in hydraulic fractures, the flow resistance of fluid on the surface of proppant is pretty large, so the standard rebound scheme without slip boundary is adopted as the boundary conditions.
For the flowing of a single-phase fluid in the proppantfilled fracture, the permeability can be obtained by Darcy's law: where K is absolute permeability, μm 2 ; μ is the viscosity of the fluid, mPa·s; v a is the average seepage velocity of fluid, m/s; and dp/dx is the average pressure gradient, kPa/m.
The fracture conductivity can be expressed as: where K d is the conductivity of hydraulic fractures, μm 2 ·cm.

| MODEL VALIDATION
In order to verify the accuracy of this model, fracture conductivity tests were conducted on the experimental equipment named ZCJ-200, and the simulation results were compared with the experimental data. The experiments were conducted according to the standards of the petroleum and natural gas industry of the People's Republic of China (SY/T 6302-2009). The sandstone rock sample and crude oil obtained from a tight oil reservoir in China were used. The static Young's modulus and Poisson's ratio of the rock were 23.4 GPa and 0.235, respectively. The viscosity and density of the oil were 293.8 mPa·s and 888 kg/ m 3 at the room temperature. The experiments were also conducted at that temperature. During the experiment, the closure pressure was loaded at the rate of 3500 kPa/min until the required closure pressure is reached, which could ensure that the proppant reached a semistable state under the closed pressure. Three different flow rates were tested, and the average of the experimental results was taken as fracture permeability. The specific physical parameters of the proppant are (11) (0,0,0) e i = 0 ( ± 1,0,0)e,(0, ± 1,0)e,(0,0, ± 1)e i= 1,2,3,4,5,6 ( ± 1, ± 1,0)e,( ± 1,0, ± 1)e,(0, ± 1, ± 1)e i= 7, 8,9,10,11,12,13,14,15,16,17,18 , 8,9,10,11,12,13,14,15,16,17,18 , Table 3, and the sphericity of all proppants was greater than 0.8.

as shown in
Three proppant packs of different particle sizes (20/40 mesh, 30/50 mesh, 40/70 mesh) were modeled in the simulation, as shown in Figure 2. The proppant particle sizes of all proppant packs in the simulation were randomly distributed. The initial and boundary conditions applied to the model are as shown in Figure 3. Under the initial conditions, there was no prestress in the rocks and proppants, and the fluid in the fracture did not flow. The closure pressures (10, 30, 50, and 70 MPa) were, respectively, loaded on both sides of the rock, and the proppant was fixed during the compression process without considering deformation. While the closure pressure was applied, the fluid entered the fractures filled with the proppants from the inlet and flowed out from the outlet. The fluid and rock properties were the same as those in the experiment. The flowing rate and pressure of fluid in the inlet were 5 cm 3 / min and 1 MPa.
The simulation and experimental conductivity data of proppants with three particle sizes are as shown in Figure 4. It can be seen that the data range and the changing trend of both data are basically consistent. Besides the ignorance of proppant deformation and breakage, proppants of the same sphericity were used in each proppant pack and the real shape of proppants may not be limited to be columnar or planar in the experiment. All these factors account for the error in the simulation. The overall error of the numerical simulation is 7.4%, which is acceptable.

| RESULTS AND DISCUSSIONS
The fluid-solid coupling model established above was used to analyze the influence of various factors such as Young's modulus of the rock, proppant particle size, proppant sphericity, and proppant shape on the width, permeability, and conductivity of the fracture filled with nonspherical proppants under different closure pressures. The control variable method was used to simulate the influence of one variable while keeping other variables unchanged. The proppants packs in the simulation are shown in Table 2.

| Fracture aperture
Under the closure pressure, the fracture surface moves toward inside of the fracture, which leads to the gradual reduction in the fracture width. The simulated displacement of fracture surface is shown in Figure 5. For the convenience of observation, this figure shows a horizontal cross-section of the  The deformation of the fracture surface is mainly effected by the rock's Young's modulus, proppant particle size, proppant sphericity, and proppant shape, as shown in Figure 6. It should be noted that Figure 6A shows the simulation result of 30/50 mesh proppant with a sphericity of 0.866 in the fracture whose width is 2 mm, and the data in the other figures are simulated under the closure pressure of 70 MPa.
It can be seen from Figure 6A that under the same closure pressure, as the Young's modulus of the rocks decreases, it is easier for the proppants to get embedded into the fracture. Also, with the increase in closure pressure, the fracture aperture decreases faster, and this trend becomes more obvious. For rocks with the same Young's modulus, the reduction in the fracture width is basically the same for increasing every 20 MPa in the closure pressure, so the relationship between closure pressure and the reduction of the fracture width (average embedment of proppants) is basically linear. Figure 6B shows the influence of proppant size on the fracture width. Under the same closure stress, the fracture width decreases with the increase in the proppant particle size, but the decreasing amplitude is getting smaller. The trend becomes more prominent with the decrease in the Young's modulus of the rocks. This is because the larger particle size of the proppant, the less proppant is filled in the fracture of the same fracture width, and then, the larger the pores in the fractures are, the worse the supporting effect is. The width changes of fractures filled with proppants of different sphericity are as shown in Figure 6C. The larger the proppant sphericity is, the worse the supporting effect is. When the sphericity increases to some value (0.866 in the simulation), it tends to be substantially stable, because when the shapes of proppants are more irregular, they are more closely packed, and the contact area with the fracture surface is larger. For the two geomorphic proppants constructed in this study, the fracture width of the fractures filled with planar proppants is larger under the same sphericity and closure pressure, and the difference between them decreases with the increase in the rock's Young's modulus. The reason is that under the same circumscribed sphere radius, the volume and minimum radius of the planar proppant are smaller than those of the cylindrical proppant. Under the same fracture width, there is more planar proppant filled in the fracture. During deformation of the fracture surface, if the contact area of the planar proppant is larger, it can effectively share the closure pressure to resist the fracture deformation.

| Fracture permeability
While the fracture surface is deformed under the closure pressure, the flowing of the fluid in the fracture pore was analyzed by the Lattice-Boltzmann method, and the influences of the rock's Young's modulus, fracture width, proppant particle size, proppant sphericity, and proppant shape on the permeability and conductivity of the fractures were studied. The flow velocity and pressure of the fluid are presented in the horizontal cross-section of the fracture, as shown in Figure 7 and 8. This fracture is filled with 30/50 mesh proppants with sphericity of 0.909. It can be seen that the flow velocity of the fluid in the pores area is ununiformly distributed, and the flow velocity in most areas is small, less than 1 mm/s, while flow velocity in few parts is pretty large, close to 7 mm/s. Fluid pressure also exhibits uneven reduction, reducing from 1000 kPa at the inlet to 991 kPa at the outlet.
Based on the analysis on the fluid pressure and flow velocity in the fractures, the changes of the fracture permeability under different conditions can be obtained, as shown in Figure 9. Figure 9A shows the influence of the rock's Young's modulus on fracture permeability. With the uniform reduction of the rock's Young's modulus, the reduction in fracture permeability continuously increases under the same closure pressure, and the nonlinear reduction in the fracture permeability along with the raising of the closure pressure is more obvious. When the rock's Young's modulus is 10 GPa, the fracture permeability is reduced by 70% under the closure pressure of 70 MPa, and it tends to be stable. The changes of the fracture permeability under different fracture widths are shown in Figure 9B. When the closure pressure is small (less than 30 MPa in the simulation), the fracture permeability of fractures with the width of 1 mm is significantly higher than that in the fractures with the width of 2 and 3 mm, but with the reduction in fracture width, its stress sensitivity of the permeability becomes stronger. When the closure pressure is greater than 40 MPa, the permeability of the fracture with the width of 1 mm begins to be lower than those of the other two types of fractures, and the difference increases with the increase in the closure pressure. This is because the wider the fractures are, the less the proppants are effected by the fracture surface during settlement, the proppant can be arranged more closely, and its stability is better during compression. The influence of the proppant particle size on the permeability of the fractures is shown in Figure 9C. It can be seen that the larger the proppant particle size is, the greater the fracture conductivity is, resulting from the larger pore between the larger proppants, which makes the fluids flow more easily. However, the permeability of the large-size proppant decreases more rapidly after the closure pressure is loaded. Under the higher closure pressure, the difference in the permeability of the fractures filled with different proppants is gradually reduced. In contrast to the fracture width changes in the fractures ( Figure 6C), the fracture permeability increases with the increase in proppant sphericity as shown in Figure 9D. However, the increased amplitude of permeability decreases with the increase in sphericity. Within the pressure in the simulation, the permeability of the fractures filled with proppants of different sphericity basically changes linearly with the closure pressure, and there is no significant F I G U R E 7 Fluid velocity in the fracture F I G U R E 8 Fluid pressure in the fracture difference in stress sensitivity. For the proppants of two different geometries designed in this simulation, the permeability of the fractures filled with cylindrical proppants is higher than that in the fractures filled with planar proppant in the whole closure pressure range, but its stability under pressure is poor, so that the permeability of both types of fractures decreases as the closing pressure increases.

| Fracture conductivity
The change of fracture conductivity under the closure pressure is the result of the variation of fracture width and fracture permeability. And the influencing factors include rock's Young's modulus, fracture width, proppant particle size, proppant sphericity, and proppant shape. The simulation results are as shown in Figure 10. Under the closure pressure, the width of fractures varies substantially linearly, so the change of the fracture conductivity is basically consistent with the change of the fracture permeability. The fracture conductivity decreases with the increase in the closure pressure, and the nonlinear variation becomes more obvious with the decrease in the Young's modulus of the rocks and the increase in the particle size of the proppant. Under the same closure pressure, the fracture conductivity decreases with the decrease in the rock's Young's modulus, fracture width, proppant particle size, and proppant sphericity. It is necessary to pay special attention to the fact that within the closure pressure range of the simulation, the conductivity of the wider fractures is larger than that of the narrower fractures. By comparing with the change of the fracture permeability ( Figure 9B), it can be seen that increasing fracture width can improve fracture conductivity more significantly than increasing fracture permeability. In addition, the conductivity of the fractures filled with cylindrical proppants is higher than that of the fractures filled with planar proppants.

| CONCLUSIONS
In this study, we investigated the effect of proppant sphericity on the fracture conductivity by combining a plurality of equal-diameter spheres to form proppants of different sphericity and using the discrete element method to construct the disordered proppant packs. Solid mechanics and Lattice-Boltzmann coupling method has the high accuracy in numerical simulation of fracture conductivity. We can draw the following conclusions: 1. Under the same closure pressure, the fracture width decreases with the decrease in the rock's Young's modulus and decreases with the increase in proppant particle size and sphericity. The main reason lies in the amount of proppants filling in fractures and their closeness. 2. Under the same closure pressure, the fracture permeability and fracture conductivity decrease with the decrease in the rock's Young's modulus, proppant particle size, and proppant sphericity, and their stress sensitivity increases with the decrease in the rock's Young's modulus and the increase in the particle size of proppants. 3. For the proppant with higher fracture permeability under a low closure pressure, its conductivity is smaller. It can be seen that increasing fracture width can improve fracture conductivity more significantly than increasing fracture permeability. 4. For the circumscribed sphere of the same radius, the volume and minimum radius of the cylindrical proppants are larger than those of the planar proppants.
Although the stress sensibility of the permeability and conductivity of the fractures filled cylindrical proppanst is stronger, they are larger than the permeability and conductivity of the fractures filled with the planar proppants.

ACKNOWLEDGMENTS
This work was supported by the National Science and Technology Major Project of China (grant numbers: 2016ZX05046-004).