A geometrically and locally adaptive remeshing method for finite difference modeling of mining-induced surface subsidence

Ziyu Zhang,Gang Mei,Nengxiong Xu

School of Engineering and Technology,China University of Geosciences(Beijing),Beijing,100083,China

Keywords:Mining-induced subsidence Numerical modeling Finite difference method(FDM)Distorted mesh Adaptive remeshing

ABSTRACT Surface subsidence induced by underground mining is a typical serious geohazard.Numerical approaches such as the discrete element method(DEM)and finite difference method(FDM)have been widely used to model and analyze mining-induced surface subsidence.However,the DEM is typically computationally expensive,and is not capable of analyzing large-scale problems,while the mesh distortion may occur in the FDM modeling of largely deformed surface subsidence.To address the above problems,this paper presents a geometrically and locally adaptive remeshing method for the FDM modeling of largely deformed surface subsidence induced by underground mining.The essential ideas behind the proposed method are as follows:(i)Geometrical features of elements(i.e.the mesh quality),rather than the calculation errors,are employed as the indicator for determining whether to conduct the remeshing;and(ii)Distorted meshes with multiple attributes,rather than those with only a single attribute,are locally regenerated.In the proposed method,the distorted meshes are first adaptively determined based on the mesh quality,and then removed from the original mesh model.The tetrahedral mesh in the distorted area is first regenerated,and then the physical field variables of old mesh are transferred to the new mesh.The numerical calculation process recovers when finishing the regeneration and transformation.To verify the effectiveness of the proposed method,the surface deformation of the Yanqianshan iron mine,Liaoning Province,China,is numerically investigated by utilizing the proposed method,and compared with the numerical results of the DEM modeling.Moreover,the proposed method is applied to predicting the surface subsidence in Anjialing No.1 Underground Mine,Shanxi Province,China.

1.Introduction

Coal seam mining could lead to deformation of surface and overburden movement,and further cause a series of geohazards,such as overburden fracturing,surface cracks,landslides and surface subsidence(He and Xu,2018).There are several methods to evaluate and predict the surface subsidence,including the pro file function(Asadi et al.,2005),in fluence function(Sheorey et al.,2000;Chen et al.,2016),and numerical modeling(Ghabraie et al.,2015).Numerical modeling methods have widely developed to investigate mining-induced surface subsidence,which can be roughly divided into two categories:(i)continuum-based methods,such as finite element method(FEM)(Sepehri et al.,2017)and finite difference method(FDM)(Alejano et al.,1999;Karmis and Agioutantis,1999;Xu et al.,2013),and(ii)discontinuum-based methods,such as discrete element method(DEM)(Xu et al.,2016;Lv et al.,2017;Regassa et al.,2018),meshless method(Wu et al.,2001;Kim et al.,2006;Fu and Jin,2015),and discontinuous deformation analysis(DDA)(Zuo et al.,2009).

The DEMis a discontinuum-based numerical method especially applicable to interpreting discontinuities in the jointed rock mass,and modeling the discontinuous rock movement and failure characteristics of rock masses.For example,Ren(2011)investigated the failure of the coal seam roof and the movement behavior of the overlying strata with the DEM and predicted the rule of surface subsidence.Xu et al.(2016)employed the DEMto reveal the failure characteristics of jointed rock mass in underground mining.Lee et al.(2016)used the DEM to study the spherical underground cavity affecting the ground surface subsidence.Regassa et al.(2018)proposed an equivalent simulation method to study the rock joint fissures in slope failure.However,when analyzing the large-scale computational model consisting of a large number of blocks,the DEM modeling is quite computationally expensive(Karrech et al.,2007;Karampinos et al.,2015;Xu et al.,2016;Regassa et al.,2018),and even cannot be well conducted in many cases.

The continuum-based numerical methods such as the FDM and FEM can also be used to evaluate and predict the mining-induced surface subsidence(Zhao et al.,2016;Zhou et al.,2018).Both the FDM and FEM are much more computationally ef ficient than the DEM.However,there is an inevitable problem when using the FDM or FEM to model the mining-induced surface subsidence,i.e.the mesh distortion(Wang et al.,2015;Shang and Ouyang,2018).

