Numerical Investigation of Top Coal Drawing Evolution in Longwall Top Coal Caving by the Coupled Finite Difference Method-Discrete Element Method

: In longwall top coal caving (LTCC), the resource recovery ratio of the working face is directly determined by the top coal recovery ratio. An investigation of the evolution of top coal drawing characteristics and revealing the evolution of top coal drawing parameters is necessary when providing guidance for caving parameter selection and improving the top coal recovery ratio. Based on in-situ measurements of the size distribution of caved top coal blocks in Wangjialing coal mine, a ﬁnite difference method (FDM)–discrete element method (DEM) coupled method was applied to establish a “continuous–discontinuous” numerical model and the process from the ﬁrst coal drawing to the common coal drawing was simulated with 17 separate working face advances. The evolution of the drawing body (DB), loose body (LB), and top coal boundary (TCB) was obtained. The results show that, the evolution of parameters of DB such as shape and size, drawing amount, length and deﬂection angle of the long axis of the proﬁle ellipsoid tended to decrease ﬁrst, then increase, decrease again, and ﬁnally stabilise; the increment of the LB advance coal wall distance and the coal pillar distance was close to 0 m in the common coal drawing stage, while width increment of the LB was close to the drawing interval (0.865 m). The TCB formed after each coal drawing round was ﬁtted based on the improved “Hook” function. The evolution of height and radius of curvature of TCB’s stagnation point was analysed. This was divided into three stages: the ﬁrst (ﬁrst to third drawing rounds) was the initial mining inﬂuence stage, the second (fourth to ninth drawing rounds) was the transitional caving stage, and the third (after tenth drawing round) was the common coal drawing stage. coal blocks may have a certain impact on the overall movement and drawing effect of top coal, its impact is limited; so we ignore the role of too small coal blocks in the overall drawing process of top coal. Three hundred and ﬁfteen, 240 and 208 sets of data were obtained at the middle (75# caving support), and two ends (16# and 134# caving supports) of the working face. The measured results are shown in Figure 2. within the quadrant the coordinate system ( 0). reduce errors caused by deﬁnition-domain limitation during the ﬁtting process, absolute values were taken the X -coordinate of each and a curve symmetrical to the Y -axis with the TCB was obtained. Fitting was performed on this curve. The ﬁtting formula was determined based on the evolution characteristics of the TCB. The curve had an asymptote and was asymmetrical about its own axes, is the Hook function (of the form f ( x ) = x + 1/ x ). method of determination


Introduction
As one of the main methods for underground thick seam mining, longwall top coal caving (LTCC) has been fully developed and applied in China [1]. The coal seam is divided into two parts in LTCC: bottom coal and top coal. Bottom coal is cut directly by the shearer, while top coal is finally caved from the drawing opening under its self-weight due to the fracture gradually developing under the influence of mining pressure [2,3]. Based on the drawing characteristics of top coal during the caving process of the LTCC, the caving process can be divided into first coal drawing and common coal drawing [4]. First coal drawing refers to the first round in the top coal caving process of the LTCC working face. The most remarkable characteristic of the common coal drawing stage is that the top coal drawing characteristics show a relatively consistent periodicity [5,6]. As most operations of the LTCC mining face fall within the common coal drawing stage, most scholars focus the numerical method. Castro et al. [23,24] used the physical simulation method to carry out the rock gravity caving test under the condition that the boundary stress and particle sizes differed; their test results showed that the application of vertical pressure and change of particle size had a significant impact on particle caving. By using the PFC software, Yang et al. [25] studied the top coal drawing mechanism under the influence of upward angle, focusing on the drawing amount, DB form and TCB. Based on the analytic hierarchy process (AHP) and fuzzy discrimination method, Chi et al. [26] put forward the prediction model for top coal caving characteristics and drawing features and verified them in Tashan Coal Mine and Tongxin Coal Mine. Wang et al. [13] established the prediction model of top coal recovery ratio by using physical simulation and theoretical analysis methods and verified its reliability through in situ measuring the block size of top coal on four LTCC faces and its recovery ratio. Zhang et al. [27] put forward a numerical method for automatic control of coal drawing based on time criterion by using the continuous-discontinuous numerical method, conducted a simulation study on the top coal drawing mechanism in extremely thick coal seam, and optimized the window switching principle for coal drawing and the coal drawing technology centering on the gangue mixture ratio, top coal recovery ratio and drawing amount.
The top coal drawing differs greatly in the first coal drawing stage and the common coal drawing stage, and a transition can be observed between the two stages. Most scholars focus their research on the common coal drawing stage and ignore this transition. Based on the in-situ measurement of the equivalent diameter of caved top coal blocks in the 12309 LTCC working face of the Wangjialing coal mine, a coupled finite difference method (FDM)discrete element method (DEM) was applied to established a "continuous-discontinuous" numerical model and the process from the first coal drawing to the common coal drawing was simulated with 17 separate working face advances. Consider that one of the most important stages of safe operation in a longwall excavation is the selection of the appropriate powered roof support [28,29]. The influence of hydraulic support on the drawing process is considered in the numerical simulation. The evolution of the DB, LB, and TCB was deduced, thus revealing the evolution of the top coal drawing and its characteristics during the process from first coal drawing to common coal drawing.

General Situation
The average burial depth of panel 12309 in Wangjialing coal mine was 360 m while the average thickness was 6.5 m. The cutting and caving height were 3.0 m and 3.5 m, respectively. A schematic representation of the coal seam and roofs is shown in Table 1. The hydraulic support and scraper conveyors are illustrated in Figure 1.

Top Coal Block Size Distribution
In the LTCC operation, top coal changes to a granular state under th sure, and the top coal migration is affected by the distribution of the top certain extent. As it is difficult to predict the distribution of granular top we measured the top coal lumpiness in the 12309 working face as the basis f the FDM-DEM numerical model. The measurement and analysis process h lished elsewhere [31]. In the process of field test, although some rocks w out and transported to the measurement location, we only measured the instead of rocks. Certainly, as you said, the measurement error is inevita process, we evenly spread the coal blocks drawn out by a single support on the stage loader and measure the coal blocks one by one. However, as so are too small (or even powder-like) and in huge quantity, we did not meas may lead to the decrease in the accuracy of the measurement results. We hough the smaller coal blocks may have a certain impact on the overall m drawing effect of top coal, its impact is limited; so we ignore the role of blocks in the overall drawing process of top coal. Three hundred and fiftee sets of data were obtained at the middle (75# caving support), and two end caving supports) of the working face. The measured results are shown in F  Figure 2 shows that the top coal block equivalent diameters were div sectors and the mass fraction of each sector was obtained in Figure 2: 0-9 c 18 cm (46.31%), 18

Top Coal Block Size Distribution
In the LTCC operation, top coal changes to a granular state under the mining pressure, and the top coal migration is affected by the distribution of the top coal block to a certain extent. As it is difficult to predict the distribution of granular top coal block size, we measured the top coal lumpiness in the 12309 working face as the basis for establishing the FDM-DEM numerical model. The measurement and analysis process have been published elsewhere [31]. In the process of field test, although some rocks were also drawn out and transported to the measurement location, we only measured the weight of coal instead of rocks. Certainly, as you said, the measurement error is inevitable. In the test process, we evenly spread the coal blocks drawn out by a single support on the scraper of the stage loader and measure the coal blocks one by one. However, as some coal blocks are too small (or even powder-like) and in huge quantity, we did not measure them. This may lead to the decrease in the accuracy of the measurement results. We think that, although the smaller coal blocks may have a certain impact on the overall movement and drawing effect of top coal, its impact is limited; so we ignore the role of too small coal blocks in the overall drawing process of top coal. Three hundred and fifteen, 240 and 208 sets of data were obtained at the middle (75# caving support), and two ends (16# and 134# caving supports) of the working face. The measured results are shown in Figure 2.

Top Coal Block Size Distribution
In the LTCC operation, top coal changes to a granular state under the mining pressure, and the top coal migration is affected by the distribution of the top coal block to a certain extent. As it is difficult to predict the distribution of granular top coal block size, we measured the top coal lumpiness in the 12309 working face as the basis for establishing the FDM-DEM numerical model. The measurement and analysis process have been published elsewhere [31]. In the process of field test, although some rocks were also drawn out and transported to the measurement location, we only measured the weight of coal instead of rocks. Certainly, as you said, the measurement error is inevitable. In the test process, we evenly spread the coal blocks drawn out by a single support on the scraper of the stage loader and measure the coal blocks one by one. However, as some coal blocks are too small (or even powder-like) and in huge quantity, we did not measure them. This may lead to the decrease in the accuracy of the measurement results. We think that, although the smaller coal blocks may have a certain impact on the overall movement and drawing effect of top coal, its impact is limited; so we ignore the role of too small coal blocks in the overall drawing process of top coal. Three hundred and fifteen, 240 and 208 sets of data were obtained at the middle (75# caving support), and two ends (16# and 134# caving supports) of the working face. The measured results are shown in Figure 2.  Figure 2 shows that the top coal block equivalent diameters were divided into four sectors and the mass fraction of each sector was obtained in Figure 2: 0-9 cm (13.83%), 9-18 cm (46.31%), 18-27 cm (20.34%), and 27-36 cm (19.52%).

Numerical Model
The FDM-DEM model includes coal seam, immediate roof, main roof, top coal caving hydraulic support and rear scraper conveyor. In the modeling process, the PFC 3D V6.0 (ITASCA Inc., Minneapolis, MN, USA) was used as the main body to establish the discontinuous rock stratum (top coal and immediate roof) with Balls, the Walls was used to establish the equipment model (top coal caving hydraulic support and rear scraper conveyor), and the FLAC 3D V6.0 (ITASCA Inc.) module was introduced to establish the con-   Figure 2 shows that the top coal block equivalent diameters were divided into four sectors and the mass fraction of each sector was obtained in Figure 2: 0-9 cm (13.83%), 9-18 cm (46.31%), 18-27 cm (20.34%), and 27-36 cm (19.52%).

Numerical Model
The FDM-DEM model includes coal seam, immediate roof, main roof, top coal caving hydraulic support and rear scraper conveyor. In the modeling process, the PFC 3D V6.0 (ITASCA Inc., Minneapolis, MN, USA) was used as the main body to establish the discontinuous rock stratum (top coal and immediate roof) with Balls, the Walls was used to establish the equipment model (top coal caving hydraulic support and rear scraper conveyor), Energies 2021, 14, 219 5 of 18 and the FLAC 3D V6.0 (ITASCA Inc.) module was introduced to establish the continuous rock stratum (bottom coal and main roof) with Zones. The particle size distribution of top coal complied with the test results as shown in Figure 2. The mechanical properties of the immediate roof rock are much stronger than the coal seam, so the crushing lumpiness of the immediate roof rock is bound to be greater than that of the top coal. In view of the unmeasurable of the crushing lumpiness of the immediate roof rock, we assumed that the size of its lumpiness was four times [8,12,31] that of the top coal, and the mass fraction was the same as that of the top coal as shown in Table 2. The numerical model is shown in Figure 3. lumpiness of the immediate roof rock is bound to be greater than that of the top coal. In view of the unmeasurable of the crushing lumpiness of the immediate roof rock, we assumed that the size of its lumpiness was four times [8,12,31] that of the top coal, and the mass fraction was the same as that of the top coal as shown in Table 2. The numerical model is shown in Figure 3.   Figure 3 shows that overall size of the numerical model has a length of 45 m and the bottom coal is divided into 3 sections in the lengthwise direction. The width of the numerical model is 1.75 m, which is equal to the width of top coal caving hydraulic support. The height of the numerical model is 16.1 m, which is the same as the thickness of the actual coal seam, the immediate roof and the main roof stratum. As a result, to inverse the pressure exerted by the 343.9 m (average depth of 360 m) rock stratum above the main roof on the main roof itself, a vertical compressive stress of 7.57 MPa is applied to the upper surface of the model. As this paper focuses on the evolution law of the loose top coal during its movement, and the top coal above the hydraulic support can be considered to be in a loose state, the contact constitutive model between Balls in the numerical simulation adopts a linear model, and the constitutive model of the continuous unit uses an elastomeric model. Numerical parameters are listed in Table 3.  The height of the numerical model is 16.1 m, which is the same as the thickness of the actual coal seam, the immediate roof and the main roof stratum. As a result, to inverse the pressure exerted by the 343.9 m (average depth of 360 m) rock stratum above the main roof on the main roof itself, a vertical compressive stress of 7.57 MPa is applied to the upper surface of the model. As this paper focuses on the evolution law of the loose top coal during its movement, and the top coal above the hydraulic support can be considered to be in a loose state, the contact constitutive model between Balls in the numerical simulation adopts a linear model, and the constitutive model of the continuous unit uses an elastomeric model. Numerical parameters are listed in Table 3.

Simulation Process
The process of numerical simulation was divided into three stages. The first stage was the mining preparation stage, which the original state and equipment set-up process of the working face were inverted to obtain each formation state before working face advance. The second stage was the drawing preparation stage, the process of only cutting bottom coal without caving top coal after working face advance was inverted. The calculation was stopped before proceeding to the next stage after the range of LB expanded to the immediate roof. The range and shape of LB before caving were derived as the primary boundary condition for caving operation. The third stage was the top coal drawing stage, the caving and recovery process of top coal was inverted to reveal the evolution of the DB, LB, and TCB. The numerical simulation process is demonstrated in Figure 4.

Simulation Process
The process of numerical simulation was divided into three stages. The first stage was the mining preparation stage, which the original state and equipment set-up process of the working face were inverted to obtain each formation state before working face ad vance. The second stage was the drawing preparation stage, the process of only cutting bottom coal without caving top coal after working face advance was inverted. The calcu lation was stopped before proceeding to the next stage after the range of LB expanded to the immediate roof. The range and shape of LB before caving were derived as the primary boundary condition for caving operation. The third stage was the top coal drawing stage the caving and recovery process of top coal was inverted to reveal the evolution of the DB LB, and TCB. The numerical simulation process is demonstrated in Figure 4. The "Fishcall" function [32,33] was inserted at the "−11" point of the calculation cycle of PFC 3D to realize real-time monitoring and control of the model. In order to automati cally delete the drawn particles and close the coal drawing opening (automatically stop the calculation), the "Fishcall" function is introduced, which is called before the start o each calculation step to identify whether it is necessary to delete the particles or stop the calculation. The advantage of doing this is that there is no need to manually control the termination of the calculation program; namely, when certain conditions are met, the pro gram will automatically stop the calculation. In this paper, when the specific gravity o the drawn rock particles reaches 5% of that of the top coal particles (or any other propor tion), the calculation will automatically stop. The biggest advantage of this setting is tha the calculation process can be controlled according to the calculation termination condi tions defined by users, avoiding the errors caused by manually. Whether the proportion of rock particles exceeds 5%?

Drawing preparation stage
After the 12th coal drawing round The "Fishcall" function [32,33] was inserted at the "−11" point of the calculation cycle of PFC 3D to realize real-time monitoring and control of the model. In order to automatically delete the drawn particles and close the coal drawing opening (automatically stop the calculation), the "Fishcall" function is introduced, which is called before the start of each calculation step to identify whether it is necessary to delete the particles or stop the calculation. The advantage of doing this is that there is no need to manually control the termination of the calculation program; namely, when certain conditions are met, the program will automatically stop the calculation. In this paper, when the specific gravity of the drawn rock particles reaches 5% of that of the top coal particles (or any other proportion), the calculation will automatically stop. The biggest advantage of this setting is that the calculation process can be controlled according to the calculation termination conditions defined by users, avoiding the errors caused by manually. The main functions include: (1) deleting particles dropping from the drawing opening to the rear scraper conveyor in real-time to invert the rear scraper conveyor hauling the caved top coal out; (2) monitoring the proportion of immediate roof particles in the caved particles in real-time and triggering the calculation termination condition when the proportion of immediate roof particles reaches 5% of the total number of particles caved in each drawing round.

Mining Preparation and Drawing Preparation Stages
After equipment installation, the top coal above the hydraulic support dropped and reached equilibrium ( Figure 5a); only then can the working face begin to cut the bottom coal and advance. After five steps involving this coal cutting and support advance, the range of the LB reached the immediate roof, which triggered the termination of the calculation. The particle displacement diagram and movement boundary of LB after each round of support advance were extracted, as shown in Figure 5.
OR PEER REVIEW 7 of 18 clude: (1) deleting particles dropping from the drawing opening to the rear scraper conveyor in real-time to invert the rear scraper conveyor hauling the caved top coal out; (2) monitoring the proportion of immediate roof particles in the caved particles in real-time and triggering the calculation termination condition when the proportion of immediate roof particles reaches 5% of the total number of particles caved in each drawing round.

Mining Preparation and Drawing Preparation Stages
After equipment installation, the top coal above the hydraulic support dropped and reached equilibrium ( Figure 5a); only then can the working face begin to cut the bottom coal and advance. After five steps involving this coal cutting and support advance, the range of the LB reached the immediate roof, which triggered the termination of the calculation. The particle displacement diagram and movement boundary of LB after each round of support advance were extracted, as shown in Figure 5. (a-f) respectively refers to the range of LB formed from when the support stays in place to after it is removed for the 5th time. Figure 5 shows that the range of LB expanded with the advance of the working face. Before the support was advanced, the maximum boundary height was 1.27 m (Figure 6a). When the support begins to advance, the maximum boundary height kept increasing albeit at a gradually decreasing rate. After the fourth support advance, the range of LB reached the immediate roof. In actual operation, as the mining pressure was not significant at this stage, top coal and the immediate roof broke to a relatively small extent, which did not meet the requirements for caving, therefore, only bottom coal cutting was conducted without top coal drawing in the caving preparation stage. As top coal and immediate roof crushing block size of these stages could not be measured in-situ, the particle diameter distribution was configured based on the measurement result in Figure 2. Whether the range of LB reached the immediate roof was the basis for judgment as to start of drawing. The reason for this was that when the LB reached the immediate roof, under the influence of working face advance, the top coal particles near the drawing opening in the vertical direction had all started to move, and it can be considered that the top coal particles were relieved from being squeezed, therefore, the top coal drawing was started after the fifth support advance. (a-f) respectively refers to the range of LB formed from when the support stays in place to after it is removed for the 5th time. Figure 5 shows that the range of LB expanded with the advance of the working face. Before the support was advanced, the maximum boundary height was 1.27 m (Figure 6a). When the support begins to advance, the maximum boundary height kept increasing albeit at a gradually decreasing rate. After the fourth support advance, the range of LB reached the immediate roof. In actual operation, as the mining pressure was not significant at this stage, top coal and the immediate roof broke to a relatively small extent, which did not meet the requirements for caving, therefore, only bottom coal cutting was conducted without top coal drawing in the caving preparation stage. As top coal and immediate roof crushing block size of these stages could not be measured in-situ, the particle diameter distribution was configured based on the measurement result in Figure 2. Whether the range of LB reached the immediate roof was the basis for judgment as to start of drawing. The reason for this was that when the LB reached the immediate roof, under the influence of working face advance, the top coal particles near the drawing opening in the vertical direction had all started to move, and it can be considered that the top coal particles were relieved from being squeezed, therefore, the top coal drawing was started after the fifth support advance.

Evolution of DB
The DB refers to all top coal blocks caved in each drawing round. The shape of DB refers to the shape formed by these blocks before drawing. The 3D images and profile images of the DB shapes of all drawing rounds were extracted, as illustrated in Figure 6, which shows that all DBs fell in an ellipsoid cut by the tail beam of support (the area marked by the green dotted line in Figure 6a). The DBs of the first coal drawing and the common coal drawing were different. For the range of DB (Figure 6a), it extended to above the top beam of support, completely covering the shield beam and the tail beam after the first coal drawing. None of the DBs formed after common coal drawing rounds extended to the top beam, and except that of the first coal drawing, the DB of the fourth round extended the farthest (to most parts of the shield beam), the second to the fourth DB ranges extended to local parts of the shield beam, and the fifth to the twelfth DBs only covered the tail beam. In terms of drawing amount (ADB), the first round was the greatest, reaching 47.9 t, the second and third rounds decreased to about 7.5 t, the fourth round involved 15.9 t, and the subsequent rounds gradually decreased and stabilised at around 6 t each (Figure 6b). The evolution of the horizontal width (WDB) and vertical height (HDB) of DB was the same as that of ADB (Figure 6c), first decreasing, then increasing, decreasing again, and finally stabilising. WDB (HDB) decreased from 8.16 m (6.66 m) in the first round to around 3.72 m (3.81 m) in the second and third rounds, then increased to 5.01 m (5.05 m), and then gradually decreased before stabilising at around 2.7 m (4.0 m).
From these results, it can be seen that, from the first coal drawing round to the common coal drawing stage, each parameter of the DB changed according to a certain trend. The fourth drawing round represented a turning point, with its ADB, WDB, and HDB being slightly increased compared with those of the second and third rounds; after the fourth drawing round, the parameters gradually decreased and tended to stabilise. To verify the conclusion, the long-axis length (LDB) and the included angle (θDB) of the DB profile ellipsoid of each drawing round were extracted, as shown in Figure 7.

Evolution of DB
The DB refers to all top coal blocks caved in each drawing round. The shape of DB refers to the shape formed by these blocks before drawing. The 3D images and profile images of the DB shapes of all drawing rounds were extracted, as illustrated in Figure 6, which shows that all DBs fell in an ellipsoid cut by the tail beam of support (the area marked by the green dotted line in Figure 6a). The DBs of the first coal drawing and the common coal drawing were different. For the range of DB (Figure 6a), it extended to above the top beam of support, completely covering the shield beam and the tail beam after the first coal drawing. None of the DBs formed after common coal drawing rounds extended to the top beam, and except that of the first coal drawing, the DB of the fourth round extended the farthest (to most parts of the shield beam), the second to the fourth DB ranges extended to local parts of the shield beam, and the fifth to the twelfth DBs only covered the tail beam. In terms of drawing amount (A DB ), the first round was the greatest, reaching 47.9 t, the second and third rounds decreased to about 7.5 t, the fourth round involved 15.9 t, and the subsequent rounds gradually decreased and stabilised at around 6 t each (Figure 6b). The evolution of the horizontal width (W DB ) and vertical height (H DB ) of DB was the same as that of A DB (Figure 6c), first decreasing, then increasing, decreasing again, and finally stabilising. W DB (H DB ) decreased from 8.16 m (6.66 m) in the first round to around 3.72 m (3.81 m) in the second and third rounds, then increased to 5.01 m (5.05 m), and then gradually decreased before stabilising at around 2.7 m (4.0 m).
From these results, it can be seen that, from the first coal drawing round to the common coal drawing stage, each parameter of the DB changed according to a certain trend. The fourth drawing round represented a turning point, with its A DB , W DB , and H DB being slightly increased compared with those of the second and third rounds; after the fourth drawing round, the parameters gradually decreased and tended to stabilise. To verify As shown in Figures 6a and 7, the evolution of LDB and θDB was consistent with that mentioned above. For the long-axis extension of DB, the intersection point between it and coal floor was located in the goaf after the first four rounds of coal drawing, but it was located at the back of the rear scraper conveyor from the fifth coal drawing round (black dotted line in Figure 7a). For LDB, it reached 7.97 m after the first drawing round and decreased to around 4.3 m after the second and third drawing rounds, increased slightly after the fourth drawing round, and then gradually decreased before stabilising at around 4 m. For θDB, it maintained a large deflection (over 45°) during the first four drawing rounds and decreased to around 10° then stabilised from the fifth drawing round. Therefore, a similar conclusion was drawn, suggesting that the fourth drawing round is the key turning point for DBs in the present numerical simulation.

Evolution of LB
The LB refers to the coal and rock loose range caused by top coal caving. Particle displacement diagrams after each drawing round were extracted, as shown in Figure 8. The shape of LB is marked on the basis of the displacement exceeding 25 mm.     As shown in Figures 6a and 7, the evolution of L DB and θ DB was consistent with that mentioned above. For the long-axis extension of DB, the intersection point between it and coal floor was located in the goaf after the first four rounds of coal drawing, but it was located at the back of the rear scraper conveyor from the fifth coal drawing round (black dotted line in Figure 7a). For L DB , it reached 7.97 m after the first drawing round and decreased to around 4.3 m after the second and third drawing rounds, increased slightly after the fourth drawing round, and then gradually decreased before stabilising at around 4 m. For θ DB, it maintained a large deflection (over 45 • ) during the first four drawing rounds and decreased to around 10 • then stabilised from the fifth drawing round. Therefore, a similar conclusion was drawn, suggesting that the fourth drawing round is the key turning point for DBs in the present numerical simulation.

Evolution of LB
The LB refers to the coal and rock loose range caused by top coal caving. Particle displacement diagrams after each drawing round were extracted, as shown in Figure 8. The shape of LB is marked on the basis of the displacement exceeding 25 mm.  As shown in Figures 6a and 7, the evolution of LDB and θDB was consistent with that mentioned above. For the long-axis extension of DB, the intersection point between it and coal floor was located in the goaf after the first four rounds of coal drawing, but it was located at the back of the rear scraper conveyor from the fifth coal drawing round (black dotted line in Figure 7a). For LDB, it reached 7.97 m after the first drawing round and decreased to around 4.3 m after the second and third drawing rounds, increased slightly after the fourth drawing round, and then gradually decreased before stabilising at around 4 m. For θDB, it maintained a large deflection (over 45°) during the first four drawing rounds and decreased to around 10° then stabilised from the fifth drawing round. Therefore, a similar conclusion was drawn, suggesting that the fourth drawing round is the key turning point for DBs in the present numerical simulation.

Evolution of LB
The LB refers to the coal and rock loose range caused by top coal caving. Particle displacement diagrams after each drawing round were extracted, as shown in Figure 8. The shape of LB is marked on the basis of the displacement exceeding 25 mm.    As shown in Figure 8, the shape of LB expanded as the working face advanced, from 16.35 m after the first drawing round to 38.92 m after the twelfth drawing round. Throughout the process, the LB expanded towards the working face side and goaf side to different degrees as the working face advanced. With the coal wall and the coal pillar as the benchmark, the horizontal distances between the edge of LB and coal wall (S w ) and coal pillar (S p ) were extracted, respectively. The curve and the total width of LB (W LB ) were plotted, as shown in Figure 9.
Energies 2021, 14, x FOR PEER REVIEW 10 of 18 As shown in Figure 8, the shape of LB expanded as the working face advanced, from 16.35 m after the first drawing round to 38.92 m after the twelfth drawing round. Throughout the process, the LB expanded towards the working face side and goaf side to different degrees as the working face advanced. With the coal wall and the coal pillar as the benchmark, the horizontal distances between the edge of LB and coal wall (Sw) and coal pillar (Sp) were extracted, respectively. The curve and the total width of LB (WLB) were plotted, as shown in Figure 9.  Figure 9 shows that Sw and Sp first increased, then stabilised as the working face advanced. For Sw, the LB did not reach the coal wall after the first and second drawing rounds, so the distances were negative (−1.41 m and −0.43 m). The LB shape exceeded the constraint of the coal wall from the third drawing round onwards and Sw gradually increased from 0.43 m to 7.93 m after the tenth drawing round and stabilised thereafter (at around 7.9 m). For Sp, it increased from 3.52 m after the first drawing round to 6.86 m after the seventh drawing round, then remained stable (at around 7.0 m). For WLB, it increased from 16.35 m to 38.92 m, but its rate of increase decreased as the working face was advanced. The increment of Sw, Sp, and WLB (ΔSw, ΔSp, and ΔWLB) between drawing rounds was calculated (Figure 10, where the abscissa "2-1" represents the change to the second drawing round from the first).  Figure 9 shows that S w and S p first increased, then stabilised as the working face advanced. For S w , the LB did not reach the coal wall after the first and second drawing rounds, so the distances were negative (−1.41 m and −0.43 m). The LB shape exceeded the constraint of the coal wall from the third drawing round onwards and S w gradually increased from 0.43 m to 7.93 m after the tenth drawing round and stabilised thereafter (at around 7.9 m). For S p , it increased from 3.52 m after the first drawing round to 6.86 m after the seventh drawing round, then remained stable (at around 7.0 m). For W LB , it increased from 16.35 m to 38.92 m, but its rate of increase decreased as the working face was advanced. The increment of S w , S p , and W LB (∆S w , ∆S p , and ∆W LB ) between drawing rounds was calculated ( Figure 10, where the abscissa "2-1" represents the change to the second drawing round from the first).
Energies 2021, 14, x FOR PEER REVIEW 10 of 18 As shown in Figure 8, the shape of LB expanded as the working face advanced, from 16.35 m after the first drawing round to 38.92 m after the twelfth drawing round. Throughout the process, the LB expanded towards the working face side and goaf side to different degrees as the working face advanced. With the coal wall and the coal pillar as the benchmark, the horizontal distances between the edge of LB and coal wall (Sw) and coal pillar (Sp) were extracted, respectively. The curve and the total width of LB (WLB) were plotted, as shown in Figure 9.   to around 0 m; ∆W LB first increased from 2.42 m to 4.01 m, then decreased, and tended to stabilise after the tenth drawing round (around 0.865 m). These simulated results show that ∆S w , ∆S p , and ∆W LB reached their maxima at the fourth drawing round, which agreed with the aforementioned conclusion pertaining to the evolution of DB. ∆S w and ∆S p stabilised during the common drawing stage, then remained unchanged (at around 0), while ∆W LB stabilised during cyclic drawing, then approached 0.865 m, which approximated to the working face advance interval and drawing interval.

Evolution of TCB
TCB refers to the boundary between the top coal and the immediate roof rock. The TCB formed by coal caving was analysed in this section. The TCB after each caving round was extracted, as shown in Figure 11.
Energies 2021, 14, x FOR PEER REVIEW 11 of 18 As shown in Figure 10, during the caving process, ΔSp increased from 0.98 m to 2.23 m, then decreased to around 0 m; ΔSw first increased from 0.57 m to 0.92 m, then decreased to around 0 m; ΔWLB first increased from 2.42 m to 4.01 m, then decreased, and tended to stabilise after the tenth drawing round (around 0.865 m). These simulated results show that ΔSw, ΔSp, and ΔWLB reached their maxima at the fourth drawing round, which agreed with the aforementioned conclusion pertaining to the evolution of DB. ΔSw and ΔSp stabilised during the common drawing stage, then remained unchanged (at around 0), while ΔWLB stabilised during cyclic drawing, then approached 0.865 m, which approximated to the working face advance interval and drawing interval.

Evolution of TCB
TCB refers to the boundary between the top coal and the immediate roof rock. The TCB formed by coal caving was analysed in this section. The TCB after each caving round was extracted, as shown in Figure 11. As shown in Figure 11, the shapes of the TCB are different on the goaf side and the working face side. On the goaf side (red dotted line), the TCB is parabolic; its shape was significantly affected by the first three support advances and drawing rounds and it remained practically unchanged since the fourth drawing round and was no longer affected by the advance of the working face or the coal caving. For the working face side (purple solid line in the figure), the evolution of the TCB was differed substantially from that of the goaf side during the advance of the working face, and its characteristics are summarised as follows: ① There is a definite horizontal asymptote and the curve is asymmetrical. Horizontally, the curve will approach the original TCB, that is, the original TCB is its definite asymptote. Vertically, the formation of the curve is affected by gravity and the drawing opening, the boundary condition differs from that applied to the horizontal face, indicating that it is asymmetrical.
② The curve approaches the horizontal asymptote at a slower rate. The curve above the hydraulic support shield beam, tail beam, and the rear scraper conveyor increases Figure 11. The evolutional process of TCB.
As shown in Figure 11, the shapes of the TCB are different on the goaf side and the working face side. On the goaf side (red dotted line), the TCB is parabolic; its shape was significantly affected by the first three support advances and drawing rounds and it remained practically unchanged since the fourth drawing round and was no longer affected by the advance of the working face or the coal caving. For the working face side (purple solid line in the figure), the evolution of the TCB was differed substantially from that of the goaf side during the advance of the working face, and its characteristics are summarised as follows: 1 There is a definite horizontal asymptote and the curve is asymmetrical. Horizontally, the curve will approach the original TCB, that is, the original TCB is its definite asymptote. Vertically, the formation of the curve is affected by gravity and the drawing opening, the boundary condition differs from that applied to the horizontal face, indicating that it is asymmetrical. plumb line after the rear scraper conveyor as the X-axis and the vertical upwards direction taken as positive, with the horizontal line of the original TCB as the Y-axis and the advance direction of the working face as the positive direction. That is, the TCB lay within the second quadrant of the coordinate system (x < 0). To reduce errors caused by definition-domain limitation during the fitting process, absolute values were taken for the X-coordinate of each point, and a curve symmetrical to the Y-axis with the TCB was obtained. Fitting was performed on this curve. The fitting formula was determined based on the evolution characteristics of the TCB. The curve had an asymptote and was asymmetrical about its own axes, as is the Hook function (of the form f (x) = x + 1/x). The method of determination of the fitting model is described in Figure 12.
each caving event, the TCB deflects toward the working face advance direction at t scraper conveyor, as shown by the purple dotted arrowhead in Figure 11.
To reveal the evolution of the TCB, a fitting model was selected to fit the TCB working face side. A rectangular coordinate system was established (Figure 11), w plumb line after the rear scraper conveyor as the X-axis and the vertical upwards di taken as positive, with the horizontal line of the original TCB as the Y-axis and the a direction of the working face as the positive direction. That is, the TCB lay wit second quadrant of the coordinate system (x < 0). To reduce errors caused by def domain limitation during the fitting process, absolute values were taken for the Xnate of each point, and a curve symmetrical to the Y-axis with the TCB was ob Fitting was performed on this curve. The fitting formula was determined based evolution characteristics of the TCB. The curve had an asymptote and was asymm about its own axes, as is the Hook function (of the form f(x) = x + 1/x). The me determination of the fitting model is described in Figure 12. The dashed lines in Figure 12 are the Hook function and its two terms (g1(x) g2(x) = 1/x). The expression of Hook function is f1(x) = g1(x) + g2(x) = x + 1/x. This fu agrees with the TCB in terms of the horizontal asymptote and its asymmetry, but of change exceeds that of the TCB. Fitting for the TCB cannot be performed throug the aforementioned two terms. To control the rate of change of Y along the nega rection of X, g1(x) = x is substituted by g3(x) = x 2 , and to control it along the positiv tion of X, g4(x) = x 1/2 was added. The fitting formula was determined to be as follow where a, b, c, and d are fitting parameters. It can be shown that a, b, and d > 0, and c < position and curvature radius of the stagnation point were calculated. x was differe to get F(x). Its first-order derivative and second-order derivative are shown in Form The dashed lines in Figure 12 are the Hook function and its two terms (g 1 (x) = x and g 2 (x) = 1/x). The expression of Hook function is f 1 (x) = g 1 (x) + g 2 (x) = x + 1/x. This function agrees with the TCB in terms of the horizontal asymptote and its asymmetry, but its rate of change exceeds that of the TCB. Fitting for the TCB cannot be performed through only the aforementioned two terms. To control the rate of change of Y along the negative direction of X, g 1 (x) = x is substituted by g 3 (x) = x 2 , and to control it along the positive direction of X, g 4 (x) = x 1/2 was added. The fitting formula was determined to be as follows: where a, b, c, and d are fitting parameters. It can be shown that a, b, and d > 0, and c < 0. The position and curvature radius of the stagnation point were calculated. x was differentiated to get F(x). Its first-order derivative and second-order derivative are shown in Formula (2): We set F'(x) = 0 to obtain the x-coordinate of the stagnation point. The radius of curvature thereat was calculated using Equation (3): The TCB after each caving round was fitted and calculated based on Equations (1)-(3). The results are collated in Table 4. The fitting curves are shown in Figure 13. Table 4. Fitting parameters.

Caving
Step We set F'(x) = 0 to obtain the x-coordinate of the stagnation point. The radius of curvature thereat was calculated using Equation (3): The TCB after each caving round was fitted and calculated based on Equations (1) to (3). The results are collated in Table 4. The fitting curves are shown in Figure 13.    As shown in Figure 13a, the evolution of the TCB can be divided into three stages: the first stage was from the first drawing round to the third drawing round. During this stage, the evolution of the TCB was relatively stable. The second stage was that from the fourth drawing round to the ninth drawing round. A change in the trend arose between each pair of groups, that is, the TCB of each two adjacent drawing rounds was put into one group (three groups in total in this stage), with the shape of each group showing stable development. The third stage was after the tenth drawing round. After this stage, the TCB had a relatively stable shape and advanced with the working face. For the stagnation point position (Figure 13b (Figure 13c), in the first stage, it was reduced gradually from 6.48 m (first drawing round) to 3.67 m (third drawing round); in the second stage, it was decreased from 4.47 m (fourth drawing round) to 3.15 m (fifth drawing round) in the first group by 1.32 m and from 5.42 m (sixth drawing round) to 3.54 m (seventh drawing round) by 1.88 m in the second group, but remained relatively stable at 5.00 m (eighth drawing round) and 5.49 m (ninth drawing round) in the third group; in the third stage, it was reduced slightly but remained relatively high, with an average radius of 7.48 m and a standard deviation of 0.94 m. Based on the above information, the evolution of the TCB curve can be divided into three stages from the first caving stage to the common caving stage from the perspective of its shape, stagnation point height and radius of curvature at the stagnation point. The first stage was influenced most by the initial mining process, the second stage represented a transitional drawing stage, and the third stage was the common drawing stage.

Discussion
(1) Numerical model Based on the FEM-DEM coupled method, a "continuous-discontinuous" numerical model was constructed with the 12309 working face of the Wangjialing coal mine as the engineering background to address the actual problems posed in this paper. 1 In this paper we focused on the drawing of the granular top coal without considering its crushing. The linear contact constitutive model was used for all "ball-ball" contact. The top coal and the immediate roof were all regarded as granular media, therefore the top coal and the immediate roof did not break during the advance of the working face. 2 During mining operations before caving, the extent of crushing of the top coal did not meet the requirement for caving: the top coal block size must exceed that measured in-situ as the mining pressure was relatively low during this stage, however, as the data cannot be measured, parameters were assigned according to the result of numerical simulation. 3 The motion process of top coal and rocks above the hydraulic power support is affected not only by one hydraulic support, but also inevitably by the support advancement and coal drawing in its neighbouring sections. In the actual condition, there are usually more than 100 hydraulic supports arranged on a working face; besides, the top coal caving mining methods vary greatly, which leads to difference in coal drawing sequence, such as single round sequential coal drawing, single round interval coal drawing, multi-round interval coal drawing, etc. Different coal drawing methods are suitable for different coal seam conditions, the influence of which on the motion of top coal and rocks varies greatly; actually, even different coal drawing time can have different impact on them. In our views, although there are too many uncertainties in the study on the above issues, obviously we can still find a certain rule to follow. We need to consider the influence of multiple supports in synergetic mining on top coal and rocks, which is also one of our top priorities in the future researches. In this paper, taking the "FDM-DEM" numerical model as the carrier, we conduct the study on the movement and drawing rule of loose top coal; especially, our research focus is on the transition process from the first coal drawing stage to the common coal drawing stage; in view of this, we only conduct the 2D analysis. 4 We had tried 2-d software for top coal caving simulation, but the result was not satisfactory, mainly because the particle gaps significantly increased and the force exerted was beyond the range of the 2-d model when the range of the particle size was relatively large, and because all simulated particles were spheres with a regular packing in the 2-d state. The DB was more diamond-shaped than ellipsoidal, producing a poor inversion result. To avoid these problems, 3-d software was used to build the numerical model.
(2) Limitations and weaknesses of the method In this paper, a coupling numerical simulation method based on "FDM-DEM" is proposed, which is realized by the bridging function between FLAC 3D and PFC 3D . The method enables the simulation of top coal caving mining process, including the simulation of hydraulic support and rear scraper conveyor; furthermore, it enables the automatic control of simulation process by using Fishcall function. However, there are obvious limitations and weaknesses in this method. The contact constitutive model between top coal particles is a linear model, that is, in the process of numerical simulation, top coal and immediate roof are considered to be in the loose state (This method has been used by most scholars in their related researches), so the crushing process of top coal and immediate roof under the action of mining pressure is ignored. Currently, our research team has attempted to establish a numerical method for the whole process, which enables the inversion of mining pressure as well as the progressive failure of top coal and immediate roof in the numerical simulation from continuous state to loose state (Particles are no longer simple "Balls") to final drawing. Moreover, it will also be one of our main research directions in the future.
(3) Advance abutment pressure The evolutional process of the TCB is shown in Figure 12. The TCB in front of the working face basically was a horizontal line in first to seventh drawing rounds. However, the line changed into the shape of a wave from the eighth drawing round. The reason for this is as follows: more and more top coal particles were deleted as the working face advanced. The immediate roof particles moved downward to fill the goaf, causing a disconnection between the immediate roof and the main roof (the blank part in Figure 12). The stress of the main roof and the part above was passed to the coal wall front, resulting in stress concentration. Loading and unloading was caused by the unloosening of the top coal and immediate particles within the loose area, particles being "squeezed out." As a result, the horizontal TCB appeared wavy.

(4) Transitional Process
In LTCC, the three parameters, the DB, LB, and TCB underwent a transition from the first drawing round to the common drawing stage. The simulation results had some similarities with the relevant researches of other scholars. But in the previous researches, the transition process from the first coal drawing round to the common coal drawing stage has been ignored. Seen from the numerical simulation results, this transitional process had a clear turning point. For the DB and the LB, the fourth drawing round was the key turning point, while for the TCB, the fourth and tenth drawing rounds were both key turning points. As some hypotheses were made in numerical model building and parameter selection, the result did not directly correspond to engineering practice but could provide certain guidance. We believe that the aforementioned transitional process must occur in practice, but due to the complexity of the in-situ conditions, the key transitional point cannot be a certain caving round and there may be multiple key transitional points.

Conclusions
A "continuous-discontinuous" numerical model was established to invert the evolution of top coal migration characteristics in LTCC. In terms of DB, LB, and TCB, the transition process from the first drawing round to the common drawing stage was analysed. The following conclusions are drawn:

•
The evolution of DB parameters such as A DB , W DB , H DB , L DB , and θ DB from first drawing round to common drawing stage showed a trend whereby they first decreased, then increased, decreased again, and finally stabilised. Specifically, the parameters were observed to reach their respective maxima after the first drawing round and then gradually decrease, only to increase again after the fourth drawing round, and then stabilise. After stabilisation, A DB was around 6 t, W DB was 2.7 m, H DB was 4.0 m, L DB was around 4.0 m, and θ DB was around 10 • . In the numerical simulation, the fourth drawing round was the key turning point in the DB evolution.

•
The evolution of LB parameters such as W LB , S w , and S p and increments thereof was analysed from the first drawing round to the common drawing round. S w and S p were close to 0 m in the common caving stage, while W LB was close to the caving interval (0.865 m) in this stage. In the numerical simulation, the fourth drawing round was the key turning point in the LB evolution.

•
The TCBs formed in each drawing round were fitted. From the perspective of the height and radius of curvature stagnation point, the evolution of the TCB was analysed.
The process was divided into three stages in the numerical simulation: the first (first to third drawing rounds) was the initial mining influence stage, the second (fourth to ninth drawing rounds) was the transitional caving stage, and the third (after tenth drawing round) represented the common drawing stage.

•
As some hypotheses were made in relation to the boundary conditions, the model, and numerical simulation process, its result did not directly correspond to actual engineering conditions, but can provide reference to some extent. It is believed in this present study that the above-mentioned transitional process must exist during mining operations, but due to the complexity of the in-situ conditions, the key transition point cannot be a certain drawing round and there may be multiple key transition points.

•
We believe that the research results of this paper are applicable to other coal mines in the world with occurrence conditions similar to those of Wangjialing coal mine with the same mining method. In the process of numerical modeling, all parameters are obtained according to the actual site conditions. Under different engineering backgrounds, the obtained parameters are different. Therefore, this numerical method has high flexibility, which is suitable for the simulation of top coal movement in any top coal caving mining.