Coal and gas outburst is one of the most serious disasters caused by mining, which is a mechanical failure process, mainly controlled by gas pressure, stress, and coal properties (Hu et al. 2008; Jin et al. 2011; Kursunoglu and Onder 2019; Lei et al. 2021). In 2003, Wu et al. (2003) proposed a novel method to prevent coal and gas outburst by gas hydration, which lied in the fact that the formation of hydrate can reduce the gas pressure and increase the failure strength of coal (Gao et al. 2015), illustrated in Fig. 1. It can be seen from Fig. 1 that gas pressure decreases by 43.88% and failure strength of GHBC increases by 43.94%, after gas hydrate formation in coal at 20 MPa.

Mechanism of coal and gas outburst prevention based on hydrate method: a Coal sample; b Formed gas hydrate; c Hydrate distribution at the coal sample; d Stress–strain curves of coal before and after gas hydrate formation; e Pressure–temperature–time curves during gas hydrate formation in coal sample

It is crucial to understand the strength behavior of coal after gas hydrate formation. Laboratory test is an effective method to gain a first insight into the strength behavior of GHBC. Up to now, there are few experimental studies on the mechanical behavior of GHBC. However, triaxial compression tests have been conducted on the methane hydrate-bearing sediment (MHBS) (Miyazaki et al. 2011; Ghiassian et al. 2013; Hyodo et al. 2013a; Song et al. 2014; Yan et al. 2017). Although the host sediments are different, there are common characteristics: MH (methane hydrate) formation can increase the strength and stiffness of the sediments (Winters et al. 2007; Makogon 2008; Waite et al. 2009; Hyodo et al. 2013a, b; Li et al. 2016; Sun et al. 2019; Wu et al. 2020). Due to the high cost and long test cycle of experimental research, as well as the limitations to revealing the mesoscopic mechanism, numerical simulations such as DEM aided laboratory tests allow for the granular material from micro–meso-macro perspectives, thus making it a promising tool in investigating the mechanical behavior of GHBC.

Some researchers had studied the effects of hydrate saturation on the mechanical behavior of MHBS using DEM methods, and they found that the mechanical behavior of MHBS mainly depended on hydrate saturation (Masui et al. 2005a, b; Yun et al. 2007; Brugada et al. 2010; Miyazaki et al. 2011; He and Jiang 2016a, b; Tang et al. 2020), hydrate distribution (Jung et al. 2012; Dai et al. 2012; Jiang et al. 2013, 2017; Yang and Zhao 2014a, b; Malinverno and Goldberg 2015; Shen and Jiang 2016; Shen et al. 2016), particle size (Yu et al. 2014, 2016), and confining pressure (Miyazaki et al. 2010; Jiang et al. 2017; Zhou et al. 2019). For instance, Brugada et al. (2010) conducted a conventional triaxial drainage simulation test on hydrate sediments with different saturations and confining pressures using PFC 3.0, and found that the contribution of hydrate to sediment strength was the friction angle. He and Jiang (2016a, b) conducted the discrete element modeling on the drained triaxial compression test of the energy soil. The results revealed that volume reduction decreased with increasing hydrate saturation, and dilation angle increases linearly. Jung et al. (2012) modeled the triaxial compression characteristics of cementation and pore-filling hydrate, and concluded that stiffness, strength and expansion trend of sand increased with increasing sediment density or hydrate saturation. Yu et al. (2016) studied the mechanical behavior of MHBS by using spherical or elongated particles to simulate soil particles and triaxial compression tests to simulate two different hydrate formation patterns: pore-filling and cementation. Zhou et al. (2019) modeled the mechanical behavior of hydrate sediments under six different saturations and confining pressures. Some studies had shown that hydrate effects are large on the mechanical properties of MHBS when MH is cemented with soil particles (Masui et al. 2005a, 2005b; Hyodo et al. 2009).

