Study on Meso-Material Parameters of Submarine Weathered Granite Based on Parallel Bond Model

In order to study the mechanical properties of submarine weathered granite under marine geological conditions, uniaxial compression tests were carried out on the original medium weathered granite of the seafloor of an offshore area in Pingtan, Fujian Province by using triaxial experimental apparatus to analyze the fracture characteristics, stress–strain characteristics, and compressive strength indexes. Based on the theory of discontinuous medium, the uniaxial compression and uniaxial tensile tests of rocks were simulated, and the microscopic mechanical parameters of discrete elements of granite samples were determined based on the indoor macroscopic mechanical tests: effective modulus Et, compressive elastic modulus Ec, macro Poisson’s ratio μ, and uniaxial compressive strength σc. The results show that the parallel bond model has good simulation results for the uniaxial compression test, but the tensile strength and tensile–compression ratio were quite different from the experimental values. When the confining pressure is large, the calibrated parameter adaptability by uniaxial compression is poor. The reason for certain errors is a large resistance of the parallel bond model to particle rotation and the influence of normal stress on shear strength is not considered. The cementation model can be modified by adding coefficients based on laboratory test results.


Introduction
As a major renewable energy source, wind energy is considered the most effective way to reduce greenhouse gas emissions and has become an increasingly important renewable energy market in the past few decades [1].
The coastal wind resource in China (e.g., Fujian, Guangdong, and other coastal areas) is abundant and the offshore wind turbine (OWT) is being developed rapidly [2][3][4]. However, due to the influence of submarine islands and reefs, the topographic and geological conditions are complex [5,6], the rock formations are shallow, and the depth varies greatly. The properties of the overburdened soil are weak. Therefore, the rock-socketed construction method is needed for offshore wind turbine pile foundations [7].
However, the existing safety design and analysis methods of the offshore wind power pile foundation are more for soil seabed [8][9][10][11] and less for rock seabed. For the rock-soil mass of rock-based seabed, it is difficult to test due to the large difference in mechanical properties caused by different [12] weathering degrees. Additionally, under the extreme ocean load, the microscopic failure mechanism of weathered rock mass is not clear and the simulation accuracy is low. So, it is urgent to study the mechanical properties of rock-soil mass of rock-based seabed and clarify the mechanical properties and failure mechanism of submarine weathered granite. Rock macroscopic mechanical properties are affected to varying degrees by the differences in internal microscopic components properties of rock. Therefore, studying the fracture law of rock under various loads from the perspective of micromechanics is an important basis for revealing the fracture failure mechanism and constitutive behavior of rock, and it is also a key to establishing the relationship between macroscopic mechanical response and micromechanics of rock. In this paper, based on the original strongly weathered granite samples from offshore wind farms of an offshore area in Fujian Province, the mechanical models of uniaxial compression, triaxial compression, and Brazilian disc splitting of weathered rock masses were established by PFC and bondedparticle model to investigate the micromechanics.

Physical Properties of Rock Samples
The test soil samples were taken from the engineering area of a wind farm in an offshore area in Pingtan Comprehensive Experimental Zone, Fujian Province where a large amount of late Yanshan granite is distributed. In the field, The rotary soil sampler with a single hole and double pipes was used, and the plant gum was supplemented as the circulating fluid for drilling and sampling several sections of weathered granite with depths of 30 to 35 m in the same borehole. The sampling method not only has a small disturbance on soil samples but also can ensure the undisturbed soil samples.
For the experiment, 27, 19, and 20 groups of samples (ZWJZK1-27, ZWJZK4-21, ZWJZK5-23) that contain residual soil, fully weathered soil, granular soil, and fragmented soil were taken from three different boreholes. The particle size distribution curve of weathered granite is shown in Figure 1. The percentage of rock mass larger than 1 mm in weathered rock particle size can reach 60%, and the percentage of rock mass smaller than 0.1 mm in weathered rock particle size is relatively small. The uneven coefficient C u = 23 and the curvature coefficient C c = 2.23 were obtained, which meet the requirements of good gradation. affected to varying degrees by the differences in internal microscopic components properties of rock. Therefore, studying the fracture law of rock under various loads from the perspective of micromechanics is an important basis for revealing the fracture failure mechanism and constitutive behavior of rock, and it is also a key to establishing the relationship between macroscopic mechanical response and micromechanics of rock. In this paper, based on the original strongly weathered granite samples from offshore wind farms of an offshore area in Fujian Province, the mechanical models of uniaxial compression, triaxial compression, and Brazilian disc splitting of weathered rock masses were established by PFC and bonded-particle model to investigate the micromechanics.

