Study on Meso-Structure Evolution in Granular Matters Based on the Contact Loop Recognition and Determination Technique

On the mesoscopic scale, granular matter is tessellated into contact loops by a contact network. The stability of granular matter is highly dependent on the evolution of contact loops, including the number and area evolutions of contact loops with different geometric shapes (which can reflect the mechanical variables in the macroscale). For the features of numerous loops with complex geometry shapes in contact network images, a contact loop recognition and determination technique was developed in this study. Then, numerical biaxial compression tests were performed by the discrete element method (DEM) to investigate how the meso-structural indexes evolve along with the macro-mechanical indexes. The results show that the proposed Q-Y algorithm is effective in determining the geometric types of contact loops from contact network images. The evolution of contact loops is most active in the hardening stage, during which the number percentages of L3 (loops with three sides) and L6+ (loops with six or more sides) show opposite evolution patterns. For the area percentage, only L6+ increases while others decrease. Considering the meso-structural indexes (number percentage and area percentage of loops) are sensitive to the change of macro-mechanical indexes (deviatoric stress, axial strain, and volumetric strain) in the hardening stage. Multivariate models were established to build a bridge between the meso-structure and the macro-mechanics.


Introduction
The analysis of microstructure (such as particle rolling, contact sliding, coordination number, and contact anisotropy, etc.) is helpful to explain the complex macroscopic behaviors (such as strength and deformation characteristics) of granular materials [1][2][3][4][5]. However, establishing the relationship between microscopic variables and macroscopic behavior requires an intermediate scale, namely the mesoscale [6,7]. It has been widely accepted that the geometric arrangement and combination of particles determine the macroscopic mechanical behaviors of particle assembly [8]. At the mesoscale, structural features are accounted for with clusters of a few interacting particles. Therefore, the structure of particle assembly can be divided into elementary "contact loops" (hereinafter referred to as loops), which are enclosed by the contact branches of interacting particles [9]. The stability of assembly is highly dependent on the evolution of contact loops, which can help understand the underlying mechanisms behind macroscopic observations [10].
The edge number of the contact loop plays an important role in the mechanical behavior of assembly. For example, the increase in edge number will increase the freedom degree of contact loops, causing a complex mechanical response of assembly [11][12][13][14]. However, how to identify contact loops and determine their edge number from complex contact networks has not been specifically described in relevant studies. For the assembly containing a large number of particles, the contact network has the characteristics of numerous nodes and complex loop shape, which requires a higher recognition accuracy of the contact loops. Based on the contact network, similar recognition methods of polygon loops include artificial recognition, classical image recognition techniques, and deep learning recognition methods (which developed rapidly in recent years) [15]. However, artificial recognition is inefficient and is not suitable for contact networks with a large number of loops. Classical image recognition techniques have been developed and have high accuracy for images with simple and good semantic results, but their algorithms need to be improved according to the image characteristics. The deep learning recognition method can recognize polygon types, but it needs a lot of samples to be labeled and trained before it can be applied [16]. Furthermore, the recognition of single or small numbers of images of contact loops with complex geometric shapes is inefficient. Therefore, based on the characteristics of the contact network, a new recognition method needs to be proposed to recognize the contact loop and judge its geometric type.
To understand the evolution of contact loops, some studies have been carried out by the discrete element method (DEM) [7,8,17,18]. For example, the shear dilatancy in the biaxial compression test can be interpreted as the transformation of contact loops from the small-dense structure to the large-loose structure [7]. However, the evolution representation of contact loops is often qualitative and lacks quantitative representation, especially the quantitative relationship between the structural indexes of contact loops at mesoscale and the mechanical indexes of assembly at the macroscale. The purpose of this study is to put forward a determination technique of contact loops with different geometric shapes and to establish the quantitative relationship between structural indexes and mechanical indexes. The structure of the paper is as follows. In Section 2, the recognition and determination technique for contact loops based on contact network images is introduced. Section 3 analyzes the evolution of meso-structural indexes in the biaxial compression test based on the above techniques. In Section 4, the quantitative relationship between meso-structural indexes and macro-mechanical indexes was established, and thus the bridge between meso-structure and macro-mechanics was also established. Finally, the conclusion and prospect of future work are given in Section 5.

Contact Loop Recognition and Determination Technique
For the DEM simulation with contact data, the contact network can be discretized into sub-domains (polygonal loops) by the "Delaunay triangulation" method and other analytical methods [7,17]. However, these methods are not suitable for cases without original contact data, such as contact network images obtained by 2D simulations and experiments ( Figure 1). It is necessary to propose a 2D contact loop recognition and determination technique from the perspective of contact network images, which lacks original contact data. Therefore, one of the purposes of this study is to recognize and determine contact loops in contact network images without contact data. numerous nodes and complex loop shape, which requires a higher recognition accuracy of the contact loops. Based on the contact network, similar recognition methods of polygon loops include artificial recognition, classical image recognition techniques, and deep learning recognition methods (which developed rapidly in recent years) [15]. However, artificial recognition is inefficient and is not suitable for contact networks with a large number of loops. Classical image recognition techniques have been developed and have high accuracy for images with simple and good semantic results, but their algorithms need to be improved according to the image characteristics. The deep learning recognition method can recognize polygon types, but it needs a lot of samples to be labeled and trained before it can be applied [16]. Furthermore, the recognition of single or small numbers of images of contact loops with complex geometric shapes is inefficient. Therefore, based on the characteristics of the contact network, a new recognition method needs to be proposed to recognize the contact loop and judge its geometric type.
To understand the evolution of contact loops, some studies have been carried out by the discrete element method (DEM) [7,8,17,18]. For example, the shear dilatancy in the biaxial compression test can be interpreted as the transformation of contact loops from the small-dense structure to the large-loose structure [7]. However, the evolution representation of contact loops is often qualitative and lacks quantitative representation, especially the quantitative relationship between the structural indexes of contact loops at mesoscale and the mechanical indexes of assembly at the macroscale. The purpose of this study is to put forward a determination technique of contact loops with different geometric shapes and to establish the quantitative relationship between structural indexes and mechanical indexes. The structure of the paper is as follows. In Section 2, the recognition and determination technique for contact loops based on contact network images is introduced. Section 3 analyzes the evolution of meso-structural indexes in the biaxial compression test based on the above techniques. In Section 4, the quantitative relationship between meso-structural indexes and macro-mechanical indexes was established, and thus the bridge between meso-structure and macro-mechanics was also established. Finally, the conclusion and prospect of future work are given in Section 5.

