Numerical Damping Calibration Study of Particle Element Method-Based Dynamic Relaxation Approach for Modeling Longwall Top-Coal Caving

: When using the explicit dynamic relaxation approach (DRA) to model the quasi-static rock breakage, fragmentation, and ﬂow problems, especially the top-coal caving question, introducing numerical damping into the solution equation is inevitable for reducing the vibration frequency and impact speed of mesh nodes, which is signiﬁcantly affect the ﬁdelity of the computation results. Although the DRA has been widely adopted to simulate top-coal caving, the reasonable value and calibration method of numerical damping are still open issues. In this study, the calibration process of reasonable numerical damping for modeling top-coal caving is investigated by comparing with the experimental results, in which several geometry parameters of the drawing funnel are selected as the calibration indexes. The result shows that the most reasonable numerical damping value is 0.07 for the numerical modeling of interval top-coal caving in extra-thick coal seams. Finally, the correlation between the numerical damping and the physical top-coal drawing process is discussed. The numerical damping indirectly reﬂects the fragmentation in multi scale of coal mass and friction interaction between coal particles during the caving process, which reduces the vibration intensity of the top-coal caving system and dissipates the kinetic energy.


Introduction
The longwall top-coal caving technology (LTCC) is the most popular mining approach in extra-thick coal seams in China, vigorously promoting the coal production of single working faces up to 10 million tons in many northwestern China coal mines [1][2][3]. Through nearly 40 years, the LTCC has replaced the conventional slicing mining method [4,5], as its great potential in improving the safety, effectiveness, and production of coal mines. According to statistics, the mining period of the fully mechanized top-coal caving face is shortened, and the recovery rate can reach 60~85% [6,7].
Nowadays, the LTCC in extra-thick coal seam has entered the stage of automation and intelligence. Therefore, the evolution characteristics of the coal-rock interface under different caving technologies need to be studied more accurately. The shape of coal rock interface after top-coal caving in the strike direction of working face was studied and verified by Wang et al. [8] through laboratory similar experiment. The coal and rock flow law, which is affected by coal drawing step, top coal lump size, coal drawing process and drawing opening size, was studied by Huang et al. [9] through using the experimental

Dynamic Relaxation Principle
The DRA was a new mechanical concept proposed by Rayleigh in the last century [19]. It was developed by Day [20,21] and successfully applied to solve plane rigid frame and plate structures in 1965. The DRA solves static problems by using the dynamic solution of the stressed system. In the DRA, it is assumed that the structure's mass was concentrated on the nodes, and each node adds virtually damping force and inertia force. Thus, the static equilibrium equation of each node was transformed into the dynamic equilibrium equation. Starting from the known initial state of the node, the linear differential equation with time as independent variable and velocity as dependent variable was integrated, which was called the dynamic relaxation process of the node. When the node reached the static state, the dynamic relaxation end of the node was defined. After the dynamic relaxation of a node was completed, it would inevitably lead to the imbalance of the related nodes to start the dynamic relaxation process until all nodes were statically balanced [22].
In DRA, the equation to be solved is a linear function about time. Using the initial state to iteratively calculate the new function value can also be applied to nonlinear problems. The algorithm of DRA is as follows: It is assumed that the system is discretized into n connected nodes, and the cndition that the whole system is static equilibrium is as follows: where {P(∆)} represents the inner force induced by displacements of nodes 1~N; F is the force on a node, given by: P(∆)} = {P 1 , P 2 , . . . , P N } T (2) where P N represents the inner force induced by displacements of the node N.
The three freedom degrees of each node is expressed by the following formula: Then the residual force of each node is equal to the difference value between the external load and the inner force, and its matrix is as follows: The residual force of each node can be decomposed by: For any node i, taking the x direction as an example, the residual force R xi is zero when x direction is in a static equilibrium state, and if it's unbalanced, the residual force R xi is not zero. Assuming that the inertia force along the x direction is M xi a xi and the damping force is C i V xi , the motion equation of the node along the x direction at t time can be expressed as follows: R t xi = M xi a t xi + C i V t xi (6) where a t xi and V t xi are the acceleration and velocity of node i in the x direction at time t respectively, and M xi is the mass of node i in the x direction, and C i is the numerical damping coefficient.
The difference equation of the formula is as follows: According to Equation (7), the velocity V xi t+∆t/2 can be calculated by: Through Equation (8), we can get the new coordinate value in x-direction of node i at t + ∆t: The velocity and displacement of node i in y or z direction are the same as that in the x direction.
From the node dynamic relaxation process, it can be seen that the model should be calculated from the initial position and initial velocity. Simultaneously, time goes from beginning to end until all nodes complete this process, and dynamic relaxation ends. In the process of operation, the selection of parameters M i , C i , ∆t have a direct impact on the operation results. The mass parameter M i is determined by the material properties, and ∆t is the time step which is generally very small. So the value of numerical damping C i has a significant influence on the application results of the DRA. Numerical damping dramatically influences the accuracy of static simulation results, and its value is generally related to the mechanical properties of material damage, simulation boundary conditions, contact, and other factors.