Physical Properties of Rock Samples
The test soil samples were taken from the engineering area of a wind farm in an offshore area in Pingtan Comprehensive Experimental Zone, Fujian Province where a large amount of late Yanshan granite is distributed. In the field, The rotary soil sampler with a single hole and double pipes was used, and the plant gum was supplemented as the circulating fluid for drilling and sampling several sections of weathered granite with depths of 30 to 35 m in the same borehole. The sampling method not only has a small disturbance on soil samples but also can ensure the undisturbed soil samples.
For the experiment, 27, 19, and 20 groups of samples (ZWJZK1-27, ZWJZK4-21, ZWJZK5-23) that contain residual soil, fully weathered soil, granular soil, and fragmented soil were taken from three different boreholes. The particle size distribution curve of weathered granite is shown in Figure 1. The percentage of rock mass larger than 1 mm in weathered rock particle size can reach 60%, and the percentage of rock mass smaller than 0.1 mm in weathered rock particle size is relatively small. The uneven coefficient Cu = 23 and the curvature coefficient Cc = 2.23 were obtained, which meet the requirements of good gradation.

Uniaxial Compression Test of Weathered Granite
The whole process uniaxial stress-strain test was carried out by the MTS815 full digital hydraulic servo rock triaxial testing machine that consists of a loading system, testing system, control system, and program control system, independently developed by Shandong University, as shown in Figure 2.

Uniaxial Compression Test of Weathered Granite
The whole process uniaxial stress-strain test was carried out by the MTS815 full digital hydraulic servo rock triaxial testing machine that consists of a loading system, testing system, control system, and program control system, independently developed by Shandong University, as shown in Figure 2. The main parameters of MTS815 include: the maximum static axial test force is 2667 kN, the maximum dynamic axial test force is 1335 kN, the strain measurement range is ±0.03 mm, and the displacement measurement range is ±50 mm. The static test force accuracy is 1%, and the dynamic test force accuracy is 2%. The uniaxial stress-strain whole process test and triaxial stress-strain whole process test can be completed in MTS815. The programmable uniaxial and triaxial tests are carried out according to the requirements of the special test process. The axial load, axial displacement, the stroke of axial large range, and the circumferential displacement can be switched without impact in the control.
The height and diameter of weathered granite are 100 mm and 50 mm, respectively. The sample was packaged with rubber film and installed on the base for fixing. The sensor was installed, and the sensitivity of the sensor was tested from the control system. It was confirmed that the sensor could operate the test instrument for loading after normal operation. After setting the confining pressure in the working condition, at the beginning of the test, the confining pressure was applied to the set value at 3 MPa/min, and then the axial pressure was applied at 30 kN/min. When the specimen entered the yield stage (exceeding the peak value), the circumferential displacement rate was adjusted to 0.04 mm/min. When the specimen reached peak stress, in order to obtain the residual strength quickly, the axial displacement rate was adjusted to 0.08 mm/min until the specimen was damaged ( Figure 3). The main parameters of MTS815 include: the maximum static axial test force is 2667 kN, the maximum dynamic axial test force is 1335 kN, the strain measurement range is ±0.03 mm, and the displacement measurement range is ±50 mm. The static test force accuracy is 1%, and the dynamic test force accuracy is 2%. The uniaxial stress-strain whole process test and triaxial stress-strain whole process test can be completed in MTS815. The programmable uniaxial and triaxial tests are carried out according to the requirements of the special test process. The axial load, axial displacement, the stroke of axial large range, and the circumferential displacement can be switched without impact in the control.
The height and diameter of weathered granite are 100 mm and 50 mm, respectively. The sample was packaged with rubber film and installed on the base for fixing. The sensor was installed, and the sensitivity of the sensor was tested from the control system. It was confirmed that the sensor could operate the test instrument for loading after normal operation. After setting the confining pressure in the working condition, at the beginning of the test, the confining pressure was applied to the set value at 3 MPa/min, and then the axial pressure was applied at 30 kN/min. When the specimen entered the yield stage (exceeding the peak value), the circumferential displacement rate was adjusted to 0.04 mm/min. When the specimen reached peak stress, in order to obtain the residual strength quickly, the axial displacement rate was adjusted to 0.08 mm/min until the specimen was damaged ( Figure 3).
In the small-strain stage, the stress-strain of the medium weathered rocks in both groups in the test showed an obvious linear relationship ( Figure 4), with the peak partial stress of 139.92 MPa, and the strains reaching the peak partial stress of 0.318% and 0.322%, respectively, which indicates stable rock properties. According to the stress-strain curve, the average deformation modulus and Poisson's ratio of rock can be calculated to be 49. 17   In the small-strain stage, the stress-strain of the medium weathered rocks in both groups in the test showed an obvious linear relationship ( Figure 4), with the peak partial stress of 139.92 MPa, and the strains reaching the peak partial stress of 0.318% and 0.322%, respectively, which indicates stable rock properties. According to the stress-strain curve, the average deformation modulus and Poisson's ratio of rock can be calculated to be 49.17 GPa and 0.266.