Contact Loop Recognition and Determination Technique
For the DEM simulation with contact data, the contact network can be discretized into sub-domains (polygonal loops) by the "Delaunay triangulation" method and other analytical methods [7,17]. However, these methods are not suitable for cases without original contact data, such as contact network images obtained by 2D simulations and experiments ( Figure 1). It is necessary to propose a 2D contact loop recognition and determination technique from the perspective of contact network images, which lacks original contact data. Therefore, one of the purposes of this study is to recognize and determine contact loops in contact network images without contact data.  There are three main reasons for developing the 2D recognition and determination technique. Firstly, contact forces in the 3D contact network are nonplanar. Such hierarchical relationships are difficult to represent in the image, which is a great challenge for the segmentation and identification of contact loops. The 2D technique in this study is an important basis for extending to future 3D researches. Secondly, much of the literature has demonstrated that 2D DEM simulations and idealized 2D experiments (such as Figure 1b) can reproduce the main mechanical features of mixed granular matters [19][20][21]. More importantly, the 2D model allows one to explore the microstructures in an effective but easier way than the 3D model [22]. For the evolution of the contact network, the effect of the 2D model is more intuitive. Therefore, the 2D recognition and determination technique is developed to investigate the meso-structure evolution of granular matters.
The principle of the 2D contact loop recognition and determination technique is that the corner number in the polygonal loop equals the side number of these polygonal loops. Therefore, the core of the technology is the recognition of corners. The existing corner detection technologies are mainly divided into gradient-based, template-based, and template-based gradient combinations. These algorithms include the FAST algorithm [23], HARRIS algorithm [24], SUSAN algorithm [25], and other improvement methods for corner detection. Based on the contact network image, we tested and compared the results of the above three corner detection algorithms.
In Figure 2, the FAST algorithm and the SUSAN algorithm recognize the points on the boundary as corners, resulting in a huge increase in the number of corners. However, the recognition effect of the HARRIS algorithm is the opposite and results in the loss of many corners. The traditional corner detection algorithm cannot obtain accurate corner information in the contact network image. Therefore, this study proposes a Q-Y algorithm to identify corners and then determine the geometric types of contact loops. The process of the recognition and determination technique based on the Q-Y algorithm is as follows. There are three main reasons for developing the 2D recognition and determination technique. Firstly, contact forces in the 3D contact network are nonplanar. Such hierarchical relationships are difficult to represent in the image, which is a great challenge for the segmentation and identification of contact loops. The 2D technique in this study is an important basis for extending to future 3D researches. Secondly, much of the literature has demonstrated that 2D DEM simulations and idealized 2D experiments (such as Figure  1b) can reproduce the main mechanical features of mixed granular matters [19][20][21]. More importantly, the 2D model allows one to explore the microstructures in an effective but easier way than the 3D model [22]. For the evolution of the contact network, the effect of the 2D model is more intuitive. Therefore, the 2D recognition and determination technique is developed to investigate the meso-structure evolution of granular matters.
The principle of the 2D contact loop recognition and determination technique is that the corner number in the polygonal loop equals the side number of these polygonal loops. Therefore, the core of the technology is the recognition of corners. The existing corner detection technologies are mainly divided into gradient-based, template-based, and template-based gradient combinations. These algorithms include the FAST algorithm [23], HARRIS algorithm [24], SUSAN algorithm [25], and other improvement methods for corner detection. Based on the contact network image, we tested and compared the results of the above three corner detection algorithms.
In Figure 2, the FAST algorithm and the SUSAN algorithm recognize the points on the boundary as corners, resulting in a huge increase in the number of corners. However, the recognition effect of the HARRIS algorithm is the opposite and results in the loss of many corners. The traditional corner detection algorithm cannot obtain accurate corner information in the contact network image. Therefore, this study proposes a Q-Y algorithm to identify corners and then determine the geometric types of contact loops. The process of the recognition and determination technique based on the Q-Y algorithm is as follows. Firstly, we segment the original contact network image and get separate connected domains in Image I, wherein each connected domain represents a contact loop; secondly, we assign different gray value to each connected domain and obtain Image II; thirdly, we perform the morphological opening operation to integrate the divided edges of the two adjacent connection domains to obtain an edge to form Image III; fourthly, we use the Q-Y corner detection algorithm to perform corner detection on Image III to determine each corner point and obtain Image IV; fifthly, we set a gray value uniformly for the adjacent corner points as an intersection, and extract Image V containing only the intersection point; sixthly, we perform morphological corrosion processing on each intersection of Image V, and each intersection intrudes into the adjacent connection domain to obtain Image VI; seventhly, the geometric type of each connection domain is determined according to the principle that the number of sides is equal to the number of corners. The area of the connected domains is calculated by the number of pixels in each connected domain, and the statistical results (number and area of each geometric type) are obtained in the final step. The above process is demonstrated in Figure 3. Firstly, we segment the original contact network image and get separate connected domains in Image I, wherein each connected domain represents a contact loop; secondly, we assign different gray value to each connected domain and obtain Image II; thirdly, we perform the morphological opening operation to integrate the divided edges of the two adjacent connection domains to obtain an edge to form Image III; fourthly, we use the Q-Y corner detection algorithm to perform corner detection on Image III to determine each corner point and obtain Image IV; fifthly, we set a gray value uniformly for the adjacent corner points as an intersection, and extract Image V containing only the intersection point; sixthly, we perform morphological corrosion processing on each intersection of Image V, and each intersection intrudes into the adjacent connection domain to obtain Image VI; seventhly, the geometric type of each connection domain is determined according to the principle that the number of sides is equal to the number of corners. The area of the connected domains is calculated by the number of pixels in each connected domain, and the statistical results (number and area of each geometric type) are obtained in the final step. The above process is demonstrated in Figure 3.

