A Shaft Pillar Mining Subsidence Calculation Using Both Probability Integral Method and Numerical Simulation

2018-12-10 02:51PeixianLiZhixiangTanandLiliYan

Peixian Li, Zhixiang Tan and Lili Yan

Abstract: In order to prolong the life cycle of the coal mine, Jinggezhuang (‘JGZ’) coal mine decided to excavate the shaft pillar.The first panel 0091 was designed near the pillar boundary as an experiment in shaft pillar mining.Both probability integral method(PIM) and FLAC3D were used to evaluate the influence on the shaft safety.PIM parameters were obtained from previous surface subsidence station.The rock property is based on the lab mechanical test.A simulated FLAC3D model containing shafts and a panel was built based on stratigraphic information.Surface subsidence results of PIM show that the 0091-panel excavation has no influence on the shafts.The simulated results show that the subsidence of the main shaft and air shaft is small and can be ignored, but it could cause the auxiliary shaft 220 mm horizontal displacement.So, the stress and displacement of the underground part of shaft were analyzed, it shows that the stress changes, subsidence and displacement are mainly located at the top part of the shafts.According to the stress and movement of the simulated shafts, 0091 was decided to be excavated and a surface monitor line was built and measured.In comparison of PIM,FLAC3D, and measured data, the PIM results fit the surface subsidence better.And the FLAC3D results have smaller maximum subsidence and greater influence area than measured.But FLAC3D can provide more details such as displacement, subsidence, stress and strain of both surface and underground.So, for a planned mining excavation, both methods should be used especially for the evaluation of deformation of underground constructions.In the future, with the development of the rock numerical computation technology, the numerical simulation method will be recommended first.The research shows compare of two methods of the coal mine subsidence calculation and provides a solution method for shaft pillar mining.

Keywords: Mining subsidence, deformation monitor, PIM, rock mechanical, FLAC3D.

1 Introduction

After excavation of the underground coal, the goaf space makes a new unbalance of the overlaying rock.The rock will collapse, results in deformation and redistribution until a new balance attained.The deformation will spread to surface if the excavation space big enough.This process and phenomenon are called ‘mining subsidence’ [Zhao and Meng(2010)].Mining subsidence can cause the surface constructions such as house, railway,water body, etc.to be destroyed and also create environment problems.People have recognized this problem very early; the history of mining subsidence research can be traced back to the 15thcentury.In the first years of the 20thcentury, a lot of surface deformation surveys were done by coal mine companies, and mining subsidence became a new area of research.Then, with the massive energy demand of industry, the research of mining subsidence increased, and a lot of mining subsidence prediction theories and methods were proposed based on both continuum and discontinuous mechanics theory[Cui and Deng (2017)].Those methods played important roles in reducing the damage caused by coal mine subsidence.In China, a stochastic medium theory from Poland was introduced and used in coal mine subsidence calculation in the 1960s [Litwiniszyn(1957,1974); Liu and Liao (1965)].The method was improved into an easier model named probability integral method (‘PIM’) [Liu and Liao (1965); He, Yang, Ling et al.(1991); Cui and Deng (2017)].PIM is an influence function to calculate the surface subsidence based on the geometry integration of excavation space.There are 8 parameters of PIM that can be obtained from back-analysis of the site measured data.It is the most used method in China [Coal Industry Bureau of People’s Republic of China(2000)].But PIM is based on the influence function, it is a semi-empirical function which depends on the parameters what used but not on the deformation mechanism of the rock mass.So, it can calculate only the ground surface displacement but cannot be used to calculate deformation and stress distribution of under lying rock.

With development of the rock simulation theories and computer performance, numerical simulation technology become an important way to deal with coal mine engineering problems [Costa, Zingano and Koppe (2000); Nuric, Nuric, Kricak et al.(2000); Pereira,Costa, Salvadoretti et al.(2012)].Simulation of mining subsidence also have been used widely since the 1980s [Chrzanowski, Monahan, Roulston et al.(1997); Szostak and Chrzanowski (1991)].Fast Lagrangian Analysis of Continua in three dimensions(FLAC3D) is a numerical modeling and simulation software for rock, groundwater,constructions, coal mine and ground support [Bock (2015); Zhao, Jiang and Lan (2015);Shibata, Zarlin, Shimada et al.(2014)].FLAC3Dcan define both linear and nonlinear constitutive models for the simulation elements.The models can calculate deformation and displacement with large deformation mode.The software uses the explicit Lagrange difference to deal with 3D, large area plastic and flowing deformation problems with small computer memories.

