Modelling heterogeneous coal-rock (HCR) failure patterns under dynamic impact loads using image-based finite element (FE) and discrete element (DE) model

Journal Pre-proof Modelling heterogeneous coal-rock (HCR) failure patterns under dynamic impact loads using image-based finite element (FE) and discrete element (DE) model Kehong Zheng, Bingjing Qiu, Zhenyu Wang, Jianping Li, Kuidong Gao PII:

S0032-5910(19)30826-5

DOI:

https://doi.org/10.1016/j.powtec.2019.10.036

Reference:

PTEC 14798

To appear in:

Powder Technology

Received Date: 3 July 2019 Revised Date:

19 August 2019

Accepted Date: 1 October 2019

Please cite this article as: K. Zheng, B. Qiu, Z. Wang, J. Li, K. Gao, Modelling heterogeneous coal-rock (HCR) failure patterns under dynamic impact loads using image-based finite element (FE) and discrete element (DE) model, Powder Technology (2019), doi: https://doi.org/10.1016/j.powtec.2019.10.036. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Modelling Heterogeneous Coal-rock (HCR) failure patterns under dynamic impact loads using image-based finite element (FE) and discrete element (DE) model Kehong Zhenga, Bingjing Qiua, Zhenyu Wanga*, Jianping Lib, Kuidong Gaoc a b

College of Civil Engineering and Architecture, Zhejiang University, 310058, China

College of Mechanical and Electrical Engineering, China University of Mining and Technology, 221116, China

c

College of Mechanical & Electrical Engineering, Shandong University of Science & Technology, Qingdao, Shandong, China;

*

Corresponding author at: College of Civil Engineering and Architecture, Zhejiang

University, Yuhangtang Road 866, Hangzhou 310058, Zhejiang, China. Email address: [email protected] Abstract: The paper described two different numerical 3D mesoscopic approaches (the meso-scale image-based FE and DE model) for quantitative and qualitative characterization of the dynamic impact mechanical behavior of HCR. The 3D meso-structure of HCR in calculations, which was directly obtained from real coal rock specimen using non-destructive XCT techniques, was modelled as a random heterogeneous two-phase material composed of coal and gangue respectively. Dynamic impact loading in three directions (x-, y-, and z-axis) were modelled to investigate the 3D meso-structure and different input energies effects on failure patterns of HCR. The very different qualitative and quantitative results demonstrated that the spatial distribution of mineral phases, loading directions and input energies have pronounced influenced on material strength as well as failure patterns of heterogeneous rock materials. The XCT image-based numerical model proved to be an effective tool that gives insights into the meso-deformation mechanisms of HCR undergoing dynamic impact failure behavior. Keywords: Heterogeneous coal-rock (HCR); Impact loads; X-Ray computed tomography (XCT); Image-based modelling; Failure pattern. 1

1. Introduction In coal mining process, coal rock is a well-known composite material wherein different phases may be separated: coal, gangue and cracks [1-4]. Each of mineral phases has a different set of physical and mechanical parameters (e.g., volume fractions, sizes, shapes and spatial distribution of mineral phases). The failure patterns of HCRs at macroscale level are strongly influenced by the key meso-structure parameters of different mineral phases and their various interactions. However, modelling fracture behavior of HCR at the heterogeneous material scale is difficult due to the damage and failure behaviors of HCR involve complex physical and mechanical behaviors at different scales. In order to properly describe the failure patterns of HCRs in detail, material mseo-structure should be taken into account since its effects on the macro-mechanical responses are pronounced [5-8]. The good understanding of 3D meso-structure effects on the failure patterns is of major importance to help reveal the failure conditions and prevention mechanism of coal rock dynamic disaster. Modelling coal rock fracture behavior at meso-scale provides a practical way to quantitative characterize their damage and failure behaviors under dynamic impact loads through computational simulations. Recently, the non-destructive X-Ray computed Tomography (XCT) technique is now widely used to characterize the internal micro/meso-structures of various materials, such as construction materials (e.g. concrete [9-11], asphalt [12-13] and sand [14-15]) and geological materials (e.g. copper oxide ore [16], iron ore [17], and coke [18]). Many XCT studies

acquire the fracture evolution characteristics of

heterogeneous materials with quasi-static loading conditions (i.e. in-situ X-ray computed tomography), or study damaged materials after loading [15,19-20]. Due to the impact loading process always occurs within a few microseconds, the continuous investigations of the entire fracture process with the current XCT techniques is not workable for technical reasons yet. To examine the meso-structure damage and fracture evolution of HCR in detail, the mainly two approaches were typically used, namely the user-defined models [21-22] and the image-based model [14]. The 3D 2

