Introduction

As an important technology for the development of low-permeability tight oil/gas reservoirs, the core goal of hydraulic fracturing is to form a high-permeability zone between formation and wellbore and maintain high conductivity for a long time to achieve favorable productivity (Xinhua 2021). The accumulation form of sand in fractures will directly affect the effectiveness of fracturing. Furthermore, the effective accumulation of proppant in fracturing fracture will directly affect the productivity of oil and gas wells (Sun et al. 2021a), which is closely correlated with proppant inherent properties. As a result, the optimization of proppant has become an important problem in the actual field practices (Liu et al. 2019).

For a long time, scholars, all over the world, have carried out a large number of numerical simulation studies regarding proppant transportation and settlement. As for the numerical simulation method (Zeng et al. 2021), direct numerical simulation method, CFD-DEM method (Tsai et al. 2012), and Euler two-fluid method (Deng 2019) are mainly used. Huang et al. (2012) established a full-scale fracture model, in which fluid–solid interactions are incorporated to describe fracturing fluid flow. Eskin (2009) established a non-Newtonian fluid migration model based on the particle dynamics theory. Tsai et al. (2012) established a three-dimensional fracture model and tracked the proppant particles by utilizing CFD-DEM method. Wang et al. (2003) explored settlement law of proppant particles in Newtonian and non-Newtonian fluids, considering that the proppant has intergranular interference during the settlement process. Tong and Mohanty (2016) captured the key characteristics of sand dike formation and migration in fractures, based on the dense discrete phase model.

To sum, the above research mainly involves the numerical simulation and proppant migration theory, focusing on proppant accumulation principle in various complex fractures, while influence of proppant physical properties is rarely to be considered. Yang et al. (2016) and Zhang et al. (2014) investigated the influence of proppant density and size on sand embankment accumulation form, but did not quantify the relationship between sand accumulation form and proppant physical properties at each pumping time node. In view of this, on the basis of the Euler two-fluid method, this paper analyzes influence of proppant physical properties, including concentration, density, and particle size, on the sand accumulation law in fractures. This article is able to achieve the purpose of optimizing proppant and provide reference for on-site fracturing design.

Particle–fluid coupling flow model in fractures

Based on the Euler two-fluid model, proppant particles are regarded as the pseudo-fluid to describe the flow behavior of proppant in fractures (Sun et al. 2021b). This method is frequently utilized to solve multiphase flow problems in various industries.

Governing equation

The continuity equations of liquid phase and solid phase are as follows (Srivastava and Sundaresan 2003):

$$\frac{\partial }{\partial t}\left( {\alpha_{l} \rho_{l} } \right) + \nabla \cdot \left( {\alpha_{l} \rho_{l} \overrightarrow {{v_{l} }} } \right) = 0$$
(1)
$$\frac{\partial }{\partial t}\left( {\alpha_{{\text{s}}} \rho_{{\text{s}}} } \right) + \nabla \cdot \left( {\alpha_{{\text{s}}} \rho_{{\text{s}}} \overrightarrow {{v_{{\text{s}}} }} } \right){ = }0$$
(2)

where \(\alpha\) is the volume fraction; \(\rho\) density, kg/m3; \(t\) time, s; \(\overrightarrow {{v_{{}} }}\) velocity vector, m/s.

The corresponding momentum conservation equations are

$$\frac{\partial }{\partial t}(\alpha_{l} \rho_{l} \overrightarrow {{\nu_{l} }} ) + \nabla \cdot (\alpha_{l} \rho_{l} \overrightarrow {{\nu_{l} }} \overrightarrow {{\nu_{l} }} ) = - \alpha_{l} \nabla p_{l} + \nabla \cdot \tau_{l} + \alpha_{l} \rho_{l} \overrightarrow {g} + \beta (\overrightarrow {{\nu_{s} }} - \overrightarrow {{\nu_{l} }} )$$
(3)
$$\frac{\partial }{\partial t}(\alpha_{s} \rho_{s} \overrightarrow {{\nu_{s} }} ) + \nabla \cdot (\alpha_{s} \rho_{s} \overrightarrow {{\nu_{s} }} \overrightarrow {{\nu_{s} }} ) = - \alpha_{s} \nabla p_{s} + \nabla \cdot \tau_{s} + \alpha_{s} \rho_{s} \overrightarrow {g} + \beta (\overrightarrow {{\nu_{l} }} - \overrightarrow {{\nu_{s} }} )$$
(4)

where \(g\) is the gravity vector, m/s2; \(p\) pressure component, Pa; \(\tau\) shear stress tensor, Pa; \(\beta\) exchange coefficient, kg/(m3·s).

The model, utilized to describe turbulent flow, including turbulent kinetic energy equation as well as turbulent diffusion equation, is provided as follows:

$$\frac{\partial }{\partial t}\left( {\alpha_{l} \rho_{l} k} \right) + \nabla \cdot \left( {\alpha_{l} \rho_{l} \overrightarrow {{v_{l} }} k} \right) = \nabla \cdot \left( {\alpha_{l} \frac{{\mu_{t} }}{{\sigma_{k} }}\nabla k} \right) + G_{k,l} + S_{k} - \alpha_{l} \rho_{l} \varepsilon$$
(5)
$$\frac{\partial }{\partial t}\left( {\alpha_{l} \rho_{l} \varepsilon } \right) + \nabla \cdot \left( {\alpha_{l} \rho_{l} \overrightarrow {{v_{l} }} \varepsilon } \right) = \nabla \cdot \left( {\alpha_{l} \frac{{\mu_{t} }}{{\sigma_{\varepsilon } }}\nabla \varepsilon } \right) + \alpha_{l} \frac{\varepsilon }{k}\left( {C_{1\varepsilon } G_{k,l} - C_{2\varepsilon } \rho_{l} \varepsilon } \right) + S_{\varepsilon }$$
(6)

where \(k\) is the turbulent kinetic energy of continuous phase, m2/s2; \(\varepsilon\) dissipation rate of turbulent kinetic energy, m2/s2; \(\mu_{{\text{t}}}\) turbulent viscosity coefficient, Pa·s; \(G_{{\text{k,l}}}\) generation term of turbulent kinetic energy, kg/(m·s3); \(\sigma_{k}\), \(\sigma_{\varepsilon }\) turbulent Prandtl number, which is 1 and 3, respectively; \(S_{k}\), \(S_{\varepsilon }\) two-phase turbulent exchange term, kg/(m·s3); \(C_{1\varepsilon }\), \(C_{2\varepsilon }\) empirical constants, considered as 1.44 and 1.92, respectively.

Constitutive equation

The simulation is carried out based on the above control equations, and the resistance correlation of Gidaspow model is adopted.