JGZ coal mine is located in the east of China; the coal resource is in a basin-like syncline.It is about 3.5 Km from North to South and 3.4 Km East to West, the area is 9.23 Km2.The coal mine is in the alluvial valley plain area and has great traffic network.The elevation is +38.9 m in the north and +23.9 m in the south, average +35.0 m.The coal mine was first built in 1958 and began production in 1978.It has 140 million tons of coal resource of which 72.68 million tons are recoverable reserves.The designed production is 1.2 million tons per year.The coal mine development method is mainly multi-level based on vertical shaft.There are three shafts located in the industrial site.With about 30 years excavation, JGZ coal mine faced the problem of resource exhausting.The coal mine manager had a plan for mining the industrial site pillar.In order to maintain the shaft safety,an experimental panel was designed first.This paper predicted the mining subsidence of the experimental panel with both PIM and FLAC3D, to evaluate excavation influence on the shaft, and then results comparison with the data measured was done in the paper.

2 Models of PIM and FLAC3D

2.1 Introduction of PIM model

PIM is an influence function of the mining surface subsidence prediction model based on stochastic medium theory [Litwiniszyn (1957); Liu and Liao (1965)].In order to calculate the surface subsidence, the coal mine influenced rock can be hypothesized as stochastic, discontinuous elements as shown in Fig.1(a).

The figure shows that if one element of the first layer is excavated, one of two elements in the second layer will drop into the excavated space, and it is random which of the two elements will drop in.With this hypothesis, suppose that the subsidence will spread to the surface after one element of the first layer excavation; the surface subsidence curvature and distribution are similar with the Gauss bell function as shown in Fig.1(b) [Yang, Liu and Wang (2004); Zhang (1998)].For the large space excavation, the subsidence can be calculated by the geometry integral.As shown in Fig.2, surface point A (x, y)’s movement can be decomposed into subsidence W (x, y) and horizontal displacement U (x,y, φ).Djis the area of jthpanel which have influence to point A.According to the theory of PIM, the surface subsidence W (x, y) and horizontal displacement U (x, y, φ) can be calculated by Eq.(1)-Eq.(2) [He, Ling, Yang et al.(1991)].In China, PIM is the most used function for coal mine subsidence prediction and plays an important role in reduction the loss of mining subsidence.

Figure 1:Theory model of particles movements as a stochastic medium [He, Yang, Ling et al.(1991)]

Figure 2:A diagram of surfer movement calculated by PIM

where (x, y) are coordinate of the surface point, W, U are mining subsidence and horizontal displacement, φ is direction angle of surface movement; W0is the maximum ground subsidence, given by W0=mq cosα; We(x, y) is subsidence of a small unit mining;H is mining depth; m is mining thickness; D is calculation mining area of the panel, n is total number of panels, j is the number of panel; α is dip angle of coal seam; q is subsidence factor; b is displacement factor; tanβ is tangent of major influence angle; r is major influence radius given by r = H/tanβ; dsdt is integration variable of double integral of area Dj.

The PIM function results depend on those parameters.The parameters can be obtained from back analysis of measured data [Li, Peng, Tan et al.(2017)] or by the experience function of overlaying rock mechanical parameters [Coal Industry Bureau of People’s Republic of China (2000)].

2.2 Introduction of FLAC3D

FLAC3Dis a simulation software developed by ITASCA company.It is a 3D finite difference solving software.It can be used for rock mechanical simulation and plastic flow in soil, rock and other material.It has a beautiful performance on large deformation [Xie,Zhou, Wang et al.(1999)].Compared with other simulation systems like ANASYS,ADINA, FLAC3Dcan solve the physical instability problems.It can be used in stability analysis of highly nonlinear geometry and physical problems like coal mine rock deformation, tunnel crack, and stress field of surrounding rock [Li, Tan, Gu et al.(2015);Rawamurthy (2014)].

FLAC3Duses the dynamic explicit method to get the elements’ time step solution, then the model’s failure and crack can be traced.It is important to research time and space effect of the rock deformation.FLAC3Dhas 11 different material constitutive models like null model,elastic model, and mohr-coulomb model, etc.Those materials can be used to simulate the complex geological and mining environment.The software also has powerful pre- and postprocessing functions.Complicated models can be built with its basic grid and a lot of results like stress, strain, and displacement can be print out with values or plots.