To address the mesh distortion in the numerical modeling when using the FEM or FDM,adaptive remeshing or re finement is commonly exploited.For example,Labergere et al.(2011)proposed a two-dimensional(2D)adaptive remeshing method to accurately simulate various metal forming processes in FDM.You et al.(2015)used an adaptive mesh regeneration scheme that can apply the mesh repartition according to the prede fined material threshold.Yang et al.(2018)employed a h-adaptive mesh re finement technique to deal with the distorted mesh in FEM modeling.

As introduced above,adaptive remeshing of the distorted mesh is an effective solution to dealing with the mesh distortion in FEM or FDMmodeling.However,to the best of the authors’knowledge,almost all related work focuses on the adaptive remeshing for the distorted mesh with only a single attribute(for example,for the distorted mesh model representing a single stratum),without consideration of the concurrent remeshing for those distorted meshes with multiple attributes(for example,for the distorted mesh model representing several strata).

In summary,when employing numerical methods such as DEM or FDM to investigate the mining-induced surface subsidence,several problems are remaining to be addressed.For complex and large geological models,the DEM modeling is quite computationally expensive,and even cannot be conducted.The adaptive remeshing is commonly used to deal with the mesh distortion in the FDMmodeling,but generally aims at the remeshing of distorted meshes with only a single attribute,without consideration of those distorted meshes with multiple attributes.

To address the above problems,this paper presents a geometrically and locally adaptive remeshing method for the FDM modeling of largely deformed surface subsidence induced by underground mining.The essential ideas behind the proposed method are as follows:(i)Geometrical features of elements(i.e.the mesh quality),rather than the calculation errors,are employed as the indicator for determining whether to conduct the remeshing;and(ii)Distorted meshes with multiple attributes,rather than those with only a single attribute,are locally regenerated.

In the proposed method,the distorted meshes are first adaptively determined based on the mesh quality,and then removed from the original mesh model.The tetrahedral mesh in the distorted area is first regenerated,and then the physical field variables of old mesh are transferred to the new mesh.Numerical calculation process recovers when finishing the regeneration and transformation.

To verify the effectiveness of the proposed method,the surface deformation of the Yanqianshan iron mine,Liaoning Province,China,is investigated using the proposed method,and compared with the numerical results of the DEM modeling.Moreover,the proposed method is applied to predicting the surface subsidence in Anjialing No.1 Underground Mine,Shanxi Province,China.

The main contributions of this paper can be summarized as follows:(i)A geometrically and locally adaptive remeshing method speci fically for distorted mesh with multiple attributes is proposed;and(ii)The proposed adaptive remeshing method is utilized in the FDM modeling to numerically investigate the largely deformed surface subsidence induced by underground mining.

2.Methodology

2.1.Overview

Surface subsidence caused by underground mining is a gradual development process.With the continuous advancement of the working face,the overburden is in a dynamic change of stabilityinstability-stability.After the coal seam is mined out,the overlying strata may cause a large displacement.DEM and FDM are usually used to predict surface subsidence.However,for the largescale and complex model,the DEMis considerable computationally costly or may fail.The FDM cannot avoid the mesh distortion problemwhen simulating the large deformation.The appearance of distorted meshes will cause the non-convergence of numerical calculations.

To address the above problems,in this paper,a geometrically and locally adaptive remeshing method is proposed for the FDM modeling of largely deformed surface subsidence induced by underground mining.The main objective of the proposed geometrically and locally adaptive remeshing method is to regenerate highquality computational meshes and avoid distorted meshes in the numerical modeling process.

The work flow of the FDM modeling method for analyzing mining-induced surface subsidence by employing the proposed geometrically and locally adaptive remeshing is illustrated in Fig.1.The proposed geometrically and locally adaptive remeshing method is mainly composed of six sub-processes:

(1)In the iterative FDM calculation process,the distorted meshes whose quality is less than the prede fined threshold are searched(see Section 2.2).

(2)The distorted area is determined by finding the third-order neighboring element of the distorted meshes(see Section 2.3.1).

(3)The distorted area is removed from the original computational model(see Section 2.3.2).

Fig.1.Work flow of the FDM modeling of mining-induced surface subsidence by employing the proposed geometrically and locally adaptive remeshing method.

(4)The surface triangular mesh of the distorted area is obtained(see Section 2.3.3).

(5)New meshes are regenerated and merged with the remaining mesh(see Section 2.3.4).

(6)The physical field variables are transferred from the old meshes to the new meshes(see Section 2.4).

2.2.Determination of the mesh quality

This section summarizes four scale-invariant quality measures for a tetrahedron in FDM modeling.There are several properties that the quality measure should have:

(1)All degenerate elements have a quality of zero.

(2)The measure is scale-invariant;in other words,two elements of different sizes and identical shapes have the same quality.

(3)The measure is normalized to achieve a maximum value of one,which means an equilateral tetrahedron has a quality of one.

(4)All inverted elements have a negative quality.

In practice,we consider that an element is distorted if mesh quality is less than a speci fic threshold.The distorted mesh is evaluated according to built-in functions in FLAC3D software based on several quality metrics such aszqualitytest(pz, 1),zqualitytest(pz,2)andzaspect(pz).In those functions,pzmeans the index of the element.

More speci fically,zqualitytest(pz,1)andzqualitytest(pz,2)indicate a test of volume over edge length and a test of skew,respectively(Eq.(1)and(2)).zaspect(pz)returns a measure of the aspect ratio of a zone(Eq.(3)).The aspect ratio of the tetrahedron also lies between 0 and 1,and the larger aspect ratio implies the better quality of the element.

whereVis the volume of a tetrahedron,Llongestis the longest edge of a tetrahedron,Voptis the same tetrahedral volume of the circumscribed sphere as a tetrahedron,R0is the radius of the inscribed sphere of a tetrahedron,andRiis the radius of the circumscribed sphere of a tetrahedron.Fig.2 simply shows the inscribed and circumscribed spheres of a tetrahedron.

2.3.Process of the proposed adaptive remeshing method

A simpli fied process of the proposed geometrically and locally adaptive remeshing method is illustrated in Fig.3.In the iterative calculation,the distorted mesh is selected by continuously judging the quality of the mesh,and then the distorted area is determined.The tetrahedral meshes in the distorted area are remeshing,and the new meshes are merged with the remaining mesh.The majorpurpose is to detect,remove and regenerate distorted meshes in the FDM modeling of mining-induced surface subsidence.

2.3.1.Obtainingdistortedarea

As explained in Section 2.2,it needs to de fine the distorted area when the distorted mesh has been identi fied.Considering the data structure type of computational models in the FDM modeling,the third-order neighboring elements of the distorted mesh are selected as the distorted area in this study.

The judging principle of the neighboring element is as follows.The element that shares the nodes with elementAis the first-order neighboring element ofA.The calculation of the 2D neighboring mesh is the same as that of the three-dimensional(3D)neighboring mesh.In order to show the selection method of the distorted area conveniently,we take a 2D triangular mesh as an example.As shown in Fig.3c,the red element is the distorted element,and the green,orange and blue elements constitute the third-order neighboring elements of the distorted mesh,which is the distorted area expected.

2.3.2.Distortedarearemoval

According to Section 2.3.1,the mesh can be divided into two parts;one is the distorted area(Fig.3c)which needs to be removed from the mesh,and the other is the remaining mesh(Fig.3g)which needs to be kept.We use the method that obtains the surface triangular mesh based on the distorted area in Section 2.3.3.

2.3.3.Acquisitionofthesurfacetriangularmeshes

After removing the distortion area,it needs to group the meshes according to the attributes of the mesh in the distorted region,as shown in Fig.3d.Then,the boundary elements are obtained according to different groups.As shown in Fig.3e,the boundary meshes of the 2D triangular model are boundary edges.While in the 3D tetrahedral model,the surface mesh is a collection of triangles.When traversing all the triangles that make up the tetrahedron,if there is a triangular face that is not shared by other elements,then this triangular face is the surface of the model.After traversing the triangular surfaces of all tetrahedral meshes,we can obtain the surface mesh for the entire remeshing region.

Fig.3.An illustration of the adaptive remeshing in two dimensions for meshes with multiple attributes.

2.3.4.Tetrahedralmeshregenerationandmerging

In the 2D remeshing procedure,the remaining mesh(Fig.3g)needs to be combined with high-quality meshes.That is,it needs to remesh the distorted area according to the surface of the distorted area(Fig.3e).As shown in Fig.3f,the quality of the regenerated mesh is higher than before by the mesh generation schemes,which also improves both the average mesh quality and the extremal mesh quality of the local mesh.

In this work,we use the open-source tetrahedral mesh generator-TetGen developed by Si(2008)to regenerate high-quality meshes.TetGen is a program to generate tetrahedral meshes from 3D polyhedral domains.Its objective is to generate high-quality tetrahedral meshes suitable for various applications in scienti fic computing.We invoke the commands in TetGen to regenerate the high-quality tetrahedral meshes and then use the meshes to replace the tetrahedral meshes in the distorted area.This could be used to improve the average quality of the tetrahedral mesh and the worst-quality element in the model.