The Principle of Progressive Failure Process of CDEM
CDEM is a method of coupling the finite element and the discrete element. And It can carry out finite element calculation inside the block and discrete element calculation at the boundary of the block. CDEM can simulate the deformation and motion characteristics of materials in continuous and discontinuous states and simulate materials' progressive failure process from continuum to in-continuity.
As shown in Figure 1, a link bar model is established between the two particles, which can simulate the fracture and crushing process of top coal under the action of external force in the initial complete state [18,19]. The link bar model can be regarded as a rectangle. The short side of the rectangle is the smaller particles' diameter, and the long side is the sum of the radii of the two particles. It can be used to calculate the contact force or cohesion between two particle elements. Based on the connecting bar model, the surface contact relationship between two particle elements is established. The equivalent contact area A c is the projected area of smaller particles, and A c is mainly used to calculate the contact stiffness between particle elements [17]. The contact stiffness parameters can be obtained by Equation (10): where Ki and Ei are the contact stiffness and modulus tensors, respectively; the subscript i = 1 represents the normal component, and it represents tangential components if i = 2. R1 The contact stiffness parameters can be obtained by Equation (10): where K i and E i are the contact stiffness and modulus tensors, respectively; the subscript i = 1 represents the normal component, and it represents tangential components if i = 2. R 1 and R 2 are the radii of the two particle elements, respectively. The contact force and momentum between two contact particle elements are calculated by: where F i (t + ∆t) and F i (t) are respectively the contact force at the time of t + ∆t and t. M i (t + ∆t) and M i (t) are respectively the moment at the time of t + ∆t and t; J i and ∆θ i are respectively the moment of inertia and increment angle difference between two contact particles; ∆u i is the increment of contact displacement. According to Mohr-Coulomb criterion and maximum tensile stress criterion, the calculation formula of contact force is shown as below: Equation (13) is the criterion of contact between two particle elements. If any inequality in the formula is satisfied, the contact between particles will no longer transfer moment: where T, C and ϕ are tensile strength, cohesion, and internal friction angle, respectively; I is the moment of inertia; R ave = (R 1 + R 2 )/2. In the two-dimensional numerical model, the contact state between the rigid boundary and the particle element is determined by the relative position between the particle center and the rigid wall edge. When the distance between the center of the particle and the rigid wall is less than or equal to the radius of the particle, it is considered that the particle and the rigid boundary are in contact, given by: where d is the distance from the particle center to the rigid wall, V kl is the relative position vector from the particle center to the endpoint in which k and l are the two points of the vector, respectively (see Figure 2); n is the unit normal vector of the rigid wall, and R is the radius of the particle element. The method to determine the contact between particles and rigid wall is shown in Figure 2. Once the contact between the particle and the rigid wall is established, the normal and tangential contact behaviors are generated automatically.
where d is the distance from the particle center to the rigid wall, Vkl is the relative position vector from the particle center to the endpoint in which k and l are the two points of the vector, respectively (see Figure 2); n is the unit normal vector of the rigid wall, and R is the radius of the particle element. The method to determine the contact between particles and rigid wall is shown in Figure 2. Once the contact between the particle and the rigid wall is established, the normal and tangential contact behaviors are generated automatically.

Numerical Model
(1) Basic assumptions In the process of top-coal caving, the top coal and immediate roof experienced the process of 'damaged body → broken block → loose coal flow' under the action of overlying strata advanced supporting pressure and repeated support pressure. The numerical