numerical image-based FE and DE model provides a feasible way to study the fracture evolution process of HCR under dynamic impact loads. The major advantage of the image-based model is that it has the ability to represent accurate meso-structure of HCR. The following sections will introduce the image-based FE and DE models that were applied for modelling geological materials. In regard to image-based FE model, Xu et al. introduced the technique of digital image processing (DIP) based on the finite element method (FEM) is to study the mechanical properties of soil–rock mixture (SRM) [23]. Yang et al. developed a novel 3D XCT-image based FE fracture modelling method, and attempted to validate 3D damage initiation evolution until failure for concrete [24-25]. Ren et al. used the image-based FE model to study 3D homogenization of elastic properties of concrete [26-27]. Djordjevic et al. validated that the image based numerical modelling is a useful tool for investigation of the rock breakage process[28].Wang et al. used the numerical modelling and image analysis technologies to investigate the behavior of actual rocks in the breakage process[29]. Yue et al. presented a DIP-FE method for the two-dimensional mechanical analysis of geo-materials by actually taking into account their material in-homogeneities and microstructures [30-31]. Zhang et al. applied the FEM to compute the static effective elastic properties from 3D micro-tomographic images of Berea sandstone saturated with different fluids [32]. From above analysis, many researches have validated that FE model has the ability to simulate materials discontinuous, additional assumptions are still needed. Since DEM was developed to simulate material discontinuous, Refs [33-36] proved validated that a DE model is a very intuitive approach in modeling HCR mixture mechanics. As to meso-scale modelling of geological materials in discrete methods, the parallel-bond model (PBM) [37] in PFC3D was used to simulate complicated damage initiation and evolution in geological rock materials under dynamic impact loads. The PBM has been used to investigate the mechanical behavior of rocks [38-39] and characterize the breakage properties of rocks in crushing process [40-41]. Lisjak A et al. described a new acoustic emission (AE) modelling technique based on the combined finite-discrete element method (FDEM) 3

[42]. Beckmann et al. presented a two-dimensional simulation of concrete behavior using the DEM [43]. Nitka M et al. showed the 3D numerical mesoscopic results on progressive concrete fracture at aggregate level in a notched concrete beam under bending [44]. From above analysis, the use of DEM to simulate the fracture behavior of heterogeneous rock is very powerful, including meso-quantitative analyses which are difficult to be obtained directly from laboratory experiments. In this study, we illustrated the comparative calculation results on HCR fracture under impact loads at meso-scale level using two different approaches: the meso-scale image-based FE and DE model. In this model, the HCR was characterized by two phases, namely, coal and gangue respectively. The geometry of HCR meso-structure was obtained from real coal rock specimen using non-destructive XCT techniques. The paper was organized as follows: in section 2, the techniques of XCT image segmentation and image-based FE modelling were illustrated; in section 3, the XCT image-based DE modelling were carried out; in section 4, the comparison results of the non-linear fracture process in 3D models were discussed, the advantages and disadvantages of the selected two approaches to analyze fracture mechanism in HCR were outlined. 2. Image-based FE model generation and set up 2.1 XCT image acquisition, processing and analysis A coal cube specimen of size 40mm (Figure 1(a)) was first cut from raw coal sample in the lab, and then scanned using a Micro XMT-400 scanner with 150Kv energy level. Finally, the raw 8 GB 3D dataset of 16 bit float with 1372 slices were obtained. For each slice of images from the scanning test, there are 1720×1771 pixels of 32µm in x×y direction in the sampled image dataset (Figure 1(b)). To reduce the data processing time and possible edge effects, a 32 mm cube region of interest (ROI) was cropped from the 40 mm cube volume (Figure 1(c)). The 3D data size was reduced to 500MB for the 3D original dataset with the resolution also compressed from 32µm to 100µm (Figure 1(d)) using Avizo software. These raw compressed CT images are first processed using the non-local mean 4

filters to reduce noise. This is followed by segmentation of material phases and internal cracks using the multi-threshold segmentation algorithm. However, to model complicated internal structure in HCR of this study so as to obtain the accurate morphology of gangue phase, further operation, such as the marker controlled watershed algorithm [5], on the 2D binary images are carried out to get the better results as a double check to avoid cases such as some noise similar to the insufficient contrast boundaries. A 3D reconstructed model can be converted from the segmented slices. After reconstructed, the volume fraction for coal and gangue phases are 84.42% and 15.58% respectively. Figure 2 shows the 3D meso-morphologies of HCR after 3D reconstruction. 2.2 Computational Details of Volumetric Mesh Generation After different mineral phases are segmented, the 3D surface mesh model are generated based on the marching cube algorithm. However, the original surface mesh model may always contain topological deficiencies (such as non-equilateral shapes and overlapped triangular elements), which would be great barrier to generate volumetric tetrahedral mesh. Some triangular meshes are overlapped or intersected in the original surface mesh model, especially the interfaces between coal and gangue phases. This is often the case and especially in HCR because of its complex components. So it is quite important to a FE model that the meshes for different mineral phases have share nodes at the interface, the intersecting edges shown in Figure 3(a) were deleted to make the assembled HCR have shared surfaces (Figure 3(b)). The 3-Matic software [45] can solve this problem by creating a non-manifold assembly based on the Boolean operation. This tool allows assembling different parts into one combined part. Since non-conforming surface generally tends to create high stress concentration, thus resulting in unrealistic loads transferred from coal phase to gangue phase. It is well known that the higher regularity of triangular element may produce better and more reliable FE results. In this regard, the mesh smoothing procedure is performed to the original surface mesh model using 3-Matic software. It is worth note that the major features of the original geometry (e.g., volumes, sizes, and shapes of 5

gangue phases) should be preserved in smoothing process. The smoothed surface mesh shown in Figure 4(a) for HCR is often a very large dataset, too many triangle elements can sharply increase the computation cost to an unacceptable level. Thus, surface simplification [46] is performed to the smoothed surface model to reduce the triangle density (i.e. triangles per unit area), and finally the simplified surface model shown in Figure 4(b) is generated. As shown in Figure 4, there are 215266 and 36386 triangles in original surface model and simplified surface model respectively. Using the improved surface mesh, the 3-Matic software is used to generate 3D four-node tetrahedral meshes at a specified mesh density (i.e. elements per volume) while preserving the various inclusions in HCR. The elements along the interfaces form conformal interfaces between coal phase and gangue phase (i.e. two tetrahedral elements share the same facet across the coal-gangue interfaces). 938,034 tetrahedral meshes, which are suitable for a finite element computation, were then generated in the specific geometric regions (Figure 4(c)). In Figure 4(c), the 3D tetrahedral mesh of coal and gangue is shown in green and blue respectively. 2.3 Constitutive model of HCR The continuous surface cap model (CSCM), validated in the previous studies[47-48], and has been found to be suitable for modelling brittle materials, such as coal rock and concrete. The principal strain is selected as the failure criterion due to the strength of HCR varies with strain rate under dynamic impact loads. The failure of the CSCM model of HCR is controlled by the erode parameter of CSCM. That is, when the maximum principal strain exceeds the set threshold, the failure of the element is no longer involved in the following analysis. The complete theoretical analysis of CSCM is described in the LS-DYNA keyword user’s manual [49]. Here, the yield equation is expressed as Y ( I1 , J 2 , J 3 ) = J 2 − ℜ2 ( J 3 ) Ff2 ( I1 ) Fc ( I1 , κ )