After the locally adaptive remeshing,the new meshes need to be merged with the remaining mesh based on the topology relationship(Fig.3g).

We use a generic,simple and practical algorithm,MeshCleaner(Mei et al.,2018),to merge the meshes.The algorithm is composed of the stage of compacting and reordering nodes and the stage of updating mesh topology.More details about the mesh merge can be found in Mei et al.(2018).

2.4.Transfer of the physical field variables

After merging the mesh,the position and topological relationship of the computational mesh have changed.As a consequence,the original physical field variables are no longer applicable to the new mesh.The physical field variables must be transferred fromthe old mesh to the new mesh.In the FLAC3D software,the nodal variables include velocity and displacement variables,and element variables include stress variables.Field transfer means that the velocity and displacement variables of the old nodes are transferred to the new nodes,and the stress variables of the old elements are transferred to the new elements.

Physical field transformation has been widely investigated(Peric et al.,1996,1999;Mediavilla et al.,2006;De Boer et al.,2007).In our work,we use the field transfer method by speci fically utilizing interpolation algorithms,and cooperate with some functions for output and input physical field variables in the FLAC3D software.

2.4.1.Transferofthenodeandelementvariables

The nodal and element variables of the new mesh are created by the moving least squares(MLS)interpolation algorithm(Ding et al.,2018)using the nodal and element variables of the old mesh as information data.

In the original MLS method(Lancaster and Salkauskas,1981),the following test function is used as the approximation function of a functionf(x):

wherefh(x)is the global approximation function of a functionf(x),PT(x)=[p1(x),p2(x),…,pm(x)],mis the number of basis,pi(x)is the basis function,andai(x)is the corresponding coef ficient for the basis functions.

Corresponding to the global approximation of Eq.(5),the local approximation in the neighborhood of a pointxis de fined as

In Eq.(6),xis the point to be fitted,is the point in the neighborhood of the pointx,and the coef ficientsai(x)are determined using the weighted least squares method by minimizing the local approximation error of the functionf(x):

whereJis the local approximation functionat all nodes(x=xk)weighted square sum of the error,w(x-xk)is the compact-supported weight function,andxkis the node in the compact-supported region of pointx(Ding et al.,2018).

The transfer of element physical field variables is composed of three steps:

(1)Step 1:The volume-weighted averaging method is used to transfer the stress of the old meshes to the old nodes,and the equivalent physical field variables of the old node are obtained.The equivalent physics field variables of the old nodes can be obtained by

(2)Step 2:The equivalent element physical field variables of the new node are obtained by interpolation(Ding et al.,2018)of the known equivalent element physical field variables of the old node.

(3)Step 3:By calculating the physical field variables of the equivalent elements of the new node,the variables are assigned to the new mesh.The average of the equivalent element physical field variables of the four nodes is the new physical field variables of the tetrahedral element,which can be written by

The node field transfer method uses the second and third steps of the element field variables transfer method.

2.4.2.Transferofthestatevariables

In the FLAC3Dsoftware,the zone state is represented by the sum of multiple integers.For example,“Failure in shear now”is represented by the indicator of“1”;“Failure in tension now”is represented by the indicator of“2”;“Failure in shear in the past”is represented by the indicator of“4”;“Failure in tension in the past”is represented by the indicator of“8”.Taking into account the zoning of the state of the plastic zone,in most cases,the zone states are similar in a certain local region.In the process of state variables interpolation,if a new node is inside a certain old element,the equivalent element state variable of the new node is the same as the state variable of the old element.If there are two or more new nodes in any new element with the same equivalent element state variable,then the equivalent element state variable of the element is the state variable of the new element composed of these four new nodes.If the equivalent element state variables of the four nodes are different,the instability state field variables of the four nodes are selected as the state variables of the new element for safety.

3.Numerical examples for veri fication

To verify the effectiveness of the proposed method,the surface deformation of Yanqianshan iron mining is investigated and compared.More speci fically,the roof subsidence and failure of the Yanqianshan iron mining is calculated first used the geometrically and locally adaptive remeshing method,and the calculated results are compared with those of Regassa et al.(2018)using the 3DEC software.

Fig.4.A simple illustration of the equivalent physical field variables of the node ni.

Fig.5.A simple illustration of the equivalent physical field variables of the element ei.

3.1.Computational models and parameters

In this work,the geological modeling software RockModel(Xu and Tian,2009)is used to establish a geological model of the Yanqianshan iron mine(Fig.6).

The geological model is partitioned into a tetrahedral mesh composed of tetrahedrons with the size of 25 m for the north and south subareas and 15 m for the middle subarea in size.The tetrahedral geological model is imported into FLAC3D software as the numerical computational model.The numbers of nodes and elements in the computational model are 30,963 and 165,804,respectively.Computational parameters are listed in Table 1.

Table 1Mechanical parameters of the rock mass for Yanqianshan iron mine.

3.2.Numerical results

In the FDM modeling,the adaptive remeshing method can be used to determine the distorted mesh based on the geometrical feature of elements in the iterative calculation process,and regenerate the new tetrahedral mesh in the distorted area,and then transfer the new physical field variables to new meshes.For the computational model,the quality of the mesh is checked every 2000 iterations.The appearance of distorted mesh will lead to the cessation of the numerical calculation and the implementation of adaptive remeshing.Numerical calculation recovers when finishing the regeneration and transformation.

Fig.6.Computational models of Yanqianshan iron mine:(a)Model in 3DEC(Regassa et al.,2018),and(b)Model in FLAC3D.

In the numerical calculation process of the mining of the first layer,the quality of the mesh is maintained well,and there is no distorted mesh.When the mining of the first layer is completed,there is no large displacement or deformation at this level of the mining(Fig.7a),and the rock mass overlying the mining area is relatively stable.However,the roof tends to deform towards the center of the goaf under the action of gravity(Fig.7b).

During the mining of the second layer,the overlying strata further collapse,which lead to the extension of part of the computational mesh.Thus,the mesh quality decrease signi ficantly,and a certain number of distorted meshes appear,as shown in Fig.8.Then,the locally and geometrically adaptive remeshing method have been used to improve the mesh quality.In Fig.9,the average element quality and the worst mesh quality are improved after the remeshing method.After the remeshing step is executed,the physical field variables at the nodes and elements of the new mesh are then obtained by the interpolation algorithm(Ding et al.,2018)(Fig.10).Although the partial mesh topology is changed,the new displacement field is still identical to the original displacement field which verify the effectiveness of the field transfer method.After the new mesh has the physical field variables,the numerical calculation can be resumed until convergence.

When modeling the mining of the third layer,the locally and geometrically remeshing method is also used to calculate the displacement and stress fields.After mining the third layer,the overlying rock mass in the goaf fails progressively in an upward direction.At the same time,the shear stress is relatively large at the boundary of the mining area,which is easy to produce cracks and dislocation,and the mining landslides may exist at the mining area boundary.

In some related studies,the researchers used shear strain variables to assess surrounding rock stability and damage during underground excavation(Hajiabdolmajid et al.,2002;Zhu et al.,2005;Chang et al.,2007;Corkum and Martin,2007;Huang and Xiao,2010).For example,Zhang et al.(2019)and Darvishi et al.(2020)analyzed the deformation of the goaf roof,surrounding rocks and surface by calculating the distribution of the plastic zone and the shear strain field in FLAC3D software.In this work,we analyze the surface subsidence after the mining of the third layer in the above two aspects.It can be seen from Fig.11a that when the mining of the third layer is completed,the surface cover is affected by gravity,which cause tensile and shear damage at the mining boundary,making it susceptible to the formation of surface cracks and irregular subsidence.Although the FDM is a continuous medium method and cannot produce discontinuous deformation,we can predict that large differential damage such as landslides or sinkholes and other undesirable geological phenomena will occur in the red area in Fig.11b.By comparing the numerical results of FDMand DEM,i.e.the location of the plastic zone in Fig.11a,the region with a large increment of shear strain in Fig.11b(red part)and the location of the collapse crater in Fig.11c is the same.Meanwhile,the field observation reveals a sinkhole at this location(Fig.11c).Therefore,it can be concluded that the proposed method for solving the large deformation problem caused by underground mining in FDM is correct and applicable.

4.Application case

In this section,we present an application case of the proposed geometric adaptive local remeshing method.More speci fically,thismethod is used to analyze the subsidence coef ficient and surface displacement of the mining area under the condition of mesh distortion caused by underground mining.

Fig.7.(a)Stratum movement(in m)and(b)stress contour(in Pa)after mining of the first layer.

Fig.8.Appearance of the distorted meshes during mining the second layer.

4.1.Geology background

Fig.9.Comparison of the mesh quality before and after remeshing when modeling the mining of the third layer.

The study area of the Anjialing No.1 underground mine is located in Shuozhou City,Shanxi Province,China(Fig.12).In the coal mining area,the average thickness of the coal seam is 12 m,with the maximum and minimum thicknesses of 18.6 m and 5.4 m,respectively,and the mining depth is approximately 120-240 m.

4.2.Computational model and parameters

4.2.1.Computationalmodel

By exploiting the contour map of the mining area of Panel 4106,a 3Dmodel is established.As shown in Fig.13,the positive direction of theX-axis of the model indicates the east direction,the positive direction of theY-axis indicates the north direction,and theZ-axis indicates the altitude.The plane boundary of the computational model can be simply obtained by extending the boundary of the mining area 1000 m along the west,east,south and north directions,respectively.

If the computational mesh model is established by equal-sized tetrahedral mesh,the number of elements will be large and the calculation will have a lower ef ficiency.To reduce the number of elements,smaller tetrahedral meshes are used near the mining area and larger tetrahedral meshes are used elsewhere.The total number of elements is 771,340 and the number of nodes is 137,581.

In this case,the displacement boundary condition is employed.The boundary conditions of the model are set as follows:(i)Applying the fixed constraint of theX-direction atX=11,639 m and 15,865 m to the eastern and western boundaries of the model,(ii)Applying the fixed constraint of theY-direction atY=7170 m and 9470 m,and(iii)Applying the fixed constraint of theZ-direction atZ=812 m.In addition,no constraints are added to the top surface of the model.The acceleration of gravity is set to 9.8 m/s2,and the maximum unbalanced force ratio is 1×10-5.The initial stress field of the model is calculated under the condition of self-weight,and then displacement,velocity fields and state zone are reset.

Fig.10.Comparison of the field transformation results(a)before and(b)after remeshing(displacement in m).

Fig.12.Location of the Anjialing No.1 underground mine in China.

As described in the User’s Guide of FLAC3D software(Itasca,2007),the interface elements are the connection elements between sub-grids,which can be used to simulate the separation during the calculation process.An interface can represent a physical discontinuity such as a fault(Kim and Larson,2019),contact plane(Zhao and Zhang,2017),or interface between two different materials(Bergado and Teerawattanasuk,2008;Zhang et al.,2015;Yang et al.,2017).In this study,we set the interface elements between the roof and the coal seam,and the interface elements belong to the coal seam roof(Fig.14).The above method is used to model the contact between roof and floor as the roof collapses after excavation.

Fig.11.Numerically calculated failure area and deformation of the roof:(a)Plasticity state of zones in the FDMmodeling when using the adaptive remeshing method,(b)Maximum shear strain increment in the FDM modeling when using the adaptive remeshing method,(c)The uniform vertical subsidence(in m)of the roof simulated by 3DEC(Regassa et al.,2018),and(d)The observed sinkhole(Regassa et al.,2018).

Fig.13.Computational tetrahedral mesh model of the study area.

4.2.2.Computationalparameters

In the FDM modeling of underground mining,the Mohr-Coulomb failure criterion is adopted.Based on a detailed in situ investigation and laboratory mechanical experiments on the intact rocks,the physico-mechanical parameters of the rock masses are determined,which are summarized in Table 2.

Table 2Physico-mechanical parameters of soil and rock masses used in the FDMmodeling.

4.3.Process of numerical modeling

Before the mining activities,the rock mass is in a state of equilibrium under the effect of the geostatic stress field.After the coal seam is extracted,the goaf is formed in the rock mass,which causes the stress state of surrounding rock mass to change and redistribute,thus causing the rock mass to move and deform until reaching a new balance.With the advance of working face,this process is repeated.It is a complex process of physico-mechanical change,as well as the process of stratum movement and failure,which is called stratum movement.

In the initial mining process,starting from the direct roof,along the normal direction of the bedding surface,the overburden of the goaf bends towards the goaf direction until the surface.In the bending range,the strata do not fall off and basically maintain the layered structure.When the working face is advanced to～200 m,the coal seam begins to thicken.The overburden of the goaf deforms and hence moves downward under the effect of gravity.When the applied stress exceeds the strength of the overlying rock strata,the immediate roof will fall off and fill the goaf.The overburden experiences a cyclic stable-unstable-stable-unstable process.With the advance of working face,the above process occurs repeatedly.

Fig.14.A simple illustration of the interface between the roof and coal seam.

In the numerical calculation of underground mining,we use the mesh quality metric to search for distorted meshes.Overall,the numerical calculation process can be divided into two categories:mesh distortion and non-distortion.We apply the geometrically and locally adaptive remeshing method to settle distorted mesh problem in the numerical calculation.When the mesh quality is less than the threshold,the mesh should be considered as a distorted mesh.After the remeshing procedure is completed,the next step of excavation is continued.In the remeshing procedure,the command of computational mesh quality needs to be set in the iterative calculation.In this study,it is set to repeat every 2000 steps to compute the mesh quality.If there is a distorted mesh,the calculation needs to be stopped for the remeshing;otherwise the calculation is continued until the excavation is completed.

4.3.1.Modelingoftheexcavationof0-100m

The deformation of the calculation mesh is strongly relevant to the deformation of the overlying rock mass in the goaf.We set the cross-sectionI-Ito observe the caving process of the coal seam roof(Fig.15).Fig.15a illustrates that starting from the initial mining face,the pressure arch is formed in the overlying rock which is supported by the surrounding rocks at the goaf boundary.In addition,the structure of the overburden stays stable and the pressure arch can bear the weight of the overlying rocks,and the deformation of the overburden has little effect on the mesh.Therefore,there is no distorted mesh and all computational mesh quality satis fies the metric requirements during the calculation process.

4.3.2.Modelingoftheexcavationof100-200m

Because of the thick coal seam,the tensile failure of mesh occurs in the calculation model when the working face being advanced to 200 m,and the tetrahedron elements are greatly deformed and stretched to be slender.In Fig.16,the roof caving causes deformation of the mesh,and then several distorted meshes occur at the same time.As a result,we magnify those distorted meshes and then observe that these meshes become slender,which are not suitable for numerical computation.

Fig.17a shows the third-order neighbor meshes of the distorted meshes,and the distorted area contains two topologically disconnected areas.In Fig.17b,blue meshes represent the distorted meshes,and the green and red meshes represent attributesAandB,respectively.As a result,such distorted meshes are impossible to simulate the subsequent mining,and the numerical calculation will be terminated.

It can be seen from Fig.18 that the remeshing method eliminates the distorted meshes in the zone where requires remeshing,and improves the average mesh quality as well.The lowest mesh quality is 0.272 before remeshing,and 0.308 after remeshing.The average mesh quality increases from 0.594 to 0.605 after remeshing.

Then,we use the regenerated mesh to resume the calculation.Fig.19 shows the magnitude of displacement on cross-sectionI-Ibefore and after remeshing.Comparing the black frame areas in Fig.19b and c,the local remeshing process of the mesh changes the local topology of the original mesh,but it does not alter the physical field variables of the original mesh.

Fig.15.Illustration of the collapsing of coal seam on cross-section I-I when the working face advanced to(a)100 m and(b)200 m.

4.4.Prediction of mining-induced surface subsidence and ground movements

When the mining is completed,we observe that the maximum surface subsidence is～10.47 m(Fig.20),and the subsidence coef ficient is～0.87,which is similar to that of the mining area in Shanxi Province,China.According to Pillar Design and Mining Regulations under Buildings,Water,Rails and Major Roadways(in China)(National Bureau of Coal Industry,2000),the subsidence coef ficients are 0.8 for Panel 22108 in Xiqu coal mine,0.85 for Panel 12101 in Zhenchengdi coal mine,0.8 for Panel 32903 in Ximing coal mine,and 0.83 for Yangquan coal mine.Therefore,it is clear that the geometrically and locally adaptive remeshing method proposed in the study can deal with the problem of mesh distortion caused by underground mining.

5.Discussion

5.1.Advantages

In this study,the severe subsidence damage caused by the mining of underground coal seams can be effectively simulated by the proposed method.The DEM based on block theory is usually used to simulate discontinuous movements and large deformation processes of rock masses.However,for complex and large-scale geological models,DEM is time-consuming.Instead,the method proposed in this study is based on FDM,which is more ef ficient than DEM.