Numerical Model
The PFC program based on the discontinuous medium theory is one of the methods that can effectively analyze the macroscopic nonlinear mechanical behavior response of rock from the mesoscopic perspective. The analysis results are in good agreement with the actual results [13][14][15][16]. Therefore, the PFC program has been recognized and affirmed by researchers in the field of geotechnical engineering. In the detailed observation and study of rock mechanical properties using the PFC program, the selected model has an important influence on the accuracy of the final research results. At present, the commonly  In the small-strain stage, the stress-strain of the medium weathered rocks in both groups in the test showed an obvious linear relationship ( Figure 4), with the peak partial stress of 139.92 MPa, and the strains reaching the peak partial stress of 0.318% and 0.322%, respectively, which indicates stable rock properties. According to the stress-strain curve, the average deformation modulus and Poisson's ratio of rock can be calculated to be 49.17 GPa and 0.266.

Numerical Model
The PFC program based on the discontinuous medium theory is one of the methods that can effectively analyze the macroscopic nonlinear mechanical behavior response of rock from the mesoscopic perspective. The analysis results are in good agreement with the actual results [13][14][15][16]. Therefore, the PFC program has been recognized and affirmed by researchers in the field of geotechnical engineering. In the detailed observation and study of rock mechanical properties using the PFC program, the selected model has an important influence on the accuracy of the final research results. At present, the commonly

Numerical Model
The PFC program based on the discontinuous medium theory is one of the methods that can effectively analyze the macroscopic nonlinear mechanical behavior response of rock from the mesoscopic perspective. The analysis results are in good agreement with the actual results [13][14][15][16]. Therefore, the PFC program has been recognized and affirmed by researchers in the field of geotechnical engineering. In the detailed observation and study of rock mechanical properties using the PFC program, the selected model has an important influence on the accuracy of the final research results. At present, the commonly used models in the PFC program mainly include the traditional bonded particle model (BPM) [17,18], the flat-plate joint model (FJM), and the aggregation particle model (CPM). BPM has a low computational cost and does not require complex constitutive relations to describe the mechanical behavior of complex rock. Zhou [19][20][21] studied the fracture mechanism of prefabricated fractured rock mass under a series of different loads (uniaxial compression, Brazilian splitting, and triaxial compression) using traditional BPM. Park [22] introduced a rock structural plane based on the BPM model and conducted a large number of rock direct shear tests. To study the deformation and strength characteristics of granite under true triaxial loading, Zhang [23] performed 3D PFC numerical simulations using BPM. Zhang [24] discussed the acoustic emission characteristics of rocks under different compression loading rates based on BPM.
In this paper, the discrete element model of weathered granite was established by BPM model. Through the parameter calibration of the uniaxial compression test and uniaxial tensile test, the fitting formula of BPM model was obtained. Subsequently, Brazilian splitting and biaxial compression tests were simulated by the calibrated parameters, and the calculation results were analyzed according to the calculation principle of BPM model. Finally, the advantages, disadvantages, and reasons for the calculation results of BPM model were discussed, which provided reference and guidance for revealing the micromechanical properties of weathered granite.