(1)

where, I1, J2 and J3 are the first, second and third invariants of the stress tensor, respectively; κ is the cap hardening parameter, and ℜ is the Rubin three-invariant reduction factor; Fc is the hardening cap, and Ff is the shear failure surface. 6

The impact simulations for HCR have been used to examine how the breakage behaviors are affected by different loading directions(x-, y-, and z-axis). The following parameters were assigned to each individual phase of HCR: coal phase (ρc =1540 kg/m3, Ec=2.3e9 Pa, vc= 0.21) and gangue phase (ρg=2950 kg/m3, Eg = 6.0e9 Pa, vg = 0.3) respectively. The Ls-dyna program was used to simulate the loading process of HCR, where the displacement-controlled loading scheme is applied to the loaded surface node of the HCR sample, while the opposite face is fixed vertically, and lateral faces are allowed freely. The explicit solver with multiple cores was used to speed up computation. 3. 3D DE model and input parameters 3.1 Linear parallel bond model The contact model adopted in PFC3D is linear parallel bond model [34], which provides an effective tool to describe the mechanical behavior of crushable materials. The parallel-bond component, which can transmit both force and moment between the balls, can be simplified as a set of elastic springs with constant normal and shear stiffness. The relative motions between the two contact balls cause increments of contact forces and moments due to the contact stiffness. The force–displacement law for the parallel bond force and moment is shown in Figure 5. The particle movements, the resultant forces and moments follow Newton's law of motion. The parallel bond model can be mainly described by five meso-parameters: the normal stiffness k n , the shear stiffness k s , the normal bond strength σ c , the shear bond strength τ c , and the radius multiplier λ .The change of contact forces and moments due to the relative particle movements are determined respectively by:

F n = −k n A∆δ n , F s = −k s A∆δ s , M t = −k s J ∆θt , M b = −k n I ∆θb

(2)

where, F n and F s are the parallel-bond force of the contact zone in the normal (n) and shear (s) directions; M t and M b are the parallel-bond moment of the contact zone in the normal (n) and shear (s) directions; ∆δ n and ∆δ s are the relative 7

displacement between the two bonded particles in the normal (n) and shear (s) directions respectively; ∆θt and ∆θb are the relative rotations between the two bonded particles in the normal (n) and shear (s) directions respectively; A , I, J are the cross-section area, moment of inertia of the parallel bond cross-section, and inertia moment of the parallel bond cross-section respectively, as expressed in the following equations 2

A =π R ,I =

4 4 1 1 πR ,J = πR 4 2

(3)

where, R = λ min( R A , RB ) is the bond radius; RA , RB is the radius of the bonded particles respectively. If either of the maximum stresses acting on the two contacting balls exceeds its corresponding bond strength, the parallel bond breaks along with its forces, moments and stiffness is then removed from the model. The parallel bond fails if one of the equations below is satisfied: σ max =

Mb R Fn +β ≥σc I A

(4)

τ max =

Mt R Fs +β ≥τ c A J

(5)

where, σ max and τ max are the maximum tensile and shear stresses; β is the moment contribution factor having a range of 0~1. 3.2 3D DE model generation and set up In HCR DE modelling process, the total 91672 discrete elements (balls) with the radius 0.433mm were putted into the reconstructed surface model up corresponding to the specimen size (32mm in x-, y- and z-axis) until the required compactness was obtained. 77387 and 14285 balls were identified as coal and gangue clusters respectively based on the simplified surface model. The followed procedures are the settlement of the balls until their total kinetic energy became insignificant, and then all forces between the balls were deleted and friction coefficient was set on the target value. 8

