Numerical Investigation of Top-Coal Migration in the First Coal-Drawing Process by an FDM–DEM Coupling Method

: In Longwall Top-Coal Caving (LTCC), the shape ofthe loose body (LB) and top-coalboundary (TCB) formed after the ﬁrst coal-drawing are the initial conditions for the common coal-drawing period. Taking the Panel 12309 in Wangjialing coal mine as the research object, the weight of the blocks of caved top coal was measured on-site, and the distribution of their equivalent diameter was calculated. By using a coupled numerical method, the “Finite di ﬀ erence method (FDM)–Discrete element method (DEM)” numerical model was established, and the evolutions of the drawing body (DB), LB, and TCB were obtained. The results show that, in the initial stage of ﬁrst coal-drawing (0–8.09 s), the amount of DB reached its maximum of 7.18 m 3 (0.89 m 3 / s) and then decreased to a stable value of 0.44 m 3 / s. The relationships between the characteristic parameters of DB, LB, and TCB and the drawing time were ﬁtted. Taking the second derivative of each parameter with respect to time as its sensitivity ( η ), it was concluded that due to the large coal-drawing volume in the initial stage (0–8.09 s), the values of the above parameters increased, and the sensitivities reached 6.02 × 10 − 3 , 3.09 × 10 − 3 , and 6.99 × 10 − 3 , respectively. Here, the top-coal migration rule in the ﬁrst coal-drawing process was revealed from the perspectives of DB, LB, and TCB, thus providing a theoretical basis for the further study of common coal-drawing processes.


Introduction
Several longwall mining methods such as Multi-Slice Longwall (MSL), High Reach Single Pass Longwall (SPL), and Longwall Top-coal Caving (LTCC) have been mainly applied to extract thick coal seams. Among them, LTCC can provide substantially less economic sensitivity and a reduced level of technical risk over others. Therefore, LTCC has been widely used to extract thick coal seams in China since the 1980s [1,2]. In LTCC, the coal seam is divided into two parts (the bottom coal and the top coal). Bottom coal is undercut directly by the shearer. Top-coal transitions from a continuous state to a discontinuous crushing state under front abutment pressure and finally caves from the coal-drawing opening under the action of gravity [3]. Most research studies about LTCC focus on the mechanism of stress distribution, top-coal crushing, and top-coal migration with the aim of improving recovery ratio [4,5]. However, a rational method of selecting and optimising coal-drawing process parameters can be obtained only on the basis of an in-depth understanding of top-coal migration. The shape of the loose body (LB) and top-coal boundary (TCB) formed after the first coal-drawing represent the initial conditions of the later coal-drawing period. Therefore, it is necessary and meaningful to study the first coal-drawing process, which is the basis of the whole LTCC mining operation. roof is siltstone with an average thickness of 4.2 m. The schematic of the coal seam and roofs [24] is shown in Figure 1.

Caving Processes and Equipment
In the mining operation of Panel 12309, the cutting (bottom coal) and caving (top coal) height are 3.0 m and 3.5 m, respectively, the cutting-caving ratio is 1:1.167, and the drawing interval is equal to web cut of shearer (0.865 m). One hundred fifty-two sets of ZFY12000/23/34D two-column top-coal hydraulic supports were arranged on the working face, and the length of the top beam and cover beam of the supports was 7.19 m. When the draw opening was closed, the angle between the tail beam and the horizontal plane was 23°, while when the draw opening was open, the angle between the tail beam and the horizontal plane was 58°. The SGZ-1000/2000 rear-scraper conveyor with the middle width of 1m was adopted. The hydraulic support and scraper conveyors are shown in Figure  2.

Test Purpose
In the process of LTCC, under the influence of mining pressure, the top-coal gradually breaks up, transitioning from a quasi-continuous state to a granular state [25,26]. According to previous research results, the distribution of top-coal lumpiness affects its migration to a certain extent. The migration of top-coal is not the same under different lumpiness distribution conditions. Moreover, due to different geological conditions and mining-production technical conditions, it is difficult to predict the lumpiness distribution of top coal after crushing [27]; therefore, to improve the reliability of the numerical model, it is necessary to install field instrumentation to measure the distribution of top-coal lumpiness in the 12,309 working face. 3.0 Siltstone Fine sandstone Coal · · · · · · · · · · · · · · · · · · · · · · · · · · ·· ·· ·· ·· ·· ·· ·· ·· ·· ·· ·· ·· ·· ·· ·· ·· ·· ··

Caving Processes and Equipment
In the mining operation of Panel 12309, the cutting (bottom coal) and caving (top coal) height are 3.0 m and 3.5 m, respectively, the cutting-caving ratio is 1:1.167, and the drawing interval is equal to web cut of shearer (0.865 m). One hundred fifty-two sets of ZFY12000/23/34D two-column top-coal hydraulic supports were arranged on the working face, and the length of the top beam and cover beam of the supports was 7.19 m. When the draw opening was closed, the angle between the tail beam and the horizontal plane was 23 • , while when the draw opening was open, the angle between the tail beam and the horizontal plane was 58 • . The SGZ-1000/2000 rear-scraper conveyor with the middle width of 1m was adopted. The hydraulic support and scraper conveyors are shown in Figure 2.
Energies 2020, 13, x 3 of 19 roof is siltstone with an average thickness of 4.2 m. The schematic of the coal seam and roofs [24] is shown in Figure 1.

Caving Processes and Equipment
In the mining operation of Panel 12309, the cutting (bottom coal) and caving (top coal) height are 3.0 m and 3.5 m, respectively, the cutting-caving ratio is 1:1.167, and the drawing interval is equal to web cut of shearer (0.865 m). One hundred fifty-two sets of ZFY12000/23/34D two-column top-coal hydraulic supports were arranged on the working face, and the length of the top beam and cover beam of the supports was 7.19 m. When the draw opening was closed, the angle between the tail beam and the horizontal plane was 23°, while when the draw opening was open, the angle between the tail beam and the horizontal plane was 58°. The SGZ-1000/2000 rear-scraper conveyor with the middle width of 1m was adopted. The hydraulic support and scraper conveyors are shown in Figure  2.