Fundamentals of the Standard BPM
The parallel bond can be expressed as two sets of springs, one parallel to the contact surface between two contacting particles, and the other perpendicular to the surface. Each set of springs has its stiffness and strength parameters. When the stress acting on the bond exceeds the corresponding strength, the bond breaks. The maximum tensile and shear stresses ( Figure 5) are calculated according to beam theory as follows.
where F n i and F s i are the normal and shear forces, respectively. The normal and shear forces of particle B; M n i and M s i are the bending and twisting moments acting at the center of the parallel bond; A, I, and J are the area of the bond cross section, the area of the moment of inertia and the polar inertia, the moment of inertia and the polar moment of inertia, respectively; and R is the average radius of particle A and particle B. If the maximum stress is greater than the corresponding strength σ max > σ c or τ max > τ c , the parallel bond will fracture as a tensile or shear crack, respectively, and it will disappear or become a zero-length bond. At the same time, the bond stress will be zero and redistributed to the adjacent bond [25].

Calibration of Uniaxial Compression Test Parameters for Weathered Granite
The two-dimensional model was used for calculations. According to the experimental data (Table 1), the specimen model was 50 mm wide and 100 mm high, and the particle radius was randomly distributed in the range of 0.5 to 5 mm. The void ratio 3 Figure 5. Force-displacement behavior of bonded-particle model. Modified after Potion and Cundall (2004) [13].

Calibration of Uniaxial Compression Test Parameters for Weathered Granite
The two-dimensional model was used for calculations. According to the experimental data (Table 1), the specimen model was 50 mm wide and 100 mm high, and the particle radius was randomly distributed in the range of 0.5 to 5 mm. The void ratio was 0.296, and the density was 2280 kg/m 3 . According to the characteristics of the bonded-particle model, the tensile elastic modulus is first fitted by the stress-strain curve under direct tension. The linear stiffness (emod) was kept to a relatively small value (1 × 10 5 ), and the parallel bonded effective modulus (Pb_deform) was adjusted to 1 × 10 9 , 5 × 10 9 , 10 × 10 9 , 20 × 10 9 , 40 × 10 9 , and 50 × 10 9 , and other parameters were taken to higher values. Stress-strain curves under direct tension were obtained ( Figure 6).