In DEM analysis, the macroscale mechanical parameter, such as elastic modulus E, Poisson ratio µ and compressive strength σ, cannot directly imposed into the numerical model. Thus, it is critical to select appropriate meso-parameters for the numerical model that aims to represent the accurate macroscale mechanical behavior of HCR. The uniaxial compression simulations (UCS) on a given DE model are carried to reproduce the experimental behaviors of the selected HCR sample accurately. The calibration approach is used to determine these appropriate meso-parameters, and the detailed calibration procedures [48] have demonstrated in Figure 6. The following six main meso-parameters are needed for our DE simulations: Ball density , ball stiffness ratio k n k s , effective modulus E*, tensile strength , cohesion c and friction angle ϕ , which might be successfully calibrated with real laboratory uniaxial compression experiments of HCR specimens. In addition, the radius multiplier λ , normal damping ratio  , installation gap gr, and friction coefficient f are also required. It is also noting that material softening is not assumed in the numerical model. The final meso-parameters of the simulation model are shown in Table1. With the assumed meso-parameters of the HCR, the DEM calculations provided the uniaxial compressive strength of 8.9 MPa, the compressive elastic modulus of 3.9 GPa and Poisson ratio of 0.21, similarly as in the standard uniaxial compression experiments of coal rocks (σ = 9.6MPa, E = 4.2GPa and µ = 0.23). In this study, two rigid and frictionless walls are created on the DE model of HCR in three directions (x-, y- and z-axis) respectively, which serves as the loading wall and fixed wall respectively. Dynamic impact loading in three directions (x-, y- and z-axis) are then modelled to investigate the 3D meso-structure effects on failure patterns and fracture mechanism of HCR under a certain impact velocity. A series of dynamic impact loading simulations in negative y-axis direction are carried to study the effects of different input energies on failure patterns and fracture mechanism of HCR. The DE model of HCR under positive x-, y- and z-load are shown in Figure 7 (a), (b) and (c) respectively. As shown in Figure 7(b), the top and bottom walls are created to apply 9

the load in positive y- direction, where a certain impact velocity is applied on the top wall and the bottom wall remains fixed, and lateral deformations are allowed freely. 4. Result and discussion In this section, the effects of loading impact modes (x-load, y-load and z-load), the heterogeneous 3D internal meso-structures and impact velocities (2-6m/s) on failure patterns of the HCR are investigated, and then the comparison results of fracture mechanism and failure patterns are also analyzed. 4.1 Effects of different loading directions 4.1.1 Finite-element prediction of failure patterns The simulated total nodal reaction forces for HCR as function of the time at 6m/s under different loading directions are shown in Figure 8. As shown in Figure 8, the force increases first and then decreases, reflecting that the HCR sample becomes more compact under impact loads and then expanded after cracks propagate at the following stage. It can also be seen that the peak loads and the post-peak softening responses exhibit high similarities for all six curves. The maximum vertical force in negative and positive load modes (x-, y-, z-axis) are 19.40, 19.00, 19.50, 19.60, 19.30 and 18.30kN respectively, with 6% maximum difference. These differences are caused by the spatial distribution of mineral phases under different loading directions. Figure 9 shows the damage evolution state of HCR sample (at peak loading points shown in Figure 8) under different loading directions with impact velocity 6m/s. As shown in Figure 9, if the loading direction is perpendicular to the growth direction of gangue interface (in x-and y-load modes), the failure zone tends propagate along the gangue interface. Otherwise, if the loading direction deviates from gangue growth direction, the failure zone tends to form failure zone at the corner of HCR sample, which can be validated in z-load mode. It can be concluded that the damage bands tend to propagate from the surfaces to the interior of HCR, and the contacting area of gangue phase on the loading surface has a significant influence on failure patterns of HCR. Qualitative results showed that the stress concentration appears near the gangue interface, which might be one of the reasons that cause damage failure zone tend to 10

propagate from gangue interface. The very different pictures reflect different failure patterns caused by the spatial distribution of gangue phases and different loading directions. 4.1.2 Discrete-element prediction of failure patterns The wall forces for HCR as functions of the time at 6m/s under different loading directions are shown in Figure 10. It should be noted that the numerical simulation process was ended with the condition of the wall force reaches 40% of the maximum wall force. The maximum vertical force in negative and positive load modes (x-, y-, z-axis) are 20.61, 17.72, 20.37, 18.68, 19.66 and 20.35kN respectively, with 14% maximum difference. Compared with the calculated forces shown in Figure 8 and 10, the calculated FEM and DEM curves are all increased first and then decreased in the entire loading time. However, the results in DEM are slightly higher than these results in FEM, and the calculated wall forces shown in Figure 10 fluctuate obviously, and have multi-peak characteristics under different loading directions, which indicated that the damage results in DEM are more sensitive to material meso-parameters, and the differences in load-carrying capacities are caused by different loading directions and the heterogeneous internal meso-structures. Figure 11 shows the damage evolution state of HCR (at peak loading points shown in Figure 10) under different loading directions. It is clear from Figure 11 that different failure patterns can be seen in different loading mode. It should be also noticed that the fragment color, which is determined by the built-in language FISH provided by PFC3D, mainly caused by the evolution state of HCR. As shown in the failure patterns in x- and y- load, the damage bands tend to propagate from gangue interface to form several main fragments due to the relative movement between HCR and the wall. The failure patterns in z-load are quite different with the results shown in x- and y-load. The damage bands initiate from the upper-right corner of HCR, and then expand to several main fragments that can be seen from Figure 11 (e) and (f). Compared with the numerical results of DEM, the differences in both failure patterns and load-carrying capacities are caused by different loading directions and the heterogeneous internal meso-structures. Qualitative understanding of the failure 11