Contact Network Image Pre-Processing
The first step of the recognition and determination technique is image segmentation. That is, each contact loop is separated from the other in the contact network image. Each separated contact ring is regarded as a connected domain. The segmentation method based on the OTSU algorithm [26] has always been regarded as the optimal method for automatic image segmentation. The basic idea of this algorithm is to divide image pixels into two groups by a threshold, and then determine the optimal threshold by the maximum interclass variance between the pixels of two groups. Suppose where 0 E μ and 1 E μ are the expectations of G0 and G1, respectively; 0 α and 1 α are the probabilities of G0 and G1, respectively. Therefore, the interclass variance of the two groups can be expressed as ( ) , then * t is the optimal threshold. If the value * t is not unique, the average value of all * t is used as the optimal threshold. For the contact network image, the OTSU segmentation method gives a more satisfactory segmentation result, as shown in Figure 4.

Contact Network Image Pre-Processing
The first step of the recognition and determination technique is image segmentation. That is, each contact loop is separated from the other in the contact network image. Each separated contact ring is regarded as a connected domain. The segmentation method based on the OTSU algorithm [26] has always been regarded as the optimal method for automatic image segmentation. The basic idea of this algorithm is to divide image pixels into two groups by a threshold, and then determine the optimal threshold by the maximum interclass variance between the pixels of two groups.
Suppose the grey levels of the contact network image is G = [0, L − 1] and the probability of each grey level is P i . The threshold t divides the image into two groups G 0 = [0, t] and G 1 = [t + 1, L − 1]. The probabilities of the two groups are where µ E 0 and µ E 1 are the expectations of G 0 and G 1 , respectively; α 0 and α 1 are the probabilities of G 0 and G 1 , respectively. Therefore, the interclass variance of the two groups can be expressed as If η 2 (t * )= max η 2 (t) , then t * is the optimal threshold. If the value t * is not unique, the average value of all t * is used as the optimal threshold. For the contact network image, the OTSU segmentation method gives a more satisfactory segmentation result, as shown in Figure 4.
In the segmented image, different grayscales are assigned to the segmented connected domains. Since the existence of boundary lines of connected domains affects the effect of corner detection, the image needs to be processed using the algorithm of binary open operation to remove the boundary lines. The binary open operation includes corrosion calculation and expansion calculation, which is a multiple-point pattern-based unconditional simulation algorithm using morphological image processing tools [27,28].  In the segmented image, different grayscales are assigned to the segmented connected domains. Since the existence of boundary lines of connected domains affects the effect of corner detection, the image needs to be processed using the algorithm of binary open operation to remove the boundary lines. The binary open operation includes corrosion calculation and expansion calculation, which is a multiple-point pattern-based unconditional simulation algorithm using morphological image processing tools [27,28].
The corrosion calculation can cause a shrinkage erosion of the image's boundaries and is used to eliminate the small and meaningless areas. The definition of the corrosion calculation is the probing of an image B with a probe S to find a region E inside the image [29], which can be expressed as The expansion calculation is a pairwise operation of the corrosion calculation, which can be used to fill certain voids in the target area as well as to eliminate small particles of noise contained in the target area. The expansion calculation can be expressed as The corrosion calculation can corrode boundary lines, and the expansion calculation can be used to fill the cavity in the target area and eliminate the noise contained in the target area. The processing effect of the binary open operation is shown in Figure 5.

Corner Detection
For the greyscale assigned image without boundary lines, there is a greyscale inside the connected domain, two greyscales around the boundary (Figure 6a), and three or more grey values around the corner (Figure 6b). The Q-Y algorithm for corner detection is proposed based on this characteristic. The corrosion calculation can cause a shrinkage erosion of the image's boundaries and is used to eliminate the small and meaningless areas. The definition of the corrosion calculation is the probing of an image B with a probe S to find a region E inside the image [29], which can be expressed as The expansion calculation is a pairwise operation of the corrosion calculation, which can be used to fill certain voids in the target area as well as to eliminate small particles of noise contained in the target area. The expansion calculation can be expressed as The corrosion calculation can corrode boundary lines, and the expansion calculation can be used to fill the cavity in the target area and eliminate the noise contained in the target area. The processing effect of the binary open operation is shown in Figure 5.  In the segmented image, different grayscales are assigned to the segmented connected domains. Since the existence of boundary lines of connected domains affects the effect of corner detection, the image needs to be processed using the algorithm of binary open operation to remove the boundary lines. The binary open operation includes corrosion calculation and expansion calculation, which is a multiple-point pattern-based unconditional simulation algorithm using morphological image processing tools [27,28].
The corrosion calculation can cause a shrinkage erosion of the image's boundaries and is used to eliminate the small and meaningless areas. The definition of the corrosion calculation is the probing of an image B with a probe S to find a region E inside the image [29], which can be expressed as The expansion calculation is a pairwise operation of the corrosion calculation, which can be used to fill certain voids in the target area as well as to eliminate small particles of noise contained in the target area. The expansion calculation can be expressed as The corrosion calculation can corrode boundary lines, and the expansion calculation can be used to fill the cavity in the target area and eliminate the noise contained in the target area. The processing effect of the binary open operation is shown in Figure 5.

Corner Detection
For the greyscale assigned image without boundary lines, there is a greyscale inside the connected domain, two greyscales around the boundary (Figure 6a), and three or more grey values around the corner (Figure 6b). The Q-Y algorithm for corner detection is proposed based on this characteristic.