Numerical Model
(1) Basic assumptions In the process of top-coal caving, the top coal and immediate roof experienced the process of 'damaged body → broken block → loose coal flow' under the action of overlying strata advanced supporting pressure and repeated support pressure. The numerical simulation assumed that the top coal was discharged through the coal caving opening as uniform particles. The particle size increased gradually from bottom to top. In coal particle migration, it was only affected by its gravity and mutual friction. If there was a little gangue in the top coal, it was simplified as coal particles.
(2) Engineering background and numerical model In this paper, the 8222 working face of Tashan coal mine is taken as the engineering background. The average thickness of coal seam is 14.0 m, the mining height is 4.0 m, the caving height is about 10.0 m, the average dip angle of coal seam is 2 degrees. The working face adopts the one-round interval mining technology with 'large, medium, small and micro' four-stage automatic top-coal caving method, as shown in Figure 3.

Simulation Scheme
The numbers of the top-coal caving openings of the numerical model are 1 divided into four levels in advance: large, medium, small and micro. There are sev The top-coal caving sequence is as follows: the coal in 'large' position is drawn first, the coal in 'medium' position second, the 'small' third, and the 'micro' last. In order to simplify the calculation and accurately reflect the coal rock shape after top-coal caving in the working face, the numerical model is set with an angle of 0 • and a length of 58 m. There are 25 top-coal caving openings in total, and the width of each coal caving opening is 1.75 m. The top-coal drawing openings are numbered as 1#, 2#, 24#, 25#, with 7.05 m and 7.2 m reserved on both sides. Rigid plates are used as displacement constraints on the left and right sides and the top of the model to limit the movement range of particles; velocity constraints are used on the bottom boundary surface. The velocity and displacement are zero in the initial state, as shown in Figure 4. In the model, the top coal is divided into coal 1 layer, coal 2 layer and coal 3 layer, and the particle size is 0.05 m, 0.08 m and 0.15 m, respectively. The physical and mechanical parameters of coal and rock are shown in Table 1.

Simulation Scheme
The numbers of the top-coal caving openings of the numerical model are 1#~2 divided into four levels in advance: large, medium, small and micro. There are seven tervals between two 'large' top-coal caving openings, seven intervals between 'mediu four intervals between 'small' and one interval between 'micro'. The support classifi tion is shown in Table 2. The sequence of the numerical model is as follows: large → m dium → small → micro. This paper studies the influence of different numerical dampi on the numerical simulation of top-coal caving and puts forward the damping va method, which is reasonably determined. Combined with experience, numerical dam ing is set as 0.01, 0.03, 0.07, 0.11 and 0.15, respectively, to simulate the coal-rock interfa rules of top coal caving.

Simulation Results
(1) Results of the level 1 caving under different damping values The top-coal drawing openings which are set as level 1 are 1#, 9#, 17#, 25#, and coal-rock boundaries are formed as a symmetrical 'V' type with the top-coal drawi

Simulation Scheme
The numbers of the top-coal caving openings of the numerical model are 1#~25#, divided into four levels in advance: large, medium, small and micro. There are seven intervals between two 'large' top-coal caving openings, seven intervals between 'medium', four intervals between 'small' and one interval between 'micro'. The support classification is shown in Table 2. The sequence of the numerical model is as follows: large → medium → small → micro. This paper studies the influence of different numerical damping on the numerical simulation of top-coal caving and puts forward the damping value method, which is reasonably determined. Combined with experience, numerical damping is set as 0.01, 0.03, 0.07, 0.11 and 0.15, respectively, to simulate the coal-rock interface rules of top coal caving.

