Quantitative Study on the Contact Force and Force Chain Characteristics of Ore Particle Systems during Ore Drawing from a Single Drawpoint under the Influence of a Flexible Barrier

It is of great significance to carry out quantitative research on the contact force chain characteristics of ore particle systems during ore drawing to reveal the microscopic and mesoscopic characteristics of ore particle systems during implementation of the synchronous filling shrinkage stoping method. Based on the particle discrete element method, combined with the relevant knowledge of contact mechanics and statistical mechanics, microscopic properties of the ore particle system were studied quantitatively. (1) The probability distributions of the normal and tangential contact forces during ore drawing from a single drawpoint under a flexible barrier are similar, and both show exponential decay. (2) In the regions on both sides of the model, the ore particles will not be released because they are far from the drawpoint; there, the coordination number is stable, and the density of contact network is large. In the upper part of the drawpoint, ore particles flow, the coordination number fluctuates violently, and the density of contact network is small. (3) In the initial stage, the stress distribution of the ore particle system is uniform, and the strong force chain does not show obvious directivity. After that, the directions of the strong force chains in the ore particle system are inclined to the vertical direction, and the strong force chains mainly bear the load in the vertical direction. (4) In the initial stage of drawing, the normal contact force is mainly concentrated in the vertical direction. With the progression of drawing, the normal contact force at an angle of 30° to the horizontal direction increases gradually, and the number of the main direction of the average contact force distribution changes from one to three (the vertical direction and the direction with a ±30° angle to the horizontal direction).


Introduction
In 2010, the authors proposed the idea of "synchronous filling" and invented a representative mining method, namely, "the mass draw and synchronous filling with no-top-pillar shrinkage stoping method" (synchronous filling shrinkage stoping method) [1]. This method is based on the traditional shrinkage method; before drawing the ore, workers can stand on the ore pile (6 in Figure 1) to lay the flexible barrier. In the process of drawing, the filling materials are sent into the stope room through the return airway (1 in Figure 1) in a timely fashion, and by means of vibration and gravity, the filling materials and ore particles sink synchronously and uniformly [2]. The ore particles are drawn under the influence of a flexible barrier, which does not come into direct contact with the waste rocks. Accordingly, there is no ore loss or dilution in the ore particles drawn. At the same time, surrounding rock pressure can be controlled and surface subsidence can be reduced. However, the flow behaviour of the ore particles in the new mining method is considerably improved over that using traditional draw theories. In view of this problem, the authors carried out a series of studies [3][4][5][6][7]. During the research process, it was found that the caved ore in a stope is a typical soft condensed matter system (i.e., granular matter system) composed of the combination of solid, liquid, and gas particles of different shapes. In the granular matter system, particles squeeze and contact each other under the action of an external load and their own gravity, forming a chain path transferring the external load, which is called the force chain [8]. As a bridge connecting the macroscopic and microscopic properties of granular materials, the study of force chains is the foundation of the multiscale mechanical research framework of granular matter [9]. Therefore, to further reveal the internal action mechanism of the flow process of ore particles under the influence of a flexible barrier, it is of great significance to study the contact forces and the force chain formed by them.
Research on the contact forces and resulting force chain in a granular matter system has been conducted for a long time. Initially, researchers mainly studied the micro-and mesoscopic characteristics of the granular matter system through physical experiments. For example, Oda [10] observed the distribution of the particle contact angle in the process of triaxial compression via thin section of the sample and found that the contact angles of the particles in the granular matter system follow the direction of the major stress. Matsuoka [11] studied the evolution law of the particle contact angle during the direct shear process of aluminium rods through a photoelastic test to derive the stress-dilatancy relation equation. Since 1979, with the introduction of the discrete element method (DEM) [12], the DEM has gradually become a major tool for studying the microscale rules of granular matter systems, and the DEM has previously been demonstrated to be able to reproduce certain key features of granular materials [13][14][15][16]. Matuttis et al. [17] studied the characteristics of the arch and force chain inside an accumulation body formed by spherical particles and polygonal particles by using DEM and found that the distribution characteristics of the arch and force chain inside the accumulation body formed by the two types of particles were not very different, but the accumulation history process, such as point-source-type pouring and the one-by-one accumulation of monolayer particles, has a great influence on the accumulation result and force chain distribution. Sun et al. [18][19][20] described the research framework of multiscale mechanics of granular matter and used the DEM to generate a particle system to simulate the evolution law of the force chain and the force transfer characteristics of particles during uniaxial compression and biaxial compression tests.
The above studies have been carried out on the microand mesoscopic characteristics of the granular matter system under shear, compression, and static accumulation; however, these studies are all in a preliminary stage of exploration, and there is no unified standard for quantitative research on contact force and force chain characteristics; meanwhile, none of them have been involved in the field of ore drawing. In the field of ore drawing, the three classical ore drawing theories [21] were put forward without defining the force chain, let alone the related problems of the representation of the force chain. Therefore, on the basis of previous studies and with the background of the synchronous filling retention method, this paper quantified the contact force and force chain characteristics in the ore particle system under the influence of a flexible barrier to provide a new approach to the study of metal ore drawing.

Model Construction and Statistical
Parameter Description According to the dimensions of the structure outline of the physical experimental model in the synchronous filling mining method [1][2][3][4], the dimensions of the draw model in this 2D numerical simulation were 168 × 128 cm, and the spaces between drawpoints were 24 cm, as shown in Figure 2. The contact model of this ore drawing test is the rolling resistance linear model, and the particle generation is performed by the gravity accumulation method. To realize the synchronous process of filling and ore drawing under the flexible isolation layer, particle generation is carried out in the following three steps: (1) Generation of initial particles: through the ball generate command, several ore particles were generated in the range of 8 cm to 130 cm along the y-axis of the wall model, the gravity acceleration of these ore top-pillar, (3) raise, (4) cross heading, (5) barrier pillar, (6) remaining ore, (7) drift pillar, (8) drawpoint, (9) haulage drive, (10) unmined ore, (11) filling material, and (12) barrier. 2 Geofluids particles is g = −9:81 m · s −2 , and the microscopic mechanical parameters are shown in Table 1. To compact the ore particles in the granular matter system as soon as possible, the contact model of the initial particles was set as a linear contact model, and the friction coefficient between the particles was set to 0.3; at the same time, to facilitate the observation of the ore particle flow phenomenon during ore drawing, after the model was balanced, the particles were assigned different colours at 10 cm intervals, and ore particles above the 128 cm level were removed (2) Generation of the flexible barrier: to simulate the flexible barrier, the Cubic command was used to generate a 250 cm row of fine particles with a radius of 0.0015 m above the resulting ore particles. These particles were bonded in parallel, and the barrier was given the microscopic mechanical parameters shown in Table 2 (3) Generation of the calculating particles: to simulate the influence of ore particles with different sizes and shapes, the contact model of particles was changed from the linear contact model to the rolling resistance linear model. The resulting micromechanical calculation parameters of particles in the granular matter system are shown in Table 3. After deleting the bottom wall at drawpoint no. 4, the drawpoint was opened, ore particles were released from the drawpoint, and the flow of the ore particles started immediately. Using Fish of PFC 2D to compile loop statements in the process of ore drawing while computing several time steps, drawpoint no. 4 was closed, and an appropriate amount of filling waste rock particles were generated on the surface of the ore particles (the micromechanical calculation parameters of the waste rock particles are the same as those of the ore particles). To achieve the synchronous filling effect, after the model is bal-anced under the action of dead weight, the excess filling waste rock particles were deleted, and the no. 4 ore drawpoint was opened again to start the calculation of the next cycle; when no more ore particles are released, the ore drawing is stopped 2.2. Description of Statistical Parameters. The main statistical parameters describing the characteristics of the contact forces and force chain at the microscale are as follows: the contact force probability distribution density, average particle coordination number, force chain component ratio, and contact angle. These parameters can be used for quantitative characterization of the contact force magnitude, distribution concentration, force chain shape, and direction.

Probability Distribution Density of the Contact Force.
The contact forces of granular media show certain distribution characteristics under shear, static stacking, uniaxial compression, and other conditions, and many researchers have explored this in depth [22,23]. Based on the abovementioned research, this paper intends to study the intensity characteristics of the contact force between particles in the flow process of granular matter under a flexible barrier through the distribution probability density of the contact force. Taking the calculation of the normalized normal contact force as an example, the calculation process of the probability density of the contact force distribution is as follows [24]: where superscript n is the normal contact force, subscript i represents the contact number, f n i is the normalized magnitude of the normal contact force of contact i, F n i is the magnitude of the normal contact force of contact i, and N is the total number of similar contacts. Similarly, the normalized tangential contact force and total contact force of contact i are f t i and f i , respectively. Equation (2) is used to normalize the normal contact forces of all the contacts in the granular matter system. The normal contact force is divided into 10 [3.5, 4], [4, 4.5], and [4.5, 5]. The total number of normal contact forces after normalization in each interval was calculated,   3 Geofluids and the probability density of normal contact force in each interval can be obtained: where subscript k is the kth interval, and N k is the total number of normal contact forces in the kth interval. Similarly, the distribution probability densities of the tangential contact force and total contact force in each interval can be obtained as f s k and f k , respectively. 2.2.2. Coordination Number. The coordination number Z, also known as the average number of contacts, is an important index that can be used to evaluate whether the contacts of the granular matter system are good and compact. Generally, the larger the coordination number, the denser the particle extrusion is, the more stable the system is, the less likely it is to break and recombine, and the denser the distribution of the contact network is. Therefore, the coordination number is an important indicator that reflects the density of the contact network. The coordination number is defined as [25]: where N c is the number of contacts in the granular matter system (the contacts with a normal contact force is greater than 0), and N p is the number of particles in the granular matter system. The coordination number can be understood as the average number of particles in contact with each particle, which to some extent reflects the state of the granular matter in the measured area (in terms of factors such as the stress level and contact force distribution intensity). Furthermore, Thornton [26] proposed the concept of mechanical coordination number, pointing out that if there is only one or no particles in contact with the particles around a particle, then the contact does not contribute to the bearing capacity of the whole granular matter system; therefore, theoretically, these particles should be removed. Therefore, the coordination number is an important parameter to evaluate the state of the granular matter system, and this paper intends to use it to characterize the distribution density of contact networks.

Force Chain Component
Ratio. The force chain component ratio quantitatively represents the overall extension direction of the force chain morphology through the relationship between the contact force components in the x-axis direction and the y-axis direction [27]. The calculation method is as follows: where f x i represents the magnitude of the component of the contact force of contact i on the x-axis, f y i represents the magnitude of the component of the contact force of contact i on the y-axis, n represents the total number of contacts, μ c represents the force chain component ratio (when μ c > 0, the force chain is inclined to the vertical direction, that is, the y-axis direction; when μ c < 0, the force chain is inclined to the horizontal direction, that is, the x-axis direction), and the greater the difference between μ c and 0 is, the greater the deflection of the force chain. In this paper, the force chain component ratio is compared between the total force chain and the strong force chain, in which the strong chain refers to the force chain formed by the contact force greater than the average contact force, and the total force chain refers to the force chain formed by all the contacts in the granular system.

Contact
Angle-Average Contact Force. The directivity and local mechanical properties of the force chain can be reflected by the statistical distribution of the particle contact angles and average contact force. Dividing 360°into k intervals (in this paper, k is 36), by counting the number of contacts and contact forces within each interval, the average contact forces within each interval are calculated. The calculation formula of the average contact force is as follows: where θ c is the included angle between the normal and horizontal direction of the contact plane, F k i is the magnitude of the contact force of contact i in the kth interval, N k c is the total number of contacts in the kth interval, and f ðθ c Þ is the average contact force (the larger the value of f ðθ c Þ is, the larger the average value of the contact force in the direction of θ c , and the greater the impact of the contact force on the stability of the system in the direction of θ c ).

Macroscopic Flow Characteristics of Ore Particles
The whole ore drawing process is completed in 15 iterations. The macroscopic flow phenomenon of ore particles is shown in Figure 3, and the corresponding contact force distribution information is shown in Figure 4 (due to space constraints, only a few representative drawing snapshots were selected). It can be seen from Figures 3 and 4 that after opening the drawpoint, with the continuous release of ore particles, the barrier gradually drops, and the flow range of particles expands upward in a funnel shape from the drawpoint; the labelled particles in each layer descend successively from the bottom to the barrier, and the shape of the labelled particles presents a Gaussian distribution curve on the whole; after the isolation layer sinks, its morphology is similar to that of the labelled particles, as it also presents a Gaussian distribution, but the morphology at the bottom of the barrier was relatively smooth, and a cavity appeared near the bottom. With the increase in ore particles drawn, the volume of the cavity becomes increasingly obvious.
According to the traditional ore drawing theory, the labelled particles at each layer sink in a funnel shape during ore drawing from a single drawpoint. When the labelled particles in the funnel reach the drawpoint, the labelled particles on the plane begin to form a broken funnel. After the labelled 4 Geofluids particles at the highest level form a broken funnel, all the labelled particle layers form a broken funnel. If the ore particles continue to be drawn, the waste rock will mix in. However, if a barrier is set on the top layer of the ore surface, the waste rock will be quarantined, and ore drawing can continue. To explore the flow characteristics of ore particles under the influence of a flexible barrier, the authors' team conducted a detailed study on the draw column morphology of ore particles, the relationship between the height of the draw column and the accumulated mass drawn, the interface movement law of the barrier, and the evolution law of the cavity under the condition of mass draw and synchronous filling [3][4][5][6][7]. To further explore the internal action mechanism of granular matter in the process of ore drawing under the influence of a flexible barrier from the microscale, the following part describes a statistical analysis of the evolution laws of the magnitudes and distribution density of the contact forces and the shape and direction of the force chain.

Distribution Law of the Contact Force Magnitude.
According to the macroscopic flow characteristics of ore particles in the draw process with a single drawpoint under the influence of a flexible barrier, four representative ore drawing iterations (the 1st drawing, the 5th drawing, the 10th drawing, and the 15th drawing) were selected, and the contact force data of each node model were derived to obtain the probability distribution of the normal contact forces, tangen-tial contact forces, and total contact forces between the ore particles in each node, as shown in Figure 5.
To describe the distribution characteristics of the contact forces quantitatively, the fitting function was used. As shown in Figure 5, the changes in the normal contact force, tangential contact force, and total contact force are similar, and they all decay exponentially. Therefore, the exponential decay function can be used to fit the normal contact force, tangential contact force, and total contact force. The exponential decay function y = A 1 exp ð−x/t 1 Þ + A 2 exp ð−x/t 1 Þ + y 0 has the best fitting effect on the above points, and the fitting coefficient R 2 is above 0.99. The fitting curves of each point under different iterations of ore drawing are shown in Figure 5.
According to the distributions of the normal contact forces, tangential contact forces, and total contact forces as well as the corresponding probability density function, during the ore drawing process with a single drawpoint under the influence of a flexible barrier, the distributions of contact forces of different types are similar and show exponential attenuation, while the distributions of contact force intensity at different iterations of ore drawing are not significantly different. The distribution of the contact forces shows exponential attenuation, which indicates that there are few strong force chains in the granular matter system and that most of the force chains are weak. There are fewer strong force chains and weaker force chains in the ore particle system, and the strong and weak force chains are interwoven to form a complete force chain network that maintains the stability of the whole system. There is no significant difference in the distribution laws of contact force intensity under different drawing

Geofluids
iterations, which indicates that the ore particle system in the whole drawing process is controlled by the continuous fracturing and reorganization of the internal force chain; there is no great change on the whole, and the interaction with the surrounding rock does not change the stress state considerably from the original stress state.

Density Distribution Law of the Contact Force Network in Different Positions.
To describe the variation in the density of the contact force network at different positions during ore drawing, the ore particle system was evenly divided into 9 parts, and a measuring circle with a radius of 0.2 m was arranged at the centre of each part (as shown in Figure 6). The distribution density of the contact force in the ore particle system during the ore drawing process with a single drawpoint under the influence of a flexible barrier was quantitatively characterized by using the measuring circle to monitor the change law of the coordination number in each part of the ore particle system. The whole ore drawing process was carried out in 15 iterations, and the change law of the coordination number in the ore particle system in the first ore drawing process is shown in Figure 7.
It can be seen from Figure 7 that in the process of ore drawing with a single drawpoint under the influence of a flexible barrier, the ore particles above the ore drawing port (i.e., the region including measurement circles 2, 5, and 8) are greatly affected by the drawpoint; therefore, the force chain of the ore particles in this region is constantly broken and recombined, and the coordination number fluctuates violently. The ore particles in the rest of the model are far from the drawpoint and are less affected by the drawpoint; therefore, their ore particle coordination numbers are stable near a certain fixed value without much fluctuation. According to this phenomenon, the ore particle system is divided into three regions in the horizontal direction, and different heights of the measurement circles are selected, as shown in Figure 8. The whole ore drawing process is monitored in the measuring circles, and the monitoring results are shown in Figure 9.
As shown in Figure 9, the coordination number in measuring circle 1 (coordination number 1) fluctuates at approximately 3.5 throughout the ore drawing process, with a small fluctuation range; coordination number 3 fluctuates approximately 3.0, also with a small fluctuation range. The coordination number in measuring circle 2 (coordination number 3) fluctuated sharply at approximately 1.5 in the early stage of ore drawing. With the continuous release of ore particles, coordination number 2 gradually increased in the late stage of ore drawing, and the fluctuation range gradually decreased and finally stabilized at approximately 3.0. This is mainly because the measure circle corresponding to the coordination number was eventually filled by the waste rock particles.
It can be seen from the above phenomenon that the ore particles on both sides of the model during the drawing  7 Geofluids process with a single drawpoint are far from the drawpoint; therefore, they do not participate in the release of ore particles. The ore particles in these areas are stable; thus, the coordination numbers are large, and the densities of the corresponding contact networks are also high. Moreover, the further down the ore particles move, the closer their contact becomes, and the particle coordination number and contact network density increase. Affected by the flow of ore particles, the internal force chain of ore particles in the upper part of the drawpoint is constantly broken and recombined, and the coordination number fluctuates violently. With the continuous release of ore particles, more filling waste rocks are added, the ore drawing area decreases, and fewer particles flow. Finally, the ore particle system above the drawpoint gradually becomes stable.

Analysis of Force Chain Morphology and
Contact Directivity 5.1. Morphology Evolution Law of the Force Chain. The forms of the force chain network of the ore particle system under different iterations of ore drawing are shown in Figure 10 (due to space limitations, only a few representative ore drawing snapshots were selected). The thickness of the lines in the figure indicates the strength of the force chain: the thicker a line is, the stronger the force chain, and the thinner a line is, the weaker the force chain. Figure 10 shows that in the initial iteration, the number of strong force chains was small, mainly distributed around the drawpoint on both sides of the model, the whole force chain network was relatively uniform, and the network gap was small; with the continuous release of ore particles, the number of strong force chains gradually increased, and their positions gradually spread upward along the sidewall from both sides of the drawpoint, and the force chain network showed obvious directionality, mainly along the direction of the vertical barrier.
The x-axis and y-axis components of the force chain strength in the ore particle system in the global coordinate system after the end of each ore drawing were counted, and the morphological evolution characteristics of the force chain network were further quantified. The variation law of the force chain component ratio for the different iterations of ore drawing is shown in Figure 11. Figure 11 shows that the total force chain and strong force chain component ratios in the ore particle system during the ore drawing process are greater than zero, and the strong force chain component ratio is approximately 2% larger than the total force chain component ratio (except for the first ore drawing). This indicates that the direction of the force chain network in the ore particle system is always inclined to the y direction, and the strong force chains mainly bear the load in the vertical direction. At the end of the first ore drawing, the total force chain component ratio is greater than the strong force chain component ratio, which shows that the stress distribution in the ore particle system is uniform in the initial stage, and the strong chain does not show obvious directivity. The change rules of the component ratios of the strong force chain and the total force chain are consistent. The strong force chain and total force chain component ratios fluctuate with the increase in ore drawing iterations: the fluctuation range in the early and middle stages is small, while the fluctuation range in the later stage is large, which shows that the force chain network in the ore particle system is constantly  8 Geofluids breaking and restructuring in the process of ore drawing. Due to the influence of the shape of the barrier, the direction of the overlying load acting on the ore particle system changed greatly in the late period of ore drawing; therefore, the fluctuation range of the force chain component ratio also increased significantly.

Analysis of the Change in the Contact Force Angle
Distribution. The whole ore drawing process is completed in 15 iterations, and the contact coordinates and size of the ore particles at the end of each drawing are recorded and output. After dividing 360°into 36 intervals, counting the number of contacts and contact forces within each interval, and then calculating the average contact forces within each interval, the distribution of the contact angle and normal contact force of particles in the ore drawing process is statistically obtained, as shown in Figure 12.
To quantitatively describe the directional distribution characteristics of the average contact force inside the ore particle system during synchronous filling and ore drawing, based on the previous description function of the distribution of interparticle contact force [28,29], Equation (6) was used to compare the statistical results of the distributions of interparticle contact force direction: where f n ðθÞ is the distribution function of the normal contact force between particles, f 0 is the average value of all the contact forces in the ore particle system, a n is the Fourier coefficient, whose value represents the anisotropy of the contact force in the ore particle system, θ n is the main direction of the distribution of the force chain in the ore particle   9 Geofluids system, w is the periodic control coefficient of the trigonometric function, and the minimum positive period of the trigonometric function is T = 2π/ðπ/wÞ = 2w. Equation (6) is used to compare the statistical results of the distributions of contact force and contact angle in the ore particle system under different ore drawing times. The blue line in Figure 12 is the fitting line, and the fitting parameters are shown in Table 4.
It can be seen from Figure 12 and Table 4 that the normal contact force shows obvious anisotropy at the initial stage of ore drawing, that is, the average contact force along the y-axis is relatively large, while the average contact force along the x-axis is relatively small, and the normal contact force is mainly concentrated in the vertical direction; with the progress of ore drawing, the normal contact force at an angle of ±30°from the horizontal direction gradually increased, the number of the main distribution direction of the average contact force changed from one to three, and the anisotropy degree increases gradually. According to the fitting results, it can be known that from the first ore drawing to the end of the fifth ore drawing, the fitting figure presented a peanut-shaped distribution, the anisotropy coefficient remained at 0.34, and the main angle changed from 94.03°to 89.93°, that is, the distribution of the contact force was slightly deflected to the vertical direction; from the end of the 5th ore drawing to the end of the 15th ore drawing, the fitting figure gradually changed from a peanut shape to a petal shape, the anisotropic coefficient experienced a variation of 0.34~0.37~0.53, and the main angle experienced a variation of 89.93°~12.64°~14.65°, which shows that with the progress of ore drawing, the anisotropy of the contact force distribution in the ore particle system gradually increases. The above-mentioned analysis shows that, during the initial drawing, the normal contact   force mainly bears the load in the vertical direction. Then, as the ore particles continue to release, the barrier gradually bends and sinks, the load in the direction of the vertical barrier gradually increases, and the normal contact force is initially distributed along the direction of the vertical barrier, that is, the contact force strength at an angle of ±30°from the horizontal direction gradually increases.

Conclusions
(1) During the ore drawing process with a single drawpoint, the probability distributions of the normal and tangential contact forces are similar and show exponential attenuation, while the distribution of contact force intensities under different iterations of ore drawing is not significantly different (2) The ore particles on both sides of the model during the drawing process with a single drawpoint are far from the drawpoint; therefore, they do not participate in the release of ore particles. The ore particles in these areas are stable, the coordination number is large, and the density of the corresponding contact network is also large. Moreover, the further down the ore particles are, the more closely they come into contact, increasing the particle coordination number   and contact network density. The internal force chains of ore particles in the upper part of the drawpoint are constantly broken and recombined, and the coordination number fluctuates violently (3) The total force chain and strong force chain component ratios in the ore particle system during the ore drawing process are greater than zero, and the strong force chain component ratio is approximately 2% larger than the total force chain component ratio (except for the first ore drawing). This indicates that the direction of the force chain network in the ore particle system is always inclined to the y direction, and the strong force chains mainly bear the load in the vertical direction. At the end of the first ore drawing, the total force chain component ratio is greater than the strong force chain component ratio, which shows that the stress distribution in the ore particle system is uniform in the initial stage, and the strong chain does not show obvious directivity. The change rules of the component ratios of the strong force chain and the total force chain are consistent (4) In the early stage of ore drawing, the normal contact force is mainly concentrated in the vertical direction; with the increase in ore drawing iterations, the normal contact force at an angle of ±30°f rom the horizontal direction gradually increased, the main distribution direction of the average contact force changed from one to three (the horizontal direction and the direction at an angle of ±30°from the horizontal direction), and the anisotropy degree increased gradually

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.