There are 6 basic steps of simulation with FLAC3Das shown in Fig.3.

Figure 3:Flow chart of simulation with FLAC3D

Step 1.First of all, finite difference grids similar to the initial geological should be built.The software provides 13 different grids which can simulate arbitrarily topography.

Step 2.Then, the grid should be given a constitutive model and material parameters.And the model also needs a boundary and initial condition.

Step 3.Next step is to get the initial rock stress state.

Step 4.The model excavated or changed according to the actual engineer conditions.

Step 5.Calculation to get a new balance for the engineering changes.

Step 6.The results can be printed out and analyzed after the rock is re-balanced through calculation.

3 Mining subsidence calculation and its parameters

3.1 Basic geological of 0091 panel

There are four coal seams in the research area; the 9#coal seam is which investigated.0091 panel is the first experimental excavation of industrial site pillar.The panel is located at the south borderline of the industrial site pillar.The average elevation of the panel is -295 m and the surface elevation are 35.5 m, mining thickness is 7 m, coal seam dip angle is 4o~8o, average dip angle is about 7o.The panel’s length of strike is 1050 m,length of tendency is 112 m, strike azimuth 317o17’, the geometry and location are shown in Fig.4, and basic information of 0091 is shown in Tab.1.As shown in Fig.4, there are three vertical shafts in the industrial site.

Figure 4:Geometry of the 0091 panel

Table 1:Basic information of panel 0091

The strata of this area from bottom to top are Carboniferous upper series (C3), Permian lower series (P1), and Quaternary loose deposit (Q).The thickness of the Quaternary loose layer is about 160 m.The geology of investigated panel is simple.The geometry and geology of section 1-1 are shown in Fig.5.The direct roof of the 9#coal is mud stone,shale, sandstone, etc.And direct floor of 9#coal is clay stone.The 0091 panel is the first panel in industrial site pillar; an observation station was set in order to get the displacement law of the ground and protect the constructions in the industrial site.

Figure 5:Geology of Section 1-1

3.2 Parameters used to predict mining subsidence based on PIM

To use the PIM method, there are 8 parameters that need to be determined first.Those parameters are often obtained based on previous measured mining subsidence data [Zha,Feng and Zhu (2001)].The JGZ coal mine and other collieries in Tangshan area have built a lot of sites monitoring stations to study the ground movement, and PIM parameters are shown in Tab.2.

Table 2:Some of measured PIM parameters

In this research, the average parameters are used to calculate the 0091-mining surface subsidence.

3.3 Rock mechanical parameters of FLAC3D model

3.3.1 Boundary condition and initial stress

A reliable boundary condition is the basic requirement to get correct results of the mechanical simulation.In this shaft pillar mining simulation, the model boundary is displacement boundary.It meets the following conditions:

1.Fix the horizontal displacement of the periphery.As shown in Fig.6, fix the X-direction displacement in right and left side plane, and fix the Y-direction displacement in front and back side plane.

2.Fix the vertical displacement of the bottom;

3.Make the top surface of the model move free.

The model is simplified for ignoring tectonic stress, it just in the hydrostatic state.Gravity stress field is the only thing to be considered.The initial stress of model depends on the weight and property of the rock.The model’s boundary conditions are shown in Fig.6.In addition, to reduce the error of the boundary effect, distances between the model geometry boundary and the edge of the ground subsidence basin are more than 200 m.The model geometry is big enough to simulation for 0091 panel excavation deformation.

Figure 6:Boundary conditions of the model

3.3.2 The geometry

The simulation geometry of the model depends on the actual geological conditions and its boundary conditions.The 3D model was established by layers according to exploration results, and the 0091-excavation area is represented as the null elements, as shown in Fig.7.The model elements were built with wedge-shaped mesh which is a built-in grid of FLAC3D.The null model was used for 0091panel excavation.

Figure 7:The geometry and 0091 panel of the model

3.3.3 Rock material properties