Test Purpose
In the process of LTCC, under the influence of mining pressure, the top-coal gradually breaks up, transitioning from a quasi-continuous state to a granular state [25,26]. According to previous research results, the distribution of top-coal lumpiness affects its migration to a certain extent. The migration of top-coal is not the same under different lumpiness distribution conditions. Moreover, due to different geological conditions and mining-production technical conditions, it is difficult to predict the lumpiness distribution of top coal after crushing [27]; therefore, to improve the reliability of the numerical model, it is necessary to install field instrumentation to measure the distribution of top-coal lumpiness in the 12,309 working face. 3.0 Siltstone Fine sandstone Coal · · · · · · · · · · · · · · · · · · · · · · · · · ·

Test Purpose
In the process of LTCC, under the influence of mining pressure, the top-coal gradually breaks up, transitioning from a quasi-continuous state to a granular state [25,26]. According to previous research results, the distribution of top-coal lumpiness affects its migration to a certain extent. The migration of top-coal is not the same under different lumpiness distribution conditions. Moreover, due to different geological conditions and mining-production technical conditions, it is difficult to predict the lumpiness distribution of top coal after crushing [27]; therefore, to improve the reliability of the numerical model, it is necessary to install field instrumentation to measure the distribution of top-coal lumpiness in the 12,309 working face.

Test Method
In order to protect the gateways, the top coal near the two ends of the working face will not be caved in LTCC mining. Therefore, 75#, 16#, and 134# hydraulic supports were selected from the middle and two ends of the working face, respectively. We opened the designated coal-drawing opening to draw the top coal and transport it to the stage loader, where the tester takes samples of, and measurements on, the top coal using an explosion-proof electronic scale ( Figure 3).

Test Method
In order to protect the gateways, the top coal near the two ends of the working face will not be caved in LTCC mining. Therefore, 75#, 16#, and 134# hydraulic supports were selected from the middle and two ends of the working face, respectively. We opened the designated coal-drawing opening to draw the top coal and transport it to the stage loader, where the tester takes samples of, and measurements on, the top coal using an explosion-proof electronic scale ( Figure 3). No less than 200 samples of top-coal were randomly selected for each support. First, an explosion-proof electronic scale is used to weigh each block. The equivalent diameter of each block is then calculated using Equation (1).
where d is equivalent particle diameter, M is the mass of a top-coal block, and ρc is the density of the top-coal block.

Test Results
Three hundred fifteen, 240, and 208 groups of data were measured from the 16#, 75#, and 124# hydraulic support respectively, and a total of 763 values were obtained. According to Equation (1), the data were counted, and the results of the number and mass distribution of blocks in different lumpiness ranges were obtained ( Figure 4). No less than 200 samples of top-coal were randomly selected for each support. First, an explosion-proof electronic scale is used to weigh each block. The equivalent diameter of each block is then calculated using Equation (1).
where d is equivalent particle diameter, M is the mass of a top-coal block, and ρ c is the density of the top-coal block.

Test Results
Three hundred fifteen, 240, and 208 groups of data were measured from the 16#, 75#, and 124# hydraulic support respectively, and a total of 763 values were obtained. According to Equation (1), the data were counted, and the results of the number and mass distribution of blocks in different lumpiness ranges were obtained ( Figure 4).  It can be seen from Figure 4a that the maximum lumpiness of top-coal caved by support No. 75 is 24-28 cm, while that from supports No. 16 and No.134 is 32-36 cm. The top-coal crushing in the middle of the working face is more complete. In Figure 4b, the equivalent diameter of caved top-coal blocks was divided into four intervals (0-9 cm, 9-18 cm, 18-27 cm, and 27-36 cm), and the mass fractions were 13.83%, 46.31%, 20.34%, and 27.36%, respectively. In the range of 9-18 cm, the mass fraction is the largest, while in the range of 0-9 cm, the mass fraction is the smallest. In the FDM-DEM coupled numerical model, the "Ball" unit in PFC 3D was used to simulate the top coal, and the model was established based on the results in Figure 4b.

Numerical Model
The FLAC 3D module is called in the PFC 3D software to establish the FDM-DEM numerical model of top-caving coal, including the coal seam, immediate roof, main roof, and a set of hydraulic supports and the rear-scraper conveyor. The length, width, and height of the numerical model were 45 m, 1.75 m and 16.1 m, respectively. The top coal and immediate roof were modelled by discontinuous unit (Ball), and the size distribution of top-coal particles conforms to the measured results summarised in Section 2.2. In view of the unmeasurable distribution of the lumpiness of the immediate roof, it was assumed that the equivalent diameter of the immediate roof was four times that of the top coal [7,13,20], with the same probability distribution ( Table 1). The bottom coal and the main roof were modelled by continuous unit (Zone). The main roof was located at the top of the whole model, and the unit size was 1 m. Hydraulic supports and the rear-scraper conveyor were simulated by "Wall" elements in PFC 3D . In the modelling process, the bottom coal model was first established, and then the top-coal and the immediate roof particles were successively established by means of random particle generation, automatic gravity field balance, and gradual servo-controlled compaction under vertical stress until the actual thickness was reached. The numerical model is shown in Figure 5.  It can be seen from Figure 4a that the maximum lumpiness of top-coal caved by support No. 75 is 24-28 cm, while that from supports No. 16 and No.134 is 32-36 cm. The top-coal crushing in the middle of the working face is more complete. In Figure 4b, the equivalent diameter of caved top-coal blocks was divided into four intervals (0-9 cm, 9-18 cm, 18-27 cm, and 27-36 cm), and the mass fractions were 13.83%, 46.31%, 20.34%, and 27.36%, respectively. In the range of 9-18 cm, the mass fraction is the largest, while in the range of 0-9 cm, the mass fraction is the smallest. In the FDM-DEM coupled numerical model, the "Ball" unit in PFC 3D was used to simulate the top coal, and the model was established based on the results in Figure 4b.