Simulation Results
(1) Results of the level 1 caving under different damping values The top-coal drawing openings which are set as level 1 are 1#, 9#, 17#, 25#, and the coal-rock boundaries are formed as a symmetrical 'V' type with the top-coal drawing opening as the center. The 'V' type is divided into upper and lower parts. The upper coal-rock boundary is similar to a curve, and the lower part is identical to a straight line, as shown in Figure 5. The simulation results show that under the same conditions of other parameters, the larger the numerical damping value, the smaller the top-coal drawing funnel's angle, and the narrower the 'V' shape. opening as the center. The 'V' type is divided into upper and lower parts. The upper coal-rock boundary is similar to a curve, and the lower part is identical to a straight line, as shown in Figure 5. The simulation results show that under the same conditions of other parameters, the larger the numerical damping value, the smaller the top-coal drawing funnel's angle, and the narrower the 'V' shape.  Table 3.  Table 3.
(2) Results of the level 2 caving under different damping value After the hydraulic support coal caving of level 2 (5#, 13#, 21#), the coal-rock boundaries also form a 'V' shaped funnel, as shown in Figure 6. With the increase of numerical damping value, the width of 'V' type top-coal drawing funnel gradually decreases and the drawing funnel becomes smaller. When the numerical damping value is 0.01 and 0.03, the funnel formed by the secondary caving is irregular 'V' shape, and the angle is about 64 • and 56.7 • , respectively. When the numerical damping value is 0.01 and 0.03, the funnel is irregular 'V' shape, and the angle is about 64 • and 57 • respectively; When the numerical damping value is 0.07, a regular 'V' type caving funnel is formed, and the top angle is about 43.3 • . When the numerical damping value is 0.11, the drawing funnel is irregular 'V' shape, and the top angle of 'V' is about 30.7 • . When the numerical damping is 0.15, a 'V' shape is formed with an average angle of 10 • . The maximum width of the upper part of 'V' shape formed by different numerical damping is also different. From small to large, Energies 2021, 14, 2348 9 of 17 the average maximum width of 'V' shape is 6.88 m, 6.72 m, 7.15 m, 7.74 m and 7.88 m, respectively, as shown in Table 4. After the hydraulic support coal caving of level 2 (5#, 13#, 21#), the coal-rock boundaries also form a 'V' shaped funnel, as shown in Figure 6. With the increase of numerical damping value, the width of 'V' type top-coal drawing funnel gradually decreases and the drawing funnel becomes smaller. When the numerical damping value is 0.01 and 0.03, the funnel formed by the secondary caving is irregular 'V' shape, and the angle is about 64° and 56.7°, respectively. When the numerical damping value is 0.01 and 0.03, the funnel is irregular 'V' shape, and the angle is about 64° and 57° respectively; When the numerical damping value is 0.07, a regular 'V' type caving funnel is formed, and the top angle is about 43.3°. When the numerical damping value is 0.11, the drawing funnel is irregular 'V' shape, and the top angle of 'V' is about 30.7°. When the numerical damping is 0.15, a 'V' shape is formed with an average angle of 10°. The maximum width of the upper part of 'V' shape formed by different numerical damping is also different. From small to large, the average maximum width of 'V' shape is 6.88 m, 6.72 m, 7.15 m, 7.74 m and 7.88 m, respectively, as shown in Table 4.    (

3) Results of the level 3 and level 4 caving under different damping value
As shown in Figure 7, the morphology of coal-rock is further affected after the level 3 caving, and the top coal sinks gradually as a whole. Due to different numerical damping, coal-rock movement and boundary distribution are complex and irregular. With the increase of damping value, the more coal left, the higher the 'V' shaped funnel. Among all the horns formed by hydraulic support coal caving, the funnels' height created by 3#, 11# and 19# support top-coal caving is higher. When the numerical damping is 0.01 and 0.03, the funnel's maximum height is 4.41 m and 4.07 m, respectively. When the numerical damping is 0.07, the average height of the funnel is about 4.09 m. When the numerical damping is 0.11 and 0.15, the height of the funnel gradually increases to 5.06 m and 5.33 m. After the end of level 4 ( Figure 8) micro caving, the top coal is further released, and the remaining coal is less. The caving funnel is irregular, and the coal-rock boundary is also very rough.  After the end of level 4 ( Figure 8) micro caving, the top coal is further released, and the remaining coal is less. The caving funnel is irregular, and the coal-rock boundary is also very rough.
(e) Damping value = 0.15 After the end of level 4 ( Figure 8) micro caving, the top coal is further released, and the remaining coal is less. The caving funnel is irregular, and the coal-rock boundary is also very rough.

Verification by Laboratory Experiments
The purpose of the laboratory experiments is to compare with the numerical simulation results under different numerical damping values, and put forward the method by which the reasonable damping value can be determined. The engineering background of the laboratory mining experiment is the working face no. 8222 of Tashan coal mine, and the similarity ratio is 1:10.