Previous studies had shown that contact models, including contact bond (CB) and parallel bond (PB) models, were suitable for describing the mesoscopic characteristics of soil particles and hydrate in sediment, the effect of hydrate on the mechanical behavior was significant for cementation-type MHBS. Moreover, hydrates are simulated using ball or contact model. However, it should be noted that coal is usually buried in underground less than 2000 m in depth with larger ground stress, resulting in mechanical behavior differ from that of MHBS. Additionally, the different pore characteristics of coal and MHBS will lead to the divergence of gas hydrate formation and distribution in coal and MHBS. Therefore, it is necessary to study different hydrate saturations on the mesoscopic mechanical characteristics of GHBC. The relationship between the macro–meso parameters are of great importance for the numerical model establishment of triaxial compression behavior of GHBC.

This paper will focus on this cementation-type GHBC, and conduct the numerical simulation of triaxial compression tests at confining pressure of 16 MPa and saturations of 20% and 80%, respectively. It is organized as follows: Firstly, the basic assumptions and key points of 3D particle model of GHBC are introduced. Then, normalization discusses the mathematical relations between macro and meso parameters. Next, the numerical simulation results and sample results are compared and analyzed. Finally, the internal relationship between macro and meso mechanical failure is discussed to reveal the strengthening mechanism of MH on coal. Some preliminary conclusions are summarized.

2.1 Model assumption

In the numerical simulation tests of GHBC, we model the mechanical behavior of GHBC from our previous tests of Yu et al. (2019). For the experimental tests, first, the excessive gas method is used to form hydrates in coal samples and then triaxial compression tests are conducted on GHBC with hydrate saturation S_h of 20% and S_h of 80% and at confining pressure \({\sigma }_{3}\) of 16 MPa. According to the test results of Chen et al. (2018), cemented hydrate is formed, when assuming that water and gas fully participate in the reaction using the excessive gas method, as shown in Fig. 2. It can also be seen that hydrate is mainly distributed in coal as cemented type. This paper focuses on cemented hydrate and makes the following three hypotheses for the model:

1. (1)

   Gas is not considered.

2. (2)

   Gas hydrate and coal are simplified as spherical particles.

3. (3)

   Hydrate exists in the form of cementation.

Schematic flow of the 3D image processing and modeling: a Actual GHBC (scanning by X-CT); b Extract the shapes of GHBC; c Determination of the hydrate distribution mode of GHBC—cementation; d Numerical model

2.2 Numerical modelling of triaxial compression

2.2.1 Specimen specification

In the triaxial compression test, the sample size has a great influence on the test results (Tao et al. 1981). When the ratio of the height and the diameter of the sample is about 2 (Yin et al. 2011), the stress in the specimen is evenly distributed and the compressive strength keeps stable (Zhou et al. 2015). It is found that the influence of size effect on calculation results can be ignored, if the numerical sample size is 30 – 40 times the average particle size (Jensen et al. 1999), or when the total number of simulated particles is greater than 2000 (Zhou et al. 2000). To improve computation efficiency, some scholars set the simulated sample size as b × h = 2 × 4 mm to simulate the physical triaxial compression tests with the sample size of 50 × 100 mm (Yang and Zhao 2014a, b; Xu et al. 2010). Hence, the sample size is set as b × h = 2 × 4 mm with the total number of particles being 3477 (S_h of 20%), which meets the requirements, shown in Fig. 3. Particle expansion method (O’Sullivan 2011) is introduced to generate the specimen. Coal particles are firstly shrunk to 25% of their experimental size. Next, all the spheres are expanded ten times to increase the computation efficiency. The particle size distribution in DEM modeling ranges from 0.072 mm to 0.1 mm. Gas hydrate particle size is 0.06 mm (Yang and Zhao 2014a, b). The sample density is 1220 kg/m³. In this study, the displacement control mode (usually used in laboratory tests) is adopted for conducting the numerical test (Wang et al. 2016, 2019). As the simulation is carried out under quasi-static conditions, the loading speed can be ignored (Zhao et al. 2021). Therefore, the loading rate is set as 0.1 mm/s, which is higher than the experimental axial loading speed of 0.01 mm/s, to improve the calculation efficiency.