$$\beta = \left\{ {\begin{array}{*{20}l} {\frac{3}{4}C_{D} \frac{{\rho_{l} \alpha_{l} \alpha_{s} \left| {\overrightarrow {{\nu_{l} }} - \overrightarrow {{\nu_{s} }} } \right|}}{{d_{s} }}\alpha_{1}^{ - 2.65} {\kern 1pt} } \hfill & {\alpha_{l} \ge 0.8} \hfill \\ {\frac{{150\alpha_{s} \left( {1 - \alpha_{l} } \right)\mu_{1} }}{{\alpha_{1} d_{s}^{2} }} + \frac{{1.75\rho_{1} \alpha_{s} \left| {\overrightarrow {{\nu_{l} }} - \overrightarrow {{\nu_{s} }} } \right|}}{{d_{s} }}} \hfill & {\alpha_{l} < 0.8} \hfill \\ \end{array} } \right.$$
(7)

Here

$$C_{D} = \left\{ {\begin{array}{*{20}l} {\frac{24}{{{\text{Re}}_{s} }}\left[ {1 + 0.15\left( {{\text{Re}}_{s} } \right)^{0.687} } \right]} \hfill & {{\text{Re}}_{s} < 1000} \hfill \\ {0.44} \hfill & {{\text{Re}}_{s} \ge 1000} \hfill \\ \end{array} } \right.$$
(8)

where \({\text{C}}_{D}\) is the drag coefficient of solid particles; \({\text{Re}}_{{\text{s}}}\) Reynolds number of two-phase slip velocity.

Boundary condition

Liquid–solid two-phase velocity and volume concentration are given as the inlet boundary condition, and the outlet boundary condition is set as the standard atmospheric pressure.

Also, non-slip boundary condition (Sun et al. 2017) is used to describe the dynamic solid phase on the wall.

$$\left( {\mu_{s} + \mu_{f} } \right)\frac{{\partial \overrightarrow {{\nu_{sl} }} }}{\partial x} = - \frac{{\phi \pi \rho_{s} \alpha_{s} g_{0} \sqrt {{\kern 1pt} \Theta_{s} } }}{{2\sqrt 3 \alpha_{s}^{\max } }}\overrightarrow {{\nu_{sl} }}$$
(9)
$$\kappa_{s} \frac{{\partial \Theta_{s} }}{\partial x} = \frac{{\phi \pi \overrightarrow {{\nu_{sl} }} \rho_{s} \alpha_{s} g_{0} \sqrt {{\kern 1pt} \Theta_{s} } }}{{2\sqrt 3 \alpha_{s}^{\max } }} - \frac{{\sqrt 3 \pi \rho_{s} \alpha_{s} g_{0} \left( {1 - e_{w}^{2} } \right)\sqrt {{\kern 1pt} \Theta_{s} } }}{{4\alpha_{s}^{\max } }}\Theta_{s}$$
(10)

where \(\overrightarrow {{v_{sl} }}\) is the slip velocity vector between solid wall, m/s; \(\phi\) reflection coefficient, considered as 0.001; \(\Theta_{s}\) particle pseudo-temperature; \({\text{e}}_{{\text{w}}}\) recovery coefficient of particle–wall collision, considered as 0.9.

Simulation model and solution method

The numerical simulation is carried out with the help of a flat fracture model. Unlike molecular simulation (Sun et al. 2020, 2017), requiring some complex molecular system construction procedures, the numerical simulation can yield prediction results rapidly. The simulated fracture size is the approximate proportional scaling result (2 m), in terms of the average fracture length. Figure 1 shows the width and height of hydraulic fractures of tight oil well (GD1701h) in Cangdong Sag (2 m × 0.2 m × 0.01 m) (Pu et al. 2019). The middle position, on the left side of the fracture, is set as the proppant pumping inlet to simulate the perforation position, with the height of 0.02 m. And, five outlets, with the same size, are set as the proppant inlets, evenly arranged in the right section of the crack. In the numerical simulation, the fracturing fluid is set as clean water, and the geometry of proppant particles is assumed to be spherical.

Fig. 1
figure 1

Schematic diagram and simulation physical model of pressure crack

Full size image

The governing equations are discretized by the finite volume method and solved by the coupled simple algorithm. The convection term is approximated by an upwind difference scheme. In order to avoid divergence, underrelaxation technology is adopted, and the convergence standard is that each residual error is less than 10−4.

Results and discussion

Based on the above model, this article focuses on the influence of proppant physical properties on the stacking law of sand embankment. Related key influential factors include proppant concentration (10%, 15%, 20%, 25%), proppant density (2000 kg/m3, 2250 kg/m3, 2500 kg/m3, 2750 kg/m3), and proppant particle size (0.6 mm, 0.7 mm, 0.8 mm, 0.9 mm). Other parameters are shown in Table 1.

Table1 Calculation parameters of sand bank accumulation form
Full size table

Influence of proppant concentration on sand embankment accumulation law

Three stages, in terms of sand dike formation, are experimentally reported: (1) Under the influence of particle gravity, a small sand dike is formed by settlement or accumulation near the inlet; (2) with the further pumping of fracturing fluid, the sand dike grows vertically; and (3) when it grows vertically to a certain height, the sand dike grows horizontally under the influence of particle migration.

In this article, the accumulation principle of sand dike at each stage is recorded by numerical simulation. Based on the numerical simulation results, the change law, regarding sand embankment accumulation form with time under different proppant concentrations, is presented in Fig. 2. As illustrated, the migration, settlement, and accumulation processes of the sand embankment are similar. The greater the proppant concentration, the less time it takes to reach the accumulation equilibrium state. At the initial stage of accumulation, due to the difference between particle density and water density, particles tend to deposit to the bottom of the crack at the inlet. After the accumulation of sand dike at the inlet, the sand dike is formed again, between 0.25 and 0.5 m away from the inlet, due to the pumping speed. With the further pumping, the valley bottom, formed between the two sand embankments, gradually rises, reflecting the vertical growth stage of the sand embankment. Meanwhile, the higher proppant concentration corresponds to the higher valley bottom height. Thus, although the increasing proppant concentration is beneficial to rapid accumulation to achieve equilibrium, it weakens the fluidization capacity of the proppant to a certain extent. As a result, the greater pumping speed is required to allow the proppants to flow to the deep part of fractures. When a sand dike grows vertically to a certain height, due to the reduction in the gap in fractures, the proppant flow rate increases, the sand carrying capacity improves, and the sand dike moves to the deep part of fractures and gradually reaches the equilibrium state.

Fig. 2
figure 2

Sand embankment accumulation form under different proppant concentrations

Full size image

The space ratio of proppant in fractures is directly related to the conductivity, particularly when the fractures are closed. The change of proppant space ratio with time and the final equilibrium height under different proppant concentrations are shown in Fig. 3. It can be observed that the high proppant concentration helps to accelerate migration of proppant in fractures. At the same time, the settlement of proppant becomes more evident, and the comprehensive equilibrium height of sand dike is higher. Also, it should be noted that the viscosity in this simulation assumes unchanged, which is actually not completely consistent with realistic situations. The increase in fluid concentration is often accompanied with the increase in viscosity, which has an impact on the migration of proppant (Liang et al. 2018).

Fig. 3
figure 3

The change of proppant space ratio under different proppant concentrations

Full size image

Influence of proppant density on sand embankment accumulation

The variation law, in terms of sand embankment accumulation form with time under different proppant densities, is shown in Fig. 4. It can be seen that the sand dike behavior in Fig. 4B–D has the same characteristics, compared with that in Sect. 3.1. However, when the proppant density is small, the density difference between particles and fluid decreases, the particle fluidization capacity strengthens, and proppants are easier to be brought to the deep fractures. Therefore, a small sand dike at the inlet is not formed in Fig. 4a, but the accumulation height gradually rises from the inlet to the deep fractures. Taking the pumping time of 5 s as the observation point, it can be seen that with the increase in proppant density, the front edge of proppant is located at 1.54 m, 1.41 m, 1.31 m, and 1.26 m, respectively. It further shows that the lower the proppant density, the stronger the sand carrying capacity, and the easier the proppant transporting to the deep fracture. Taking the time of pumping of 10 s as the observation point, the position of the highest proppant accumulation in fractures is 0.137 m, 0.150 m, 0.168 m, and 0.173 m, respectively. It can be seen that the higher the proppant density, the faster the settlement speed, and the higher the stacking height. As a result, pumping high-density proppants can produce a better stacking form. High-density proppants are conducive to ensure the effectiveness of early fracturing, by forming a good accumulation form, while the low-density proppants are easier to fluidize and migrate to the deep part of fractures, improving conductivity when fractures are closed. Therefore, density selection should be seriously considered based on actual working conditions.

Fig. 4
figure 4

Change of sand embankment accumulation form under proppant density impact

Full size image

The change of proppant space ratio with time and the final equilibrium height under different proppant densities are shown in Fig. 5. It can be seen from the figure that the high density of proppant results in settlement of the larger proppant and the higher comprehensive equilibrium height of sand embankment at the same time, both of which are conducive to the formation of a fracture zone with high permeability. However, within the simulation range, between 2000 kg/m3 and 2250 kg/m3, comprehensive equilibrium height changes rapidly. Although the comprehensive equilibrium height rises with the further increase in proppant density, the key lifting degree decreases gradually, and the fluidization capacity of proppant cannot be properly guaranteed.

Fig. 5
figure 5

Change of proppant space ratio under proppant density impact

Full size image

Influence of proppant particle size on sand embankment accumulation law

The variation law of sand embankment accumulation with time under different proppant particle sizes is shown in Fig. 4. It can be seen from the figure that the sand embankment accumulation process is similar under different proppant particle sizes. Taking the pumping time of 10 s as an observation point, the proppant with smaller particle size as shown in Fig. 4a, b can migrate to the outlet end of the fracture, while the proppant with larger particle size as shown in Fig. 4c, d stagnates before reaching the outlet end, indicating that the smaller the particle size of the proppant, the stronger the fluidization capacity, and the easier the proppant migrating to the deep part of fractures. Taking the pumping time of 20 s as an observation point, it can be demonstrated that the larger the particle size of proppant, the longer the top of the sand embankment formed by accumulation extends. Under the same proppant concentration, the smaller the particle size, the higher the number of particles pumped per unit time, and the collision interactions between particles affect the settlement effect. As a result, the larger the particle size, the faster the settlement speed, and the easier the sand embankment accumulation (Fig. 6).

Fig. 6
figure 6

Change of sand embankment accumulation form under different proppant particle sizes

Full size image

The change of proppant space ratio with time and the final equilibrium height under different proppant densities are presented in Fig. 7. It can be seen from the figure that particles, with large particle size, can accumulate to achieve a higher comprehensive equilibrium height. Also, at the same time, the space proportion of large particle-size proppant in fractures is larger, maintaining a certain high conductivity after the fracture closure.

Fig. 7
figure 7

Change of proppant space ratio with time under different proppant particle sizes

Full size image

Conclusions

  1. 1.

    The proppant accumulation process in fractures can be roughly divided into four stages: a) Two sand banks are formed at the inlet of fractures and a position, with a certain distance away from the inlet; b) the second sand dike grows vertically; c) after the sand dike reaches a certain height, it grows horizontally; d) when the balance state is reached, comprehensive balance height of the sand embankment remains unchanged.

  2. 2.

    Increasing the proppant concentration is conducive to improving sedimentation rate, as well as accumulation of proppant particles. However, with the increase in sand concentration, fluid viscosity will generally increase, which will mitigate the migration of proppant in fractures.

  3. 3.

    Reducing proppant density and particle size will enhance fluidization capacity of particles, which is conducive to increasing the migration depth of particles in fractures at the initial pumping stage. Also, the higher the density and particle size, the higher the comprehensive equilibrium height, which is helpful to improve the conductivity after the closure of high-pressure fracture.