patterns of HCR further confirmed that the growth direction of strong gangue interfaces and loading directions have greatly influenced on failure patterns of HCR. It can be concluded that the damage and fracture appears to propagate from the surfaces to the interior of HCR, and the number of discrete elements of gangue phase contacting the wall has a significant influence on the failure patterns of HCR. It should be noted that the failure patterns in DEM are similar with these results in FEM, and the DEM better follows the fracture process, including the micro-crack initiation, macro-crack formation and the fragment propagation. 4.2 Effects of different impact velocities 4.2.1 Finite-element prediction of failure patterns The Ls-dyna solver is employed with total time 1ms to ensure the dynamic impact loading for all the simulations in this paper. The impact simulations in negative y-load mode for HCR have been used to examine how the breakage behaviors are affected by different impact velocities (2-6m/s). Figure 12 shows the calculated force for HCR as function of the time under different impact velocities. As shown in Figure 12, the maximum nodal reaction force is 17.80, 18.10, 18.40, 18.90 and 19.50kN with impact velocities of 2, 3, 4, 5 and 6 m/s respectively. This together with different pre-peak curves demonstrated that the input energy for the HCR affects the macro-scale mechanical response. As for 3D models under dynamic impact loading, the differences in both failure patterns and load-carrying capacities are greatly influenced by different input energies, or the heterogeneous internal meso-structures. Figure 13 shows the evolution state of HCR sample (at peak loading points shown in Figure 12) under different impact velocities. The evolution state of HCR at 6m/s can be seen in Figure 9(c). The very different pictures reflect different failure patterns caused by different input energies. It can be observed that the stress concentration appears near the gangue interface, which indicated that the failure zone tends to propagate along the interface from inside to outside surface. It can also be noted that the failure pattern propagation process depends on both the models and parameters used for the HCR. 12

4.2.2 Discrete-element prediction of failure patterns Figure 14 shows the relationship of the wall force as function of the time under negative y-load mode under different impact velocities. The maximum wall forces of the HCR are calculated to be 17.79, 17.87, 19.13, 19.24 and 20.37kN in negative y-load with the impact velocity of 2, 3, 4, 5 and 6 m/s respectively, with 12% maximum difference. As shown in Figure 14, the maximum wall forces increased as impact velocities, which further indicated that more input energies can lead to greater difference in damage characteristics of HCR. The maximum wall forces in DEM are slightly higher than those in FEM under difference impact velocities. Although the overall change trend of the wall force is similar to the results in FEM, the calculated wall forces seem more sensitive, and the curves in DE model need shorter time to reach its peak load than FE model. It mainly due to the mesoscale parameters have greatly influenced on the DE model. The internal damage evolution states of HCR (at peak loading points shown in Figure 14) in negative y-load direction are shown in Figure 15. It is clear from the evolution images of 2m/s and 3m/s that only the upper corner of coal phase has been changed to failure state, which indicated that HCR is not sufficiently damaged under the relatively lower impact velocities. As the impact velocity increases, the damaged bands tend to propagate along the gangue phase, which can be validated by the evolution image (c) and (d). For the results shown in negative y-load, gangue phases are not crushed sufficiently, which further confirmed that the propagation of damage bands are hindered by strong gangue phase under the condition of the loading direction is perpendicular to the growth direction of gangue interface. Combined with the analysis of FEM and DEM, it can be concluded that the approximate location of gangue interfaces and the loading directions are both play an important role in the breakage behavior of HCR. The effect of impact velocity on coal rock damage characteristics has also been reported in reference [50]. 5. Conclusion This study presented two different numerical 3D mesoscopic approaches based on XCT images to simulate damage processes of HCR under dynamic impact loading. 13

The XCT image-based numerical model proved to be an effective tool that gives insights into the meso-deformation mechanisms of HCR undergoing dynamic impact failure behavior. The comparative calculation results of 3D meso-structure of HCR, loading directions and impact velocities on failure patterns and fracture mechanism of HCR are demonstrated. The results showed that the loading directions, input energies and the spatial distribution of different mineral phases, have pronounced influenced on material strength as well as failure patterns of HCR. It can be concluded that the damage and fracture appears to propagate from the surfaces to the interior of HCR, and the area of gangue phase contacting the loading wall has a significant influence on the failure patterns of HCR. The calculated forces in FEM and DEM are all increased first and then decreased in the entire loading time, and the results in DEM are slightly higher than these results in FEM. DEM better follows the fracture process, including the micro-crack initiation, macro-crack formation and the fragment propagation. It was also found that the calculated wall forces in DEM seem more sensitive, and the curves in DE model need shorter time to reach its peak load than FE model. It mainly due to the damage results in DEM are more sensitive to material meso-parameters, and the differences in failure patterns are caused by different loading directions and internal meso-structures. Nomenclature ρ

the density of the sample,[ kg/m3]

E

the elastic modulus of the sample,[Pa]

v

the Poisson’s ratio of the sample

ρc

the density of coal phase,[ kg/m3]

ρg

the density of gangue phase,[ kg/m3]

Ec

the elastic modulus of coal phase ,[Pa]

Eg

the elastic modulus of gangue phase,[Pa]

vc

the Poisson’s ratio of coal phase

vg