DEM model of GHBC: a Formation of gas-hydrate; b Test GHBC; c Mesoscopic characteristics of GHBC (scanning by X-CT); d Actual particle shape of coal; e Extraction of coal particle shape; f Simplification of coal particles and gas hydrate particles; g Simulated GHBC

2.2.2 Contact model

This study adopts two basic bond models: the parallel bonding model (PB) and the linear model (LB). The contact between all particles is characterized by PB model (Han et al. 2019). This model is used due to the reason that it can more accurately characterize the meso-structure of rock materials and has better applicability (Liu et al. 2015a, b; Cao et al. 2016; Zhang 2017; Jiang et al. 2021; Yang et al. 2021). The interaction between the wall and the particle is depicted using the LB. The triaxial compression test simulation system and parallel bonding model, as shown in Fig. 4. Where \(\overline{{M }_{i}^{s}}\), \(\overline{{M }_{i}^{n}}\) are the shear and normal moment, respectively; \(\overline{{F }_{i}^{s}}\), \(\overline{{F }_{i}^{n}}\) the shear and normal force, respectively; \(\overline{{k }_{s}}\), \(\overline{{k }_{n}}\) the shear and normal stiffness, respectively; \({g}_{s}\), \(\mu\), \({\sigma }_{c}\), \(\overline{c }\) and \(\overline{\varphi }\) the parallel-bond surface gap, friction coefficient, tensile strength, cohesion and friction angle, respectively.

Simulation of random GHBC particles under the triaxial compression test and mechanical response of the linear parallel bond model [modified from Itasca (2002)]

The force-displacement law for the parallel-bond force and moment consists of the following steps, as shown in Fig. 5:

Force–displacement law for the parallel bond force and moment: a Normal force versus parallel-bond surface gap; b Shear force versus relative shear displacement; c Twisting moment versus relative twist rotation; d Bending moment versus relative bend rotation; e Failure envelope for the parallel bond (Itasca 2016)

1. (1)

   Update the cross-sectional bond properties:

   $$\begin{aligned} \overline{R} & = \overline{\lambda }\min (R^{(1)} ,R^{(2)} ),\;\;ball - ball \\ \overline{R} & = \overline{\lambda }R^{(1)} ,\;\;ball - {\text{facet}} \\ \end{aligned}$$

   (1)
   $$\overline{A} = \pi \overline{R}^{2} ;\;\overline{I} = \frac{1}{4}\pi \overline{R}^{4} ;\;\overline{J} = \frac{1}{2}\pi \overline{R}^{4}$$

   (2)

   where, \(\overline{A }\) is the cross-sectional area; \(\overline{I }\) is the moment of inertia of the parallel bond cross-section; \(\overline{J }\) is the polar moment of inertia of the parallel bond cross-section.

2. (2)

   Update normal contact force \(\overline{{F }_{n}}\) and shear contact force \(\overline{{F }_{s}}\).

   $$\overline{F}_{{\text{n}}} = \overline{F}_{{\text{n}}} + \overline{k}_{n} \overline{A}\Delta \delta_{n}$$

   (3)
   $$\overline{F}_{{\text{s}}} = \overline{F}_{{\text{s}}} + \overline{k}_{s} \overline{A}\Delta \delta_{s}$$

   (4)

   where, \(\Delta {\delta }_{n}\) is the relative normal-displacement increment; \(\Delta {\delta }_{s}\) is the relative shear-displacement increment.

