Bonded-Particle Model with Nonlinear Elastic Tensile Stiffness for Rock-Like Materials

: The bonded-particle model (BPM) is a very efﬁcient numerical method in dealing with initiation and propagation of cracks in rocks and can model the fracture processes and most of macro parameters of rocks well. However, typical discrete element method (DEM) underestimates the ratio of the uniaxial compressive strength to the tensile strength (UCS/TS). In this paper, a new DEM method with a nonlinear elastic tensile model embedded in BPM is proposed, which is named as nonlinear elastic tensile bonded particle model (NET-BPM). The relationships between micro parameters in NET-BPM and macro parameters of specimens are investigated by simulating uniaxial compression tests and direct tension tests. The results show that both the shape coefﬁcient of the nonlinear elastic model and the bond width coefﬁcient are important in predicting the value of UCS/TS, whose value ranging from 5 to 45 was obtained in our simulations. It is shown that the NET-BPM model is able to reproduce the nonlinear behavior of hard rocks such as Lac du Bonnet (LDB) granite and the quartzite under tension and the ratio of compressive Young’s modulus to tensile Young’s modulus higher than 1.0. Furthermore, the stress-strain curves in the simulations of LDB granite and the quartzite with NET-BPM model are in good agreement with the experimental results. NET-BPM is proved to be a very suitable method for modelling the deformation and fracture of rock-like materials.