The proposed method can solve the problem of mesh distortion in the FDM calculation process.In view of the mesh distortion phenomenon that occurs in the FDM,a geometrically and locally adaptive remeshing method is proposed to settle the distorted mesh problem to update the mesh that satis fies the quality requirements of FDM modeling.

Fig.16.Distorted tetrahedral meshes occurred when working face is advanced to 200 m.There are six distorted tetrahedral meshes in this figure.

Fig.17.(a)Third-order neighboring elements of the distorted meshes,and(b)A view of internal element when cutting the distorted area.

In this method,the distorted area is determined by the mesh quality metric and then is removed from the original mesh.In the locally remeshing method,the mesh is updated after transferring the physical field variables of the original mesh to the regenerated mesh.Compared with updating the entire computational model,the proposed method reduces the frequency and scale of remeshing and retains the vast majority of the original nodes,thereby reducing the errors caused by interpolation.

Our method can accurately determine the distorted mesh by the mesh quality,remove the distorted area from the original mesh,and regenerate the updated mesh,which can avoid the irreversible numerical calculation termination caused by mesh distortion.

5.2.Disadvantages

Fig.18.Comparison of the mesh quality before and after remeshing when working face is advanced to 200 m.

Fig.19.Comparison of the magnitude of displacements before and after remeshing:(a)The magnitude of displacements on cross-section I-I,(b)Displacements calculated before remeshing,and(c)Displacements interpolated after remeshing.

In the proposed method,the selection of the distorted mesh should be based on the actual deformation of the mesh during the excavation calculation process,and the current mesh quality cannot re flect the dynamic change of the mesh.However,we used a tentative threshold based on experience.We did not calculate the actual value of deformation in mesh quality based on other reliable geometry theories.In addition,we did not consider the complicated hydrogeological conditions and fault development extent in the numerical simulation process of underground mining.As a result,the model was relatively simpli fied.Although we used a high-precision interpolation method,the interpolation results of the plastic zone were different from the actual results,interpolating multiple times and choosing the more precise interpolation algorithm can reduce the error.The 3D adaptive methodology is currently used in FDM with good results,but it is uncertain about the effect in other simulation methods.

5.3.Future work

In the proposed method,the error mainly comes from the transfer process of the physical field variables.The following solutions can be conducted for this problem.First,interpolating multiple times and choosing a more precise interpolation algorithm can reduce the error caused by interpolation.Second,combining the continuous medium method and discontinuous medium method,the meshless method is used at the large deformation region,and the FDMis used for the stable part.This method not only guarantees the accuracy of the calculation but also reduces the amount of calculation.

Fig.20.Investigation of the mining-induced surface subsidence in the study area.

Nowadays,various new techniques are still being proposed to solve the mesh distortion,such as the new spline element method(Chen et al.,2010),the smoothed FEM(Liu et al.,2007;Leonetti et al.,2016),and the overlapping FEM(Bathe and Zhang,2017;Zhang and Bathe,2017).We can learn from these new techniques to avoid mesh distortion during large deformation simulation.

6.Conclusions

In this study,a geometrically and locally adaptive remeshing method is proposed for FDMmodeling of largely deformed surface subsidence induced by underground mining.The essential ideas behind the proposed method are as follows:(i)The geometric characteristics of the elements(i.e.the mesh quality),rather than the calculation errors,are employed as an indicator to determine whether to conduct the remeshing;and(ii)Distorted meshes with multiple attributes,rather than those with only a single attribute,are locally regenerated.

In the proposed method,the distorted meshes are adaptively determined and removed.The distorted mesh is regenerated,and the physical field variables of the old mesh are transferred to the new mesh.The numerical calculation recovers after the regeneration and transformation.

The effectiveness of the proposed method has been veri fied by comparing the surface deformation for the Yanqianshan iron mine using the proposed method and DEM.The results indicate that the geometrically and locally adaptive remeshing method can improve both the average mesh quality and the extremal mesh quality of the distorted mesh during the calculation process,and the proposed method has been successfully applied to a surface stability analysis induced by deep underground mining and is thus suitable for addressing the mesh distortion problem induced by large displacement during the mining process.Future work is planned to achieve more accurate results by combining the continuum-and discontinuum-based methods.

Declarationofcompetinginterest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to in fluence the work reported in this paper.

Acknowledgments

This work was jointly supported by the National Natural Science Foundation of China(Grant Nos.11602235 and 41772326)and the Fundamental Research Funds for the Central Universities of China(Grant No.2652018091).