the Poisson’s ratio of gangue phase 14

I1

the first invariant of the stress tensor

J2

the second invariant of the stress tensor

J3

the third invariant of the stress tensor

κ

the cap hardening parameter

ℜ

the Rubin three-invariant reduction factor

Fc

the hardening cap

Ff

the shear failure surface.

R

the bond radius,[m]

σ max

the maximum tensile stresses,[Pa]

τ max

the maximum shear stresses,[Pa]

β

the moment contribution factor

kn ks

Ball stiffness ratio



Tensile strength,[Pa]



Cohesion,[Pa]



Friction angle[°]



Normal damping ratio

Acknowledgments The authors gratefully acknowledge the financial support received from the China Postdoctoral Science Foundation (2018M630676), National Nature Science Foundation of China (Nos. 51675521 and 51779224), Zhejiang Basic Public Welfare Research Program (No. LHZ19E090002) and Open Foundation of Shandong Province Key Laboratory of Mine Mechanical Engineering (No.2019KLMM105). References [1] F. Shi, W. Zuo, Coal breakage characterization — part 1: breakage testing with the JKFBC, Fuel. 117 (2014) 1148–1155. [2] K.H. Zheng, C.L. Du, D.L. Yang. Orthogonal Test and Support Vector Machine Applied to Influence Factors of Coal and Gangue Separation, International Journal of 15

Coal Preparation & Utilization. 34 (2) (2014) 75-84. [3] J.P. Li, D.L. Yang, C.L. Du, Evaluation of an underground separation device of coal and gangue, International Journal of Coal Preparation & Utilization, 33 (4) (2013) 188–193. [4] D.L. Yang, J.P. Li, C.L. Du, K.H. Zheng, S.Y. Liu, Particle size distribution of coal and gangue after impact-crush separation, Journal of Central South University, 24 (6) (2017) 1252-1262. [5] C.L. Lin, J.D. Miller, T. Nguyen, A. Nguyen, Characterization of Breakage and Washability of ROM Coal using X-ray Computed Tomography, International Journal of Coal Preparation and Utilization. (2017) 1-4. [6] J. Ali, J.K. Farooqi, D. Buckthorpe, A. Cheyne, P. Mummery, Comparative study of predicti ve FE methods for mechanical properties of nuclear Composites, J. Nucl. Mater. 383 (2009) 247–253. [7] D. Asahina, E.N. Landis, J.E. Bolander, Modelling of phase interfaces during pre-critical crack growth in concrete. Cement Concr. Compos. 33 (9) (2011) 966–977. [8] L. Babout, T.J. Marrow, D. Engelberg, P.J. Withers, X-ray micro-tomographic observation of intergranular stress corrosion cracking in sensitised austenitic stainless steel, Mater. Sci. Technol. 22 (2006) 1068–1075. [9] Chen G , Hao Y , Hao H . 3D meso-scale modelling of concrete material in spall tests[J]. Materials and Structures. 48(6) (2015) 1887-1899. [10] S.M. Kim, R.K. Abu Al-Rub. Meso-scale computational modelling of the plastic-damage response of cementitious composites. Cem. Concr. Res. 41 (3) (2011) 339–358. [11] S.J. Hollister, N. Kikuchi, Homogenization theory and digital imaging: A basis for studying the mechanics and design principles of bone tissue, Biotechnology and bioengineering. 43 (7) (1994), 586-596. [12] H.K. Man, J.G.M.V. Mier, Size effect on strength and fracture energy for numerical concrete with realistic aggregate shapes, International Journal of Fracture. 154 (1-2) (2008) 61-72. [13] S.A. McDonald, G. Dedreuil-Monet, Y.T. Yao, A. Alderson, P.J. Withers, In situ 16

3D X-ray microtomography study comparing auxetic and non-auxetic polymeric foams under tension, Phys. Status Solidi. 248 (2011) 45-51. [14] R. Sharma, P. Mahajan, R.K. Mittal, Elastic modulus of 3D carbon/carbon composite using image-based finite element simulations and experiments,Compos. Struct. 98 (2013) 69-78. [15] T. Hirono, M. Takahashi, S. Nakashima, In situ visualization of fluid flow image within deformed rock by X-ray CT, Eng. Geol. 70 (2003) 37–46. [16] P. Kodali, N. Dhawan, T. Depci, C.L. Lin, J.D. Miller, Particle damage and exposure analysis in HPGR crushing of selected copper ores for column leaching, Minerals Engineering. 24 (13) (2011) 1478-1487. [17] L.C. Bam, J.A. Miller, M. Becker, F.C. De Beer, I. Basson. X-Ray Computed Tomography-Determination of Rapid Scanning Parameters for Geometallurgical Analysis of Iron Ore, The Third AusIMM International Geometallurgy Conference, (2016) 209-220. [18] K. Hiraki, Y. Yamazaki, T. Kanai, Numerical Analysis of 3-D Microstructure of Coke Using Micro X-ray CT, Isij International, 52(11) (2012) 1966-1972. [19] L. Jiang, O. Nishizawa, Y. Zhang, H. Park, Z. Xue, A novel highpressure vessel for simultaneous observations of seismic velocity and in situ CO2 distribution in a porous rock using a medical X-ray CT scanner, J. Appl. Geophys., 135 (2016) 67–76. [20] H. Kang, X. Zhang, L. Si, Y. Wu, F. Gao, In-situ stress measurements and stress distribution characteristics in underground coal mines in China, Eng. Geol., 116 (2010) 333–345. [21] J. Lubliner, J. Oliver, S. Oller, E. Oñate, A plastic-damage model for concrete, Int. J. Solids Struct. 25 (3) (1989) 299–326. [22] G.H. Mahmud, Z.J. Yang, A. Hassan, Experimental and numerical studies of size effects of ultra-high performance steel fibre reinforced concrete (UHPFRC) beams, Constr. Build. Mater. 48 (2013) 1027–1034. [23] W.J Xu, Z.Q. Yue, R.L. Hu, Study on the meso-structure and meso-mechanical characteristics of the soil–rock mixture using digital image processing based finite element method, International Journal of Rock Mechanics & Mining Sciences. 45(5) 17