Corner Detection
For the greyscale assigned image without boundary lines, there is a greyscale inside the connected domain, two greyscales around the boundary (Figure 6a), and three or more grey values around the corner (Figure 6b). The Q-Y algorithm for corner detection is proposed based on this characteristic. If there are two or more grayscales within a certain range of the target pixel and they are different from its grayscale, the point is considered as a corner pixel. The specific process for corner pixel detection is as follows.
Step 1: Four pixels k i N around the target pixel i N are selected to check the grayscales, and the spatial positions of N and k N are shown in Figure 7. If at least two If there are two or more grayscales within a certain range of the target pixel and they are different from its grayscale, the point is considered as a corner pixel. The specific process for corner pixel detection is as follows.
Step 1: Four pixels N k i around the target pixel N i are selected to check the grayscales, and the spatial positions of N i and N k i are shown in Figure 7. If at least two pixels have different grayscale from the target pixel, then the target pixel may be a corner pixel and proceed to the next step; if not, then the next target pixel is checked.  Figure 7. If at least two pixels have different grayscale from the target pixel, then the target pixel may be a corner pixel and proceed to the next step; if not, then the next target pixel is checked.
Step 2: The corner response function is used to check whether the pixels i N screened in Step 1 are corner pixels. If the corner response function is satisfied, the pixel i N is regarded as a corner pixel. The corner response function is expressed as: where i g is the grayscale of i N , and x i g is the grayscale of the pixel Step 3: The pixels meeting the corner response function constitute corner pixels. It is worth noting that each corner may contain multiple corner pixels, as shown in Figure 8.   Step 2: The corner response function is used to check whether the pixels N i screened in Step 1 are corner pixels. If the corner response function is satisfied, the pixel N i is regarded as a corner pixel. The corner response function is expressed as: where g i is the grayscale of N i , and g x i is the grayscale of the pixel N x i around N i . Step 3: The pixels meeting the corner response function constitute corner pixels. It is worth noting that each corner may contain multiple corner pixels, as shown in Figure 8. If there are two or more grayscales within a certain range of the target pixel and they are different from its grayscale, the point is considered as a corner pixel. The specific process for corner pixel detection is as follows.
Step 1: Four pixels k i N around the target pixel i N are selected to check the grayscales, and the spatial positions of i N and k i N are shown in Figure 7. If at least two pixels have different grayscale from the target pixel, then the target pixel may be a corner pixel and proceed to the next step; if not, then the next target pixel is checked.
Step 2: The corner response function is used to check whether the pixels i N screened in Step 1 are corner pixels. If the corner response function is satisfied, the pixel i N is regarded as a corner pixel. The corner response function is expressed as: where i g is the grayscale of i N , and x i g is the grayscale of the pixel Step 3: The pixels meeting the corner response function constitute corner pixels. It is worth noting that each corner may contain multiple corner pixels, as shown in Figure 8.   In Figure 8, the Q-Y algorithm identifies more than 99% of corners, which has a perfect recognition effect. For comparison, three common corner detection algorithms were selected again, whose results are shown in Figure 9. The FAST algorithm [23] identified only a few corners, the HARRIS algorithm [24] identified 80% of corners, and the SUSAN algorithm [25] identified corners including numerous non-corners and lost a large number of corners. The accuracy of these three algorithms is quite inaccurate compared to the Q-Y algorithm.

Recognition and Statistics
The adjacent corner pixels are assigned to the same value and considered as a corner.
The above binary open operation affected the range of corner pixels, resulting in a decrease in the accuracy of corner allocation. Therefore, the corrosion calculation is required for equals the number of its sides, the geometric types of contact loops can be determined. Additionally, the area of a contact loop is represented by the number of pixels. Finally, the number and area of contact loops of each geometric type can be calculated. The above process can be represented in Figure 10. In Figure 8, the Q-Y algorithm identifies more than 99% of corners, which has a perfect recognition effect. For comparison, three common corner detection algorithms were selected again, whose results are shown in Figure 9. The FAST algorithm [23] identified only a few corners, the HARRIS algorithm [24] identified 80% of corners, and the SUSAN algorithm [25] identified corners including numerous non-corners and lost a large number of corners. The accuracy of these three algorithms is quite inaccurate compared to the Q-Y algorithm. Figure 9. The corner recognition results of (a) FAST algorithm, (b) HARRIS algorithm, and (c) SUSAN algorithm for the greyscale assigned image without boundary lines (the bule, green, and red points are the recognized corners by the three algorithms respectively).

Recognition and Statistics
The adjacent corner pixels are assigned to the same value and considered as a corner.
The above binary open operation affected the range of corner pixels, resulting in a decrease in the accuracy of corner allocation. Therefore, the corrosion calculation is required for Figure 8 so that every corner point can be allocated to the surrounding connected domains. After all of the corners are allocated, the number of corners in all connected domains can be obtained. Based on the principle that the number of corners in the connected domain equals the number of its sides, the geometric types of contact loops can be determined. Additionally, the area of a contact loop is represented by the number of pixels. Finally, the number and area of contact loops of each geometric type can be calculated. The above process can be represented in Figure 10.

Recognition and Statistics
The adjacent corner pixels are assigned to the same value and considered as a corner.
The above binary open operation affected the range of corner pixels, resulting in a decrease in the accuracy of corner allocation. Therefore, the corrosion calculation is required for Figure 8 so that every corner point can be allocated to the surrounding connected domains. After all of the corners are allocated, the number of corners in all connected domains can be obtained. Based on the principle that the number of corners in the connected domain equals the number of its sides, the geometric types of contact loops can be determined. Additionally, the area of a contact loop is represented by the number of pixels. Finally, the number and area of contact loops of each geometric type can be calculated. The above process can be represented in Figure 10.  The 2D contact loop recognition and determination technique can accurately determine and count the geometric types and their areas of contact loops, which creates conditions for quantitative analysis of contact network images.