Numerical Model
The FLAC 3D module is called in the PFC 3D software to establish the FDM-DEM numerical model of top-caving coal, including the coal seam, immediate roof, main roof, and a set of hydraulic supports and the rear-scraper conveyor. The length, width, and height of the numerical model were 45 m, 1.75 m and 16.1 m, respectively. The top coal and immediate roof were modelled by discontinuous unit (Ball), and the size distribution of top-coal particles conforms to the measured results summarised in Section 2.2. In view of the unmeasurable distribution of the lumpiness of the immediate roof, it was assumed that the equivalent diameter of the immediate roof was four times that of the top coal [7,13,20], with the same probability distribution ( Table 1). The bottom coal and the main roof were modelled by continuous unit (Zone). The main roof was located at the top of the whole model, and the unit size was 1 m. Hydraulic supports and the rear-scraper conveyor were simulated by "Wall" elements in PFC 3D . In the modelling process, the bottom coal model was first established, and then the top-coal and the immediate roof particles were successively established by means of random particle generation, automatic gravity field balance, and gradual servo-controlled compaction under vertical stress until the actual thickness was reached. The numerical model is shown in Figure 5.  The actual average buried depth of coal seam was 360 m, and the height of the numerical model was 16.1 m; therefore, during the establishment of the initial balance in the model, vertical stress should be applied to the surface of the main roof to simulate the weight of the rock layer with a thickness of 345.86 m lying above. Calculated using an average bulk unit weight for this rock stratum of 22 kN/m 3 [24], the main roof should be applied to a vertical stress of 7.57 MPa. Elastic models were used for bottom coal and main roof, and linear models were used for contact between Ball-Ball and Ball-Facet components ( Table 2 lists unit parameters; Table 3 lists contact parameters). The linear model provides the behaviour of an infinitesimal interface that does not resist relative rotation so that the contact moment (Mc) equals zero. The contact force is resolved into linear and dashpot components. The linear component provides linear elastic (no-tension) frictional behaviour, while the dashpot component provides viscous behaviour (Figure 6a). The linear force is produced by linear springs with constant normal and shear stiffnesses, kn and ks. The dashpot force is produced by dashpots with viscosity given in terms of the normal and shear critical-damping ratios, βn and βs. The linear springs act in parallel with the dashpots. A surface gap, gs, is defined as the difference between the contact gap gc and the reference gap gr, so that when the reference gap is zero, the notional surfaces coincide with the piece surfaces ( Figure 6b). The contact is active if and only if the surface gap is less than or equal to zero; the force-displacement law is skipped for inactive contacts. The actual average buried depth of coal seam was 360 m, and the height of the numerical model was 16.1 m; therefore, during the establishment of the initial balance in the model, vertical stress should be applied to the surface of the main roof to simulate the weight of the rock layer with a thickness of 345.86 m lying above. Calculated using an average bulk unit weight for this rock stratum of 22 kN/m 3 [24], the main roof should be applied to a vertical stress of 7.57 MPa. Elastic models were used for bottom coal and main roof, and linear models were used for contact between Ball-Ball and Ball-Facet components ( Table 2 lists unit parameters; Table 3 lists contact parameters). The linear model provides the behaviour of an infinitesimal interface that does not resist relative rotation so that the contact moment (M c ) equals zero. The contact force is resolved into linear and dashpot components. The linear component provides linear elastic (no-tension) frictional behaviour, while the dashpot component provides viscous behaviour (Figure 6a). The linear force is produced by linear springs with constant normal and shear stiffnesses, k n and k s . The dashpot force is produced by dashpots with viscosity given in terms of the normal and shear critical-damping ratios, β n and β s . The linear springs act in parallel with the dashpots. A surface gap, g s , is defined as the difference between the contact gap g c and the reference gap g r , so that when the reference gap is zero, the notional surfaces coincide with the piece surfaces ( Figure 6b). The contact is active if and only if the surface gap is less than or equal to zero; the force-displacement law is skipped for inactive contacts.

Numerical Methods
In the numerical simulation, the actual production process used at the LTCC face was fully considered, and the whole numerical calculation process was divided into two stages: a preparation stage and the first coal-drawing stage. Among them, the preparation stage was used to invert the process whereby the bottom coal was cut off without top-coal caving. The first coal-drawing stage was used to invert the first coal-drawing process of the LTCC working face.
In the preparation stage, the excavation of the working face before coal-drawing is simulated, including the extraction of open-off cut and the installation of the hydraulic support and the rearscraper conveyor (Figure 7). With the coal-drawing opening closed, the operation of cutting bottom coal, shifting supports, and the rear-scraper conveyor shall be carried out in succession. When the moving boundary of the LB reaches the immediate roof, the calculation of the first stage shall be terminated. The second stage is the first coal-drawing stage. On the basis of the first stage, the coal-drawing opening is opened, and the calculation is terminated when the immediate roof particle reaches the rear-scraper conveyer. In this stage, the monitoring of each indicator is realised by the FISHCALL function which allows for the FISH function to be executed at a specific point in the cycle sequence during cycling. In the present work, the FISH function is inserted at point −11.0 (before validate data structures in PFC) of each calculation cycle to achieve two ends: (1) delete the particles that fall from the coal opening to the rear-scraper conveyor; (2) when immediate roof particles fall into the scraper conveyor, the calculation will be terminated. The process of numerical simulation is shown in Figure  8.

Numerical Methods
In the numerical simulation, the actual production process used at the LTCC face was fully considered, and the whole numerical calculation process was divided into two stages: a preparation stage and the first coal-drawing stage. Among them, the preparation stage was used to invert the process whereby the bottom coal was cut off without top-coal caving. The first coal-drawing stage was used to invert the first coal-drawing process of the LTCC working face.
In the preparation stage, the excavation of the working face before coal-drawing is simulated, including the extraction of open-off cut and the installation of the hydraulic support and the rear-scraper conveyor ( Figure 7). With the coal-drawing opening closed, the operation of cutting bottom coal, shifting supports, and the rear-scraper conveyor shall be carried out in succession. When the moving boundary of the LB reaches the immediate roof, the calculation of the first stage shall be terminated.