(2008) 749-762. [24] Y.J. Huang, Z.J. Yang, W.Y. Ren, 3D meso-scale fracture modelling and validation of concrete based on in-situ X-ray Computed Tomography images using damage plasticity model, International Journal of Solids and Structures. 68 (2015) 340-352. [25] Y.J. Huang, Z.J. Yang, X.W. Chen, Monte carlo simulations of meso-scale dynamic compressive behavior of concrete based on X-ray computed tomography images, International Journal of Impact Engineering, (2016) S0734743X1630389X. [26] W.Y. Ren, Z.J. Yang, R. Sharma, 3D meso-scale image-based fracture modelling of concrete using cohesive elements, 22nd UK National Conference on Computational Mechanics in Engineering. Exeter. 2014. [27] W.Y. Ren, Z.J. Yang, R. Sharma, C.H. Zhang, P.J. Withers, Two-dimensional Xray CT image based meso-scale fracture modelling of concrete. Eng. Fract. Mech. 133 (2015) 24–39. [28] N. Djordjevic, Image based modeling of rock fragmentation, Minerals Engineering 47 (2013) 68-75. [29] Y.C. Wang, Numerical modelling of heterogeneous rock breakage behavior based on texture images, Minerals Engineering. 74 (2015) 130-141. [30] Z.Q. Yue, S. Chen, L.G. Tham, Finite element modeling of geo-materials using digital image processing, Computers and Geotechnics. 30 (5) (2003) 375-397. [31] W. Wang, P.L. He, D.H. Zhang. Finite element simulation based on soil meso-structures extracted from digital image, Soil Mechanics and Foundation Engineering. 51 (1) (2014) 17-22. [32] Y. Zhang, Z.M.N. Toks, Impact of the cracks lost in the imaging process on computing linear elastic properties from 3D micro-tomographic images of Berea sandstone, Geophysics, 77 (2) (2012) 95-104. [33] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1979) 47–65. [34] P.W. Cleary, C.S. Campbell, Self-lubrication for long runout landslides: examination by computer simulation, J. Geophys. Res. Solid Earth. 98 (1993) 21911– 18

21924. [35] C. Thornton, L. Liu. How do agglomerates break?, Powder Technology. 143(26) (2004) 110-116. [36] P.W. Cleary, Predicting charge motion power draw, segregation and wear in ball mills using discrete element methods, Minerals Engineering. 11 (1998) 1061–1080. [37] P.W. Cleary, R. Morrisson, S. Morrell, Comparison of DEM and experiment for a scale model SAG mill, International Journal of Mineral Processing. 68 (1) (2003) 129–165. [38] P.W. Cleary, M.D. Sinnott, Simulation of particle flows and breakage in crushers using DEM: Part 1–Compression crushers, Minerals Engineering, 74 (2015) 178–197. [39] G.W. Delaney, P.W. Cleary, R.D. Morrison, S. Cummins, B. Loveday, Predicting breakage and the evolution of rock size and shape distributions in AG and SAG mills using DEM, Minerals Engineering. 50 (2013) 132–139. [40] X. Ding, L. Zhang, A new contact model to improve the simulated ratio of unconfined compressive strength to tensile strength in bonded particle models, International Journal of Rock Mechanics and Mining Sciences. 69 (2014) 111–119. [41] N. Djordjevic, Influence of charge size distribution on net-power draw of tumbling mill based on DEM modelling, Minerals Engineering, 18(3) (2005) 375– 378. [42] A. Lisja, Q. Liu , Q. Zhao, Numerical simulation of acoustic emission in brittle rocks by two-dimensional finite-discrete element analysis, Geophysical Journal International. 195(1) (2013) 423-443. [43] D.I.B Beckmann, D.I.K. Schicktanz, D.M.D. Reischl, DEM simulation of concrete fracture and crack evolution, Structural Concrete, 13(4) (2012) 213-220. [44] M. Nitka, J. Tejchman, A three-dimensional meso-scale approach to concrete fracture based on combined DEM with X-ray µCT images, Cement and Concrete Research. 107 (2018) 11-29. [45] J. Li, H. Zhang, Y. Peng, A New Measurement Technique of the Characteristics of Nutrient Artery Canals in Tibias Using Materialise's Interactive Medical Image Control System Software, Biomed Research International. (3) (2015) 1-7. 19