Structure Design of Test Device
The top-coal caving test platform in the dip direction of the extra thick coal seam is shown in Figure 9. The main body is a cuboid frame with a length of 3600 mm, a width of 250 mm, and a height of 2000 mm. It is mainly composed of the test frame, the top-coal caving openings, the observable baffle, the tilt adjustment cylinder, and other supporting devices. (1) The oil cylinder can adjust the tilt angle of the test frame at the bottom of the frame. The range of tilt angle adjustment is 0~30 • . (2) The top-coal caving opening simulation device is composed of box body, tail beam simulation inclined plane, coal chute inclined plane, coal drawing inserting plate, glass plate groove, and coal chute. And the height of the box is 300 mm, the width is 175 mm, and the length is 250 mm. (3) The observable baffle, whose thickness is 10 mm, is used to fix the front and rear sides of the experimental frame. We can observe the experiment process through the glass plate. (4) Other supporting devices mainly include a three-layer vibrating machine, electronic platform scale, camera, etc. The three-layer vibrating machine has three layers of sieves with a diameter of 12, 9 and 6 mm from top to bottom, which can automatically divide the mixed bulk materials into different particle sizes. The electronic platform scale adopts K-FINE high-precision electronic platform scale with a weighing range of 200 g~300 kg and an error value of 50 g~200 g. The camera adopted was a D7000 high-resolution SLR camera (Nikon, Bangkok, Thailand) which is used for real-time video recording of observed changes in the whole experimental platform during coal caving. experimental frame. We can observe the experiment process through the glass plate. (4) Other supporting devices mainly include a three-layer vibrating machine, electronic platform scale, camera, etc. The three-layer vibrating machine has three layers of sieves with a diameter of 12, 9 and 6 mm from top to bottom, which can automatically divide the mixed bulk materials into different particle sizes. The electronic platform scale adopts K-FINE high-precision electronic platform scale with a weighing range of 200 g~300 kg and an error value of 50 g~200 g. The camera adopted was a D7000 high-resolution SLR camera (Nikon, Bangkok, Thailand) which is used for real-time video recording of observed changes in the whole experimental platform during coal caving.

Scheme of Laboratory Similarity Experiment
(1) Material laying scheme. In the actual fully mechanized top coal caving mining process, the coal above the hydraulic support has been broken under the overburden pressure and repeated support pressure. The coal block size gradually increases from bottom to top. The top coal is laid in three layers. The first layer is yellow terrazzo, whose particle size is 3-6 mm, and the height is 334 mm. The second layer is red terrazzo with a height of 333 mm, and its particle size is 6-9 mm. The third layer is black terrazzo with a

Scheme of Laboratory Similarity Experiment
(1) Material laying scheme. In the actual fully mechanized top coal caving mining process, the coal above the hydraulic support has been broken under the overburden pressure and repeated support pressure. The coal block size gradually increases from bottom to top. The top coal is laid in three layers. The first layer is yellow terrazzo, whose particle size is 3-6 mm, and the height is 334 mm. The second layer is red terrazzo with a height of 333 mm, and its particle size is 6-9 mm. The third layer is black terrazzo with a height of 333 mm, and its particle size is 9-12 mm. The immediate roof is replaced by white terrazzo which its particle size is 10-20 mm and its laying height is 400 mm. The laying scheme is shown in Figure 10. height of 333 mm, and its particle size is 9-12 mm. The immediate roof is replaced by white terrazzo which its particle size is 10-20 mm and its laying height is 400 mm. The laying scheme is shown in Figure 10. (2) Experiment scheme of coal drawing. The simulation device of coal drawing opening is numbered 1#~19#. The laboratory similarity experiment scheme is the same as the numerical simulation scheme, and the intervals of 'large, medium, small and micro' are the same. The sequence of coal drawing is shown in Table 5.  (2) Experiment scheme of coal drawing. The simulation device of coal drawing opening is numbered 1#~19#. The laboratory similarity experiment scheme is the same as the numerical simulation scheme, and the intervals of 'large, medium, small and micro' are the same. The sequence of coal drawing is shown in Table 5.