Numerical Methods
In the numerical simulation, the actual production process used at the LTCC face was fully considered, and the whole numerical calculation process was divided into two stages: a preparation stage and the first coal-drawing stage. Among them, the preparation stage was used to invert the process whereby the bottom coal was cut off without top-coal caving. The first coal-drawing stage was used to invert the first coal-drawing process of the LTCC working face.
In the preparation stage, the excavation of the working face before coal-drawing is simulated, including the extraction of open-off cut and the installation of the hydraulic support and the rearscraper conveyor (Figure 7). With the coal-drawing opening closed, the operation of cutting bottom coal, shifting supports, and the rear-scraper conveyor shall be carried out in succession. When the moving boundary of the LB reaches the immediate roof, the calculation of the first stage shall be terminated. The second stage is the first coal-drawing stage. On the basis of the first stage, the coal-drawing opening is opened, and the calculation is terminated when the immediate roof particle reaches the rear-scraper conveyer. In this stage, the monitoring of each indicator is realised by the FISHCALL function which allows for the FISH function to be executed at a specific point in the cycle sequence during cycling. In the present work, the FISH function is inserted at point −11.0 (before validate data structures in PFC) of each calculation cycle to achieve two ends: (1) delete the particles that fall from the coal opening to the rear-scraper conveyor; (2) when immediate roof particles fall into the scraper conveyor, the calculation will be terminated. The process of numerical simulation is shown in Figure  8. The second stage is the first coal-drawing stage. On the basis of the first stage, the coal-drawing opening is opened, and the calculation is terminated when the immediate roof particle reaches the rear-scraper conveyer. In this stage, the monitoring of each indicator is realised by the FISHCALL function which allows for the FISH function to be executed at a specific point in the cycle sequence during cycling. In the present work, the FISH function is inserted at point −11.0 (before validate data structures in PFC) of each calculation cycle to achieve two ends: (1) delete the particles that fall from the coal opening to the rear-scraper conveyor; (2) when immediate roof particles fall into the scraper conveyor, the calculation will be terminated. The process of numerical simulation is shown in Figure 8.

Preparation Stage
After the equipment is installed, the top-coal begins to fall ( Figure 9A), and the termination condition of the first stage is reached after cutting the bottom coal and advance the support five times. The boundary of the LB after each extraction was extracted based on the particle displacement of up to 25 mm, as shown in Figure 9.  Figure 9 shows that the boundary of LB expands with the advance of the working face during the preparation stage. The endpoint position of the boundary of the LB is relatively stable and is located at the transition point between the top beam and the cover beam of the support on the side of working face and at the apex of the bottom coal pillar on the side of the goaf. Before moving the support, the maximum height of the LB is 1.27 m ( Figure 9A). After moving the support, the maximum height of LB increases continuously, and the rate of increase thereof decreases over time.

Preparation Stage
After the equipment is installed, the top-coal begins to fall ( Figure 9A), and the termination condition of the first stage is reached after cutting the bottom coal and advance the support five times. The boundary of the LB after each extraction was extracted based on the particle displacement of up to 25 mm, as shown in Figure 9.

Preparation Stage
After the equipment is installed, the top-coal begins to fall ( Figure 9A), and the termination condition of the first stage is reached after cutting the bottom coal and advance the support five times. The boundary of the LB after each extraction was extracted based on the particle displacement of up to 25 mm, as shown in Figure 9.   Figure 9A). After moving the support, the maximum height of LB increases continuously, and the rate of increase thereof decreases over time. Whether the immediate roof particles have been caved? Figure 9. Displacement of particles in the first stage. Figure 9 shows that the boundary of LB expands with the advance of the working face during the preparation stage. The endpoint position of the boundary of the LB is relatively stable and is located at the transition point between the top beam and the cover beam of the support on the side of working face and at the apex of the bottom coal pillar on the side of the goaf. Before moving the support, the maximum height of the LB is 1.27 m ( Figure 9A). After moving the support, the maximum height of LB increases continuously, and the rate of increase thereof decreases over time. The boundary Energies 2020, 13, 5493 9 of 18 of the LB reaches the immediate roof after the fourth support-shifting operation, and the particles in the immediate roof move locally after the fifth support-shift.

•
Total DB The total DB refers to all top-coal particles caved from the drawing opening during the first coal-drawing stage. In the numerical simulation, the first coal-drawing stage lasted for 64.72 s, with an error of 10.1% from the actual average first coal-drawing time of 72 s, which was divided into eight time intervals, and the caved particles in each interval were grouped ( Figure 10). The boundary of the LB reaches the immediate roof after the fourth support-shifting operation, and the particles in the immediate roof move locally after the fifth support-shift.