3. (3)

   Update twist moments \(\overline{{M }_{t}}\) and bend moments \(\overline{{M }_{b}}\) (Crandall et al. 1987).

   $$\overline{M}_{{\text{t}}} = \overline{M}_{{\text{t}}} + \overline{k}_{s} \overline{J}\Delta \overline{\theta }_{{\text{t}}}$$

   (5)
   $$\overline{M}_{{\text{b}}} = \overline{M}_{{\text{b}}} + \overline{k}_{n} \overline{J}\Delta \overline{\theta }_{{\text{b}}}$$

   (6)

   where, \(\Delta \overline{{\theta }_{t}}\) is the relative twist-rotation increment; \(\Delta \overline{{\theta }_{b}}\) is the relative bend-rotation increment.

4. (4)

   For three-dimensional discrete element simulation, within the scope of cementation, the maximum tensile stress and maximum shear stress are:

   $$\overline{\sigma } = \frac{{\overline{F}_{{\text{n}}} }}{{\overline{A}}} + \overline{\beta }\frac{{\left\| {\overline{M}_{b} } \right\|\overline{R}}}{{\overline{I}}}$$

   (7)
   $$\begin{aligned} \overline{\tau } & = \frac{{\left\| {\overline{F}_{s} } \right\|}}{{\overline{A}}} + \overline{\beta }\frac{{\left\| {\overline{M}t} \right\|\overline{R}}}{{\overline{J}}} \\ \overline{\beta } & < (0,1] \\ \end{aligned}$$

   (8)

   The moment-contribution factor (\(\overline{\beta }\)) is discussed in (Potyondy 2011).

5. (5)

   Fig. 5e shows the failure envelope of the cement. If it is greater than tensile strength or shear strength, the cement will fail:

   $$\begin{gathered} \overline{\tau }_{c} = \overline{c} - \sigma {\text{tan}}\varphi = \overline{c} - \frac{{\overline{F}_{n} }}{{\overline{A}}}{\text{tan}}\varphi \hfill \\ \hfill \\ \end{gathered}$$

   (9)

2.2.3 Sample generation

Figure 6 shows the flow chart of the sample generation of triaxial compression tests of GHBC. Firstly, the wall size was set as 2  mm × 4 mm, with two infinitely large loading plates (top wall and bottom wall) and a cylindrical side wall created.

Simulation process for sample generation of GHBC

Secondly, sufficient particles were randomly generated in the closed cylindrical region. After the initial "coal particles" were prepared, isotropic consolidation of the sample was carried out, so that the effective stress was 1 MPa and the porosity was 0.4. The number of hydrates that need to be filled in the sediment sample of cemented hydrate is determined to reach the target saturation. Accordingly, coal and gas hydrate particles that meet the requirements are generated at the same time. After that, 1000 steps are cycled to bring the particles to equilibrium. When the hydrate saturation is 20%, the porosity is 0.4 and there are 946 hydrate particles. When the hydrate saturation is 80%, the porosity is 0.4 and there are 3775 hydrate particles.

In addition, the origin of X, Y and Z axes was specified to ensure particles moving in the inter of the specimen.

Finally, “solve” was entered to achieve self-equilibrium by the particle's own gravity.

To calibrate the mesoscopic parameters, the single factor sensibility analysis was carried out to quantify the mathematical relationship between macroscopic and mesoscopic parameters. Secondly, the sensitivity of meso-parameters to macro-parameters was studied by multifactor sensibility analysis to further obtain the fine calibrated meso-parameters. Furthermore, the mesoscopic parameters were adjusted, and the numerical models were verified using the physical tests of GHBC.

3.1 Initial calibration of mesoscopic parameters

Through the literature research (He and Jiang 2016a, b; Yang and Zhao 2014a, b; Li et al. 2018; Duan et al. 2015), the mesoscopic parameters of GHBC with the saturation of 20% and 80% at confining pressure of 16 MPa were preliminarily determined, as shown in Table 1. They were calibrated using the triaxial test of GHBC with the saturation of 20% and 80% at confining pressure of 16 MPa. The numerical comparison with the experimental tests was shown in Fig. 7. It is clear that the variation trend of the numerical test almost reproduces that of the experiment at the initial elastic stage, while varying at the hardening stage. Therefore, it is necessary to refine the mesoscopic parameters. According to our previous research literature (Zhang et al. 2021), the friction angle and cohesion are 18.26° and 4 MPa (S_h of 20%), 14.57° and 6.81 MPa (S_h of 80%), respectively. To make more simulated parameter values close to the real parameter values of laboratory tests, the friction angle and cohesion values are set as 18° and 4 MPa (S_h of 20%), 15° and 7 MPa (S_h of 80%) in this study, respectively.

Comparison of stress–strain curves for numerical and experimental tests

As shown in Fig. 7. By comparing the test results with the simulation results, it can be seen that the initial stages of the two curves were similar to each other. However, the variation trend of the two curves was quite different, which indicates that the model's mesoscopic parameters should be recalibrated.

3.2 Correlation of macro-parameters with meso-parameters

The determination of meso-parameters was an important step of the numerical simulation of DEM. Scholars often used trial-and-error method (Yang et al. 2016; Yang et al. 2018; Wang and Tian 2018; Huang et al. 2019), which was time-consuming, empirical and parameter selection was very random. However, the researches had shown that there was a certain relation between meso-parameters and macro-mechanical parameters (Liu et al. 2015a, b; Xing et al. 2017; Xiao et al. 2020; Xiao et al. 2021; Li and Rao 2021). The influence of meso-parameters, such as bond stiffness ratio \(\overline{{k }_{n}}/\overline{{k }_{s}}\), friction coefficient \(\overline{\mu }\), normal bond strength \(\overline{{\sigma }_{c}}\) and bond radius coefficient \(\overline{\lambda }\), were discussed on macroscopic parameters, such as elastic modulus E and failure strength \({\sigma }_{c}\) to provide a reference for calibrating the meso-parameters of GHBC.

For the numerical simulation scheme of the triaxial tests of GHBC, 20 groups of simulation tests were carried out in this paper, as shown in Table 2.

In this paper, the normalization method is adopted and the "normalization" equation is as follows:

where, x is macro-parameter; x₀ is the initial macro parameter.

According to the mesoscopic parameter values set in Table 2, stress–strain curves under different mesoscopic parameter conditions are drawn, as shown in Fig. 8.

Stress–strain curves under different meso-parameters

Figure 8a shows that under different bond stiffness ratios, stress–strain curves all show strain softening characteristics, and with the increase of bond stiffness ratio, a large strain softening phenomenon appears in the curves. The failure strength decreases and the elastic model increases with the bond stiffness ratio increasing.

Figure 8b shows that the stress–strain curve changes from strain-hardening to strain-softening as the friction coefficient increases. The failure strength and elastic modulus increase with the friction coefficient increasing.

Figure 8c shows that stress–strain curves show strain softening characteristics under different normal bond strengths. The influence of particle normal bond strength on elastic modulus is small and the failure strength increases with the normal bond strength increasing.

Figure 8d shows that the stress–strain curve changes from strain-hardening to strain-softening with the increase of the bonding radius coefficient. The material exhibits greater failure strength and will show greater strain softening.

And the normalized elastic modulus and failure strength corresponding to twenty groups of simulated tests were listed in Table 3. For the stress–strain curve exhibiting strain-hardening type, the failure strength \({\sigma }_{c}\) was the determined as the deviator stress corresponding to the axial strain of 12%. If the curve was strain-softening type, the failure strength was the failure of the curve. The elastic modulus E was adopted as the gradient of the stress–strain curve corresponding strain ranges from 0.5% to 4%. Effects of meso-parameters on \(E\) and \({\sigma }_{\mathrm{c}}\), as shown in Fig. 9.

Effects of meso-parameters on \(E\) and \({\sigma }_{c}\)

The relationship between macro and meso parameters were fitted in Fig. 9, and the fitting formulas with the correlation coefficients were listed in Table 4. It was clear that the goodness of fit \({R}^{2}\) were all greater than 0.95, indicating a good relationship between macro and meso parameters.

3.3 Sensitivity analysis of mesoscopic parameters

3.3.1 Influence of meso-parameters on elastic modulus

Table 5 lists the orthogonal design results of the elasticity modulus to evaluate the sensitivity of the meso-parameters. It can be seen from Table 5 that the friction coefficient was the main factor affecting the elastic modulus, followed by the bond radius coefficient, the normal bond strength and the bond stiffness ratio. Figure 10 shows that the elastic modulus significantly increases with the increase of the friction coefficient, the bond radius coefficient and the bond stiffness ratio. However, the elastic modulus significantly decreases with normal bond strength. With the increase of the normal bond strength, the elastic modulus decreases first and then increases.

Trend diagram of elastic modulus and index

3.3.2 Failure strength

It can be seen from Table 6 that the friction coefficient was the main factor affecting the failure strength, The first is the friction coefficient, followed by the bond radius coefficient, normal bond strength and bond stiffness ratio. Figure 11 shows that the failure strength significantly increases with the increase of the friction coefficient and the normal bond strength. The failure strength increases first and then decreases with the increase of the bond stiffness ratio and the bond radius coefficient. It can be concluded that the friction coefficient was the most sensitive parameter to the elastic modulus and the failure strength of GHBC.

Variation of failure strength with the meso-parameters

3.4 Meso-parameter calibration approach

1. (1)

   Firstly, the internal friction angle and cohesion can be determined by the triaixial compression tests and should be constant.

2. (2)

   Then, both the bond stiffness ratio and the normal bond strength can be equal to 1, which will promote the calibration efficiency.

3. (3)

   Moreover, the friction coefficient and the bond radius coefficient should be calibrated first due to they have a significant effect on the macro parameters. Keep the coefficient of the bond radius constant, and increase the friction coefficient at high saturation while decreasing it at low saturation with the interval of 0.01.

4. (4)

   Finally, the friction coefficient and the bond radius coefficient should be calibrated first due to they have a significant effect on the macro parameters. Keep the coefficient of the bond radius constant, and increase the friction coefficient at high saturation while decreasing it at low saturation with the interval of 0.01.

3.5 Calibration of the meso-parameters

Basis on the above-mentioned relationship between macro and meso-parameters, the meso-parameters were re-calibrated by trial-and-error method. The calibration process was conducted as follows, as shown in Fig. 12. First, the geometric and physical parameters of numerical sample size and density were determined. Then, the contact types of particles in the numerical simulation were specified and the initial meso-parameters were set according to previous results. Furthermore, the simulated stress–strain curves were compared with the experimental stress–strain curves. If the experimental and numerical curves basically coincide with each other, with the error rate of failure strength and elastic modulus within 15% (Huang and Yang 2014; Liu et al. 2015a, b; Han et al. 2019; Zhenhua et al. 2019), then the meso-parameters were adopted. Otherwise, each meso-parameter was recalibrated. The value of numerical meso-parameters was shown in Table 7.

Simulation process for meso-parameters calibration of PB model

Numerical triaxial compression tests on GHBC with different saturations were conducted, and the results were compared with laboratory tests. As shown in Fig. 13a, it can be obviously observed that the simulation stress–strain curves agree well with the test curves. The failure strength and elastic modulus from laboratory tests were compared with the numerical results, shown in Fig. 13b. It was clear that the error rates of elastic modulus and failure strength were both less than 10%.

Stress–strain–volumetric responses of GHBC from DEM simulations and experiments at 20% and 80% of saturation: a Stress–strain relationships; b Elastic modulus and failure strength; c Failure pattern

In addition, as can be seen from the Fig. 13c, the arrow is the direction of motion, and its length (color depth) represents the velocity vector size (relative). The velocity direction of the particles inside the sample is chaotic, and the size is different. The reason is that the particles inside the sample occur mismovement, relative slip and tumbling, and the macro manifestation is volume dilatancy deformation and failure of the sample. This is consistent with the fact that annular expansion failure is the main failure mode in laboratory test. The failure patterns of the numerical almost reproduce those of the experimental.

Mesoscopic properties such as average pore ratio, pore distribution, displacement field, velocity field, contact force chain and contact force direction were important indicators to explain the macroscopic mechanical behavior of geomaterials (Jiang et al. 2009, 2008, 2006). In this study, meso properties such as contact force chains will be discussed.

Figure 14 shows the force chains distribution of GHBC at axial strain stages of 3%, 6%, 9% and 12% under different saturations. According to the transfer external load, it can be divided into strong chain and weak chain (Sun and Wang 2008). The thickness of the force chains represents the magnitude of the force (the thicker the line was, the greater the contact force was). In this Fig. 14, the green was the strong chain and the blue was the weak chain. Take the saturation of 20% for example, it was shown from Fig. 14a that there was a large force chain gap between the particles at the initial stage of loading (ε₁ = 3%), during the loading process, and the particles with large contact force were mainly distributed near the loading plate. In the later compression stage (ε₁ = 6% – 12%), more and more strong chains formed and transferred loading vertically in the sample, indicating that the sample can resist more strength.

Distributions of force chains observed in DEM GHBC samples of different saturations at different axial strains

Figure 15 shows that both the contact force and contact number of GHBC with the high saturation were higher than those of GHBC with saturation of 20%. For GHBC with saturation of 80%, more hydrate particles can fill the pore, which increases the contact number and the contact force. Thus, the friction between particles were increased. Therefore, the bonding strength will increase with the increasing of the hydrate saturation.

Comparison of the contact force and contact number at 20% and 80% of saturations

In this paper, the calibration method of meso-parameters for the DEM models of GHBC was presented, based on the parallel bond model using PFC3D. The influence law of meso-parameters on the macro-parameters was first studied by the single factor sensitivity method. Then, the sensitivity of the meso-parameters was evaluated using the multi-factor sensitivity method. Based on the presented numerical models, the meso-mechanism was discussed in terms of force chain to explain the effect of saturation on the mechanical behavior of GHBC.

1. (1)

   The gas hydrate was made to be characterized by parallel bond, and meso-mechanical parameters were characterized by six indicators such as friction coefficient, the normal bond strength, the bond radius coefficient, the bond stiffness ratio, the cohesion and the internal friction angle. The macroscopic mechanics property parameters were characterized by elastic modulus and failure strength.

2. (2)

   According to the initial calibration results of meso-mechanical parameters, twenty groups of numerical tests were conducted to establish the macro-parameter models using the meso-parameters. The elastic modulus linearly increases with the bonding stiffness ratio and the friction coefficient while exponentially increasing with the normal bonding strength and the bonding radius coefficient. The failure strength increases exponentially with the increase of the friction coefficient, the normal bonding strength and the bonding radius coefficient, and remain constant with the increase of bond stiffness ratio. Additionally, four factor and three levels of nine orthogonal simulation tests were conducted to investigate the influence of the meso-parameters, with the friction coefficient the most sensitive parameter.

3. (3)

   The numerical results were compared with the laboratory triaxial compression tests. The profile of the deviator stress–strain of the numerical and the failure pattern almost reproduces those of the laboratory results. Additionally, the error rates of elastic modulus and failure strength were both less than 10%. It was found that the proposed DEM model can predict the mechanical properties of GHBC materials.

4. (4)

   The deviator stress–strain curves exhibit strain hardening behavior for different GHBC with different saturations. The higher the saturation, the larger the failure strength. The contact force and contact number increase with the saturation increase, which increases the friction coefficient between the particles. The effect of friction increases the higher bonding strength of GHBC.