Calibration of Uniaxial Compression Test Parameters for Weathered Granite
The two-dimensional model was used for calculations. According to the experimental data (Table 1), the specimen model was 50 mm wide and 100 mm high, and the particle radius was randomly distributed in the range of 0.5 to 5 mm. The void ratio was 0.296, and the density was 2280 kg/m 3 . According to the characteristics of the bondedparticle model, the tensile elastic modulus is first fitted by the stress-strain curve under direct tension. The linear stiffness (emod) was kept to a relatively small value (1 × 10 5 ), and the parallel bonded effective modulus (Pb_deform) was adjusted to 1 × 10 9 , 5 × 10 9 , 10 × 10 9 , 20 × 10 9 , 40 × 10 9 , and 50 × 10 9 , and other parameters were taken to higher values. Stress-strain curves under direct tension were obtained ( Figure 6).   The green means shear crack and the blue means tension crack. The relationship between tensile elastic modulus and Pb_deform is obtained by fitting the peak stress and strain (Figure 7). It can be seen from the figure that the fitting formula is E t , E t = 1.1161x + 0.005, where E t is the effective modulus and the tensile elastic modulus to be calibrated is 49.17 GPa. Then, the tensile elastic modulus was substituted into the fitting formula, x should be approximately taken as 44.1 GPa, that is, the Pb_deform emod is 44.1 GPa. A new fitting formula is obtained by reducing the number of data sets. The uncertainty of a variable E t with respect to parameter x is d = ∆k/k = 10.12%, k is the fitting coefficient of x.
For the linear stiffness (emod) of the linear contact modulus, the size of the parallel bond modulus is kept unchanged (Pb_deform emod is 58.7 GPa), and the linear elastic modulus is changed to 1 × 10 9 , 10 × 10 9 , 20 × 10 9 , 40 × 10 9 and 50 × 10 9 , respectively. The stress-strain curve and the failure curve of the specimen under uniaxial compression are obtained (Figure 8).
The red means shear crack and the purple means tension crack. The relationship between contact linear modulus and elastic modulus is shown in Figure 9. The dot line in the Figure 9 is the original data, and the solid line is the curve image after fitting. It can be seen from the figure that the fitting formula is E c = 1.1804x + 36.407, where E c is the compressive elastic modulus and x is the actual linear deform emod/1 × 10 9 . At this time, the tensile modulus to be calibrated is 49.17 GPa. Then, the tensile modulus was substituted into the fitting formula, x should be approximated as 10.8 GPa; that is, the linear deform emod is 10.8 GPa. A new fitting formula is obtained by reducing the number of data sets. The uncertainty of a variable E c with respect to parameter x is d = ∆k/k = 7.89%.
The green means shear crack and the blue means tension crack. The relationship between tensile elastic modulus and Pb_deform is obtained by fitting the peak stress and strain (Figure 7). It can be seen from the figure that the fitting formula is Et, Et = 1.1161x + 0.005, where Et is the effective modulus and the tensile elastic modulus to be calibrated is 49.17 GPa. Then, the tensile elastic modulus was substituted into the fitting formula, x should be approximately taken as 44.1 GPa, that is, the Pb_deform emod is 44.1 GPa. A new fitting formula is obtained by reducing the number of data sets. The uncertainty of a variable Et with respect to parameter x is d = Δk/k = 10.12%, k is the fitting coefficient of x. For the linear stiffness (emod) of the linear contact modulus, the size of the parallel bond modulus is kept unchanged (Pb_deform emod is 58.7 GPa), and the linear elastic modulus is changed to 1 × 10 9 , 10 × 10 9 , 20 × 10 9 , 40 × 10 9 and 50 × 10 9 , respectively. The stress-strain curve and the failure curve of the specimen under uniaxial compression are obtained (Figure 8). The red means shear crack and the purple means tension crack. The relationship between contact linear modulus and elastic modulus is shown in Figure 9. The dot line in the Figure 9 is the original data, and the solid line is the curve image after fitting. It can be seen from the figure that the fitting formula is Ec = 1.1804x + 36.407, where Ec is the compressive elastic modulus and x is the actual linear deform emod/1 × 10 9 . At this time, 0.005, where Et is the effective modulus and the tensile elastic modulus to be calibrated is 49.17 GPa. Then, the tensile elastic modulus was substituted into the fitting formula, x should be approximately taken as 44.1 GPa, that is, the Pb_deform emod is 44.1 GPa. A new fitting formula is obtained by reducing the number of data sets. The uncertainty of a variable Et with respect to parameter x is d = Δk/k = 10.12%, k is the fitting coefficient of x. For the linear stiffness (emod) of the linear contact modulus, the size of the parallel bond modulus is kept unchanged (Pb_deform emod is 58.7 GPa), and the linear elastic modulus is changed to 1 × 10 9 , 10 × 10 9 , 20 × 10 9 , 40 × 10 9 and 50 × 10 9 , respectively. The stress-strain curve and the failure curve of the specimen under uniaxial compression are obtained (Figure 8). The red means shear crack and the purple means tension crack. The relationship between contact linear modulus and elastic modulus is shown in Figure 9. The dot line in the Figure 9 is the original data, and the solid line is the curve image after fitting. It can be seen from the figure that the fitting formula is Ec = 1.1804x + 36.407, where Ec is the compressive elastic modulus and x is the actual linear deform emod/1 × 10 9 . At this time, the tensile modulus to be calibrated is 49.17 GPa. Then, the tensile modulus was substituted into the fitting formula, x should be approximated as 10.8 GPa; that is, the linear deform emod is 10.8 GPa. A new fitting formula is obtained by reducing the number of data sets. The uncertainty of a variable Ec with respect to parameter x is d = Δk/k = 7.89%.  For Poisson's ratio, linear deform emod and Pb_deform emod were kept unchanged, assuming that the stiffness ratio (kratio) of parallel bond component and linear contact component were equal to study the corresponding relationship between stiffness ratio and macroscopic Poisson's ratio. Setting stiffness ratios of 0.5, 1.0, 1.5, 2.0 and 3.0, uniaxial compression was carried out, and the failure results and stress curves of uniaxial compression are shown in Figure 10, where red means shear crack and purple means tension crack. In the figure, gray color represents particles, and blank color is pores. For Poisson's ratio, linear deform emod and Pb_deform emod were kept unchanged, assuming that the stiffness ratio (kratio) of parallel bond component and linear contact component were equal to study the corresponding relationship between stiffness ratio and macroscopic Poisson's ratio. Setting stiffness ratios of 0.5, 1.0, 1.5, 2.0 and 3.0, uniaxial compression was carried out, and the failure results and stress curves of uniaxial compression are shown in Figure 10, where red means shear crack and purple means tension crack. In the figure, gray color represents particles, and blank color is pores. Poisson's ratio was calculated by taking the transverse and longitudinal nominal strains that correspond to half of the peak stress intensity, respectively, and the relationship between Poisson's ratio and kratio is fitted as shown in Figure 11. The dot line in the Figure 11 are the original data, and the blue circles are the curve image after fitting From the figure, the formula is μ = 0.131x − 0.1094, where μ is macro Poisson's ratio, and x is the stiffness ratio. The Poisson's ratio value of 0.266 was substituted into the formula; the stiffness ratio can be obtained as 2.82. A new fitting formula is obtained by reducing the number of data sets. The uncertainty of a variable μ with respect to parameter x is d = Δk/k = 12.04%. Poisson's ratio was calculated by taking the transverse and longitudinal nominal strains that correspond to half of the peak stress intensity, respectively, and the relationship between Poisson's ratio and kratio is fitted as shown in Figure 11. The dot line in the Figure 11 are the original data, and the blue circles are the curve image after fitting From the figure, the formula is µ = 0.131x − 0.1094, where µ is macro Poisson's ratio, and x is the stiffness ratio. The Poisson's ratio value of 0.266 was substituted into the formula; the stiffness ratio can be obtained as 2.82. A new fitting formula is obtained by reducing the number of data sets. The uncertainty of a variable µ with respect to parameter x is d = ∆k/k = 12.04%. For the bonding ratio, after determining the effective contact modulus and stiffness ratio of linear bond and parallel bond, the normal bond or tangential bond strength is defined as the bonding ratio (ten_coh). To study the failure modes of specimens under different bond combinations, the bonding ratio was defined as 0.1, 0.5, 1.0, 1.2, 1.5, and 2.0, and the damage results of the specimens are shown in Figure 12, where green means shear crack and red means tension crack.
The results show that, when the bonding ratio is less than 1, the specimen is prone to tensile damage, and when the bonding ratio is larger, the specimen is prone to shear damage. According to the uniaxial compression damage of the rock in the actual specimen, the bonding ratio can be taken between 1.0 and 1.5, and the cohesion ratio is taken as 1.2 in this model.
For the amplification factor of bonding strength, the bonding ratio is fixed, assuming that the tangential bonding strength pb_coh = 1 × 10 7 Pa, the normal bonding strength pb_ten = 1.2 × 10 7 Pa can be obtained according to the bonding strength ratio 1.2, which will be positioned as the base bonding strength. On the basis of the reference bond strength, the normal bond strength and the tangential bond strength were multiplied by the magnification factor k, which was set to 0.5, 1.0, 2.0, 5.0, and 10.0, in turn. The failure modes and stress curves of different samples were obtained (Figure 13), where red means shear crack and green means crack-mix. In the figure, gray color represents particles, and blank color is pores.. The fitting relationship between the peak stress strength and the magnification factor is shown in Figure 14. The dot line in the Figure 14 are the original data, and the blue circles are the curve image after fitting. The fitting function σ c = 24.682x − 1.5931 can be obtained, where σ c is the uniaxial compressive strength of the rock and x is the cohesive strength ratio. The σ c that value is 139.92 was substituted into the function, the x can be obtained as 5.74, the normal cohesive strength is 68.8 MPa, and the tangential cohesive strength is 57.4 MPa. A new fitting formula is obtained by reducing the number of data sets. The uncertainty of a variable σ c with respect to parameter x is d = ∆k/k = 9.61%.
After the preliminary calibration, the uniaxial compression test and uniaxial tensile test ( Figure 15) were carried out by using a series of mesoscopic deformation and bonding parameters determined previously.
For the uniaxial compression fracture of the specimen, the simulation results are basically consistent with the experimental results ( Figure 3) that the specimen presents one main crack. As can be seen from the simulation results, the cracking type is a tensile-shear composite failure, which is consistent with the research results of Einstein [26]. The elastic modulus under compression was 47.01 GPa, Poisson's ratio 0.241, compressive strength 136.5 MPa, and tensile elastic modulus 42.94 GPa. Since each parameter will also interact with each other, the compressive elastic modulus and tensile elastic modulus obtained by the combination were slightly smaller than those under 65.5 GPa, and the compressive strength was slightly smaller than that under 185 MPa. The results show that, when the bonding ratio is less than 1, the specimen is prone to tensile damage, and when the bonding ratio is larger, the specimen is prone to shear damage. According to the uniaxial compression damage of the rock in the actual specimen, the bonding ratio can be taken between 1.0 and 1.5, and the cohesion ratio is taken as 1.2 in this model.
For the amplification factor of bonding strength, the bonding ratio is fixed, assuming data, and the blue circles are the curve image after fitting. The fitting function σc = 24.682x − 1.5931 can be obtained, where σc is the uniaxial compressive strength of the rock and x is the cohesive strength ratio. The σc that value is 139.92 was substituted into the function, the x can be obtained as 5.74, the normal cohesive strength is 68.8 MPa, and the tangential cohesive strength is 57.4 MPa. A new fitting formula is obtained by reducing the number of data sets. The uncertainty of a variable σc with respect to parameter x is d = Δk/k = 9.61%.  After the preliminary calibration, the uniaxial compression test and uniaxial tensile test ( Figure 15) were carried out by using a series of mesoscopic deformation and bonding parameters determined previously. After the preliminary calibration, the uniaxial compression test and uniaxial tensile test ( Figure 15) were carried out by using a series of mesoscopic deformation and bonding parameters determined previously.