Introduction
Discrete element method (DEM) is proved to be a very effective method in modelling rock-like materials [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. The discontinuous nature of DEM makes DEM very suitable for dealing with the initiation and propagation of micro cracks in rock-like materials. Within DEMs, the bonded-particle model (BPM) is widely used to simulate the deformation and fracture of rock-like materials [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Since the BPM was first proposed [1], it has been rapidly developed. A clumped particle model was proposed in [2] based on the BPM, which improves the predictive capability of the DEM by successfully predicting the whole failure envelope and increasing the ratio of uniaxial compressive strength to tensile strength (UCS/TS) making it closer to that of rocks. A damage-rate law was first included by Potyondy [3] to extend the BPM and to make it more suitable for simulating the static-fatigue behavior of granite. Utili [19] introduced the Mohr-Coulomb failure criterion into the BPM without considering the bending stiffness, and successfully reproduced the mechanical behavior of frictional cohesive materials. A progressive failure model was applied later in [21], which makes the BPM more

Contact Model
The contact particles are bonded together by the bonds, which can resist normal force, shear force and bending moment. As shown in Figure 1a, particle p and particle q are in contact with each other at point A at the original states. X 0 p and X 0 q are their coordinates respectively. dir 0 pq and dir 0 qp are the vectors pointing from their centers to the contact point A respectively. At time t, the two particles move to new positions X t p and X t q , as shown in Figure 1b. The contact point A on particle p move to point B and the point A on particle q move to point C. The vector dir 0 pq and the vector dir 0 qp rotate to dir t pq and dir t qp , respectively. In addition, θ p and θ q are rotation angles of particles p and q in time t respectively, defining anticlockwise as positive direction. The contact forces and moment applied to particle p from particle q via the bond are normal force F nr pq , shear force F tn pq and moment M pq , as shown in Figure 1c. Appl. Sci. 2017, 7, 686 3 of 20

Contact Model
The contact particles are bonded together by the bonds, which can resist normal force, shear force and bending moment. As shown in Figure 1a, particle p and particle q are in contact with each other at point A at the original states.
where n k is normal stiffness, n pq ε is normal strain of the bond between p and q, the b is shape coefficient of the nonlinear elastic model and 0 n ε is the ultimate normal tensile strain. The shear force and moment are determined by the following equations: where k s and k m are shear stiffness and bending stiffness respectively, dis 0 pq is the original distance between p and q, ε s pq is shear strain of the bond between p and q which is the ratio of deformation in tangential direction to dis 0 pq . A nonlinear elastic tensile stiffness model is proposed to describe the relationship between normal force and normal strain undergoing tension. Then the force in normal direction is expressed as: where k n is normal stiffness, ε n pq is normal strain of the bond between p and q, the b is shape coefficient of the nonlinear elastic model and ε n 0 is the ultimate normal tensile strain. Figure 2 presents the relationship between normal force and normal strain with different value of b. The value of b has great effect on the normal force versus normal strain curve and the ultimate normal tensile force, F nr m . ε n pq and ∆ε s pq which is the increment of ε s pq in each time step are given by: where dis t pq is the distance between the centers of p and q at time t, r p and r q are the radii of particles p and q respectively, u p and u q are the centroid displacements of particles p and q respectively, U pq is an Appl. Sci. 2017, 7, 686 4 of 20 unit vector obtained by rotating the unit vector which points to the center of q from the center of p 90 degrees in the anticlockwise direction, ∆t is time step. Appl. Sci. 2017, 7, 686 4 of 20 where t pq dis is the distance between the centers of p and q at time t, p r and q r are the radii of particles p and q respectively, p u and q u are the centroid displacements of particles p and q respectively, pq U is an unit vector obtained by rotating the unit vector which points to the center of q from the center of p 90 degrees in the anticlockwise direction, t Δ is time step. Another modification to the BPM is that the bonds are assumed as many tiny elastic beams, which distribute between two contact particles in a region wide 2wrm and are parallel to the center line of the two particles, as presented in Figure 3, where the w is bond width coefficient, rm is the smaller one of rp and rq. The range of w is from 0 to 1 and that w equals 0 means all the tiny beams gather in the center line of the two contact particles. Bonds with different bending stiffness but with the same normal stiffness and shear stiffness can be modelled by changing the distribution of the tiny beams. is the upper bound of shear strain, 0 s ε is the ultimate pure shear strain of the bond, ϕ is the internal friction angle and f ϕ is the friction angle between particles when the bond between them fails. In this paper, the stresses are replaced by the strains in order to obtain better simulation results by cooperating with the nonlinear tensile model and the effect of the bend on the tensile strain is taken into consideration as follows.
If the bond is broken or no bond exists between particles p and q which contact with each other, then only normal force and shear force exist between them. The normal force nr pq F is calculated by Equations (3) and (4).  Another modification to the BPM is that the bonds are assumed as many tiny elastic beams, which distribute between two contact particles in a region wide 2wr m and are parallel to the center line of the two particles, as presented in Figure 3, where the w is bond width coefficient, r m is the smaller one of r p and r q . The range of w is from 0 to 1 and that w equals 0 means all the tiny beams gather in the center line of the two contact particles. Bonds with different bending stiffness but with the same normal stiffness and shear stiffness can be modelled by changing the distribution of the tiny beams. Appl. Sci. 2017, 7, 686 4 of 20 where t pq dis is the distance between the centers of p and q at time t, p r and q r are the radii of particles p and q respectively, p u and q u are the centroid displacements of particles p and q respectively, pq U is an unit vector obtained by rotating the unit vector which points to the center of q from the center of p 90 degrees in the anticlockwise direction, t Δ is time step. Another modification to the BPM is that the bonds are assumed as many tiny elastic beams, which distribute between two contact particles in a region wide 2wrm and are parallel to the center line of the two particles, as presented in Figure 3, where the w is bond width coefficient, rm is the smaller one of rp and rq. The range of w is from 0 to 1 and that w equals 0 means all the tiny beams gather in the center line of the two contact particles. Bonds with different bending stiffness but with the same normal stiffness and shear stiffness can be modelled by changing the distribution of the tiny beams. is the upper bound of shear strain, 0 s ε is the ultimate pure shear strain of the bond, ϕ is the internal friction angle and f ϕ is the friction angle between particles when the bond between them fails. In this paper, the stresses are replaced by the strains in order to obtain better simulation results by cooperating with the nonlinear tensile model and the effect of the bend on the tensile strain is taken into consideration as follows.   The modified Mohr-Coulomb criterion with a tension cut-off and an upper bound is used for the shear strength of the bonds, as shown in Figure 4, where ε s max is the upper bound of shear strain, ε s 0 is the ultimate pure shear strain of the bond, ϕ is the internal friction angle and ϕ f is the friction angle between particles when the bond between them fails. In this paper, the stresses are replaced by the strains in order to obtain better simulation results by cooperating with the nonlinear tensile model and the effect of the bend on the tensile strain is taken into consideration as follows.
If the bond is broken or no bond exists between particles p and q which contact with each other, then only normal force and shear force exist between them. The normal force F nr pq is calculated by Equations (3) and (4). dis 0 pq in Equation (4) should be replaced by the sum of the radii of p and q. The shear force F tn pq is calculated as follows: where µ is the coulomb friction coefficient. The resultant contact force vector F p and the resultant contact moment M p on particle p are calculated as follows: where V pq is an unit vector pointing from the center of particle p to the center of particle q. G p is a set of particles which are in contact with particle p.
where pq V is an unit vector pointing from the center of particle p to the center of particle q. p G is a set of particles which are in contact with particle p.

Numerical Simulation Model
The disc-shaped particles are considered as rigid bodies and their translational motions and rotational motions are governed by the standard equations of rigid body dynamics: are damping force and damping moment respectively, which are introduced to dissipate kinetic energy and gain a quasi-static state of equilibrium of the particles and are expressed as [34]: where t α and r α are translational and rotational damping constant respectively. The centroid displacement and rotation angle of particle p can be obtained by using an explicit central difference method to integrate Equations (10) and (11). The algorithm is described as follows:

Numerical Simulation Model
The disc-shaped particles are considered as rigid bodies and their translational motions and rotational motions are governed by the standard equations of rigid body dynamics: where F ext p and M ext p are the external force vector and external moment applied to particle p respectively. m p and J p are mass and moment of inertia of particle p respectively. F d p and M d p are damping force and damping moment respectively, which are introduced to dissipate kinetic energy and gain a quasi-static state of equilibrium of the particles and are expressed as [34]: where α t and α r are translational and rotational damping constant respectively. The centroid displacement and rotation angle of particle p can be obtained by using an explicit central difference method to integrate Equations (10) and (11). The algorithm is described as follows: where the superscript i indicates ith step. ∆t should not exceed the critical time step ∆t c to remain numerical stability as introduced in [34]. The forces and moments are transmitted through the contacts between particles or the contacts between particles and plates. There are two kinds of particle contacts according to whether an effective bond exists between the contact particles: the bonded particle contacts and the unbonded particle contacts. Two particles p and q are bonded together in the original rock specimens if their relative position satisfies the Equation (18): where δ is a nonnegative constant defined by user to eliminate the numerical errors caused during the particle generating. The errors may result in separation of two particles which should touch each other. A bond is created between every bonded particles. If δ is too large, it would lead to forming bond between two particles far away from each other. In our model, δ = 0.001 is suggested. When the bond is broken, the bonded particle contact becomes unbonded particle contact. Furthermore, two particles p and q which are not close to each other may approach each other during the simulation process, and form an unbonded particle contact when their relative position satisfies the Equation (19): The forces and moments on the bonded particle contacts are calculated by Equations (1)-(3). While forces on the unbonded particle contacts, where no moment exists are calculated by Equations (3) and (7). The bonded particle contacts are identified at the beginning of the simulation and updated in each time steps if any bond is broken. The unbonded contacts should be identified in every time steps.
In order to improve the computational efficiency, background grids and cells are introduced in searching particle contacts. In each step, all the particles are stored in the cells formed by the grids according to their centroid position coordinates. The cell size is equal to 2(1 + λ)R max , where R max is the maximum radius of the particles and λ is a user-defined parameter, which is positive and very small to obtain high searching efficiency. The suggested value for λ is 0.01 in our model. A particle can only form a particle contact with particles located in the same cell with it and with particles located in the neighbor cells. To avoid duplicate calculation, particles in cell j are only used to try to form particle contacts with particles in cell j or in neighbor cells whose numbers are larger than j.

Dimensionless Parameters in NET-BPM
In the conventional BPM, the micro parameters are normal stiffness k n , shear stiffness k s , bending stiffness k m , ultimate normal strain ε n 0 , ultimate pure shear strain ε s 0 , upper bound of shear strain ε s max , ultimate relative rotation angle θ 0 , friction angles ϕ and ϕ f , translational damping constant α t and rotational damping constant α r , etc. The macro parameters are UCS, TS, the ratio of UCS/TS, Young's modulus E and Poisson's ratio v, etc. A few more micro parameters, which are w and b are included in the NET-BPM. In addition, because Young's modulus and Poisson's ratio in compression and tension are different, they are denoted as E c and E t , v c and v t , respectively. In this paper, the effects of particle size, particle arrangement, porosity and density of the specimens are not investigated. Young's modulus and Poisson's ratio discussed here are the apparent ones. By using dimensionless analysis similar to [28,34], the dimensionless micro parameters are k s /k n , k m /(k n R), ε s 0 /ε n 0 , ε s max /ε n 0 , ϕ, ϕ f , w, b, and the dimensionless macro parameters are E c R 2 /k n , v c , σ c R 2 /k n ε n 0 , σ t R 2 /k n ε n 0 , E t /E c , v t /v c in our model. R is the characteristic length which is equal to R max . As the nonlinear elastic model is applied, the ratio of E c to E t and the ratio of v c to v t are not equal to 1 anymore, so E t /E c and v t /v c should be taken into consideration.

Specimen Preparation
Uniaxial compression tests and direct tension tests are simulated in this study. The rock specimens used for simulations are created by the particle generating method proposed in [44], which can create dense and stable particle packings with disks once the geometric contour of the specimen and the minimum and maximum particle size are given. The effect of the size of particles and specimens on our method is not concerned in this study. Thus, square specimens with length L = 50 mm created by the method from [44] with R max = 0.9 mm and minimum particle radius R min = 0.3 mm are used for our simulations. Five specimens are created by using different seeds in the initiation of the particle generation, as shown in Figure 5. They are used to study the effects of the particle arrangement on the simulations. Some main parameters of the specimens are listed in Table 1, where R a is the average radius of the particles and overlap ratio is the ratio of the overlap length to the sum of the radii of the two overlap particles.

Specimen Preparation
Uniaxial compression tests and direct tension tests are simulated in this study. The rock specimens used for simulations are created by the particle generating method proposed in [44], which can create dense and stable particle packings with disks once the geometric contour of the specimen and the minimum and maximum particle size are given. The effect of the size of particles and specimens on our method is not concerned in this study. Thus, square specimens with length L = 50 mm created by the method from [44] with Rmax = 0.9 mm and minimum particle radius Rmin = 0.3 mm are used for our simulations. Five specimens are created by using different seeds in the initiation of the particle generation, as shown in Figure 5. They are used to study the effects of the particle arrangement on the simulations. Some main parameters of the specimens are listed in Table 1, where Ra is the average radius of the particles and overlap ratio is the ratio of the overlap length to the sum of the radii of the two overlap particles.  As presented in Figure 6, when the uniaxial compression tests are simulated, the specimens are placed between two parallel plates which are considered as rigid bodies. The bottom one is fixed and the top one moves towards the bottom one with a given constant velocity, V. When the direct tension tests are simulated, the particles in the specimens near the top and bottom boundaries are fixed to the top and bottom plates respectively and have the same velocity with the plate. In order to eliminate the effect of the boundaries on the simulation results, a layer with a height of 0.6L in tensile direction and a width of L perpendicular to the tensile direction is defined, which located in the center of the specimens. Only the bonds existing between the particles located in this layer are allowed to fail, similar to that in [12]. The bottom plate is fixed and the top plate moves away from the bottom one with a constant velocity V. The lateral strain of the specimens is determined by the displacement of two particle groups. One particle group consists of four particles, which are closer to the middle point  As presented in Figure 6, when the uniaxial compression tests are simulated, the specimens are placed between two parallel plates which are considered as rigid bodies. The bottom one is fixed and the top one moves towards the bottom one with a given constant velocity, V. When the direct tension tests are simulated, the particles in the specimens near the top and bottom boundaries are fixed to the top and bottom plates respectively and have the same velocity with the plate. In order to eliminate the effect of the boundaries on the simulation results, a layer with a height of 0.6L in tensile direction and a width of L perpendicular to the tensile direction is defined, which located in the center of the specimens. Only the bonds existing between the particles located in this layer are allowed to fail, similar to that in [12]. The bottom plate is fixed and the top plate moves away from the bottom one with a constant velocity V. The lateral strain of the specimens is determined by the displacement of two particle groups. One particle group consists of four particles, which are closer to the middle point of the left side boundary of the specimen than other particles at the initial time. The other particle group consists of four particles, which are closer to the middle point of the right side boundary than other particles. The average distance of the two groups of particles in the direction perpendicular to the loading direction at the initial time and its variation during simulations are used to calculate the lateral strain. The axial strain is determined by the displacement of the two rigid plates. The loads are calculated by averaging the forces on the top and bottom plates. UCS and TS can be obtained by dividing the ultimate loads by the width of the specimen. The thickness of the specimen is assume to be 1 m in 2-D. E c , E t , v c and v t are determined at the state when the load is equal to half the ultimate load and the secant lines are used to estimate them as introduced in [ of the left side boundary of the specimen than other particles at the initial time. The other particle group consists of four particles, which are closer to the middle point of the right side boundary than other particles. The average distance of the two groups of particles in the direction perpendicular to the loading direction at the initial time and its variation during simulations are used to calculate the lateral strain. The axial strain is determined by the displacement of the two rigid plates. The loads are calculated by averaging the forces on the top and bottom plates. UCS and TS can be obtained by dividing the ultimate loads by the width of the specimen. The thickness of the specimen is assume to be 1 m in 2-D. Ec, Et, vc and vt are determined at the state when the load is equal to half the ultimate load and the secant lines are used to estimate them as introduced in [8,45,46]. The NET-BPM model is compiled in Fortran 90 language for simulating uniaxial compression tests and direct tension tests. The values of the micro parameters in the model are presented in Table 2.
The density of the particles is set to be 2900 kg/m 3 . The critical time step α and r α are set to be 0.02 which are proved to be large enough to guarantee the quasistatic state of equilibrium of the particles. All five specimens presented in Figure 5 are used for simulations with the micro parameters in Table 2.  A simulation case costs about two hours on a machine with an Intel Core 5 3.2 GHz Processor. The compressive loads and tensile loads of specimen A at different displacements are presented in Figure 7. The creak patterns of specimen A under compression and tension are shown in Figure 8. The tcm and ttm are the time when the loads reach their peak values in compression and tension, respectively. The color of the particle varies from black to red gradually according to the ratio of the present number to the original number of the bonded contacts on the particle, which varies from 1 to 0. Tensile failure is the dominant failure pattern of the bonded contacts in both uniaxial compression tests and direct tension tests. The number of the failed bonded contacts was 420 by the time t = tcm, which increased sharply to 931 after 10,000 time steps in compression test. The number of the failed bonded contacts was 7 by the time t = ttm, which increased sharply to 88 after 4000 time steps in tension test. Only one main crack appeared in tension test, which is almost perpendicular to the tensile direction, but several oblique cracks appeared in compression tests. All cracks observed in our simulation are similar to those observed in experiments presented in [47,48].  Table 2. The density of the particles is set to be 2900 kg/m 3 . The critical time step ∆t c = 2.5 × 10 −8 s is obtained with the method introduced in [34] if the safety factor is set to 0.1. Thus, the time step ∆t = 10 −8 s is proper in our model. The loading speed V = 0.1 m/s in both compression and tension tests. Damping constant α t and α r are set to be 0.02 which are proved to be large enough to guarantee the quasi-static state of equilibrium of the particles. All five specimens presented in Figure 5 are used for simulations with the micro parameters in Table 2.  A simulation case costs about two hours on a machine with an Intel Core 5 3.2 GHz Processor. The compressive loads and tensile loads of specimen A at different displacements are presented in Figure 7. The creak patterns of specimen A under compression and tension are shown in Figure 8. The t cm and t tm are the time when the loads reach their peak values in compression and tension, respectively. The color of the particle varies from black to red gradually according to the ratio of the present number to the original number of the bonded contacts on the particle, which varies from 1 to 0. Tensile failure is the dominant failure pattern of the bonded contacts in both uniaxial compression tests and direct tension tests. The number of the failed bonded contacts was 420 by the time t = t cm , which increased sharply to 931 after 10,000 time steps in compression test. The number of the failed bonded contacts was 7 by the time t = t tm , which increased sharply to 88 after 4000 time steps in tension test. Only one main crack appeared in tension test, which is almost perpendicular to the tensile direction, but several oblique cracks appeared in compression tests. All cracks observed in our simulation are similar to those observed in experiments presented in [47,48].  The properties of specimen A-E from the simulations are listed in Table 3. It shows that the difference of the properties between different specimens is small. The maximum ratio of the stand deviation value to the average value is 0.15. This indicates that the specimens created by the method in [44] with the same Rmax and Rmin share almost the same properties.

Relationship between Micro and Macro Parameter
In this Section, the relationships between micro and macro parameters in the NET-BPM introduced in Section 2.3 are investigated, which are revealed by simulating the uniaxial compression tests and direct tension tests with different values of one or two micro parameters. If not specified particularly, the values of the micro parameters are the same as those in Table 2. In order to eliminate the effects of the particle arrangement, specimen size and particle size on the simulations, specimen A introduced in Section 2.4 is used for all the simulations.

Effect of Bond Stiffness on Macro Parameters
Changing the value of ks/kn from 0.1 to 1.0 with km/(knR) equal to 0.064, 0.160 and 0.256, different values of the macro parameters can be observed. The relationships between the macro parameters  The properties of specimen A-E from the simulations are listed in Table 3. It shows that the difference of the properties between different specimens is small. The maximum ratio of the stand deviation value to the average value is 0.15. This indicates that the specimens created by the method in [44] with the same Rmax and Rmin share almost the same properties.

Relationship between Micro and Macro Parameter
In this Section, the relationships between micro and macro parameters in the NET-BPM introduced in Section 2.3 are investigated, which are revealed by simulating the uniaxial compression tests and direct tension tests with different values of one or two micro parameters. If not specified particularly, the values of the micro parameters are the same as those in Table 2. In order to eliminate the effects of the particle arrangement, specimen size and particle size on the simulations, specimen A introduced in Section 2.4 is used for all the simulations.

Effect of Bond Stiffness on Macro Parameters
Changing the value of ks/kn from 0.1 to 1.0 with km/(knR) equal to 0.064, 0.160 and 0.256, different values of the macro parameters can be observed. The relationships between the macro parameters The properties of specimen A-E from the simulations are listed in Table 3. It shows that the difference of the properties between different specimens is small. The maximum ratio of the stand deviation value to the average value is 0.15. This indicates that the specimens created by the method in [44] with the same R max and R min share almost the same properties.

Relationship between Micro and Macro Parameter
In this Section, the relationships between micro and macro parameters in the NET-BPM introduced in Section 2.3 are investigated, which are revealed by simulating the uniaxial compression tests and direct tension tests with different values of one or two micro parameters. If not specified particularly, the values of the micro parameters are the same as those in Table 2. In order to eliminate the effects of the particle arrangement, specimen size and particle size on the simulations, specimen A introduced in Section 2.4 is used for all the simulations.

Effect of Bond Stiffness on Macro Parameters
Changing the value of k s /k n from 0.1 to 1.0 with k m /(k n R) equal to 0.064, 0.160 and 0.256, different values of the macro parameters can be observed. The relationships between the macro parameters and k s /k n at different k m /(k n R) are presented in Figure 9. The results show that E c increases with the increasing of k s /k n at different values of k m /(k n R), and v c decreases with the increasing of k s /k n . The trends are in accord with those in [28,29,34]. With the values of the k s /k n changing from 0.1 to 1.0, Poisson's ratio decreases from 0.45 to 0. The effects of the bending stiffness on E c and v c are very small and can be ignored. However, increasing the bending stiffness can enhance the strength (UCS and TS) of the specimens and increase UCS/TS distinctly when k s /k n is higher than 0.4, as shown in Figure 9c-e. The values of k s /k n and k m /(k n R) have very small effects on the value of E c /E t , as shown in Figure 9f, and small effects on the value of v c /v t .
Appl. Sci. 2017, 7, 686 10 of 20 and ks/kn at different km/(knR) are presented in Figure 9. The results show that Ec increases with the increasing of ks/kn at different values of km/(knR), and vc decreases with the increasing of ks/kn. The trends are in accord with those in [28,29,34]. With the values of the ks/kn changing from 0.1 to 1.0, Poisson's ratio decreases from 0.45 to 0. The effects of the bending stiffness on Ec and vc are very small and can be ignored. However, increasing the bending stiffness can enhance the strength (UCS and TS) of the specimens and increase UCS/TS distinctly when ks/kn is higher than 0.4, as shown in Figure 9c-e. The values of ks/kn and km/(knR) have very small effects on the value of Ec/Et, as shown in Figure 9f, and small effects on the value of vc/vt.

Effect of Shape Coefficient on Macro Parameters
As presented in Figure 2b is a very important micro parameter whose effects on macro parameters are investigated by changing its value from 0.05 to 1.0 with ks/kn equal to 0.2, 0.5 and 0.8 in the simulations. The results are presented in Figures 10 and 11.

Effect of Shape Coefficient on Macro Parameters
As presented in Figure 2b is a very important micro parameter whose effects on macro parameters are investigated by changing its value from 0.05 to 1.0 with k s /k n equal to 0.2, 0.5 and 0.8 in the simulations. The results are presented in Figures 10 and 11.     The results show that with the increasing of b, TS increases but UCS almost remains unchanged and the value of UCS/TS decreases. The effects of b on TS and UCS/TS decrease with the increasing of b and when b is larger than 0.5, the effects become very small, as shown in Figure 10. By choosing a small b, e.g., b = 0.05, a value of UCS/TS as large as 45 can be obtained to meet the practical values of many hard rocks, which is difficult to achieve by the conventional BPM [14]. The value of b has distinct effects on the value of E t , while small effects on the value of E c , as shown in Figure 11a,b. Decreasing b from 1.0 to 0.05 results in increasing of E c /E t from 1.05 to 1.6 in our study at different values of k s /k n , as shown in Figure 11c. Similar to the Young's modulus, the Poisson's ratio are also affected by the value of b. Increasing of b leads to decrease of v c , increasing of v t and sharply decreasing of v c /v t , as shown in Figure 11d-f. The value of v c /v t can reach 4.5 when b is 0.05. When b is larger than 0.7, the variation of b has small effect on the macro parameters.

Effect of Bond Width Coefficient on Macro Parameters
The effects of w on macro parameters are investigated by changing its value from 0.4 to 0.8 with k s /k n = 0.2, 0.5 and 0.8. In order to guarantee the identical value of k m /(k n R) in all the cases, the distributions of the tiny beams in the bonds are changed according to the value of w. The simulation results are presented in Figure 12. The results show that the width of the bond has little effects on Young's modulus and Poisson's ratio, and has almost no effects on E c and v c , as shown in Figure 12a-d. While larger w leads to smaller UCS, as shown in Figure 12e, and has almost no effect on TS, thus results in smaller UCS/TS, as shown in Figure 12f. Noticing Equation (6), if the relative rotation angle is not zero, larger w indicates smaller ε n when the bond fails, that means smaller tensile forces in the bonds. However, this relationship is not shown in the simulated tension tests because the relative rotation angles between contact particles are very small in the tension tests. While in the compression tests, the relative rotation angles are large enough to affect the strength of the specimens and the major failure type of the bonds is tensile failure.
Appl. Sci. 2017, 7, 686 12 of 20 of b and when b is larger than 0.5, the effects become very small, as shown in Figure 10. By choosing a small b, e.g., b = 0.05, a value of UCS/TS as large as 45 can be obtained to meet the practical values of many hard rocks, which is difficult to achieve by the conventional BPM [14]. The value of b has distinct effects on the value of Et, while small effects on the value of Ec, as shown in Figure 11a

Effect of Bond Width Coefficient on Macro Parameters
The effects of w on macro parameters are investigated by changing its value from 0.4 to 0.8 with ks/kn = 0.2, 0.5 and 0.8. In order to guarantee the identical value of km/(knR) in all the cases, the distributions of the tiny beams in the bonds are changed according to the value of w. The simulation results are presented in Figure 12. The results show that the width of the bond has little effects on Young's modulus and Poisson's ratio, and has almost no effects on Ec and vc, as shown in Figure 12a-d. While larger w leads to smaller UCS, as shown in Figure 12e, and has almost no effect on TS, thus results in smaller UCS/TS, as shown in Figure 12f. Noticing Equation (6), if the relative rotation angle is not zero, larger w indicates smaller n ε when the bond fails, that means smaller tensile forces in the bonds. However, this relationship is not shown in the simulated tension tests because the relative rotation angles between contact particles are very small in the tension tests. While in the compression tests, the relative rotation angles are large enough to affect the strength of the specimens and the major failure type of the bonds is tensile failure.

Effect of Ultimate Pure Shear Strain on Macro Parameters
The ultimate pure shear strain has great effects on UCS, as shown in Figure 13a. Changing

Effect of Ultimate Pure Shear Strain on Macro Parameters
The ultimate pure shear strain has great effects on UCS, as shown in Figure 13a. Changing ε s 0 /ε n 0 from 2.0 to 4.5 with k s /k n equal to 0.2, 0.5 and 0.8, UCS increases by more than 50%, but TS, E t and v t remain the same, as shown in Figure 13b,d,f. Increasing ε s 0 /ε n 0 can decrease E c and increase v c , but the effects are very small and can be ignored, as shown in Figure 13c,e.

Effect of Ultimate Pure Shear Strain on Macro Parameters
The ultimate pure shear strain has great effects on UCS, as shown in Figure 13a. Changing

Application of the NET-BPM
In this section, the Lac du Bonnet (LDB) granite which are widely used for model calibration and a quartzite studied in [39] under uniaxial compression and direct tension are modelled with the NET-BPM. Firstly, two specimens are created and the micro parameters in our model are determined for the two rocks. Then the simulation results of LDB granite and the quartzite are compared with the experimental results.

Modelling LDB Granite
A rock specimen with the size same as that in [42] is created using the method in [44] with Rmax = 0.9 mm and Rmin = 0.3 mm. The effect of particle size is not discussed in this study, so the particle size same as that in Section 3 is adopted. The parameters of the specimen for LDB granite are listed in Table 4. The properties of the LDB granite are presented in Table 5, where Young's modulus are from [49], and the rest of the data are from [42]. It should be noticed that Young's modulus and Poisson's ratio are converted to the apparent values.
From the investigation in Section 3, we know that b is the major factor affecting the values of Ec/Et, and that b and ks/kn is the major factor affecting the values of vc. Combining Figures 9b and 11c Table 6. The macro parameters obtained are presented in Table 5.
As the ultimate axial strain in tensile test is small, a lower loading speed of 0.01 m/s is applied in direct tensile tests. The results show that not only Ec, vc, UCS, TS and UCS/TS from the NET-BPM model are in good agreement with the experimental values, which is achieved in [14], but also Et and Ec/Et are in good agreement with the experimental values. The stress-strain curves of the LDB granite under uniaxial compression test from our method are similar to those from experiments in [42], as shown in Figure 14. The crack patterns of LDB granite under compression and tension tests are shown in Figure 15. Same method introduced in Section 2.4 is used to present the cracks. Figure 15a shows the cracks at the time when the compressive load reaches its peak. After that, the cracks develop quickly. Most of them are oblique. Figure 15d shows the cracks at time when the tensile load reaches

Application of the NET-BPM
In this section, the Lac du Bonnet (LDB) granite which are widely used for model calibration and a quartzite studied in [39] under uniaxial compression and direct tension are modelled with the NET-BPM. Firstly, two specimens are created and the micro parameters in our model are determined for the two rocks. Then the simulation results of LDB granite and the quartzite are compared with the experimental results.

Modelling LDB Granite
A rock specimen with the size same as that in [42] is created using the method in [44] with R max = 0.9 mm and R min = 0.3 mm. The effect of particle size is not discussed in this study, so the particle size same as that in Section 3 is adopted. The parameters of the specimen for LDB granite are listed in Table 4. The properties of the LDB granite are presented in Table 5, where Young's modulus are from [49], and the rest of the data are from [42]. It should be noticed that Young's modulus and Poisson's ratio are converted to the apparent values.
From the investigation in Section 3, we know that b is the major factor affecting the values of E c /E t , and that b and k s /k n is the major factor affecting the values of v c . Combining Figures 9b and 11c,d, the values of b and k s /k n can be roughly estimated to satisfy that v c = 0.35 and E c /E t = 1.43. Then k n can be estimated by amplifying initial value in Section 3 to obtain that E c = 68.1 GPa. The next step is to get an appropriate value of ε n 0 so that TS = 9.3 MPa is obtained. The last step is to guarantee that UCS = 200 MPa, which can be achieved by choosing proper values of k m /(k n R), w and ε s 0 /ε n 0 . Proper micro parameters in NET-BPM for LDB granite can be obtained through several simulations and modifications, as shown in Table 6. The macro parameters obtained are presented in Table 5.
As the ultimate axial strain in tensile test is small, a lower loading speed of 0.01 m/s is applied in direct tensile tests. The results show that not only E c , v c , UCS, TS and UCS/TS from the NET-BPM model are in good agreement with the experimental values, which is achieved in [14], but also E t and E c /E t are in good agreement with the experimental values. The stress-strain curves of the LDB granite under uniaxial compression test from our method are similar to those from experiments in [42], as shown in Figure 14. The crack patterns of LDB granite under compression and tension tests are shown in Figure 15. Same method introduced in Section 2.4 is used to present the cracks. Figure 15a shows the cracks at the time when the compressive load reaches its peak. After that, the cracks develop quickly. Most of them are oblique. Figure 15d shows the cracks at time when the tensile load reaches its peak. The cracks extend perpendicular to the loading direction quickly, which results in the fracture of the specimen.  Table 6. Micro parameters in NET-BPM (nonlinear elastic tensile bonded particle model) for LDB granite and the quartzite.
Parameters k n (10 6 N) k s /k n k m /(k n R) its peak. The cracks extend perpendicular to the loading direction quickly, which results in the fracture of the specimen.

Modelling of Quartzite
A rock specimen with the size same as that in [39] is created using the method in [44]. Since the width of the specimen is only 31 mm, Rmax = 0.6 mm and Rmin = 0.2 mm are used for the specimen preparation. The parameters of specimen for the quartzite are listed in Table 4. The properties of the quartzite from experiments in [39] are listed in Table 5. Similar to that in Section 4.1, a set of micro parameters is obtained in order to get a response close to the quartzite, as listed in Table 6.
A loading speed V = 0.01 m/s is used in the direct tension tests. The properties of the specimen from our simulations are listed in Table 5. All the properties from simulations are in good agreement with those from experiments. The stress-strain curves of the quartzite under compression and tension are presented in Figure 16. The stress-strain curves of sample C-Q1 and sample T-Q1 in [39] are presented for comparison. The results show that the stress-strain curves of the specimen modelled by the NET-BPM are in good agreement with those of the quartzite. The stress-strain curve of the quartzite under direct tension tests is nonlinear, which is reproduced well by our model. The crack patterns of the quartzite in the simulations are presented in Figure 17. Most of the cracks under compression are parallel to the loading direction, which is different from that in the simulation of LDB granite. Even though failure pattern of sample C-Q1 is shear failure found in [39], split failure pattern is observed in another sample, C-Q2 under uniaxial compression loading. The cracks of the quartzite form a line perpendicular to the loading direction, which results in failure of the specimen. This is in accord to experimental results in [39].

Modelling of Quartzite
A rock specimen with the size same as that in [39] is created using the method in [44]. Since the width of the specimen is only 31 mm, R max = 0.6 mm and R min = 0.2 mm are used for the specimen preparation. The parameters of specimen for the quartzite are listed in Table 4. The properties of the quartzite from experiments in [39] are listed in Table 5. Similar to that in Section 4.1, a set of micro parameters is obtained in order to get a response close to the quartzite, as listed in Table 6.
A loading speed V = 0.01 m/s is used in the direct tension tests. The properties of the specimen from our simulations are listed in Table 5. All the properties from simulations are in good agreement with those from experiments. The stress-strain curves of the quartzite under compression and tension are presented in Figure 16. The stress-strain curves of sample C-Q1 and sample T-Q1 in [39] are presented for comparison. The results show that the stress-strain curves of the specimen modelled by the NET-BPM are in good agreement with those of the quartzite. The stress-strain curve of the quartzite under direct tension tests is nonlinear, which is reproduced well by our model. The crack patterns of the quartzite in the simulations are presented in Figure 17. Most of the cracks under compression are parallel to the loading direction, which is different from that in the simulation of LDB granite. Even though failure pattern of sample C-Q1 is shear failure found in [39], split failure pattern is observed in another sample, C-Q2 under uniaxial compression loading. The cracks of the quartzite form a line perpendicular to the loading direction, which results in failure of the specimen. This is in accord to experimental results in [39].

Conclusions
In this paper, a new model is proposed for rock-like materials, which is named as NET-BPM model. The NET-BPM model, which adopts nonlinear tensile stiffness and a new form of bonds in its contact model, has the following advantages. (1) It is very difficult to gain a large ratio of UCS/TS more than 20 in conventional BPM model [1,2,20]. However, compared to our model, a large range of UCS/TS ratio (from 5 to 45 in the presented simulations) can be obtained by using different b and ε s 0 /ε n 0 ; (2) The nonlinear behavior of some rocks under direct tension, which is observed in many experiments [8,[36][37][38][39][40], can be reproduced well in our model with a proper b; (3) The differences between Young's modulus, Poisson's ratio of specimens, when they are under compression and tension, which are observed in experiments [38,41,49,50], are captured well in the NET-BPM model. e.g., the value of E c /E t ranging from 1.0 to 1.6 can be obtained; (4) It is possible to change the width of the bond without changing the stiffness of the bond in the NET-BPM model, which helps to obtain different UCS without affecting TS and E c .
The relationships between micro and macro parameters in our model were investigated through simulating uniaxial compression tests and direct tension tests. And the tests were performed upon specimens sharing identical particle arrangement. The results show that there is no significant difference for the effect of k s /k n on rock's strength, Young's modulus and Poisson's ratio between our model and typical DEM model [28]. Larger k s /k n will results in larger E c , UCS, TS and UCS/TS, and smaller v c . Besides k s /k n , k m /(k n R) also has great effects on UCS, TS and UCS/TS, and the effects will become larger with the increasing of k s /k n . It should be noted that the larger k m /(k n R) will results in greater UCS, TS and UCS/TS, while it still has ignorable effects on Young's modulus and Poisson's ratio. In present investigation, the parameter b has great effects on strength, Young's modulus and Poisson's ratio of specimens in tensile tests but small effects in uniaxial compressive tests. In conclusion to this, parameter b has great influences on TS, UCS/TS, E c /E t and v c /v t . However, as the increasing of b, those effects will become smaller. In addition, the effect of the bond width coefficient w on UCS is obvious in the NET-BPM. When w increases, UCS will decrease, but TS, E c and v c will remain unchanged. Moreover, the ultimate shear strain has great effects on UCS, but very small effect on TS, E c , E t , v c and v t . Increasing ultimate shear strain will increase UCS and results in larger UCS/TS.
The LDB granite and a quartzite under uniaxial compression tests and direct tension tests were modelled by the NET-BPM model. All the properties, including E c , E t , v c , UCS, TS and UCS/TS from simulations are in good agreement with those from experiments. In addition, the stress-strain curves of LDB granite and the quartzite are also reproduced correctly. All those results indicate that NET-BPM is a very suitable method for modelling the deformation and fracture of rock-like materials. Some work in the subsequent studies are listed as follows. (1) The effects of particle size and particle size ratio on the simulation results, which are not concerned in this work, need to be investigated. (2) More simulations like biaxial and Brazilian tests need to be performed to validate our model. (3) The 3-D NET-BPM model needs to be established to broaden the application scope of the NET-BPM model.