1. Introduction
Currently, the mainstream idea in research on the physical properties of proppants is physical experiments or a combination of physical experiments and numerical simulation. The simulation of proppant particles is generally implemented using auxiliary software such as DEM (Altair EDEM v2.7) and CFD (Autodesk CFD 2021). However, there is almost no one focused on the computer simulation of hydraulic fracturing proppant particles.
Fan et al. [
1] used DEM software to simulate the compaction state of proppant particles under pressure. The pore structure of the proppant filling layer simulated by DEM was extracted and inputted into the LB (Lattice Boltzmann) simulator as an internal boundary condition for fluid flow. Therefore, the transport and precipitation mechanisms of proppant particles within the numerical framework of DEM-LB were studied. A continuous medium-fine particle deposition model was used to fit the numerical results of pore-scale particle tracking deposition, thereby determining the macro deposition coefficient. The effects of fine particle size, proppant particle size heterogeneity, and effective stress on the macro deposition coefficient were studied. Similarly, scholars have conducted research on the properties of particles such as heat conduction, force conduction, and particle breakage under effective stress and microstructure through DEM [
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12].
As early as 1611, Johannes Kepler began to study the arrangement of spheres and proposed the densest accumulation of spheres, known as the Kepler conjecture [
13]. This conjecture is about the best sphere-loading method in three-dimensional Euclidean space. That is, the sphere-loading method with the smallest gap left. This conjecture holds that under the condition that each sphere has the same size, the density of the sphere without any loading method is greater than face-centered cubic and hexagonal densest packing method with a density slightly greater than 74%. The microstructure of the proppants studied in this paper is similar to the Kepler conjecture, but there are also significant differences.
In order to simulate the placement of particles, it is assumed that the particles are circular and are not intersecting, overlapping, or suspended. There are three main paving methods shown in
Figure 1.
Figure 1a shows the entire two-dimensional space that starts from the center circle and gradually diffuses outward by placing new circles tangent to it [
14]. This method can avoid the intersection, overlap, and suspension of particles. However, the effect simulated by this method is not suitable for the given two-dimensional rectangular space.
Figure 1b shows placing a circle inside the boundary of the two-dimensional space and gradually advancing towards the middle until the entire space is filled [
15,
16,
17]. This method can obtain a two-dimensional rectangular space with good filling effect. However, there will be large pores in the middle of the space, making the particles in this part appear linear in suspension. The stacking method shown in
Figure 1c is to place particles at the bottom of the two-dimensional space, filling the bottom space from left to right. In addition to the bottom layer, particles are gradually stacked in the same manner from bottom to top until they fill the entire two-dimensional space [
17,
18,
19,
20]. However, this method does not take into account the principle that the placement of particles is based on the minimum distance from the bottom edge. Although it is possible to ensure that there is no intersection, overlap, or suspension between particles, there will be large residual space at the top. The intuitive result is that the entire two-dimensional space is not fully filled with particles.
Due to the complexity of particle accumulation, most scholars mainly use simulation software for the two-dimensional simulation of particles. There are very few studies on the three-dimensional simulation of proppant particles. To study the properties of proppant particles, this paper will start from the microscopic arrangement of proppant particles. Considering the randomness and certainty of three-dimensional particle stacking and arrangement, the Monte Carlo stochastic method and an optimization model were used to conduct a three-dimensional computer simulation of proppant particles to reveal the randomness and regularity of the microscopic arrangement of proppant particles.
2. Physical Model of Proppant Particles
There are mainly six forms, regarding the three-dimensional arrangement structure of proppant particles, each of which corresponds to a different included angle and porosity, as shown in
Figure 2 and
Figure 3 [
21,
22,
23]. There is a correspondence between
Figure 1 and
Figure 2. The first three arrangements belong to the square system (the bottom four particle structures are square), while the last three arrangements belong to the hexagonal system (the bottom four particle structures are diamond-shaped).
To facilitate three-dimensional computer simulation, a physical model of proppant particle placement as shown in
Figure 4 has been established. It is assumed that the proppant particles are spherical with the same particle size. Under closed pressure, the fracture changes in the direction of force parallel to the
x-axis (yOz plane), while the force on the other six faces of the plane remains unchanged. The stacking rules of proppant particles are as follows:
(1) On the bottom layer, based on the randomness of particle arrangement, the four boundaries of the given three-dimensional space are laid first. Based on the balls on four boundaries, the position of each new proppant particle is determined to be the shortest distance between the center of the proppant particle and the four boundaries, which is a deterministic placement and covers the entire bottom space one by one.
(2) Starting from the bottom layer, the accumulation of proppant particles is completed layer by layer until to the top. The position where each proppant particle is placed is determined by the shortest distance between the proppant particle’s center and the bottom of the fracture.
3. Contact Model of Proppant Particles
As shown in
Figure 5, the three-dimensional coordinates of the underlying particles are established and discretized for three-dimensional processing. The spherical center coordinate of any proppant particle is recorded as (
xi,
yi,
zi). The maximum number of particles laid in various directions in three-dimensional space can be represented by
nx,
ny,
nz, where
In the formula, nx is the maximum number of proppant particles laid in the x-axis direction, which is also the width direction of the fracture, dimensionless. ny is the maximum number of proppant particles laid in the y-axis direction, which is also the height direction of the fracture, dimensionless. nz is the maximum number of proppant particles laid in the z-axis direction, dimensionless. int (x0) represents the downward rounding of any real number x0. wf is the fracture width, mm. Hf is the fracture height, mm. Lf is the fracture length, mm. R is the diameter of the proppant particle, mm.
According to the arrangement rule of proppant particles, the three-dimensional simulation of proppant particles is processed according to random phenomena and deterministic phenomena.
3.1. Random Phenomenon and Treatment of Proppant Particle Arrangement
The random phenomenon of proppant particle arrangement is mainly manifested in the bottom layer of three-dimensional space (
k = 1). As shown in
Figure 5, in
Figure 5a, first determine the balls at the four corner positions on the bottom space. Determine the maximum number of balls between adjacent balls in
Figure 5a, and randomly place them in
Figure 5b using a random weighted distribution coefficient. Use a search algorithm to find the position of a new ball between two ad-jacent balls in
Figure 5b, which is the shortest distance from any edge interface among the four face interfaces. Continue to select the position with the shortest distance from the edge interface and place the new ball until the entire bottom space cannot accommodate the new ball, as shown in
Figure 5d.
(1)
Figure 5 shows the random phenomenon of proppant particles contacting the (0,0)-(0,
wf) boundary. Random phenomena generate random numbers
ξi(
i = 1, 2, …,
nx–1). The remaining space is allocated according to the weighted average method of random numbers; that is, the random weighted distribution coefficient is introduced as
When (0,0)-(0,
wf) boundary contacts, the new spherical center coordinates (
xi,
yi,
zi) of the proppant particles are calculated by the following formula:
(2) Another case is the random phenomenon caused by the contact of proppant particles with the (0, 0)-(0,
Hf) boundary. Random phenomena generate random numbers
ξi(
i = 1, 2, …,
ny–1). The remaining space is allocated according to the weighted average method of random numbers; that is, the random weighted distribution coefficient is introduced as
When (0,0)-(0,
Hf) boundary contacts, the new spherical center coordinates (
xi,
yi,
zi) of the proppant particles are calculated by the following formula:
(3) Another case is the random phenomenon caused by the contact of proppant particles with the (0,
Hf)-(
wf,
Hf) boundary. Random phenomena generate random numbers
ξm (
m = 1, 2, …,
nx–1). The remaining space is allocated according to the weighted average method of random numbers; that is, the random weighted distribution coefficient is introduced as
When (0,
Hf)-(
wf,
Hf) boundary contacts, the new spherical center coordinates (
xijk,
yijk,
zijk) of the proppant particles are calculated by the following formula:
(4) Last, a random phenomenon is caused by the contact of proppant particles with the (
wf, 0)-(
wf,
Hf) boundary. Random phenomena generate random numbers
ξn (
n = 1, 2, …,
ny–1). The remaining space is allocated according to the weighted average method of random numbers; that is, the random weighted distribution coefficient is introduced as
When (
wf, 0)-(
wf,
Hf) boundary contacts, the new spherical center coordinates (
xi,
yi,
zi) of the proppant particles are calculated by the following formula:
The new spherical center coordinate set
Q for the bottom layer (
k = 1) is
In
Figure 5b, a temporary dynamic matrix T0 is set as the set of the innermost spherical centers of the bottom layer (
k = 1) in three-dimensional space. When the total number of proppant particles np is less than or equal to 2
nx + 2
ny − 4, the temporary dynamic matrix
T0 equals
Q.
When the total number of proppant particles
np is greater than 2
nx + 2
ny − 4, the new sphere may be generated by two or more spheres together, and it is not possible to express the new sphere centers using a certain mathematical formula. Therefore, consider using arrays to express the set of new sphere centers. Take the shortest distance between the new sphere centers and the four boundaries (0, 0)-(0,
wf), (0, 0)-(0,
Hf), (0, Hf)-(
wf,
Hf), and (
wf, 0)-(
wf,
Hf) as the target, and recursively arrange the spheres one by one in a continuous manner towards the middle. Then, the spherical center coordinates
Q of all particles at the bottom layer can be expressed as
When the bottom space has been filled with spheres, a new temporary dynamic matrix
T is set as the set of all the sphere centers at the bottom of the three-dimensional space. Based on the sphere centers in set
T, the new sphere centers are recursively arranged one by one along the direction
z =
Lf, with the shortest distance from
z = 0 as the target. When the entire three-dimensional space is fully filled with spheres, set
Q representing the spherical centers of all proppant particles in the three-dimensional space is
where
n is the total number of proppant particles in three-dimensional space, dimensionless.
m is the total number of proppant particles in the underlying space, and
m <
n, dimensionless. In the three-dimensional space,
xi is the coordinate of the spherical center in the
x-axis direction.
yi is the coordinate of the spherical center in the
y-axis direction.
zi is the coordinate of the spherical center in the
z-axis direction.
3.2. Determination Phenomenon and Treatment of Proppant Particle Arrangement
The phenomenon of determining the three-dimensional simulation of proppant particles needs to be described and expressed through three parts: bottom layer, even layer in the middle, and odd layer in the middle.
(1) Proppant particles at the bottom layer (k = 1):
A new sphere center is generated based on two or more existing sphere centers. Therefore, the new sphere center can be represented as
where (
xijk,
yijk,
zijk) is the coordinate of the new spherical center of the proppant particle. (
x(i1)(j1)1,
y(i1)(j1)1,
z(i1)(j1)1), (
x(i2)(j2)1,
y(i2)(j2)1,
z(i2)(j2)1), and (
x(i3)(j3)1,
y(i3)(j3)1,
z(i3)(j3)1) all represent the existing sphere center of the first layer.
According to the above three conditions, search for the location of the new spherical center with the shortest distance from the four boundaries of (0, 0)-(0,
wf), (0, 0)-(0,
Hf), (0,
Hf)-(
wf,
Hf), and (
wf, 0)-(
wf,
Hf) through the traversal search method
(2) Proppant particles in the intermediate layer (1 < k < n):
a. There are two situations when the proppant particles are tangent to the boundary (0, 0)-(
wf, 0) and not tangent to both the boundary (0, 0)-(0,
Hf) and the boundary (
wf, 0)-(
wf,
Hf). A proppant particle is located on top of one sphere, or on top of two spheres. The new spherical center coordinates of the proppant particles that may be placed are calculated by the following formula:
b. When the proppant particles are tangent to the boundary (0, 0)-(
wf, 0) and tangent to the boundary (0, 0)-(0,
Hf), the new spherical center coordinates of the proppant particles that may be placed are
c. When the proppant particles are tangent to the boundary (0, 0)-(
wf, 0) and the boundary (
wf, 0)-(
wf,
Hf), the new spherical center coordinates of the proppant particles that may be placed are
d. When the proppant particles are not tangent to all four boundaries, the new spherical center coordinates of the proppant particles that may be placed are
e. When the proppant particles are not tangent to the boundary (0, 0)-(
wf, 0) or the boundary (0,
Hf)-(
wf,
Hf), and are tangent to the boundary (0, 0)-(0,
Hf), the new spherical center coordinates of the proppant particles that may be placed are
When the proppant particles are not tangent to the boundary (0, 0)-(
wf, 0) or the boundary (0,
Hf)-(
wf,
Hf), and are tangent to the boundary (
wf, 0)-(wf,
Hf), the new spherical center coordinates of the proppant particles that may be placed are
g. When the proppant particles are tangent to the boundary (0,
Hf)-(
wf,
Hf), and are not tangent to either the boundary (0, 0)-(0,
Hf) or the boundary (
wf, 0)-(
wf,
Hf), the new spherical center coordinates of the proppant particles that may be placed are
h. When the proppant particles are tangent to the boundary (0,
Hf)-(
wf,
Hf) and the boundary (0, 0)-(0,
Hf), the new spherical center coordinates of the proppant particles that may be placed are
i. When the proppant particles are tangent to the boundary (0,
Hf)-(
wf,
Hf) and the boundary (
wf, 0)-(
wf,
Hf), the new spherical center coordinates of the proppant particles that may be placed are
where (
xijk,
yijk,
zijk) is the spherical center coordinate of the new proppant particle; (
x(i1)(j1)(k1),
y(i1)(j1)(k1),
z(i1)(j1)(k1)), (
x(i2)(j2)(k2),
y(i2)(j2)(k2),
z(i2)(j2)(k2)), and (
x(i3)(j3)(k3),
y(i3)(j3)(k3),
z(i3)(j3)(k3)) all represent the spherical center coordinates of the existing proppant particles.
The objective function is to find the lowest position of the new sphere center, where the objective function is
3.3. D simulation Steps of Proppant
The 3D simulation of proppant particles includes the simulation of bottom-layer proppant particles and the simulation of other-layer proppant particles, and eight criteria for proppant particle simulation. The specific steps are as follows:
3.4. Bottom Proppant Particles
The 3D simulation process of the bottom proppant particles is shown in
Figure 6.
a. The maximum boundary values corresponding to the x, y, and z axes of a three-dimensional space are given as wf, Hf, and Lf, giving the specified sphere diameter R.
b. Four spheres ((R/2, R/2, R/2), (wf–R/2, R/2, R/2), (R/2, Hf–R/2, R/2), (wf–R/2, Hf–R/2, R/2)) are placed in the four corners of the bottom layer. The Monte Carlo method is used to generate random numbers near four boundaries((0, 0)-(0, wf), (0, 0)-(0, Hf), (0, Hf)-(wf, Hf), (wf, 0)-(wf, Hf)) and determine the random interval between the spheres and the position of the new sphere.
c. The bottom layer is based on the position of the existing sphere and uses a traversal search algorithm to find the shortest position from the four boundaries, namely
d. Determine whether the new sphere intersects, overlaps, or floats with other filled spheres. If so, go to step c and readjust the sphere center or search for a new sphere center. If this phenomenon does not exist, determine the new sphere center (xi, yi, zi) and proceed to the next step.
e. Determine whether the remaining space satisfies the volume occupied by a sphere. If new spheres can still be placed in the bottom space, proceed to step c. Otherwise, it can be determined that the underlying space has been filled and the operation is over.
The flow chart of the three-dimensional simulation for the single-particle-size proppant at the bottom layer is shown in
Figure 6.
3.5. Other Proppant Particles
The 3D simulation flowchart of the other layers of proppant particles is shown in
Figure 7.
a. Assign the values of all the bottom spherical center matrices Bf to the top spherical center matrix Tf.
b. Search for the minimum value min(zi) in the z-axis direction in the top-level spherical center matrix Tf, and determine the minimum spherical center (xi, yi, min(zi)) in the z-axis direction.
c. Search for other spherical centers near the spherical center (xi, yi, min(zi)). When the distance between the spherical center (xi, yi, min(zi)) and the other spherical center subtracted by 2 × R is less than (3/2) R, it is the spherical center near the spherical center (xi, yi, min(zi)).
d. The set of possible new spherical centers can be determined by the smallest spherical center (xi, yi, min(zi)) in the z-axis direction and the nearby spherical center or edge interface.
e. Determine whether the possible new sphere intersects, overlaps, or floats with any other sphere. If the above phenomenon exists, readjust the possible new spherical center and proceed to step d. If the above phenomenon does not exist, proceed to the next step.
f. Delete the possible new spherical center that intersects, overlaps, or floats with any sphere in the possible new spherical center collection.
g. After determining whether the set of possible new spherical centers is empty, if the set of possible new spherical centers is empty, proceed to step b to re-find the smallest ball center (xi, yi, min(zi)) in the fixed z-axis direction. If the collection of “possible new spherical centers “ is not empty, proceed to the next step.
h. Determine the possible new spherical center collection again.
i. Determine whether the possible new spherical center exceeds [0, wf] in the x-axis direction, or exceeds [0, Hf] in the y-axis direction. If the boundary of the x-axis or y-axis is exceeded, readjust the possible new spherical centers through the eight judgment conditions and skip to step d. If it does not exceed the boundaries of the x-axis or y-axis, proceed to the next step.
j. Determine the set of possible new spherical centers again, and find the smallest new spherical center in the z-axis direction in the set, which is the most suitable new spherical center.
k. Determine whether the most suitable new spherical center zi in the z-axis direction is greater than the upper limit Hf in the z-axis direction. If it is not greater than Hf, skip to step b to continue running. If the best fit for the new spherical center zi is greater than the upper limit Hf in the z-axis direction, the operation stops and ends.
3.6. Eight Cases where Proppant Particles Exceed the Boundary
There are eight cases where the proppant particles exceed the boundary, as shown in
Figure 8 and
Figure 9.
a. If xi < R/2 and R/2 < yi < Hf–R/2 are in a possible new spherical center (xi, yi, zi), then the possible new spherical center coordinates are adjusted to (R/2, yi, zi).
b. If xi < R/2 and yi < R/2 are in a possible new spherical center (xi, yi, zi), then the possible new spherical center coordinates are adjusted to (R/2, yi, zi).
c. If xi < R/2 and yi > Hf–R/2 are in a possible new spherical center (xi, yi, zi), then the possible new spherical center coordinates are adjusted to (R/2, Hf–R/2, zi).
d. If in a possible new spherical center (xi, yi, zi), xi > wf–R/2 and R/2 < yi < Hf–R/2, then the possible new spherical center coordinates are adjusted to (wf–R/2, yi, zi).
e. If xi > wf–R/2 and yi < R/2 are in a possible new spherical center (xi, yi, zi), then the possible new spherical center coordinates are adjusted to (wf–R/2, R/2, zi).
f. If in a possible new spherical center (xi, yi, zi), xi>wf–R/2 and yi > Hf–/2R, then the possible new spherical center coordinates are adjusted to (wf–R/2, Hf–R/2, zi).
g. If in a possible new spherical center (xi, yi, zi), yi < R/2 and R/2 < xi < wf–R/2, then the possible new spherical center coordinates are adjusted to (xi, R/2, zi).
h. If in a possible new spherical center (xi, yi, zi), yi > Hf–R/2 and R/2 < xi < wf -/2R, then the possible new spherical center coordinates are adjusted to (xi, Hf–R/2, zi).
According to the above steps, a three-dimensional computer simulation of proppant particles was conducted, and the simulation results are shown in
Figure 10.