Simulation of Biaxial Compression Test
The biaxial test parameters are shown in Table 2. Based on the above parameters, the biaxial test was carried out. Firstly, the confining pressure is set to be small (1 × 10 6 ), which was compared with the uniaxial test results as shown in Figure 16.

Sample
No.  It can be concluded that the compressive strength is 132.37 MPa, and the simulation results are slightly larger than the uniaxial test by 4%. According to the elastic modulus and confining pressure under various working conditions, the fracture results are shown in Figure 17. The stress-strain curve under confining pressure of 20 × 10 6 was shown in Figure 18 and the calculation results were shown in Table 3.  It can be concluded that the compressive strength is 132.37 MPa, and the simulation results are slightly larger than the uniaxial test by 4%. According to the elastic modulus and confining pressure under various working conditions, the fracture results are shown in Figure 17. The stress-strain curve under confining pressure of 20 × 10 6 was shown in Figure 18 and the calculation results were shown in Table 3.   In the biaxial test, the error of compressive strength increases with the increase in confining pressure using the parameters calibrated by uniaxial compressive strength, and the compressive strength also increases with the increase in the internal friction angle. When the confining pressure is 30 × 10 6 , the compressive strength is 302.7 MPa, and the error between the results and the experimental value is 29.5%. When the friction angle increases from 50.6 to 70.6, the compressive strength increases to 340.1 MPa, and the error between the results and the experimental value is 21%, there is still a large error. If the  In the biaxial test, the error of compressive strength increases with the increase in confining pressure using the parameters calibrated by uniaxial compressive strength, and the compressive strength also increases with the increase in the internal friction angle. When the confining pressure is 30 × 10 6 , the compressive strength is 302.7 MPa, and the error between the results and the experimental value is 29.5%. When the friction angle increases from 50.6 to 70.6, the compressive strength increases to 340.1 MPa, and the error between the results and the experimental value is 21%, there is still a large error. If the In the biaxial test, the error of compressive strength increases with the increase in confining pressure using the parameters calibrated by uniaxial compressive strength, and the compressive strength also increases with the increase in the internal friction angle. When the confining pressure is 30 × 10 6 , the compressive strength is 302.7 MPa, and the error between the results and the experimental value is 29.5%. When the friction angle increases from 50.6 to 70.6, the compressive strength increases to 340.1 MPa, and the error between the results and the experimental value is 21%, there is still a large error. If the friction angle continues to increase, the compressive strength increases slightly. The main reason for this phenomenon is that the BPM model uses a finite length bond to simulate the contact behavior between the particles, the effect of which is equivalent to two sets of springs, as shown in Figure 5. One set of springs is parallel to the contact surface of the particles, and the other is perpendicular to the contact surface, with each set of springs being the carrier of the force applied so that it can simultaneously withstand the clutch torque. However, the slip model cannot exist at the same time as the parallel bond model, and the slip model only works after the parallel bond failure. The BPM model calculates the maximum tensile stress and shear stress according to the beam theory. In this process, bending moment and torque contribute the most to the maximum tensile stress and shear stress, respectively. The parallel bond will cause the excessive rotation resistance of spherical particles (spherical particles will produce excessive rolling after the fracture of the parallel bond), and the shear strength of the BPM model depends on the bond strength of parallel keys, which is a constant, because the shear strength is not related to the normal stress, and the friction strength only plays a role after the parallel bond failure, which is different from the actual mechanical behavior. These all lead to some differences between the calculated results of BPM model and the actual results in the tension-compression ratio and biaxial calculation. However, for the uniaxial compression test, the simulation results are in good agreement with the experimental results, which can reveal the meso-failure mechanism.

1.
The parallel bonding model was calibrated by uniaxial testing, which found the elastic modulus, Poisson's ratio, and compressive strength were in good agreement with the experimental values. However, the tensile strength and tensile-compressive ratio differed significantly from the experimental values, and fine-tuning was still unable to achieve an error within 10%.

2.
To simplify the experimental process, The parameters obtained by uniaxial compression were used for biaxial tests. By changing the confining pressure, the calculated compressive strength is compared with the results of uniaxial compression. When the enclosing pressure was very small, the calculated compressive strength was closer to the strength of uniaxial compression. However, when the enclosing pressure increased, the simulation results are quite different from the experimental results and are smaller than the experimental values. By increasing the friction angle, the calculated compressive strength increased, but the error was still large.

3.
The reason for a certain gap between the calculation results and the experimental results is the large resistance of the parallel bonding model to particle rotation, and to consider the effect of normal stress on shear strength. The bonding model can be modified by adding coefficients based on laboratory test results.