First Coal-Drawing Stage
The total DB refers to all top-coal particles caved from the drawing opening during the first coaldrawing stage. In the numerical simulation, the first coal-drawing stage lasted for 64.72 s, with an error of 10.1% from the actual average first coal-drawing time of 72 s, which was divided into eight time intervals, and the caved particles in each interval were grouped ( Figure 10). In Figure 10, the shape of the total DB is shown in the model after initial equilibrium (Figure 10a) and the model before coal-drawing (Figure 10b). In the model after initial equilibrium, the total DB protrudes at the side of the goaf (yellow line in Figure 10a), while in the model before coal-drawing, it is a smooth curve (yellow line in Figure 10b); on the side of the working face, the curves delimiting the shape of the total DB in the two models coincide. In the model after initial equilibrium and the model before coal-drawing, the position of the immediate roof particles caved at the end of coaldrawing all deflected to goaf relative to the centreline of the rear-scraper conveyor, and the deflection angle by the highest particle that moves was consistent (4°). The maximum horizontal width and vertical height of the total DB at each time point were extracted and the data plotted (Figure 11a). The vertical line at the mid-point of the rear-scraper conveyor was taken as the axis for extraction of the range of the total DB at each time point on the working face and goaf side, and the histogram was thereof plotted (Figure 11b).   In Figure 10, the shape of the total DB is shown in the model after initial equilibrium (Figure 10a) and the model before coal-drawing (Figure 10b). In the model after initial equilibrium, the total DB protrudes at the side of the goaf (yellow line in Figure 10a), while in the model before coal-drawing, it is a smooth curve (yellow line in Figure 10b); on the side of the working face, the curves delimiting the shape of the total DB in the two models coincide. In the model after initial equilibrium and the model before coal-drawing, the position of the immediate roof particles caved at the end of coal-drawing all deflected to goaf relative to the centreline of the rear-scraper conveyor, and the deflection angle by the highest particle that moves was consistent (4 • ). The maximum horizontal width and vertical height of the total DB at each time point were extracted and the data plotted (Figure 11a). The vertical line at the mid-point of the rear-scraper conveyor was taken as the axis for extraction of the range of the total DB at each time point on the working face and goaf side, and the histogram was thereof plotted (Figure 11b).
It can be seen from Figure 11a that the height (h TDB ) and width (L TDB ) of the total DB increase continuously with the increase of the coal-drawing time (t). Among them, h TDB and t conform to an exponential relationship, while L TDB and t conform to a linear relationship as follows: In the first coal-drawing process, the rate of change of L TDB is constant, while that of h TDB gradually decreases. It can be seen from Figure 11b  from 1.67 m to 3.89 m, increasing by 133%. As affected by the hydraulic support, the top coal is more likely to be caved on the working face side with lower frictional resistance, so L TDB on the working face side is always larger than that on the goaf side. On the working face side, the scope of the DB is limited by the hydraulic support, so the rate and amount of expansion on the goaf side are slightly higher than those on the working face side. and the model before coal-drawing (Figure 10b). In the model after initial equilibrium, the total DB protrudes at the side of the goaf (yellow line in Figure 10a), while in the model before coal-drawing, it is a smooth curve (yellow line in Figure 10b); on the side of the working face, the curves delimiting the shape of the total DB in the two models coincide. In the model after initial equilibrium and the model before coal-drawing, the position of the immediate roof particles caved at the end of coaldrawing all deflected to goaf relative to the centreline of the rear-scraper conveyor, and the deflection angle by the highest particle that moves was consistent (4°). The maximum horizontal width and vertical height of the total DB at each time point were extracted and the data plotted (Figure 11a). The vertical line at the mid-point of the rear-scraper conveyor was taken as the axis for extraction of the range of the total DB at each time point on the working face and goaf side, and the histogram was thereof plotted (Figure 11b).   •

Interval DB
The interval DB refers to the top-coal blocks caved by the coal-drawing opening within a specific time interval. The first coal-drawing stage was divided into eight time intervals, and the top-coal particle information relating to caved-in material in each time interval was extracted and marked in the initial model ( Figure 12). It can be seen from Figure 11a that the height (hTDB) and width (LTDB) of the total DB increase continuously with the increase of the coal-drawing time (t). Among them, hTDB and t conform to an exponential relationship, while LTDB and t conform to a linear relationship as follows: In the first coal-drawing process, the rate of change of LTDB is constant, while that of hTDB gradually decreases. It can be seen from Figure 11b that, after 8.09 s of coal-drawing, LTDB in the goaf side (1.67 m) is smaller than the width of the working face (2.06 m). With the increase of coal-drawing time, the LTDB towards both the working face and the goaf side expand constantly. On the side of the working face, LTDB increases from 2.06 m to 4.13 m (an increase of 100.49%). On the side of the goaf, LTDB increased from 1.67 m to 3.89 m, increasing by 133%. As affected by the hydraulic support, the top coal is more likely to be caved on the working face side with lower frictional resistance, so LTDB on the working face side is always larger than that on the goaf side. On the working face side, the scope of the DB is limited by the hydraulic support, so the rate and amount of expansion on the goaf side are slightly higher than those on the working face side.

•
Interval DB The interval DB refers to the top-coal blocks caved by the coal-drawing opening within a specific time interval. The first coal-drawing stage was divided into eight time intervals, and the top-coal particle information relating to caved-in material in each time interval was extracted and marked in the initial model ( Figure 12).  Figure 12 shows that, except for the first interval (0-8.09 s), the interval DB was an ellipsoidal cut inclined to the working face. In addition, at the initial stage of coal-drawing, the horizontal width (LIDB) and vertical height (hIDB) of the interval DB both exceeded those in the later stage. The width (LIDB) and height (hIDB) of the shape of interval DB in each interval were extracted and plotted ( Figure  13a). The volume of particles caved in each interval was extracted and the histogram plotted ( Figure  13b).  Figure 12 shows that, except for the first interval (0-8.09 s), the interval DB was an ellipsoidal cut inclined to the working face. In addition, at the initial stage of coal-drawing, the horizontal width (L IDB ) and vertical height (h IDB ) of the interval DB both exceeded those in the later stage. The width (L IDB ) and height (h IDB ) of the shape of interval DB in each interval were extracted and plotted (Figure 13a). The volume of particles caved in each interval was extracted and the histogram plotted (Figure 13b). Energies 2020, 13 Figure 13a shows that both LIDB and hIDB first decrease and then tend to be stable with increasing coal-drawing time, and they follow exponential relationships therewith as follows: In the first coal-drawing process, LIDB decreased faster and by more than hIDB. It can be seen from Figure 13b that, in the initial stage of coal-drawing (0-8.09 s), the amount of coal in the DB was at its maximum, reaching 7.18 m 3 (0.89 m 3 /s), and it then decreased to a stable value of 0.44 m 3 /s. That is, when the coal-drawing opening was just opened, the amount of top-coal caved was maximised, gradually decreased, and then tended to be stable.

