Next Article in Journal
On the Lock-In Phenomena near the Transonic Buffet Onset of a Prescribed Pitching Airfoil
Previous Article in Journal
The Yeast-Based Probiotic Encapsulation Scenario: A Systematic Review and Meta-Analysis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Computer 3D Simulation of Proppant Particles

1
School of Sciences, Southwest Petroleum University, Chengdu 610500, China
2
College of Management Science, Sichuan University of Science and Engineering, Zigong 643002, China
3
Department of Mathematical Sciences, Tsinghua University, Beijing 100084, China
4
School of Sciences, Civil Aviation Flight University of China, Deyang 618307, China
*
Author to whom correspondence should be addressed.
Appl. Sci. 2024, 14(13), 5462; https://doi.org/10.3390/app14135462
Submission received: 11 April 2024 / Revised: 14 June 2024 / Accepted: 17 June 2024 / Published: 24 June 2024

Abstract

:
Proppants are one of the key materials for hydraulic fracturing, whose main role is to support fractures and create a channel through which oil and gas can flow. The nature of proppants is the most talked about feature besides their cost, for example, their sphericity, turbidity, particle size, or strength. The porosity, permeability, and fracture conductivity of proppants in fractures are also the main indicators to measure the performance of them. These indicators are usually obtained through physical experiments. However, experimental results often differ depending on the experimental scheme. Different stacking methods of proppant particles lead to this phenomenon. The nature of proppant particles in fractures varies with the way they accumulate. This paper will start with 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 three-dimensional computer simulation of proppant particles. This lays an important foundation for revealing the randomness and regularity of the micro arrangement of proppant particles.

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
n x = int ( w f R ) n y = int ( H f R ) n z = int ( L f R )
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
λ i = ξ i 1 k = 1 n x 1 ξ k ( w f R n x ) , ( i = 1 , 2 , , n x )
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:
{ x i 11 = R 2 ( λ i + 2 ( i 1 ) + 1 ) y i 11 = R 2 z i 11 = R 2
(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
λ j = ξ j 1 k = 1 n y 1 ξ k ( H f R n y ) , ( j = 1 , 2 , , n y )
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:
{ x 1 j 1 = R 2 y 1 j 1 = R 2 ( λ j + 2 ( j 1 ) + 1 ) z 1 j 1 = R 2
(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
λ m = ξ m 1 k = 1 n x 1 ξ k ( w f R n x ) , ( i = 1 , 2 , , n x )
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:
{ x i n y 1 = R 2 ( λ m + 2 ( m 1 ) + 1 ) y i n y 1 = H f R 2 z i n y 1 = R 2
(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
λ n = ξ n 1 k = 1 n y 1 ξ k ( H f R n y ) , ( n = 1 , 2 , , n y )
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:
{ x n x j 1 = w f R 2 y n x j 1 = R 2 ( λ n + 2 ( n 1 ) + 1 ) z n x j 1 = R 2
The new spherical center coordinate set Q for the bottom layer (k = 1) is
Q = [ x 1 x 2 x n y 1 y 2 y n z 1 z 2 z n ] .
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 2nx + 2ny − 4, the temporary dynamic matrix T0 equals Q.
When the total number of proppant particles np is greater than 2nx + 2ny − 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
Q = [ x 1 x 2 x m y 1 y 2 y m z 1 z 2 z m ] .
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
Q = [ x 1 x 2 x m x i x n y 1 y 2 y m y i y n z 1 z 2 z m z i z n ]
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 is generated between two adjacent spheres in matrix T0, as shown in subgraph (c) in Figure 5.
  • A new sphere is generated between two non-adjacent spheres in matrix T0, as shown in subgraph (d) in Figure 5.
  • When the distance between the centers of any two spheres in matrix T0 is greater than 2R, no new spheres are generated between the two spheres.
A new sphere center is generated based on two or more existing sphere centers. Therefore, the new sphere center can be represented as
{ ( x i j 1 x ( i 1 ) ( j 1 ) 1 ) 2 + ( y i j 1 y ( i 1 ) ( j 1 ) 1 ) 2 + ( z i j 1 z ( i 1 ) ( j 1 ) 1 ) 2 = R ( x i j 1 x ( i 2 ) ( j 2 ) 1 ) 2 + ( y i j 1 y ( i 2 ) ( j 2 ) 1 ) 2 + ( z i j 1 z ( i 2 ) ( j 2 ) 1 ) 2 = R z i j 1 = z ( i 1 ) ( j 1 ) 1 = z ( i 2 ) ( j 2 ) 1 = R 2 ,
{ ( x i j 1 x ( i 1 ) ( j 1 ) 1 ) 2 + ( y i j 1 y ( i 1 ) ( j 1 ) 1 ) 2 + ( z i j 1 z ( i 1 ) ( j 1 ) 1 ) 2 = R ( x i j 1 x ( i 2 ) ( j 2 ) 1 ) 2 + ( y i j 1 y ( i 2 ) ( j 2 ) 1 ) 2 + ( z i j 1 z ( i 2 ) ( j 2 ) 1 ) 2 = R ( x i j 1 x ( i 3 ) ( j 3 ) 1 ) 2 + ( y i j 1 y ( i 3 ) ( j 3 ) 1 ) 2 + ( z i j 1 z ( i 3 ) ( j 3 ) 1 ) 2 = R z i j 1 = z ( i 1 ) ( j 1 ) 1 = z ( i 2 ) ( j 2 ) 1 = z ( i 3 ) ( j 3 ) 1 = R 2 ,
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
min ( min j x j , min j y j , min j ( w f x j ) , min j ( H f y j ) )
(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:
{ ( x i j k x ( i 1 ) ( j 1 ) ( z 1 ) ) 2 + ( y i j k y ( i 1 ) ( j 1 ) ( z 1 ) ) 2 + ( z i j k z ( i 1 ) ( j 1 ) ( z 1 ) ) 2 = R y i j k = R 2 ,   z i j k = z ( i 1 ) ( j 1 ) ( z 1 ) + R
{ ( x i j k x ( i 1 ) ( j 1 ) ( z 1 ) ) 2 + ( y i j k y ( i 1 ) ( j 1 ) ( z 1 ) ) 2 + ( z i j k z ( i 1 ) ( j 1 ) ( z 1 ) ) 2 = R ( x i j k x ( i 2 ) ( j 2 ) ( z 2 ) ) 2 + ( y i j k y ( i 2 ) ( j 2 ) ( z 2 ) ) 2 + ( z i j k z ( i 2 ) ( j 2 ) ( z 2 ) ) 2 = R y i j k = R 2
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
{ ( x i j k x ( i 1 ) ( j 1 ) ( z 1 ) ) 2 + ( y i j k y ( i 1 ) ( j 1 ) ( z 1 ) ) 2 + ( z i j k z ( i 1 ) ( j 1 ) ( z 1 ) ) 2 = R x i j k = y i j k = R 2
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
{ ( x i j k x ( i 1 ) ( j 1 ) ( z 1 ) ) 2 + ( y i j k y ( i 1 ) ( j 1 ) ( z 1 ) ) 2 + ( z i j k z ( i 1 ) ( j 1 ) ( z 1 ) ) 2 = R x i j k = w f R 2 ,   y i j k = R 2
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
{ ( x i j k x ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( y i j k y ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( z i j k z ( i 1 ) ( j 1 ) ( k 1 ) ) 2 = R ( x i j k x ( i 2 ) ( j 2 ) ( k 2 ) ) 2 + ( y i j k y ( i 2 ) ( j 2 ) ( k 2 ) ) 2 + ( z i j k z ( i 2 ) ( j 2 ) ( k 2 ) ) 2 = R ( x i j k x ( i 3 ) ( j 3 ) ( k 3 ) ) 2 + ( y i j k y ( i 3 ) ( j 3 ) ( k 3 ) ) 2 + ( z i j k z ( i 3 ) ( j 3 ) ( k 3 ) ) 2 = R
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
{ ( x i j k x ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( y i j k y ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( z i j k z ( i 1 ) ( j 1 ) ( k 1 ) ) 2 = R ( x i j k x ( i 2 ) ( j 2 ) ( z 2 ) ) 2 + ( y i j k y ( i 2 ) ( j 2 ) ( z 2 ) ) 2 + ( z i j k z ( i 2 ) ( j 2 ) ( z 2 ) ) 2 = R x i j k = R 2
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
{ ( x i j k x ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( y i j k y ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( z i j k z ( i 1 ) ( j 1 ) ( k 1 ) ) 2 = R ( x i j k x ( i 2 ) ( j 2 ) ( z 2 ) ) 2 + ( y i j k y ( i 2 ) ( j 2 ) ( z 2 ) ) 2 + ( z i j k z ( i 2 ) ( j 2 ) ( z 2 ) ) 2 = R x i j k = w f R 2
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
{ ( x i j k x ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( y i j k y ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( z i j k z ( i 1 ) ( j 1 ) ( k 1 ) ) 2 = R ( x i j k x ( i 2 ) ( j 2 ) ( z 2 ) ) 2 + ( y i j k y ( i 2 ) ( j 2 ) ( z 2 ) ) 2 + ( z i j k z ( i 2 ) ( j 2 ) ( z 2 ) ) 2 = R y i j k = H f R 2
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
{ ( x i j k x ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( y i j k y ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( z i j k z ( i 1 ) ( j 1 ) ( k 1 ) ) 2 = R x i j k = R 2 , y i j k = R 2
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
{ ( x i j k x ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( y i j k y ( i 1 ) ( j 1 ) ( k 1 ) ) 2 + ( z i j k z ( i 1 ) ( j 1 ) ( k 1 ) ) 2 = R x i j k = w f R 2 , y i j k = H f R 2
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
min i z i j 1

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), (wfR/2, R/2, R/2), (R/2, HfR/2, R/2), (wfR/2, HfR/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
min ( min j x j , min j y j , min j ( w f x j ) , min j ( H f y j ) )
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 < HfR/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 > HfR/2 are in a possible new spherical center (xi, yi, zi), then the possible new spherical center coordinates are adjusted to (R/2, HfR/2, zi).
d. If in a possible new spherical center (xi, yi, zi), xi > wfR/2 and R/2 < yi < HfR/2, then the possible new spherical center coordinates are adjusted to (wfR/2, yi, zi).
e. If xi > wfR/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 (wfR/2, R/2, zi).
f. If in a possible new spherical center (xi, yi, zi), xi>wfR/2 and yi > Hf/2R, then the possible new spherical center coordinates are adjusted to (wfR/2, HfR/2, zi).
g. If in a possible new spherical center (xi, yi, zi), yi < R/2 and R/2 < xi < wfR/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 > HfR/2 and R/2 < xi < wf -/2R, then the possible new spherical center coordinates are adjusted to (xi, HfR/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.

4. Conclusions

(1) The arrangement of proppant particles in the bottom layer of three-dimensional space has a certain randomness. In this paper, a Monte Carlo method is used to generate random numbers, and the remaining space at the bottom boundary is distributed in a weighted average manner of random numbers; that is, a random weighted distribution coefficient is introduced to randomly determine the position of particles in three-dimensional space.
(2) The arrangement of proppant particles in three-dimensional space has certain certainty. Based on existing proppant particles, an optimization model is adopted to ensure that each particle is the shortest distance from the target boundary, and proppant particles are placed one by one to achieve the stacking and arrangement of proppant particles.
(3) The stacking arrangement of proppant particles is very complex. During the simulation process, there may be intersection, overlap, and suspension between the new particles and other existing particles, as well as eight types of boundary violations between the new particles and the boundary. Computer simulation of proppant particles should be implemented according to the simulation rules of particles.

Author Contributions

Conceptualization, D.G.; writing—original draft preparation, K.L.; writing—review and editing, Y.Z.; validation, Z.G. All authors have read and agreed to the published version of the manuscript.

Funding

This paper was funded by the National Science and Technology Major Project of China (no. 2016ZX05065 and no. 2016ZX05042-003).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Acknowledgments

The authors would like to thank the editor for handling this submission and the anonymous referees for reading the manuscript.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Fan, M.; Chen, C. Numerical simulation of the migration and deposition of fine particles in a proppant-supported fracture. J. Pet. Sci. Eng. 2020, 194, 107484. [Google Scholar] [CrossRef]
  2. Zhang, H.; Zhou, Q.; Xing, H.; Muhlhaus, H. A DEM study on the effective thermal conductivity of granular assemblies. Powder Technol. 2011, 205, 172–183. [Google Scholar] [CrossRef]
  3. Foroutan, T.; Mirghasemi, A. CFD-DEM model to assess stress-induced anisotropy in undrained granular material. Comput. Geotech. 2020, 119, 103318. [Google Scholar] [CrossRef]
  4. Luo, T.; Ooi, E.; Chan, A.; Fu, S. The combined scaled boundary finite-discrete element method, Grain breakage modelling in cohesion-less granular media. Comput. Geotech. 2017, 88, 199–221. [Google Scholar] [CrossRef]
  5. Tahir, A.; Rasche, S.; Könke, C. Discrete element model development of ZTA ceramic granular powder using micro computed tomography. Adv. Powder Technol. 2018, 29, 3471–3482. [Google Scholar] [CrossRef]
  6. Qu, T.; Feng, Y.T.; Zhao, T.; Wang, M. Calibration of linear contact stiffnesses in discrete element models using a hybrid analytical-computational framework. Powder Technol. 2019, 356, 795–807. [Google Scholar] [CrossRef]
  7. Moghimi, M.; Quinlan, N. Application of background pressure with kinematic criterion for free surface extension to suppress non-physical voids in the finite volume particle method. Eng. Anal. Bound. Elem. 2019, 106, 126–138. [Google Scholar] [CrossRef]
  8. Campello, E. A computational model for the simulation of dry granular materials. Int. J. Non-Linear Mech. 2018, 106, 89–107. [Google Scholar] [CrossRef]
  9. Li, X.; Wang, Z.; Zhang, S.; Duan, Q. Multiscale modeling and characterization of coupled damage-healing-plasticity for granular materials in concurrent computational homogenization approach. Comput. Methods Appl. Mech. Eng. 2018, 342, 354–383. [Google Scholar] [CrossRef]
  10. Ren, H.; Zhuang, X.; Rabczuk, T.; Zhu, H. Dual-support smoothed particle hydrodynamics in solid, variational principle and implicit formulation. Eng. Anal. Bound. Elem. 2019, 108, 15–29. [Google Scholar] [CrossRef]
  11. Fu, P.; Dafalias, Y. Relationship between void-and contact normal-based fabric tensors for 2D idealized granular materials. Int. J. Solids Struct. 2015, 63, 68–81. [Google Scholar] [CrossRef]
  12. Poorsolhjouy, P.; Gonzalez, M. Connecting discrete particle mechanics to continuum granular micromechanics, Anisotropic continuum properties under compaction. Mech. Res. Commun. 2018, 92, 21–27. [Google Scholar] [CrossRef]
  13. Ball, P. In retrospect, On the Six-Cornered Snowflake. Nature 2011, 480, 455. [Google Scholar] [CrossRef]
  14. Wang, W.; Ming, C.; Lo, S. Generation of triangular mesh with specified size by circle packing. Adv. Eng. Softw. 2007, 38, 133–142. [Google Scholar] [CrossRef]
  15. Benabbou, A.; Borouchaki, H.; Laug, P.; Lu, J. Geometrical modeling of granular structures in two and three dimensions. application to nanostructures. Int. J. Numer. Methods Eng. 2009, 80, 425–454. [Google Scholar] [CrossRef]
  16. Benabbou, A.; Borouchaki, H.; Laug, P.; Lu, J. Numerical modeling of nanostructured materials. Finite Elem. Anal. Des. 2010, 46, 165–180. [Google Scholar] [CrossRef]
  17. Bagi, K. An algorithm to generate random dense arrangements for discrete element simulations of granular assemblies. Granul. Matter 2005, 7, 31–43. [Google Scholar] [CrossRef]
  18. Feng, Y.; Han, K.; Owen, D. Filling domains with disks, an advancing front approach. Int. J. Numer. Methods Eng. 2003, 56, 699–713. [Google Scholar] [CrossRef]
  19. Munjiza, A.; Andrews, K. NBS contact detection algorithm for bodies of similar size. Int. J. Numer. Methods Eng. 1998, 43, 131–149. [Google Scholar] [CrossRef]
  20. Fang, X.; Liu, Z.; Tang, J. Particle stacking algorithm. J. Zhejiang Univ. 2011, 45, 650–655. [Google Scholar]
  21. Bandara, K.; Ranjith, P.; Rathnaweera, T. Proppant crushing mechanisms under reservoir conditions, insights into long-term integrity of unconventional energy production. Nat. Resour. Res. 2019, 28, 1139–1161. [Google Scholar] [CrossRef]
  22. Gai, G.; Tao, Z.; Ding, M. Powder Engineering; Tsinghua University Publishing House: Tsinghua, China, 2009. [Google Scholar]
  23. Xie, C.; Ma, H.; Zhao, Y. Investigation of modeling non-spherical particles by using spherical discrete element model with rolling friction. Eng. Anal. Bound. Elem. 2019, 105, 207–220. [Google Scholar] [CrossRef]
Figure 1. Subgraph (a) shows the accumulation of particles that diffuse outward from the center. Subgraph (b) shows the packing method in which particles are filled from the outside to the inside. Subgraph (c) shows the gradual accumulation of particles from bottom to top.
Figure 1. Subgraph (a) shows the accumulation of particles that diffuse outward from the center. Subgraph (b) shows the packing method in which particles are filled from the outside to the inside. Subgraph (c) shows the gradual accumulation of particles from bottom to top.
Applsci 14 05462 g001
Figure 2. Six three-dimensional arrangement forms of proppant particles.
Figure 2. Six three-dimensional arrangement forms of proppant particles.
Applsci 14 05462 g002
Figure 3. Included angle of proppant particles under different arrangement forms (the angle θ represents the inclination of the right side of a unit body with respect to the horizontal plane).
Figure 3. Included angle of proppant particles under different arrangement forms (the angle θ represents the inclination of the right side of a unit body with respect to the horizontal plane).
Applsci 14 05462 g003
Figure 4. Three-dimensional physical model of proppant particles.
Figure 4. Three-dimensional physical model of proppant particles.
Applsci 14 05462 g004
Figure 5. Arrangement of proppant particles in the bottom layer (k = 1) of three-dimensional space.
Figure 5. Arrangement of proppant particles in the bottom layer (k = 1) of three-dimensional space.
Applsci 14 05462 g005
Figure 6. Three-dimensional bottom layer simulation flow chart of single-particle-size proppant.
Figure 6. Three-dimensional bottom layer simulation flow chart of single-particle-size proppant.
Applsci 14 05462 g006
Figure 7. Simulation of other-layer proppant particles.
Figure 7. Simulation of other-layer proppant particles.
Applsci 14 05462 g007
Figure 8. Eight cases of proppant particles crossing the boundary.
Figure 8. Eight cases of proppant particles crossing the boundary.
Applsci 14 05462 g008
Figure 9. Schematic diagram of eight situations in which proppant particles cross the boundary (top view in the z-axis direction).
Figure 9. Schematic diagram of eight situations in which proppant particles cross the boundary (top view in the z-axis direction).
Applsci 14 05462 g009
Figure 10. Proppant particle simulation results.
Figure 10. Proppant particle simulation results.
Applsci 14 05462 g010
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Li, K.; Guo, D.; Guo, Z.; Zhao, Y. Computer 3D Simulation of Proppant Particles. Appl. Sci. 2024, 14, 5462. https://doi.org/10.3390/app14135462

AMA Style

Li K, Guo D, Guo Z, Zhao Y. Computer 3D Simulation of Proppant Particles. Applied Sciences. 2024; 14(13):5462. https://doi.org/10.3390/app14135462

Chicago/Turabian Style

Li, Ke, Dali Guo, Zixi Guo, and Yunxiang Zhao. 2024. "Computer 3D Simulation of Proppant Particles" Applied Sciences 14, no. 13: 5462. https://doi.org/10.3390/app14135462

APA Style

Li, K., Guo, D., Guo, Z., & Zhao, Y. (2024). Computer 3D Simulation of Proppant Particles. Applied Sciences, 14(13), 5462. https://doi.org/10.3390/app14135462

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop