Theoretical and experimental determination of proppant‐crushing ratio and fracture conductivity with different particle sizes

Crushing ratio and fracture conductivity are two key indicators for proppant evaluation. The most common way to evaluate these two indicators is physical testing, but physical testing is costly and time‐consuming, and the results are unstable. Therefore, a method for calculating proppant‐crushing ratio and fracture conductivity is proposed to replace repetitive physical experiments. However, in past attempts, only mathematical models using the same particle size have been established, while the particle‐size distribution in actual proppant arrangement is within a certain numerical range. Therefore, in order to establish a mathematical model more consistent with the actual situation of proppant arrangement, firstly, the Monte Carlo method was used to simulate proppant particles arranged with different particle sizes, and then a mathematical model was established describing the force of proppant particles. Secondly, according to the sparse, ill‐conditioned, and asymmetric characteristics of the coefficient matrix, an improved singular value decomposition method is proposed to improve the accuracy and efficiency of calculation. Finally, the crushing ratio and fracture conductivity of proppant were calculated through the new model in this paper. It was found that the results calculated by the new model in this paper are closer to actual test results.


| INTRODUCTION
Hydraulic fracturing is an essential means of oilfield stimulation. 1,2 Proppant is a key factor in determining the success or failure of fracturing construction. 3,4 A significant amount of proppant is used in hydraulic fracturing. [5][6][7] The currently low price of oil has a serious negative impact on the production of various oil fields. 8 For example, in the exploitation of tight oil reservoirs, existing fracturing technology can produce economic benefits at high oil prices, but at low oil prices, the development cost of tight oil and gas reservoirs is greater than the benefit of the output. 7 Therefore, one of the most urgent tasks in developing oil and gas fields is to "reduce cost and increase efficiency". 8 Proppant-crushing ratio and proppant fracture conductivity are two main indicators of proppant performance, 9,10 as well as being important parameters in proppant selection, 11,12 fracture design, 13,14 and postfracture evaluation. 15,16 Proppant-crushing ratio and fracture conductivity are generally tested by physical experiments, 15,[17][18][19] but physical experiments have limitations: long testing time, high cost, limited testing conditions, and significantly different results from repeated experiments under the same conditions. 20,21 An economical approach is to solve the determination of proppant-crushing ratio and fracture conductivity numerically. In 2020, Guo et al. proposed a method to calculate proppant-crushing ratio and fracture conductivity that was supported by a simple and economical experiment of crushing ratio. Through numerical calculation, proppant-crushing ratio and fracture conductivity could be obtained under various experimental conditions. 22 In the mathematical model of this study, the same particle size was assumed in proppant arrangement. However, in actual proppant arrangement, proppant particle sizes are different and are distributed within a certain numerical range. 23,24 Taking LZ sand and XM ceramsite as examples, in Figure 1, particle sizes were shown as distributed within a certain range. This study showed that simulation of proppant particle arrangement using different particle sizes is required.
In the present work, numerical simulation of proppant particle arrangement using different particle sizes was modified based on Guo's models and study of particle placement in the fields of mathematics and materials. [25][26][27] This simulation was used to reveal the arrangement of proppant particles with different sizes. The mathematical model was established by the mechanical relationship between proppant particles. By comparing the calculation results with the original model, the advantages of the new model proposed in this paper were analyzed.
The first research work presented by this paper is simulation of proppant particles with different sizes. The Monte Carlo method and an optimization model were used to deal with the random phenomena and determination phenomena of proppant particles. [28][29][30] The traversing search algorithm was used to find the center of a circle with the smallest coordinate in order to simulate proppant particles arranged with different sizes. [31][32][33] The second research work presented is establishment of a mathematical model of the mechanical relationship between proppant particles. By introducing a contact coefficient between particles and between particles and wall surface, along with considering the three mechanical relationships of contact force, reaction force, and boundary force, a mathematical model was established to describe the mechanical relationship between proppant particles. The coefficient matrix of the mathematical model is a system of linear equations with sparse, ill-conditioned, and asymmetrical features.
The third research work introduced is the proposal of an improved numerical method for solving mathematical models. Methods for solving linear equations include direct and iterative methods. 34 Direct methods mainly include the Gaussian elimination method, 35 trigonometric decomposition method, 36 etc. and cannot be used to solve an ill-conditioned system of linear equations. Commonly used iterative algorithms include conjugate gradient (CG), 37 preconditioned conjugate gradient (PCG), 38 minimal residual (MINRES), 39 generalized minimal residual (GMRES), 40 least-squares (LS), 41 and singular value decomposition (SVD) methods, 42 as well as various transformation algorithms to solve the approximate solution of the linear equations. Among them, CG and PCG methods F I G U R E 1 Schematic diagram of fracture and proppant particles require that linear equations have symmetric positive definiteness. 43 The MINRES method requires that linear equations have symmetry. 44 The GMRES and LS methods can be used to solve asymmetric, positive definite, sparse, and ill-conditioned equations, but for severely ill-conditioned linear equations, the calculation error is high, and the calculation time is long. 45 SVD and related improvement methods are the better methods for solving sparse, ill-conditioned linear equations. 46 Common SVD methods include truncated singular value decomposition (TSVD) 47 and modified singular value decomposition (MSVD). 48 TSVD ignores the lower singular values and singular matrix, resulting in loss of calculation results. 49 Although MSVD corrects low singular values, this method distorts the determined components of linear equations, causing significant deviations in calculation results. 50 In this paper, a new, improved singular value decomposition (ISVD) is proposed that not only retains the characteristics of the calculation model, but also improves the accuracy of calculation results.
The fourth research work presented is accuracy comparison of the simulation results for proppant particle arrangement with different particle sizes and with the same particle size.
The rest of this paper is organized as follows: in the second section, a new model and simulation are established for proppant particle arrangement with different sizes; the third section presents the results and discussion; and the fourth section presents the conclusions.

FOR PROPPANT PARTICLE ARRANGEMENT WITH DIFFERENT SIZES
Firstly, a new model of proppant particle arrangement with different particle sizes was established.

| Physical model
Based on the placement pattern of proppant particles in a hydraulic fracture and the arrangement structure of the proppant particles' experimental apparatus, a physical model was established for differently sized proppant particles, as shown in Figure 1. In the x direction is the fracture width, w f , and in the y direction is the fracture length, L f . The force, F c , is on the left and right boundaries of the fracture, and the upper and lower boundaries are closed. Based on the mathematical model of mechanical relations between proppant particles, the force of each proppant particle can be calculated.

| Modeling of proppant particle arrangement with different sizes
In order to perform modeling and simulation, we assumed the proppant particles to be circular. Under boundary force F, the force of the left and right boundaries changes in the parallel direction, but the force of the vertical direction (upper and lower boundaries) remains unchanged. The rules of particle arrangement are as follows 17,20-22 : 1. Particles are laid from bottom to top. 2. The placement of each particle is determined to be the shortest distance from the proppant particle center to the bottom of the fracture. 3. Deterministic placement occurs when the location with the shortest distance between the proppant particle center and the bottom of the fracture is unique. If the position is not unique, it is placed randomly.

| Random phenomena and treatment methods
A coordinate system was established, which is shown in Figure 2. When the simulation environment is in a crushing cell, the x axis represents the length of the crushing cell. The y axis represents the diameter of the crushing cell. When the simulation environment is in a conductivity cell, the x axis represents the width of the conductivity cell. The y axis represents the length of the conductivity cell. The coordinate of each circle is defined as (x i , y i ), where i is the serial number of the circle. The width of the fracture is w f , and the length of the fracture is L f . This is different from the assumption of Guo et al. that all particles are of the same diameter. The range of diameter of proppant is [R min , R max ] and R min < R max . The Monte Carlo method was used to generate random numbers, ξ i , representing the random numbers of proppant particle positions in contact with the lower boundary. The range of random numbers is [0, 1].
For any circle, the proppant particle is centered (x i , y i ); that is, The total number of proppant particles in contact with the lower boundary is n x . The diameter of proppant is R i . As shown in Figure 2A, n x is uncertain given the uncertainty of R i for each proppant particle. But we can determine the interval of n x ; that is, Then, the center coordinates (x i , y i ) of proppant particles in contact with the lower boundary are calculated by the following formula: The set of coordinates of the center of the circle in contact with the lower boundary is

| Deterministic phenomena and treatment methods
Let Q be the set of all centers and T be a temporary dynamic matrix. When the total number of proppant particles, n p , is less than or equal to the maximum number of first layers along the x axis, n x (n p ≤ n x ), T = Q, as shown in Figure 2B. When the total number of proppant particles, n p , is greater than the maximum number of first layers along the x axis, n x (n p ≤ n x ), T is the set of the highest center of the circle (green circle) in Figure 2.
Then, the center coordinates, Q, of all particles can be expressed as where n p is the total number of proppant particles, dimensionless, R np is the diameter of particles, and R np ∈ [R min , R max ], mm.
Assuming that the coordinate of the new proppant particle center is (x i ′, y i ′), the new circle center and other circle centers can be arranged in three cases as follows.
(1) When the particle is tangent to the left boundary, as shown in Figure 2C, the center coordinates of the particles that may be placed are calculated by the following formula: (2) When the particle is not tangent to the left and right boundaries (the middle position), as shown in Figure 2D, there are some particles in the middle position. By searching for the position of the minimum distance between the center of the circle and the bottom boundary, the center The green circle represents the highest circle. Red circles represent newly generated circles. The yellow circle represents the other circles. (A-F) represent six different kinds of proppant layering y coordinates of particles can be calculated by the following formula: (3) When the particle is tangent to the right boundary, as shown in Figure 2E, the center coordinates of the particles that may be placed are calculated by the following formula: In summary, according to the three arrangements of particles, the objective function is to find the position of the minimum distance between the center of the circle and the bottom boundary: where (x k , y k ) is the center coordinates of the center set Q, and R k is the diameter of a circle in the center set Q.

| Force analysis of proppant
Different from Guo et al.'s model, proppant arrangement with different diameters is more complex. When applying mechanical analysis and establishing a mathematical model among proppant particles, it was important to consider and distinguish the forces between particles and particles and between particles and boundaries. According to the force condition, the models were built separately. Any particle and any other particles may or may not have been in contact with each other. Therefore, we started with mathematical expression of proppant particle contact and then established the model of force between proppant particles.

| Mathematical expression of proppant particle contact
For the convenience of expression, we introduce contact coefficients δ ix , δ iy , and δ ij .

| Modeling of force between proppant particles
The mechanical relationship between particles and particles and between particles and boundaries includes three aspects, described in the following paragraphs.
(1) Force condition of proppant particles. Any particle may or may not contact other particles and boundaries. F j i represents the force of contact between the ith circle and the jth circle. Any circle may be in contact with any other circle (δ ij = 1) or may not (δ ij = 0). For n circles, the sum force is ∑ n j=1 ij F j i . At the same time, the circle may or may not be in contact with four boundaries (upper and lower boundaries, left and right boundaries). Therefore, for any circle in contact with other circles and boundaries, the expression is Then, the mathematical model of n proppant particles contacting each other is as follows: (2) Force and reaction between proppant particles. The two particles in contact with each other have a relationship between force and reaction force. The relationship between the force and reaction force of each particle and other particles is as follows: where i = 1, 2, 3, ···, n p ; j = 1, 2, 3, ···, n p ; and i ≠ j.
(3) Fracture boundary forces. When the simulation environment is in a conductivity cell, the force, F c , is in the left and right boundary directions, and the forces acting on the left boundary (x = 0) and the right bound- and F x=w f c , respectively. Considering the characteristics of two-dimensional simulation, the force acting on the left and right boundaries is as follows: The force acting on the boundaries is as shown: When the simulation environment is in a crushing cell, the force, F c , is in the upper and lower boundary directions, and the forces acting on the lower boundary (y = 0) and the upper boundary (y = L f ) are F y=0 c and F y=L f c , respectively. Considering the characteristics of twodimensional simulation, the force acting on the upper and lower boundaries is as follows: The force acting on the boundaries is as shown: According to the three cases of mechanical relationship between particles and particles and between particles and boundaries, the mathematical model can be expressed by the following linear set of equations: where A is the coefficient matrix of linear equations that is sparse, ill-conditioned, and asymmetric; and X and B are the vectors composed of unknowns and right-end constant terms, respectively.

| Numerical algorithm of mathematical model
The ISVD method retains the original value of the reliable part of the singular value to participate in the calculation. After correcting the singular value of the unreliable part, it is transformed into a reliable singular value to participate in the calculation.
The calculation process is as follows.
The linear set of equations AX = B has coefficient matrix A as a sparse, ill-conditioned, asymmetric matrix of n 1 × n 2 , such that Where U is an n 1 × n 1 order matrix; is an n 1 × n 2 order matrix; V is an n 2 × n 2 order matrix; UU T = I, VV T = I; and Taking threshold , singular value σ i is divided into two parts. With 1 ~ k in σ i , σ 1~k ≥ ε, and σ 1~k , this is the reliable part of the calculation. With 1+ k ~ n 2 in σ i , σ k+1 < ε, σ k+1~n2 → 0, and σ k+1~n2 , this is the unreliable part of the calculation. The relationship expression is as follows: Singular value in diagonal matrix σ i is arranged in a monotonically decreasing manner (shown in Figure 3), and by correcting k + 1 ~ n 2 of matrix , converts to the following: Then, Equation (22) transforms into Therefore, the method for correcting unreliable singular value A is the key to calculating linear equations.
The method takes the 2 ~ k -1 value from singular values σ k to σ 1 . The LS method is used for fitting, and there are at most k -1 fitting formulas. In the fitting formula, the fitting formula with the steepest slope is found, with the purpose of making the modified singular value σ k+1 ′ ~ σ n2 ′ calculated by the fitting formula as high as 0. As shown in Figure 4, when three singular values are selected, (that is, the k -2, k -1, and k values are selected), the slope of the fitted linear equation is the steepest.
In order to ensure that calculated modified singular value σ k+1 ′ ~ σ n2 ′ was optimal, we performed a formula fitting of the first to fourth powers for the selected data points, as shown in Figure 5. The optimum selection rules are as follows: 1. Calculated modified singular values σ k+1 ′~ σ n2 ′ are arranged in a monotonically decreasing manner. 2. Calculated modified singular values σ k+1 ′ ~ σ n2 ′ are all greater than 0. 3. Modified singular value σ k+1 ′ ~ σ n2 ′ of the first to fourth power is the highest. The pseudocode of the ISVD scheme is shown in Appendix 1.

| Methods for calculating proppantcrushing ratio and fracture conductivity
For proppant particle arrangement with different particle sizes and with the same size, the experimental tests were similar, but the calculation methods were different.

Proppant-crushing ratio
The crushing experiment of proppant particles in regular arrangement refers to the crushing experiment that meets the condition n x = w f / R. Under this condition, the proppant particles were arranged in a regular diamond-shaped structure with uniform forces on the particles, excluding the influence of randomness, which can truly represent the basic properties of the proppant. Crushing cells with various diameters were required to conduct the crushing experiment on the regular arrangement, which is suitable for various proppant particle sizes. According to SY/T 5108-2006 51 and API RP56, 52 the diameter of the crushing Linear equation fitting point and slope relationship cell should be between 38.1 and 76.2 mm, which is an integral multiple of proppant particle size. This can be done easily and at low cost. Experimental materials: We used LZ sand and XY ceramsite, which are widely used in China. Both proppants were 20/40 mesh, with proppant particle sizes ranging from 0.9 to 0.45 mm and an average particle size of 0.675 mm. We selected a crushing cell with a diameter of 40.5 mm (60 times the proppant particle size), as shown in Figure 6A. According to API SY/T 5108-2006 experimental standards, the crushing experiment was carried out using independently developed testing equipment rated for proppant crushing.
Experimental setup: testing equipment rated for proppant crushing, a standard sample sieve, sieve shaker, and high-precision balance.
Experimental procedures: (I) pour proppant into the standard sample sieve and place it on the sieve shaker for 30 min to ensure that all proppant samples fall within the particle-size range; (II) determine the quality of proppant-crushing test according to the calculation formula of proppant sample quality in SY/T 5108-2006 and API RP56 standards; (III) arrange the samples in the crushing room according to the rules; (IV) apply a closure pressure of 10 to 60 MPa, take the values at 5-MPa intervals, and test for 2 min at each pressure point; (V) for the proppant particles after pressure, put the standard sample sieve on the sieve shaker for 30 min again; and (VI) weigh the proppant mass on the standard sample sieve and chassis in turn, and calculate the crushing ratio of proppant. The crushing experiment results are shown in Figure 7A.

Crushing ratio
Crushing ratio is defined as the proportion of the mass of the crushing sample to the total mass under a certain pressure. Crushing ratio η can be expressed as where m c is the mass of the crushing sample, and g. m p is the total mass, g. Several standard proppant screens of different sizes can be taken to screen the proppant particles. The specific number of proppant screens and the size of proppant screens are determined according to the actual situation, and corresponding examples are given in the subsequent discussion section. After the screening, the proppant particles were divided into several groups of different sizes. For each group, under different closure pressures, the simulation and mechanical analysis of Guo's model (the same particle size) were applied successively to obtain the force condition of each proppant particle in the crushing cell under a given closure pressure, and the average value was taken. 22 The total crushing ratio curve was obtained by weighted average of the crushing ratio of each group, as shown in Figure 7.
Then, in accordance with the sample curve and numerical algorithm, each proppant particle crushing ratio, η i , can be obtained. Finally, the proppant-crushing ratio under a given closure pressure, η, is: where ‼ is the average breakage coefficient of proppant, dimensionless, and η i is the fracture coefficient of proppant, dimensionless.

| Fracture conductivity
In Guo et al.'s model, fracture permeability was expressed as 22 In this work, because the diameter of the particles is different, the fracture permeability is expressed as where where i is the number of proppant particles, and I = 1, 2, … , n p . Fracture porosity is related to the total fracture area and the area occupied by proppant particles, so fracture porosity is related to proppant-crushing ratio. The following two cases are discussed individually.
(1) Uncrushed proppant particles (η = 0) Proppant particles do not get crushed in fractures when no fracture closure pressure is applied. In this case, fracture porosity is the ratio of the pore area in the fracture to the total fracture area; that is, (2) Crushed proppant particles (η ≠ 0) When fracture closure pressure is applied, proppant particles get crushed in the fracture. In this case, physical observations after the experiment of proppant-crushing ratio and the fracture conductivity test showed the crushed proppant particles to be mainly in the form of fine powder, filling the pores between the unbroken proppant particles and mostly being clustered at the fracture boundary. Therefore, the formula for calculating fracture porosity is where A c is the area of the crushed proppant particle; that is, By substituting Equation (31) into Equation (30), the formula for calculating fracture porosity is as follows: Meanwhile, the fracture width decreases under the closure pressure. The formula to calculate the fracture width is where

| RESULTS AND DISCUSSION
We discuss the method of linear equations, the selection of simulation times, and the comparison between the simulation model of proppant particles with different sizes and the proppant particles with the same size.

| Method of linear equations
We compared and verified three methods: TSVD, MSVD, and ISVD. The indicators of comparison and verification were accuracy, calculation time, and scope of application. The verification method and the correct standard perform as follows. Firstly, a set was taken randomly of hypothetical value X * to apply in Equation (20) to calculate B. The calculated B was set as the known data. Then, three methods were used to calculate X given A and B. Finally, the error between X and X * was calculated.

| Accuracy
The calculation results are shown in Table 1. Because TSVD neglected the lower singular values and corresponding eigenvectors, it lost part of the calculation results and caused error. The average error was 23.17% between approximate results X T and X * using TSVD. MSVD corrected the lower singular values, but also corrected the higher singular values, so the whole singular value increased, which increased the deviation of the calculation results significantly. The average error was 36.49% between approximate results X M and X * using MSVD. ISVD retained a high singular value and corrected the lower singular value so that the corrected singular value was slightly higher than the original lower singular value. The average error was 2.12% between approximate results X I and X using ISVD. ISVD not only retained the original characteristics of the calculation model, but also greatly reduced the error in calculation results.

| Calculation speed
As shown in Table 2, the calculation times of the three methods were statistically analyzed. TSVD ignored smaller  singular values and corresponding eigenvectors and, therefore, had the shortest computation time. MSVD corrected all singular values, so the calculation time was the longest. ISVD retained a large singular value and corrected a small singular value. The calculation time of ISVD was slightly longer than TSVD and much shorter than MSVD.

| Scope of application
The scope of application for these three methods was verified by different condition numbers, as shown in Table 3. When the condition number was less than 1000, all three methods were able to calculate the ill-conditioned coefficient matrix. When the number of conditions was much greater than 1000, MSVD could not calculate the results.
To sum up, TSVD could calculate the results with various condition numbers and a short calculation time, but the calculation error of TSVD was high. MSVD could not calculate the results when the number of conditions was much greater than 1000. In addition, the calculation time of MSVD was the longest. Compared with TSVD and MSVD, the calculation results of ISVD were more accurate, and the calculation speed and scope of application were acceptable. Therefore, ISVD was used to solve linear equations in this study.

| Selection of simulation times
We studied the convergence of calculation times by calculating the error between the cumulative average and the total average. As shown in Figure 8, the higher the number of simulations, the lower the error between the cumulative average and the total average. When the number of simulations reached 22 times, the error between cumulative average and total average was less than 1%. When the number of simulations continued to increase, the error changed a little. Therefore, in the simulation and calculation, when the number of simulation calculations was set to 22 times, accurate results could be calculated at a given length and width.
It can be seen from Figure 9 that the particle size of proppant was distributed and that different types of proppant had different degrees of particle-size distribution. Taking the above three proppants as an example, the particle-size distribution of XY ceramsite was the most concentrated. The particle-size distribution of XM ceramsite was more dispersed, and the particle-size distribution of LZ sand was the most dispersed.

| Simulation of proppant arrangement
The model proposed in this paper was used to simulate proppant particle arrangement with different sizes. The simulation results are shown in Figure 10.
Then, force analysis was conducted. Heat maps of proppant stress are shown in Figure 11. It can be seen from Figure 11 that (1)  In order to compare the calculation accuracy between the particle-size distribution model proposed in this paper and the model with particles of the same size, we used these two models to calculate the crushing ratio and conductivity, respectively, and compared them with the  22 The results are shown in Figure 12 and Table 4.
From Figure 12 and Table 4, it can be seen that the accuracy of the simulation model of proppant particles with different sizes was higher than that of the simulation model of proppant particles of the same size. In the case of relatively concentrated distribution, the accuracy of the two models was high. However, the more dispersed the distribution, the less accurate the same-size model. In contrast, the distributed-size model achieved an accuracy of more than 85% under any condition.

| CONCLUSION
The arrangement and particle-size distribution of proppant particles have a significant impact on oil and gas production. We studied the simulation of proppant particle arrangement with different particle sizes. Through research, the following conclusions were drawn:    1. Based on the Monte Carlo method, the randomness and certainty of the arrangement of particles were processed separately. The traversing search method was used to find the position of the minimum distance between the center of the particle and the bottom boundary so as to realize computer simulation of proppant particle arrangement and to establish a mathematical model of the force analysis of proppant particles. 2. For the coefficient matrix of a mathematical model with sparse, ill-conditioned, and asymmetrical characteristics, this paper proposes ISVD, which has a high accuracy, short calculation time, and wide application range, and it numerically solves the mathematical model. 3. The more times the simulation and calculation are performed, the smaller the error between the cumulative average and the total average. When the number of simulations reaches 22, the error between the cumulative average and the total average is less than 1%, and when the number of simulations continues to increase, the error changes a little. 4. The particle size of proppant is distributed, and different types of proppant have different degrees of particlesize distribution. Through example analysis, it was shown that because the random arrangement of particles is different, the force condition is different. The arrangement of particles has a significant effect on the force of the particles. 5. The accuracy of the simulation model of proppant particles arranged with different sizes is higher than that of the simulation model of proppant particles arranged with the same size. The more dispersed the distribution, the more obvious the advantage of calculation accuracy of size distribution model.
As a basic work, this paper simplified the model to a two-dimensional model. The next step should be to establish a three-dimensional model to better simulate the actual situation. Three-dimensional problems will encounter many challenges; for instance, a sharp increase in the amount of calculation and data requires more advanced algorithms, and changes occur in the calculation methods of porosity and specific surface in the model formula after the two-dimensional circle becomes a threedimensional ball.