Experimental Process and Results
According to Table 5, the 19 simulation devices of the test bench complete a round of coal caving process by separating into four times. The evolution process of coal and rock in front of the test platform is shown in Figure 11. The detailed coal caving process and results are as follows: 6# and 8# of hydraulic support are opened in sequence as 'large' caving opening. The migration evolution process of the coal-rock boundary is shown in Figure 11a. The drawing funnel above the coal drawing opening of 6# is approximately 'V' type with the coal drawing opening as the center. The 'V' type can be subdivided into the lower part of 'V' type, similar to a straight line. The upper part of 'V' type, which is identical to the curve. The angle of the lower part of 'V' type is about 35 • , and its height is 858 mm. The width of the uppermost part of the 'V' shape is 890 mm, and the height of the upper part of 'V' type is 142 mm. The funnel of 14# is the same as that of 6#, and it is approximately 'V' type. The angle of lower part of 'V' type is about 34 • , and its height is 848 mm. The width of the uppermost part of the 'V' shape is 871 mm, and the height of upper part of 'V' type is 152 mm. The top coal of 2#, 10# and 18# as 'medium' caving opening are drawing in turn. The migration evolution process of coal-rock boundary is shown in Figure 11b, and the coal-rock boundary is approximate 'V' type. The 'V' type formed by 2# coal caving has an angle of about 36°and a maximum width of about 685 mm in the uppermost part. The 'V' shape formed above the 10# coal drawing opening has an angle of about 37° and a length of 737 mm in the uppermost part. The top angle of the 'V' type formed above the 18# coal drawing opening is about 36° and the uppermost length is 694 mm.
The top coal of 4#, 8#, 12#, and 16# as 'small' caving are drawing in sequence as shown in Figure 6c. And the drawing funnel formed by drawing coal is irregular 'V' shape. The remaining coal drawing openings are 'micro' drawing openings. After top-coal drawing, the amount of coal left is less, and the drawing funnel is exceptionally irregular.
(a) The "large" caving (b) The "medium" caving (c) The "small" caving and "micro" caving Figure 11. Evolution process of coal-rock morphology in front of test platform.

Calibration Parameters
To accurately describe and compare the coal-rock migration rules after coal drawing Figure 11. Evolution process of coal-rock morphology in front of test platform. The top coal of 2#, 10# and 18# as 'medium' caving opening are drawing in turn. The migration evolution process of coal-rock boundary is shown in Figure 11b, and the coal-rock boundary is approximate 'V' type. The 'V' type formed by 2# coal caving has an angle of about 36 • and a maximum width of about 685 mm in the uppermost part. The 'V' shape formed above the 10# coal drawing opening has an angle of about 37 • and a length of 737 mm in the uppermost part. The top angle of the 'V' type formed above the 18# coal drawing opening is about 36 • and the uppermost length is 694 mm.
The top coal of 4#, 8#, 12#, and 16# as 'small' caving are drawing in sequence as shown in Figure 6c. And the drawing funnel formed by drawing coal is irregular 'V' shape. The remaining coal drawing openings are 'micro' drawing openings. After top-coal drawing, the amount of coal left is less, and the drawing funnel is exceptionally irregular.

Calibration Parameters
To accurately describe and compare the coal-rock migration rules after coal drawing in laboratory experiments and numerical simulation experiments, it is necessary to choose the drawing funnel's characteristic parameters. After the 'large' and 'medium' caving, the coal-rock boundary forms obvious 'V' type.
After the 'large' caving, the 'V' type can be divided into the lower part of the coal-rock boundary, which is similar to a straight line, and the upper part is similar to a curve. The parameter θ 1 can be expressed as the angle of 'V' type, and B 1 can be delivered as the maximum width of the upper part of 'V' type. And h 1 can be formulated as the height of the upper part of 'V' type. Meanwhile, The parameter θ 2 can be expressed as the angle of 'V' type, and B 2 can be formulated as the maximum width of 'V' type after the 'medium' caving. Moreover, it is difficult to extract the characteristic parameters because of the small amount of coal and the irregular caving funnel after the 'small' and 'micro' caving.
After the 'large' caving, the similarity degree of the caving funnel formed by the laboratory similarity experiment and the numerical simulation experiment can be expressed by X 1,α : where θ 1,0 is the average angle of 'V' type of simulation experiment after the 'large' caving; After the 'medium' caving, the similarity degree of the caving funnel formed by the laboratory similarity experiment t and the numerical simulation experiment can be expressed by Y 2,α : where θ 2,0 is the average angle of 'V' type of simulation experiment after the 'medium' caving; B 2,0 is the average width of 'V' type in simulation experiment after the 'medium' caving; θ 2,α is the average angle of 'V' type in numerical simulation under different numerical damping α (α = 0.01, 0.03, 0.07, 0.11, 0.15) after the 'medium' caving; B 2,α is the average maximum width of 'V' type in numerical simulation under different numerical damping α (α = 0.01, 0.03, 0.07, 0.11, 0.15) after the 'medium' caving; The sum of Equations (15) and (16) gives the comprehensive index S α , which can be expressed as follows: S α = {X 1,α + Y 2,α }/2 (17) Finally, the reasonable value of numerical damping which can be defined as α can be determined by: Figures 12 and 13 show the characteristics of the coal drawing funnel after the 'large' and 'medium' caving in numerical simulation and similar simulation. According to the numerical simulation results, the 'V' type of coal drawing funnel is different with different numerical dampings, such as the angle and the maximum influence width of 'V' type. Furthermore, it is verified that the numerical damping value has a great influence on the accuracy of top-coal caving modeling. It can be concluded that the 'V' shape formed by the laboratory similarity experiment is close when the numerical damping value is 0.07 in the numerical simulation. Through Equation (18), the similarity degree parameter S α of numerical simulation and laboratory similarity experiment can be obtained as follows:    The results show that the value of S0.07 is the minimum, and the corresponding value of α' is 0.07. In other words, the deviation between numerical simulation and similar simulation is the smallest when the numerical damping value is 0.07. And the similarity is the closest between numerical simulation and laboratory similarity experiment.

Conclusions
Although the DRA has been widely adopted to simulate top-coal caving, the reasonable value and calibration method of numerical damping are still open issues. In this study, the calibration process of reasonable numerical damping for modeling top-coal caving is investigated by comparing with the experimental results, in which several geometry parameters of the drawing funnel are selected as the calibration indexes. Finally, the correlation between the numerical damping and the physical top-coal drawing process is discussed. The numerical damping indirectly reflects the fragmentation in multi scale of coal mass and friction interaction between coal particles during the caving process, which reduces the vibration intensity of the top-coal caving system and dissipates the kinetic energy. The specific conclusions can be drawn as follows.
(1) By studying the results of different damping values in CDEM numerical simulation of top-coal caving, it is verified that the numerical damping value greatly influences the accuracy of top-coal caving modeling. (2) In the comparison of numerical simulation and similar simulation results, it is proposed that the angle, the maximum width, and the upper height of 'V' shape of coal drawing funnel are the indexes to determine the numerical damping in 'large' drawing. The angle and the maximum width of 'V' shape are the indexes to determine the numerical damping in a 'medium' drawing. By comparing the deviation of characteristic parameters between numerical simulation and similar simulation, it is concluded that 0.07 is the reasonable numerical damping value for CDEM numerical simulation of 'large, medium, small and micro' top-coal caving in Tashan coal mine.   The results show that the value of S 0.07 is the minimum, and the corresponding value of α' is 0.07. In other words, the deviation between numerical simulation and similar simulation is the smallest when the numerical damping value is 0.07. And the similarity is the closest between numerical simulation and laboratory similarity experiment.

Conclusions
Although the DRA has been widely adopted to simulate top-coal caving, the reasonable value and calibration method of numerical damping are still open issues. In this study, the calibration process of reasonable numerical damping for modeling top-coal caving is investigated by comparing with the experimental results, in which several geometry parameters of the drawing funnel are selected as the calibration indexes. Finally, the correlation between the numerical damping and the physical top-coal drawing process is discussed. The numerical damping indirectly reflects the fragmentation in multi scale of coal mass and friction interaction between coal particles during the caving process, which reduces the vibration intensity of the top-coal caving system and dissipates the kinetic energy. The specific conclusions can be drawn as follows.
(1) By studying the results of different damping values in CDEM numerical simulation of top-coal caving, it is verified that the numerical damping value greatly influences the accuracy of top-coal caving modeling. (2) In the comparison of numerical simulation and similar simulation results, it is proposed that the angle, the maximum width, and the upper height of 'V' shape of coal drawing funnel are the indexes to determine the numerical damping in 'large' drawing. The angle and the maximum width of 'V' shape are the indexes to determine the numerical damping in a 'medium' drawing. By comparing the deviation of characteristic parameters between numerical simulation and similar simulation, it is concluded that 0.07 is the reasonable numerical damping value for CDEM numerical simulation of 'large, medium, small and micro' top-coal caving in Tashan coal mine.  Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.