Evolution of the LB
The LB refers to the loose part caused by the migration of surrounding coal and rock to the mined-out space during mining. The shape of the LB at nine time points was extracted based on the particle displacement of up to 25 mm ( Figure 14).    Figure 13a shows that both L IDB and h IDB first decrease and then tend to be stable with increasing coal-drawing time, and they follow exponential relationships therewith as follows: In the first coal-drawing process, L IDB decreased faster and by more than h IDB . It can be seen from Figure 13b that, in the initial stage of coal-drawing (0-8.09 s), the amount of coal in the DB was at its maximum, reaching 7.18 m 3 (0.89 m 3 /s), and it then decreased to a stable value of 0.44 m 3 /s. That is, when the coal-drawing opening was just opened, the amount of top-coal caved was maximised, gradually decreased, and then tended to be stable.

Evolution of the LB
The LB refers to the loose part caused by the migration of surrounding coal and rock to the mined-out space during mining. The shape of the LB at nine time points was extracted based on the particle displacement of up to 25 mm ( Figure 14).  Figure 13a shows that both LIDB and hIDB first decrease and then tend to be stable with increasing coal-drawing time, and they follow exponential relationships therewith as follows: In the first coal-drawing process, LIDB decreased faster and by more than hIDB. It can be seen from Figure 13b that, in the initial stage of coal-drawing (0-8.09 s), the amount of coal in the DB was at its maximum, reaching 7.18 m 3 (0.89 m 3 /s), and it then decreased to a stable value of 0.44 m 3 /s. That is, when the coal-drawing opening was just opened, the amount of top-coal caved was maximised, gradually decreased, and then tended to be stable.

Evolution of the LB
The LB refers to the loose part caused by the migration of surrounding coal and rock to the mined-out space during mining. The shape of the LB at nine time points was extracted based on the particle displacement of up to 25 mm ( Figure 14).   Figure 14. Shape of the loose body (LB). Figure 14 shows that in the coal-drawing process, the LB expands to varying degrees in both horizontal and vertical directions. At the beginning of coal-drawing (t = 0), as affected by the preparation stage, the horizontal width of the LB (L LB ) was 8.53 m, and the vertical height of the LB (h LB ) was 6.79 m. At the end of coal-drawing, L LB was 16.90 m (an increase of 98.12%) and h LB was 9.22 m (an increase of 35.79%). The width (L LB ) and height (h LB ) of the LB at each time point were extracted (Figure 15a). The vertical line at the mid-point of the rear-scraper conveyor was taken as the axis from which to extract the range of the LB at each time point on the working face and goaf side, and the histogram was plotted (Figure 15b).
Energies 2020, 13, x 12 of 19 Figure 14 shows that in the coal-drawing process, the LB expands to varying degrees in both horizontal and vertical directions. At the beginning of coal-drawing (t = 0), as affected by the preparation stage, the horizontal width of the LB (LLB) was 8.53 m, and the vertical height of the LB (hLB) was 6.79 m. At the end of coal-drawing, LLB was 16.90 m (an increase of 98.12%) and hLB was 9.22 m (an increase of 35.79%). The width (LLB) and height (hLB) of the LB at each time point were extracted (Figure 15a). The vertical line at the mid-point of the rear-scraper conveyor was taken as the axis from which to extract the range of the LB at each time point on the working face and goaf side, and the histogram was plotted (Figure 15b). As shown in Figure 15a, the height (hLB) and width (LLB) of the LB increase continuously with the increase of the coal-drawing time (t). Among them, hLB and t conform to an exponential relationship, while LLB and t conform to a linear relationship as follows: In the first coal-drawing process, the rate of increase of LLB was constant, while that of hLB gradually decreased. It can be seen from Figure 15b that, at the beginning of the first coal-drawing, the LLB in the goaf side (5.92 m) exceeded the width of the working face (2.61 m). With increased coaldrawing time, LLB both towards the working face and the goaf side expanded constantly. On the side of the working face, LLB increased from 2.61 m to 7.69 m, increasing by 1.95 times. On the side of the goaf, LLB increased from 5.92 m to 9.21 m (an increase of 55.57%). The rate of expansion of the LB on the working face side was significantly higher than that on the goaf side.

Evolution of TCB
TCB refers to the boundary surface between coal and gangue. The immediate roof particle size distribution and boundary shape at each time point were extracted. A Cartesian coordinate system was established with the mid-point of the rear-scraper conveyor as the origin, the floor as the x-axis, and the vertical line as the y-axis (Figure 16). As shown in Figure 15a, the height (h LB ) and width (L LB ) of the LB increase continuously with the increase of the coal-drawing time (t). Among them, h LB and t conform to an exponential relationship, while L LB and t conform to a linear relationship as follows: In the first coal-drawing process, the rate of increase of L LB was constant, while that of h LB gradually decreased. It can be seen from Figure 15b that, at the beginning of the first coal-drawing, the L LB in the goaf side (5.92 m) exceeded the width of the working face (2.61 m). With increased coal-drawing time, L LB both towards the working face and the goaf side expanded constantly. On the side of the working face, L LB increased from 2.61 m to 7.69 m, increasing by 1.95 times. On the side of the goaf, L LB increased from 5.92 m to 9.21 m (an increase of 55.57%). The rate of expansion of the LB on the working face side was significantly higher than that on the goaf side.