[46] E. Béchet, J.C. Cuilliere, F. Trochu, Generation of a finite element MESH from stereo-lithography (STL) files, Comput. Aided. 34(1) (2002) 1–17. [47] M. Tao, X. Li, C. Wu. Characteristics of the unloading process of rocks under high initial stress, Computers and Geotechnics. 45 (2012) 83-92. [48] M. Tao, X. Li, D. Li, Rock failure induced by dynamic unloading under 3D stress state, Theoretical and Applied Fracture Mechanics, 65 (2013) 47-54. [49] LSTC. LS-DYNA keyword user’s manual, Livermore: Livermore Software Technology Corporation. (2007). [50] K.H. Zheng, C.L. Du, J.P. Li, Numerical simulation of the impact-breakage behavior of non-spherical agglomerates, Powder Technology. 286 (2015) 582-591.

20

Table 1 Meso-parameters of the PBM used in the DE model

Coal

Meso-parameters

Value

Meso-parameters

Value

Ball density ρc (kg/m3)

1540

Ball stiffness ratio k n k s

2.17

Effective modulus E* (Pa)

2.3e9

Tensile strength
 (Pa)

7.90e6 3.95e6

Friction coefficient f

0.3

Cohesion  (Pa)

Installation gap gr (m)

1e-4

Friction angle  (°)

1

Normal damping ratio 

0.75

Ball density ρg (kg/m3)

2950

Ball stiffness ratio k n k s

2.582

Effective modulus E* (Pa)

6.0e9

Tensile strength
 (Pa)

3.65e7 3.04e7

Radius multiplier 

Gangue

Friction coefficient f

0.3

Cohesion  (Pa)

Installation gap gr (m)

1e-4

Friction angle  (°)

Effective modulus E* (Pa)

34

1

Normal damping ratio 

0.75

4.15e9

Tensile strength
 (Pa)

3.42e7 1.72e7

Radius multiplier λ

Interfaces

28

Friction coefficient f

0.3

Cohesion  (Pa)

Installation gap gr (m)

1e-4

Friction angle  (°)

Ball stiffness ratio k n k s

2.38

Normal damping ratio 

21

31 0.75

Figure 1 Sample preparation and image compression. In this figure, (a) presents the raw coal sample; (b) presents the CT slice with 1720×1771 pixels(32µm); (c) presents the CT slice 1000× 1000 pixels(32µm); (d) presents the CT slice 320×320 pixels(100µm).

Figure 2 3D reconstruction of coal sample. In this figure, (a) presents the 3D reconstructed of raw coal sample; (b) presents the 3D reconstructed of the segmented HCR.

Figure 3 Illustration of the non-manifold assembly technology. In this figure, (a) presents a triangular element intersected in a square element; (b) presents a non-manifold assembly operation make the two selected meshes have share nodes at the interface.

22

Figure 4 3D meso-structure surface and volumetric mesh generation. In this figure, (a) presents the original surface model; (b) presents the simplified surface model; (c) presents the volumetric tetrahedral mesh generation.

Figure 5 Mechanical response laws for the parallel bond force and moment. In this figure, (a) presents the normal force versus parallel-bond surface gap; (b) presents the shear force versus shear displacement; (c) presents the twisting moment versus twist rotation; (d) presents the bending moment versus bend rotation; (e) presents the failure envelope of the parallel bond.

23

Figure 6 Calibration procedures for meso-parameters of HCR

Figure 7 3D DE model of HCR generation and set up. In this figure, (a) presents the positive x-load mode; (b) presents the positive y-load mode; (c) presents the positive z-load mode.

24

Figure 8 Simulated forces for HCR as functions of the time under different loading directions

Figure 9 Damage evolution state (at peak loading points shown in Figure 8) of HCR under different loading directions. In this figure, (a) presents the negative x-load mode; (b) presents the positive x-load mode; (c) presents the negative y-load mode; (d) presents the positive y-load mode; (e) presents the negative z-load mode; (f) presents the positive z-load mode.

25

Figure 10 Calculated wall forces for HCR as functions of the time at 6m/s under different loading directions

Figure 11 Damage evolution state (at peak loading points shown in Figure 10) under different loading directions with impact velocity 6m/s. In this figure, (a) presents the negative x-load mode; (b) presents the positive x-load mode; (c) presents the negative y-load mode; (d) presents the positive y-load mode; (e) presents the negative z-load mode; (f) presents the positive z-load mode.(a' and b' represents the failure zone of coal and gangue phases respectively, c' and d' represents the undamaged zone of coal and gangue phases respectively.)

26

Figure 12 Simulated forces for HCR as functions of the time under different impact velocities

Figure 13 Damage evolution state (at peak loading points shown in Figure 12) of HCR under different impact velocities. In this figure, (a) 2m/s; (b) 3m/s; (c) 4m/s; (d) 5m/s.

Figure 14 Wall forces for HCR as functions of the time in negative y-load modes under different impact velocities 27

Figure 15 Failure pattern of HCR (at peak loading points shown in Figure 14) under different impact velocities in negative y-load mode. In this figure, (a) 2m/s; (b) 3m/s; (c) 4m/s; (d) 5m/s; (a' and b' represents the failure zone of coal and gangue phases respectively, c' and d' represents the undamaged zone of coal and gangue phases respectively.)

28

1. Dynamic impact loading of HCR in three directions(x-, y- and z-axis) are modelled. 2. The effects of loading direction and impact velocity on failure patterns are studied. 3. The comparison results of the meso-scale image-based FE and DE model are discussed.