FLAC3Dsoftware can deal with null model, elastic model and plastic model.Mohr-coulomb model is a mathematical yield failure criterion of materials; it is used widely in rock numerical[Sainoki and Mitri (2017); Cao, Wang, Li et al.(2018)], so mohr-coulomb model was selected as the constitutive model for rock.The shafts’ linings were built by steel reinforced concrete and then the elastic model was selected.Rock mechanical property is the key point to simulation results.Some of the JGZ coal mine rock properties are shown in Tab.3.

Table 3:Mechanical property of rocks and coals (Laboratory value)

Those rock mechanical parameters are based on drill core tests.The drill core tests are based on the test pieces, but the rock mass is large scale.There are a lot of cracks in the rock mass.Some research shows that strength of rock mass is about 1/10~1/30 of test pieces [Chen, Xie, Jing et al.(2006); Kang, Hu, Sin et al.(2017); Kodama, Miyamoto,Kawasaki et al.(2013)].In the paper, the bulk modulus and shear elasticity modulus applied in the model was 1/20 of the laboratory values as in Tab.3.

4 Result, comparison and discussion

4.1 Results of surface displacement

First, PIM and FLAC3Dwere used to calculate the subsidence of the 0091-panel using the average parameters in Tab.2 and rock properties in Tab.3.The results are shown in Fig.8 in which the contour is the subsidence and the vector is the horizontal displacement.

In China, 10 mm contour line was defined as the edge of the subsidence basin [Coal Industry Bureau of People’s Republic of China (2000)].As shown in Fig.8(a), the maximum subsidence located above the center of 0091 panel and the value reduced from center to outside.The subsidence concentrated in a certain distance from panel center.The horizontal displacement is directing to the panel center.The biggest horizontal displacement is located above the panel boundary.The horizontal displacement is smaller in panel center and basin boundary.The horizontal displacement influence area is almost equal to subsidence influence area.In Fig.8(a), Shafts are not influenced by the 0091-coal mine subsidence.It is safety to excavate the pillar panel 0091.In Fig.8(b), the maximum subsidence located above the center of the panel.The horizontal displacement is directing to the panel center and the biggest horizontal displacement are located above the panel boundary.The horizontal displacement influence area is bigger than the subsidence influence area.As shown in the Fig.8(b), main shaft and auxiliary shaft are located in the mining subsidence basin area.It shows that those shafts are influenced by the 0091 excavations.It needs more analysis to decide if it is safe to excavate the panel 0091.More information can be exported from the FLAC3D.

By comparing two figures, more details can be found.First, the basic appearance of subsidence and horizontal displacement are similar.But the FLAC3Dsubsidence basin area is bigger than PIM results.The ratio between subsidence and horizontal displacement is different; FLAC3D’s horizontal displacement influence area is bigger than PIM results.The influence results on shafts are different, PIM results have no influence to shafts and FLAC3Dhave two shafts influenced.The differences of basic principles make those results difference, more information and analysis are needed for decision.

4.2 Shaft displacement and stress

As the FLAC3Dcan simulate the mechanical behavior of rock mass, more information for engineering evaluation such as stress can be exported from the model.In order to maintain the shaft safety, the subsidence (z-axis displacement), x-axis displacement, y-axis displacement, stress-ZZ(σzz), stress-YY(σyy), stress-XX(σxx) are exported and shown in Fig.9.As shown in Tab.4, the maximum subsidence of shafts is on the top of the auxiliary shaft,the subsidence is 18 mm.The subsidence of the main shaft and air shaft is less than 10 mm.

Figure 8:Results of surface displacement:PIM VS FLAC3D

Figure 9:Shaft wall deformation and stress

The horizontal displacement of the auxiliary shaft is 187 mm in x-axis and 137 mm in yaxis which is the maximum horizontal displacement of three shafts.But as shown in Fig.9, the displacement and stress of shafts are all located at the top part which near the surface.The Fig.9(d) shows that it has little stress concentration on shaft wall.The stress is mainly from gravity effect, but not from excavation of 0091 panel.So, synthetically with subsidence results and stress distribution of the shafts which are obtained from both PIM and FLAC3D, the management decided to excavate the panel 0091 and measured the surface subsidence to keep the shaft from being damaged.

Table 4:Results of shaft displacement based on FLAC3D

FLAC3Dalso can provide more information of internal rock mass, as shown in Fig.10.

Figure 10:Displacement and stress of Section 1-1