Simulation and Analysis
The biaxial compression test was taken as an example and obtains the force chain images through the DEM simulation in this section. Based on the recognition and determination technique introduced in Section 2, the meso-structural indexes were calculated, and the relationships between the meso-structural indexes and the macro-mechanical indexes are analyzed.

Biaxial Compression Numerical Model
Deluzarche and Cambou [30] indicated that volumetric contracting strains are difficult to obtain in 2D and suggested that 2D simulations should be restricted to dense materials. To better reflect the strain responses, the dense assembly is adopted to match this suggestion. The dense assembly is produced by isotropic compression. First, an initial rectangular area consisting of four rigid walls was established (Figure 11a), which contained 10,000 round particles with uniformly distributed size (the minimum and maximum diameters of particles are 0.8 mm and 1.2 mm, respectively). Then, particles in the initial rectangular area were compressed isotropically. The servo control mechanism was used to continuously adjust the positions of the rigid walls until a stable confining pressure (σ 0 = 20 kPa) was attained (Figure 11b). When the ratio of the mean unbalanced force to the mean contact force was less than 10 −5 , the assembly was considered to reach equilibrium, and a dense assembly was obtained. rectangular area consisting of four rigid walls was established (Figure 11a), which contained 10,000 round particles with uniformly distributed size (the minimum and maximum diameters of particles are 0.8mm and 1.2mm, respectively). Then, particles in the initial rectangular area were compressed isotropically. The servo control mechanism was used to continuously adjust the positions of the rigid walls until a stable confining pressure (σ0 = 20 kPa) was attained (Figure 11b). When the ratio of the mean unbalanced force to the mean contact force was less than 10 −5 , the assembly was considered to reach equilibrium, and a dense assembly was obtained.
After the dense assembly was generated, the biaxial compression was used to apply axial load. Biaxial compression means that confining pressures σ0 on the left and right walls remain constant by the servo control mechanism, and the axial load is applied by the downward movement of the upper wall (Figure 11c). The stress on the upper wall is the principal stress, which is denoted by σ1. In this study, the moving rate of the upper wall was maintained at 0.05 m/s. When the axial strain reaches 20%, it is considered that the particle has compression failure, and the loading stops. The biaxial compression simulations are conducted using the DEM program PFC 2D code [31]. The linear elastic contact model was used to describe the contact behavior between particles, whose parameters were obtained according to previous experiments and numerical simulations. Specifically, the interparticle friction coefficient (μp) is After the dense assembly was generated, the biaxial compression was used to apply axial load. Biaxial compression means that confining pressures σ 0 on the left and right walls remain constant by the servo control mechanism, and the axial load is applied by the downward movement of the upper wall (Figure 11c). The stress on the upper wall is the principal stress, which is denoted by σ 1 . In this study, the moving rate of the upper wall was maintained at 0.05 m/s. When the axial strain reaches 20%, it is considered that the particle has compression failure, and the loading stops.
The biaxial compression simulations are conducted using the DEM program PFC 2D code [31]. The linear elastic contact model was used to describe the contact behavior between particles, whose parameters were obtained according to previous experiments and numerical simulations. Specifically, the interparticle friction coefficient (µ p ) is obtained by consulting previous DEM simulations, and the wall-particle friction coefficient (µ w ) is set as 0 to ensure that the particles around walls can roll and slide without resistance. The value of normal-to-tangential stiffness ratio (k n /k s ) is adopted as 4/3, which belongs to the range of 1.0 to 1.5 of realistic granular materials [32]. Damping constant β = 0.7, as suggested in Itasca [31], has been used for effectively dissipating the kinetic energy. The values of the contact parameters are shown in Table 1.

Number Evolution of Loops
The variation of loop number N l , deviatoric stress q (q = σ 1 − σ 0 ), and volumetric strain ε v with axial strain ε a are shown in Figure 12. The deviatoric stress curve can be divided into the strain hardening stage and the strain-softening stage (Figure 12a). The curve of the loop number first decreases rapidly and then tends to be stable, and the inflection point appears with the peak stress, indicating that the change of contact loop is most active in the strain hardening stage. Figure 12b shows the assembly experienced the shear contraction stage (ε v ≤ 0) and the shear dilation stage (ε v >0). The change of loops in the shear contraction stage was more active than that in the shear dilation stage.

Number Evolution of Loops
The variation of loop number Nl, deviatoric stress q (q = σ1 − σ0), and volumetric strain v ε with axial strain a ε are shown in Figure 12. The deviatoric stress curve can be divided into the strain hardening stage and the strain-softening stage (Figure 12a). The curve of the loop number first decreases rapidly and then tends to be stable, and the inflection point appears with the peak stress, indicating that the change of contact loop is most active in the strain hardening stage. Figure 12b shows the assembly experienced the shear contraction stage ( v 0 ε ≤ ) and the shear dilation stage ( v 0 ε ＞ ). The change of loops in the shear contraction stage was more active than that in the shear dilation stage. In addition, the change of average coordination number Z was also explored. The evolution of Nl is highly similar to that of Z (Figure 13a). This phenomenon can be explained by Euler's relation of 2D topology. In the particle system, the relationship between the particle number Np, the contact number  In addition, the change of average coordination number Z was also explored. The evolution of N l is highly similar to that of Z (Figure 13a). This phenomenon can be explained by Euler's relation of 2D topology. In the particle system, the relationship between the particle number N p , the contact number N c , and the loop number N l can be expressed as N p + N l ∼ = N c [33]. Based on Z = 2N c /N p , the relationship between Z and N l can be expressed as N l ∼ = N p (Z/2 − 1). For the assembly with a constant particles number, the average coordination number is positively correlated with the loop number. In order to show the relationship between the average coordination number and the loop number better, we drew two diagrams and calculated the corresponding values of Z and N l , as shown in Figure 13b.

Evolution of the Meso-Structural Indexes of Loops
The contact network is a link between the strength and deformation of granular matters. The number and structure of contact loops are highly correlated with the strength of the material, which can be used as a medium for carrying contact forces. Additionally, the macroscopic deformation can be explained in terms of the evolution of contact loops, considering the local contact particles in a loop-like mosaic. The evolution of the contact loops consists of the evolution of the number and area of loops with different geometrical types, which can be used to characterize the changes of the whole contact network. Therefore, the number percentage and the area percentage of the contact loops with different geometrical types are defined as the meso-structural indexes in this study.

Evolution of the Meso-Structural Indexes of Loops
The contact network is a link between the strength and deformation of granular matters. The number and structure of contact loops are highly correlated with the strength of the material, which can be used as a medium for carrying contact forces. Additionally, the macroscopic deformation can be explained in terms of the evolution of contact loops, considering the local contact particles in a loop-like mosaic. The evolution of the contact loops consists of the evolution of the number and area of loops with different geometrical types, which can be used to characterize the changes of the whole contact network. Therefore, the number percentage and the area percentage of the contact loops with different geometrical types are defined as the meso-structural indexes in this study.
L i denotes the contact loop with side number i (i = 3, 4, 5, . . . ). The number percentage of L i is ω i = N i /N l , where N i is the number of L i . In this study, the loops with the side number ≥ 6 are grouped into one group based on previous studies [7], which is hereafter referred to as L 6+ . Therefore, L 3 , L 4 , L 5 , and L 6+ are counted separately in this paper and their number percentage ω 3 , ω 4 , ω 5 , and ω 6+ are calculated. The variation curves of deviatoric stress, volumetric strain, and number percentage indexes with axial strain are shown in Figure 14.

Evolution of the Meso-Structural Indexes of Loops
The contact network is a link between the strength and deformation of granular matters. The number and structure of contact loops are highly correlated with the strength of the material, which can be used as a medium for carrying contact forces. Additionally, the macroscopic deformation can be explained in terms of the evolution of contact loops, considering the local contact particles in a loop-like mosaic. The evolution of the contact loops consists of the evolution of the number and area of loops with different geometrical types, which can be used to characterize the changes of the whole contact network. Therefore, the number percentage and the area percentage of the contact loops with different geometrical types are defined as the meso-structural indexes in this study. i L denotes the contact loop with side number i (i = 3, 4, 5, ...). The number percentage In this study, the loops with the side number ≥ 6 are grouped into one group based on previous studies [7], which is hereafter referred to as L6+. Therefore, L3, L4, L5, and L6+ are counted separately in this paper and their number percentage ω3, ω4, ω5, and ω6+ are calculated. The variation curves of deviatoric stress, volumetric strain, and number percentage indexes with axial strain are shown in Figure 14. Before the peak stress, ω3 decreased significantly, while ω6+ showed an obvious opposite evolution. In the softening stage after the peak stress, the evolution rates of ω3 and ω6+ gradually stabilized. It is worth noting that the maximum curvature of ω3 and ω6+ coincide with the minimum volumetric strain. In other words, the evolution rates of ω3 and ω6+ have undergone tremendous changes in the conversion process of assembly from the dense structure to the loose structure. However, ω4 and ω5 are basically constants and are not affected by the evolution of deviatoric stress and volumetric strain. Before the peak stress, ω 3 decreased significantly, while ω 6+ showed an obvious opposite evolution. In the softening stage after the peak stress, the evolution rates of ω 3 and ω 6+ gradually stabilized. It is worth noting that the maximum curvature of ω 3 and ω 6+ coincide with the minimum volumetric strain. In other words, the evolution rates of ω 3 and ω 6+ have undergone tremendous changes in the conversion process of assembly from the dense structure to the loose structure. However, ω 4 and ω 5 are basically constants and are not affected by the evolution of deviatoric stress and volumetric strain.
The number percentage reflects the evolution of contact loops with different geometric types, but it cannot fully reflect the distribution degree of contact loops in the contact network. Therefore, it is necessary to introduce the area percentage of contact loops A i , which means the area ratio of L i in the contact network. A i is defined as A i = S i /S T , where S i and S T are the area of L i and the contact network, respectively. The variation curves of deviatoric stress, volumetric strain, and area percentage indexes (A 3 , A 4 , A 5 , and A 6+ ) with axial strain are shown in Figure 15. The evolution of A i shows that only L 6+ experiences an area increase while the other kinds of loop contract. In fact, having more sides gives rise to L 6+ the ability to transform. The evolution of contact loops actually represents the evolution of assembly volume in the 2D simulation. The evolution trend of A 6+ and ε v is similar, indicating that L 6+ is the main factor affecting volume evolution. The number percentage reflects the evolution of contact loops with different geometric types, but it cannot fully reflect the distribution degree of contact loops in the contact network. Therefore, it is necessary to introduce the area percentage of contact loops i A , which means the area ratio of i L in the contact network. i A is defined as where i S and T S are the area of i L and the contact network, respectively.
The variation curves of deviatoric stress, volumetric strain, and area percentage indexes (A3, A4, A5, and A6+) with axial strain are shown in Figure 15. The evolution of i A shows that only L6+ experiences an area increase while the other kinds of loop contract. In fact, having more sides gives rise to L6+ the ability to transform. The evolution of contact loops actually represents the evolution of assembly volume in the 2D simulation. The evolution trend of A6+ and v ε is similar, indicating that L6+ is the main factor affecting volume evolution.

Quantitative Analyses between Meso-Structure and Macro-Mechanics
Compared with the softening stage, the contact loops change more significantly in the hardening stage, which should be used as the focus of force chain evolution. The  Figure 15. The variation curves of (a) deviatoric stress, (b) volumetric strain, and area percentage indexes (A 3 , A 4 , A 5 , and A 6+ ) with axial strain.

Quantitative Analyses between Meso-Structure and Macro-Mechanics
Compared with the softening stage, the contact loops change more significantly in the hardening stage, which should be used as the focus of force chain evolution. The quantified analysis for the relationship between meso-structural and macro-mechanical indexes is established in the hardening stage (ε a ≤ 2.0%).
In this section, the macro-mechanical indexes (deviatoric stress, axial strain, and volumetric strain) are used as dependent variables and meso-structural indexes (number percentage indexes and area percentage indexes) are used as independent variables to establish multivariate models. The independent variables could be reduced in dimensionality by the principal component analysis, and obtain the principal components [34,35]. Further, the multivariate models of the meso-structural and macro-mechanical indexes can be obtained by establishing a multivariate regression equation between the principal components and the independent variables.

Principal Component Analysis of the Meso-Structural Indexes
As there are eight independent variables, multicollinearity may occur in this highdimension analysis and compromise the statistical significance of independent variables. Multicollinearity occurs when the absolute value of the Pearson correlation coefficient is higher than 0.7 [36,37]. Pearson correlation coefficient (R) is defined as where Y P is the predicted value, and Y A is the actual value. The statistical results of Pearson correlation coefficients among the eight independent variables are shown in Figure 16. where P Y is the predicted value, and A Y is the actual value. The statistical results of Pearson correlation coefficients among the eight independent variables are shown in Figure 16.  Figure 16 shows that the Pearson correlation coefficients of the diagonal can be higher than 0.7, indicating that multicollinearity can occur if all the variables are used. When multicollinearity occurs, the principal component analysis is suitable for the independent variables [38]. The principal component analysis is a multivariate statistical method that reduces multiple independent variables to a small number of principal components through dimensionality reduction techniques [39]. The principal components can reflect most information of the original variables and are linearly independent of each other. The eight meso-structural indexes in the hardening stage ( a 2.0% ε ≤ ) are shown in Table 2.   Figure 16 shows that the Pearson correlation coefficients of the diagonal can be higher than 0.7, indicating that multicollinearity can occur if all the variables are used. When multicollinearity occurs, the principal component analysis is suitable for the independent variables [38]. The principal component analysis is a multivariate statistical method that reduces multiple independent variables to a small number of principal components through dimensionality reduction techniques [39]. The principal components can reflect most information of the original variables and are linearly independent of each other. The eight meso-structural indexes in the hardening stage (ε a ≤ 2.0%) are shown in Table 2. The original data matrix X = n × p = 21 × 8 was established from the data in Table 2, where n and p represent the number of samples and variables, respectively.
According to the definition of the overall principal component, the covariance of the principal component cov(F) is a diagonal array, which is expressed as The principal components F 1 , F 2 , . . . , F p are uncorrelated with one another, which F 1 , F 2 , . . . , F p are called first, second, . . . , pth principal components, respectively. The percentage of the variance of the i principal component F i in the total variance f i / m ∑ j=1 f j (i = 1, 2, . . . , p) contribution rate is called the contribution rate of the principal component F i . The contribution rate of the principal component reflects the ability of the principal component to synthesize the original variable information, and can also be understood as the ability to interpret the original variable [40]. The sum f j of the contribution of the first m (m ≤ p) principal components is called the cumulative contribution rate of the first m principal components, which reflects the ability of the first m principal components to explain the information of the original variables [41]. X is subjected to principal component analysis. According to the total variance explained table (Table 3), the percentages of the variance of the first three principal components are all greater than 10%, and the cumulative contribution rate has reached 99.868%, so it is sufficient to extract the first three principal components. The extracted three principal components can remove implausible variables and determine the contribution of each variable to each principal component by using the component matrix. The component matrix is the coefficient of the factor expression of each variable, expressing the degree of influence of the extracted component on the mesostructural index. For the component matrix, the actual meaningful relationship between the components and the variables is not obvious. To make the coefficients more significant, the component matrix can be rotated so that the relationship between principal components and variables is redistributed and the correlation coefficients are differentiated towards 0 to 1. The relationship between principal components and meso-structural indexes can be derived from Table 3, and the rotated component matrix is shown in Table 4. Table 4. Rotated component matrix between the principal components and variables.

Variables
Principal Components −0.967 --By observing Table 4, it is found that each meso-structural index has a reasonable value of 1 (i.e., greater than 0.4), so none of the eight meso-structural indexes need to be deleted. The principal component F 1 , the highest percentage of contribution, is mainly influenced by ω 3 , ω 6+ , A 3 , A 4 , A 5 , and A 6+ indexes, which can reflect the effect of the area percentages of all loops and the number percentages of L 3 and L 6+ . The middle principal component F 2 and the third principal component F 3 are mainly influenced by ω 4 and ω 5 , respectively, which reflects the effect of the number percentages of L 4 and L 5 is weak for the macro-mechanical indexes. Based on the rotated component matrix and the standardized coefficients (4.2), we can build the relationship between meso-structural indexes, principal components, and macro-mechanical indexes are shown in Figure 17, which reflects the contribution of meso-structural indexes to principal components and the effect of principal components to macro-mechanical indexes. Additionally, the influence degree between meso-structural indexes and principal components is quantized, showing as the component score coefficient matrix in Table 5.
the macro-mechanical indexes. Based on the rotated component matrix and the standardized coefficients (4.2), we can build the relationship between meso-structural indexes, principal components, and macro-mechanical indexes are shown in Figure 17, which reflects the contribution of meso-structural indexes to principal components and the effect of principal components to macro-mechanical indexes. Additionally, the influence degree between meso-structural indexes and principal components is quantized, showing as the component score coefficient matrix in Table 5. Figure 17. The relationship between meso-structural indexes, principal components, and macromechanical indexes.   The component score matrix indicates the relationship between each meso-structural index and each component, with a high score on a component indicating the closer the relationship between that indicator and that component. Based on the component score coefficient matrix, the functions and values of the three principal components F 1 , F 2 , and F 3 can be obtained (Table 6) and used in place of the meso-structural indexes for the next step.

Establishment of Multivariate Model Based on Principal Components
The feedback of meso-structural indexes on macro-mechanics was achieved by establishing multivariate models of the three principal components F 1 , F 2 , and F 3 with axial strain ε a , volumetric strain ε v , and deviatoric stress q. Tolerance and variance inflation factor (VIF) was used to determine whether equations of the multivariate models were multicollinear, and the multivariate models were validated by variance analysis. The partial regression coefficients of the models were examined to determine the influence degree of the principal components on macro-mechanical indexes utilizing standardized coefficients [42]. The multivariate model between the axial strain ε a and the principal components F 1 , F 2 , and F 3 is shown as The variance analysis of the Equation (13) indicates an F-value of 89.912 with a p-value < 0.001, i.e., indicating that the multivariate model can be considered statistically significant at the α = 0.05 test level. Table 7 shows the results of the partial regression coefficient test. The p-values of all partial regression coefficients within the 95% confidence interval (95%CI) are less than 0.05, indicating that the significance levels of the partial regression coefficients are all in order. The standardized coefficients (β) for each principal component indicate that it can be seen that the principal component F 1 has the greatest effect on the axial strain ε a , and F 2 and F 3 have a small effect. The multivariate model between the volumetric strain ε v and the principal components F 1 , F 2 , and F 3 is shown as The variance analysis of the Equation (14) indicates an F-value of 30.83 with a p-value < 0.001. Table 8 shows the results of the partial regression coefficient test. The p-values of all partial regression coefficients within the 95% confidence interval (95%CI) is less than 0.05. The standardized coefficient (β) for each principal component shows that the principal component F 2 has the greatest effect on the volumetric strain ε v , with the second-highest influence degree of the principal component F 1 , and they are about three times the influence degree of the principal component F 3 . The multivariate model between the deviatoric stress q and principal components F 1 , F 2 , and F 3 is shown as The variance analysis of the Equation (15) indicates an F-value of 753.49 with a p-value < 0.001. Table 9 shows the results of the partial regression coefficient test. The p-values of all partial regression coefficients within the 95% confidence interval (95%CI) are less than 0.05. The standardized coefficient (β) for each principal component shows that principal component F 1 has the greatest influence on deviatoric stress q, being about five times more influential than principal component F 3 , and principal component F 2 is more than three times the degree of influence of F 3 . Table 9. Partial regression coefficient test results for Equation (15).

Principal
Components β 95%CI p-Value For Equations (13)-(15), the Tolerance = 1 > 0.2 and the VIF = 1 < 10, demonstrating that there was no multicollinearity between the independent variables and no dimensionality reduction was required. Additionally, the actual values and the calculated values of three macro-mechanical indexes are fitted in Figure 18. For the fitting lines, the closer the slope is to 1, and the closer the intercept is to 0, the better the fit is.  For Equations (13)-(15), the Tolerance = 1 > 0.2 and the VIF = 1 < 10, demonstrating that there was no multicollinearity between the independent variables and no dimensionality reduction was required. Additionally, the actual values and the calculated values of three macro-mechanical indexes are fitted in Figure 18. For the fitting lines, the closer the slope is to 1, and the closer the intercept is to 0, the better the fit is. The above analysis results show that the parameter estimation, hypothesis testing, and overall fit of the three multivariate models are good and statistically significant.

Conclusions and Outlook
A recognition and determination technique for 2D contact loops was proposed in this study. Taking the biaxial compression test as an example, the meso-structural indexes were calculated by the technique, and the relationships between the meso-structural indexes and the macro-mechanical indexes are analyzed. The main findings are summarized as follows: (1) Based on the numerous and changeable polygonal loops in contact network images, the proposed Q-Y algorithm is effective in determining the geometric types of contact loops in contact network images. The above analysis results show that the parameter estimation, hypothesis testing, and overall fit of the three multivariate models are good and statistically significant.

Conclusions and Outlook
A recognition and determination technique for 2D contact loops was proposed in this study. Taking the biaxial compression test as an example, the meso-structural indexes were calculated by the technique, and the relationships between the meso-structural indexes and the macro-mechanical indexes are analyzed. The main findings are summarized as follows: (1) Based on the numerous and changeable polygonal loops in contact network images, the proposed Q-Y algorithm is effective in determining the geometric types of contact loops in contact network images. (2) The change of contact loops is most active in the hardening stage, during which ω 3 and ω 6+ show opposite evolution patterns, while ω 4 and ω 5 are basically stable. The area evolution of contact loops represents the volume evolution of the 2D assembly and L 6+ is the main factor affecting volume evolution. (3) The variation of meso-structural indexes is active in the hardening stage, wherein the multivariate models between meso-structural indexes and macro-mechanical indexes were built.
The multivariate models in this study build a bridge between the mesoscale and macroscale of granular matters. Additionally, the contribution rate of meso-structural indexes to macroscopic mechanical indexes is quantified, which makes up for the deficiency of qualitative explanation only in existing studies. Although the multivariate models have a good verification effect, the multivariate models may have limitations due to the influence of some factors (such as stress condition and particle shape), since these factors may change the expression form of the multivariate models. Additionally, since the 2D recognition and determination technique is an image processing-based algorithm, it is difficult to extend the technique to 3D. Therefore, this study is limited to meso-structure and macro-mechanics of 2D DEM simulations.
Based on the techniques proposed in this study, it is suggested that the multivariate quantitative models can be further improved by changing the influences factors in the future in order to accurately feedback the macroscopic mechanical behavior of granular matter through the mesoscopic contact network. Of course, if the quantitative relationship between meso-structure and macro-mechanics can be verified experimentally, it would be of great significance for this area of study. Data Availability Statement: All the research data used in this manuscript will be available whenever requested.