Evolution of TCB
TCB refers to the boundary surface between coal and gangue. The immediate roof particle size distribution and boundary shape at each time point were extracted. A Cartesian coordinate system was established with the mid-point of the rear-scraper conveyor as the origin, the floor as the x-axis, and the vertical line as the y-axis (Figure 16). Figure 16 shows that in the first coal-drawing process, the range of movement of the TCB gradually increased, and its lowest point gradually dropped to the rear-scraper conveyor. The range of movement of the TCB increased 1.86 times from 5.45 m at the beginning of coal-drawing to 15.59 m at the end of coal-drawing. The vertical distance between the TCB and rear-scraper conveyor is 6.12 m at the beginning of coal-drawing. At the end of coal-drawing, the TCB intersects the rear-scraper conveyor (that is, the immediate roof particles were caved by the rear-scraper conveyor). To show the evolution of the TCB more intuitively, the morphology of the TCB at each time point is plotted (Figure 17). Energies 2020, 13, x 13 of 19 Figure 16. Shape of the top-coal boundary (TCB). Figure 16 shows that in the first coal-drawing process, the range of movement of the TCB gradually increased, and its lowest point gradually dropped to the rear-scraper conveyor. The range of movement of the TCB increased 1.86 times from 5.45 m at the beginning of coal-drawing to 15.59 m at the end of coal-drawing. The vertical distance between the TCB and rear-scraper conveyor is 6.12 m at the beginning of coal-drawing. At the end of coal-drawing, the TCB intersects the rearscraper conveyor (that is, the immediate roof particles were caved by the rear-scraper conveyor). To show the evolution of the TCB more intuitively, the morphology of the TCB at each time point is plotted ( Figure 17).     Figure 16 shows that in the first coal-drawing process, the range of movement of the TCB gradually increased, and its lowest point gradually dropped to the rear-scraper conveyor. The range of movement of the TCB increased 1.86 times from 5.45 m at the beginning of coal-drawing to 15.59 m at the end of coal-drawing. The vertical distance between the TCB and rear-scraper conveyor is 6.12 m at the beginning of coal-drawing. At the end of coal-drawing, the TCB intersects the rearscraper conveyor (that is, the immediate roof particles were caved by the rear-scraper conveyor). To show the evolution of the TCB more intuitively, the morphology of the TCB at each time point is plotted (Figure 17).     figure). The TCB on the working face side underwent the most significant evolution above the back part of shield beam and coal-drawing opening. On the side of the goaf, the TCB gradually developed downward, and its rate of growth gradually decreased. As affected by the preparation stage, the TCB as shown at t = 0 is formed, and its lowest point is inclined towards the goaf. The TCB began to evolve from the above position. The evolution track of the lowest point position of TCB is shown as the red dotted line in the figure, which is quasi-parabolic and gradually approaches the centreline of the rear-scraper conveyor (blue dotted line). We extracted the moving span (L TCB ) and minimum height (h TCB ) of the TCB (Figure 18a). The vertical line at the mid-point of the rear-scraper conveyor was taken as the axis to extract the range of the TCB at each time point on the working face and goaf side, and the histogram was plotted (Figure 18b). from the above position. The evolution track of the lowest point position of TCB is shown as the red dotted line in the figure, which is quasi-parabolic and gradually approaches the centreline of the rearscraper conveyor (blue dotted line). We extracted the moving span (LTCB) and minimum height (hTCB) of the TCB (Figure 18a). The vertical line at the mid-point of the rear-scraper conveyor was taken as the axis to extract the range of the TCB at each time point on the working face and goaf side, and the histogram was plotted (Figure 18b). As shown in Figure 18a, the moving span of the TCB (LTCB) increases with prolonged coaldrawing time (t). The minimum height of the TCB (hTCB) decreases with prolonged coal-drawing time (t). There is an exponential relationship between LTCB and t, and hTCB has a linear relationship with t.
In the first coal-drawing process, the rate of increase of the moving span of the TCB gradually decreases, while the rate of decrease of the lowest height is unchanged in the first and middle stages of coal-drawing and increases slightly at the end of coal-drawing. According to Figure 18b, as influenced by the preparation stage, at the beginning of coal-drawing, the range of LTCB on the goaf side (5.3 m) exceeds that on the working face side (0.51 m). With the beginning of the first coaldrawing operation, the range of movement of the TCB expands on both the working face side and the goaf side. On the side of the working face, LTCB increased from 0.51 m to 6.95 m, increasing by 12.63 times. On the side of the goaf, LTCB increased from 5.30 m to 8.64 m (an increase of 63.02%). The rate of expansion of the range of movement of the TCB on working face side is significantly higher than that on the goaf side.

Coal-Drawing Sensitivity
The fitting equations of parameters of total DB, interval DB, LB, and TCB and coal-drawing time are summarised in Table 4  As shown in Figure 18a, the moving span of the TCB (L TCB ) increases with prolonged coal-drawing time (t). The minimum height of the TCB (h TCB ) decreases with prolonged coal-drawing time (t). There is an exponential relationship between L TCB and t, and h TCB has a linear relationship with t.
In the first coal-drawing process, the rate of increase of the moving span of the TCB gradually decreases, while the rate of decrease of the lowest height is unchanged in the first and middle stages of coal-drawing and increases slightly at the end of coal-drawing. According to Figure 18b, as influenced by the preparation stage, at the beginning of coal-drawing, the range of L TCB on the goaf side (5.3 m) exceeds that on the working face side (0.51 m). With the beginning of the first coal-drawing operation, the range of movement of the TCB expands on both the working face side and the goaf side. On the side of the working face, L TCB increased from 0.51 m to 6.95 m, increasing by 12.63 times. On the side of the goaf, L TCB increased from 5.30 m to 8.64 m (an increase of 63.02%). The rate of expansion of the range of movement of the TCB on working face side is significantly higher than that on the goaf side.

Coal-Drawing Sensitivity
The fitting equations of parameters of total DB, interval DB, LB, and TCB and coal-drawing time are summarised in Table 4. Herein, the first coal-drawing is divided into eight equal time intervals. The fitting results for L IDB and h IDB are exponential functions; therefore, in the early stage of the first coal-drawing, due to the sudden opening of the coal-drawing opening, a large amount of top-coal is caved, exceeding that subsequently (Figure 13b). The relationship between the acceleration of each parameter in the evolution process and the drawing time is helpful to analyse the sensitivity of each parameter to time in the whole process. Therefore, the second derivative of each parameter with respect to t is taken as its coal-drawing sensitivity (η), and the results (Table 4) are discussed. The sensitive-y curves of each parameter are shown in Figure 19.
LB LLB = 0.13t + 8.864 hLB = −2.56e −0.041t + 9.349 TCB LTCB = −11.91e −0.027t + 17.67 hTCB = −0.081t + 6.317 Herein, the first coal-drawing is divided into eight equal time intervals. The fitting results for LIDB and hIDB are exponential functions; therefore, in the early stage of the first coal-drawing, due to the sudden opening of the coal-drawing opening, a large amount of top-coal is caved, exceeding that subsequently (Figure 13b).
The relationship between the acceleration of each parameter in the evolution process and the drawing time is helpful to analyse the sensitivity of each parameter to time in the whole process. Therefore, the second derivative of each parameter with respect to t is taken as its coal-drawing sensitivity (η), and the results (Table 4) are discussed. The sensitive-y curves of each parameter are shown in Figure 19. As shown in Figure 19, the height of the total DB, the height of the LB, and the moving width of the TCB are more sensitive to the coal-drawing time. The sensitivity of the DB is such that ηDB = (8.89 × 10 −3 )e −0.048t , the sensitivity of the LB is such that ηLB = (4.31 × 10 −3 )e −0.041t , and the sensitivity of the TCB is such that ηTCB = (8.69 × 10 −3 )e −0.027t . Due to the large volume of top-coal caved in the initial stage (0-8.09 s), the rate of change of the aforementioned parameters increased, and the sensitivities reached 6.02 × 10 −3 , 3.09 × 10 −3 , and 6.99 × 10 −3 , respectively.