They are contours of z-axis displacement, y-axis displacement, x-axis displacement and stress-zz (σzz) of Section 1-1.The deformation of the inner rock is symmetrical of the panel center.The subsidence decreases progressively from the roof of the panel to the surface.The floor of the panel’s z-axis displacement is positive which indicates that the floor of panel is rising.It simulated the floor unloading phenomenon after the panel excavation.Fig.10(d) shows that there are two stress concentration areas above the panel which is located at the panel boundary.It is similar to the actual excavation.

4.3 Results comparison and analysis

4.3.1 Results comparison

With the analysis above, the company decided to excavate the 0091 panel.And subsidence measurement of B line in Fig.4 was surveyed from September 26, 2008 to April 30, 2010.And subsidence of those points was obtained.In order to analyses the results, the PIM results, FLAC3Dsimulated results and measured data are plotted in Fig.11.There are some similarities and differences.

Figure 11:Results comparison

As shown in Fig.11, all results of both subsidence and horizontal displacement are symmetry over the center of the 0091 panel.The maximum subsidence is located at the center of the panel.The horizontal displacement direction is towards to the center of the goaf and the maximum displacement is located above the boundary of the panel.Both FLAC3Dand PIM results are fit for the basic law of the mining subsidence.But there are some differences, as shown in Fig.11(a); the subsidence boundary of the PIM calculated results is smaller than measured but bigger than the FLAC3Dresults.The FLAC3Dresults have larger range of influence than measured.As shown in Fig.11(a), the angles of draw of three results are 39o, 47o, 57o.The maximum subsidence result is WPIM>WMeasured>WFLAC3D.The PIM results fit better with measured data than FLAC3D.The horizontal displacement was not measured; PIM and FLAC3Dhorizontal results are compared in Fig.11(b).Both PIM and FLAC3Dinfluence area of horizontal displacement are bigger than subsidence.The maximum displacements of PIM results are greater than FLAC3D.The maximum displacement locations of PIM are above the boundary of panel, but the FLAC3Dresults are located at outside of the panel boundary.Measured data shows that the shafts were not influenced by the 0091 excavations.The shafts are still on service until 2018.The evaluation with PIM and FLAC3Dprovided reliable results to keep the shaft safety.

4.3.2 Results analysis

As shown in Fig.8 and Fig.11, there are significant differences such as influence area,value and position of maximum subsidence, horizontal displacement value.The reason of differences is originated from the basic theory of two methods.Those two methods are come from two different sects of rock mechanics.

PIM is based on discontinuous particle medium mechanics theory.The theory holds that rock is constituted by sand or small rock particle.The particles do not have any relationship with each other and can have relative movement.It assumes that the movement of a great quantity of granular medium is a stochastic process.For the center of the excavation panel, it has little stress between rock bedding.PIM theory coincides with the stochastic medium hypothesis.But at the position of panel boundary, actually the rock is subjected the force of tensile stress to the center of the panel.With this force of tensile stress, the rock above or outside the panel boundary could have horizontal displacement and subsidence movement toward the panel center.For this reason, the calculation deformation results outside the panel are smaller than measured data and the subsidence basin area is also smaller than measured area.

FLAC3Dis a fast Lagrangian analysis of continua software.The software solves rock mechanical behavior with explicit finite difference.To make a good simulation model, it needs complex geological exploration and mechanical tests.FLAC3Dis good at solving the problem of continuous medium.But rock body deformation is an extremely complicated problem.Rock body contains different material, soft surface, fault and cracks.To simplify the coal mine model to be a layered continuous rock mass need to ignore the discontinuous property of the actual rock.It will reduce the maximum subsidence value and enlarge the subsidence basin area.

As shown in Fig.8 and Fig.11, both results of PIM and FLAC3Dhave its own superiority and defects.PIM calculate the subsidence depends on the actual measured parameters,the deformation results of ground surface is better.But FLAC3Dcan calculate the dynamic change of the rock, it also can provide more information on time effect and space effect of mining rock deformation process, it is very important to safety evaluation,predict and analysis [Xie, Zhou, Wang et al.(1999)].

4.4 Discussion

There are four main categories of mining subsidence prediction method based on which theories and functions are used.They are phenomenological method, rock mechanical method, physical simulation method and numerical simulation method [Cui and Deng(2017)].PIM is one of the influence functions in phenomenological methods.PIM is mostly used in China and it has produced huge economic benefits to reduce the coal mine subsidence damage.There are a lot of excellent reasons to use PIM.

1.PIM is based on stochastic medium theory which is experimentally verified [He,Yang, Ling et al.(1991)].

2.PIM just needs 8 parameters which can be obtained from measured data easily.The parameters such as subsidence factor are very easy to be understood by engineers in the coal mine.

3.The most important is reliable surface subsidence can be obtained by using PIM.For the parameters mostly measured, the calculated subsidence can fit the actual surface motion perfectly.

4.Easily programmable with PC.This is good for promoting its application.

But PIM also has some defects.

1.The phenomenological methods are built based on the geometry of surface deformation.It cannot reveal the mechanical properties of rock mass.

2.The function can calculate the maximum subsidence area well but not so exactly at the boundary of the subsidence because of the convergence fast function, calculated results can not fit very well in the boundary of subsidence.

3.PIM is only good at calculate in the subsidence caused by horizontal or gently inclined coal seams but not so good for highly inclined coal seams.

4.PIM can not calculate the strain and stress of rock mass.It difficult to evaluate the underground constructions deformation.

Numerical simulation methods like FLAC3Dalso have some benefits.

1.The simulation model can provide information of movement, strain and stress of rock mass which are very important indicators of deformation.

2.The model is based on the basic rock mechanics.Assumptions and previous observations are not required.

3.Simulation also can provide dynamic results for the future.It is an easy way to obtain clear information about the process of rock mechanical change.

But the simulation now has some visible defects.

1.The rock and environment have a lot of complex influence factors; it difficult to get enough information to build the model.A simplified model can lead to a fatal error.It needs a lot of experience and knowledge to build a reliable model.

2.It is difficult to get reliable rock property parameters.A lot of costly rock test and experimental verification is needed in rock simulation.

3.Analysis simulation results require great familiarity with rock theory.Sometimes the results are not so reliable with simplified models and it is difficult to be analyzed by engineers.In China, a lot of coal mine engineers doubt the reliability of simulated results.

With the results and analysis above, PIM method can calculate the surface subsidence better with measured parameters.PIM is a phenomenology theory based on the geometry integral function.It cannot reflect the dynamic change process of rock mass under redisturbance of coal excavation.Numerical simulation is a prototype of actual excavation;it can simulate the whole process of initial stress of rock, coal excavation disturbance,dynamic deformation and new balance.And can provide more information like strain,stress and displacement of inner rock.But the numerical simulation also based on mechanical model, rock property parameters is complicated and difficult to be obtained.It also a problem to get a reliable results of surface displacement with FLAC3D, a unified model combines with stope-rock movement-surface subsidence is needed in the future[Zuo, Sun, Wen et al.(2018); Cui and Miao (2000)].But currently both PIM and FLAC3Dshould be used to evaluate the deformation influence of coal excavation especially for underground facilities.And we recommend that more simulation method should be used with the development of dynamic rock model and numerical simulation technology in the future.

5 Conclusion

(1) In order to excavate the shaft pillar, PIM and FLAC3Dsimulation were used to calculate the surface subsidence and influence on the shafts.The PIM results show that the panel 0091 excavation has no influence on the shafts.And the simulated results of FLAC3Dshows that the subsidence of main shaft is 9 mm and auxiliary shaft is 18 mm.FLAC3Dresults have bigger influence area than PIM calculated.With the analysis of shaft and rock displacement and stress, it shows that the shaft can keep safe after 0091 panel excavated.

(2) Comparison shows that the PIM fit better to the measured surface subsidence than FLAC3D.Mining subsidence influence area of PIM results are smaller than measured but the FLAC3Dresults are bigger.PIM is easier to be used and can predict the surface subsidence with 8 parameters; also, the surface deformation results are better than FLAC3D.But FLAC3Dcan provide more information like strain, stress, etc., of the surface and underground rock.It is useful for underground facilities evaluation.So, both PIM and numerical simulation method should be used in a coal mine subsidence evaluation especially for a project of underground deformation calculation and analysis.

(3) With development of rock mechanical theories, more simulation method should be used in the future

Acknowledgement:The work was supported by the Natural Science Foundation of China (Granted No.51404272) and China Scholarship Fund.And the authors would like to thank Kailuan (Group) Limited Liability Corporation for allowing them to write this paper and providing a lot of measured data.