Geometry of the Numerical Model
The width of the numerical model is 1.75 m, which is equal to the width of a top-coal-caving hydraulic support. Only the caving process of a single hydraulic support is simulated in our manuscript. In this work, 2D simulation is carried out by using 3D numerical software, and the reasons are as follows: (1) In most of the theoretical researches on top-coal caving, many researchers have simplified the three-dimensional problem to a plane problem. Based on reasonable simplification of the problem, the influence of various factors on the top-coal caving process has been eliminated to achieve more targeted research on the caving process.
(2) We have drawn from the experience of numerical simulation that the 2D model has drawbacks in dealing with top-coal caving. When the particles vary greatly in size, the gap between the particles will increase significantly in the 2D model, which adversely affects the contact between As shown in Figure 19, the height of the total DB, the height of the LB, and the moving width of the TCB are more sensitive to the coal-drawing time. The sensitivity of the DB is such that η DB = (8.89 × 10 −3 )e −0.048t , the sensitivity of the LB is such that η LB = (4.31 × 10 −3 )e −0.041t , and the sensitivity of the TCB is such that η TCB = (8.69 × 10 −3 )e −0.027t . Due to the large volume of top-coal caved in the initial stage (0-8.09 s), the rate of change of the aforementioned parameters increased, and the sensitivities reached 6.02 × 10 −3 , 3.09 × 10 −3 , and 6.99 × 10 −3 , respectively.

Geometry of the Numerical Model
The width of the numerical model is 1.75 m, which is equal to the width of a top-coal-caving hydraulic support. Only the caving process of a single hydraulic support is simulated in our manuscript. In this work, 2D simulation is carried out by using 3D numerical software, and the reasons are as follows: (1) In most of the theoretical researches on top-coal caving, many researchers have simplified the three-dimensional problem to a plane problem. Based on reasonable simplification of the problem, the influence of various factors on the top-coal caving process has been eliminated to achieve more targeted research on the caving process. (2) We have drawn from the experience of numerical simulation that the 2D model has drawbacks in dealing with top-coal caving. When the particles vary greatly in size, the gap between the particles will increase significantly in the 2D model, which adversely affects the contact between the particles. Moreover, the simulated particles are all spheres, and their arrangement regularity is too strong in the 2D state. Therefore, the top-coal drawing body cannot form similar ellipsoid and is closer to the diamond. As a result, the various parameters cannot be simulated well in the 2D model. For these reasons, 3D software was selected instead of 2D software in our work.

Numerical Simulation Process
The numerical model and simulation process in this paper are deficient in the following two aspects: Energies 2020, 13, 5493 16 of 18 (1) In the present work, we analysed only the first round of top-coal drawing; however, in the in situ condition, it is more important to analyse the drawing procedures after the first drawing. Based on this paper, we will further establish a common top-coal drawing model and carry out in-depth researches on the top-coal drawing process after the first top-coal drawing. (2) So far, we have modelled the top coal and immediate roof with a linear model, which indicates that the top coal is already in a granular state. The undercut of bottom coal will not influence the progressive fracturing of the top coal anymore. The emphasis of our work is on the migration rule of the discontinuous crushed top coal. Of course, we firmly believe that it is the best solution if the failure process, migration rule, and drawing characteristics of top-coal can be simulated in one numerical model. However, there is still no scientific method to realize the above process, and it will become an important work of our future research.

Conclusions
Taking the working face 12309 of Wangjialing Coal Mine as the engineering background and based on in situ measurements of top-coal lumpiness, the "continual-discontinuity" numerical model was established, and the first coal-drawing process was numerically simulated. The numerical calculation process was divided into two stages: a preparation stage and the first coal-drawing stage. The main conclusions are as follows: • The equivalent weighing method was used to measure the distribution of top-coal lumps in situ. The number of top-coal lumps was inversely proportional to their equivalent diameter. The equivalent diameter was divided into four intervals (0-9 cm, 9-18 cm, 18-27 cm, and 27-36 cm), and their mass fractions were 13.83%, 46.31%, 20.34%, and 27.36%, respectively. The relationship between the characteristic parameters of the total DB, interval DB, LB, and TCB and the coal-drawing time was fitted. The second derivative of each parameter with respect to time was taken as its sensitivity. It was concluded that the height of the total DB, the height of the LB, and the moving width of the TCB were more sensitive to the coal-drawing time.
Due to the large volume of top-coal caved in the initial stage (0-8.09 s), the rates of change of the aforementioned parameters increased, and the sensitivities reached 6.02 × 10 −3 , 3.09 × 10 −3 , and 6.99 × 10 −3 , respectively.

•
In this paper, the evolution law of DB, LB, and TCB of top-coal in the first coal drawing process was obtained, but the numerical model still needs to be improved to carry out in-depth research on the top-coal drawing process after the first top-coal drawing, as does a scientific method to realize the process that takes failure process, migration rule, and drawing characteristics of top-coal into account simultaneously.