Research on rock triaxial cyclic loading and unloading characteristics based on PFC3D

. Numerical methods provide a new way for solving complex geotechnical engineering problems. Firstly, we review the emergence and development of numerical methods, and introduce the basic principles of particle flow discrete element and its application in geotechnical engineering. After that, a method of servo combination of cylinder and loading plate is used, the numerical simulation of the three-dimensional rock in triaxial loading and unloading test is carried out by PFC3D software, and the loading and unloading mechanical characteristics of rock mass under the linear elastic contact constitutive model and parallel bonded contact constitutive model are studied respectively, and we obtain the triaxial loading and unloading stress-strain curves of rock under the elastic constitutive model and elastoplastic constitutive model. The residual strains and hysteresis phenomena of the rock under the two constitutive models are analysed, and the feasibility of the particle flow discrete element method in studying the mechanical characteristics of rocks is demonstrated.


Introduction
With the continuous improvement of computer performance, the role of numerical methods in solving geotechnical engineering problems has gradually increased. Rock mass is a kind of geological body with discontinuity, heterogeneity and anisotropy, its characteristics make it difficult to describe the constitutive relationship, and the boundary conditions are often very complicated when mechanical analysis of rock mass is carried out, which causes many difficulties to engineering analysis. The emergence of numerical methods provides a fast and accurate analysis method for rapid analysis of geotechnical engineering problems.
The more commonly used methods to solve geotechnical problems are the finite element method (FEM) and the discrete element method (DEM). The theoretical basis of FEM is the principle of minimum potential energy, which treats the research object as a continuum, that is, each unit displacement maintains a single-valued continuum. The theoretical basis of DEM is Newton's second law, which treats rock masses divided by joints and fissures as rigid bodies. For the discontinuous rock mass existing in nature, the use of DEM for simulation has greater advantages, making DEM an important means of simulating rock mass.
An important discrete method in DEM is particle flow discrete element. One of its representative software is PFC software released by ITASCA company. In recent years, a series of researches have been carried out based on PFC software at home and abroad.

Development of the discrete element method
The discrete element method was proposed by Cundall P.A. in the early 1970s and was initially used to study the mechanical characteristics of discontinuous media such as rock masses. Later, Cundall P.A. and Strack proposed a two-dimensional discrete element simulation method. Currently, experts and scholars from many countries are working to promote the development of discrete elements. The earliest discrete element systems were limited to dealing with discrete blocks of stiffness, after which the methods were expanded. In the 1980 year, ITASCA developed the discrete element UDEC to be put on the market. Later, Lorig and Brady developed a program to deal with discrete element-boundary element coupling, and Cundall P.A. developed 3DEC, a 3D numerical simulation software to simulate rock masses with joints in the 1988 year [1] . China has also conducted a series of researches on discrete element software, and in the China Rock 2019 conference, several experts and scholars gave lectures on discrete element to explain, Liu Chun et al. of Nanjing University developed the highperformance discrete element software MatDEM, which greatly improved the efficiency of the operation.

Application state of particle flow discrete element in simulating rock mechanics test
The mechanism of rock masses failure is of great importance to mechanics theory, underground engineering construction and geological disaster management. The discontinuous media theory divides rock-like material into computational models consisting of rigid particles or blocks, the discrete element method using numerical test methods has become an important breakthrough point in solving the failure mechanism of rock-like media from a mesoscopic perspective. [2] Scholars at home and abroad have carried out many researches and explorations on this issue. The research result of Xu et al. [3] show that the materials mechanical characteristics of complete specimens under compression can be simulated more accurately by using PFC. Liu et al. [4] use PFC to research the time-expansion effect and brittle-ductile-plastic transition character of marble. Huang et al. [5] use particle flow to simulate the macro and meso-mechanical behaviours of red sandstone containing two pre-existing non-coplanar fissures and analysed its strength and failure characteristics.
Based on this, a servo loading method combining loading plate and cylinder is adopted in this paper to study the simulation process of conventional triaxial loading and unloading of rock mass under elastic constructive model and elastoplastic constructive model respectively, and obtain stress-strain curves, and analyse the mechanical behaviour and mechanism of rock mass during the loading and unloading process, and study the mechanical characteristics of rock mass under cyclic loading and unloading conditions, and then verify the feasibility of particle flow discrete element in simulating triaxial test.
3. Basic principle of particle flow discrete element

Particle flow discrete element
After the development of the UDEC program by Cundall P.A., the in-depth research based on the theory of discrete elements continues with the development of the BALL program, which is used to calculate the mechanical characteristics of discrete discontinuous media bodies into particles. The PFC (Particle Flow Code) program, which currently succeeds PFC3D and PFC2D, is based on the discrete element method (DEM) framework and consists of computational engine and graphical operation interface. PFC is mainly used to simulate the movement and interaction of particles of limited size. The particles in the PFC model are rigid bodies with mass, which can be translational and rotational. The interaction between the particles has internal inertial force, moment or paired contact force.
The advantages of particle flow discrete element method are that it can get rid of the limitation of deformation, so it is easy to deal with discontinuous media problems, and can better deal with the problems related to multi-field coupling. For the non-unique deformation problems of rock cracking, separation, etc., the principles and results can be fed back by showing the dynamic process. However, the particle flow discrete element method also has its own shortcomings. Firstly, the software is similar to FLAC3D, the pre-processing of modelling is difficult, and it is not easy for beginners to get started, especially it is very difficult to establish complex models, and secondly, the calibration of software parameters is difficult, and the particle flow method lacks the test of engineering applications.

Mechanics of the discrete element method
for particle flow The movement pattern of the particles in the particle flow discrete element method is basically the same as the principles of kinematics in the DEM, that is, Newton's second law is used to analyse the movement state of a single particle.
In the above equation, F and M are the unbalance force and unbalance moment, and β is the overall damping of the system. Particle flow contact model is the core principle of PFC, and it is important to simulate the contact mode between two particles in PFC. With the development of PFC, the contact model between particles has been diversified, such as rigid body contact, surface tension generated by liquid bridge and viscous force due to the fluid's own sticky. The following is a brief introduction of the parallel bonding model among the common contact models of the particle flow discrete element method. [6] First of all, we know that rock-like materials are friction-like materials, while metals belong to artificially formed crystalline materials, their material characteristics are different, mainly in: ① The shear strength of rocks is related to the cohesion of rock particle materials and the angle of internal friction, while metal materials do not have this characteristic. ② Rock materials are multiphase materials, with the characteristics of solid, water and gas, so they are affected by hydrostatic pressure. ③ Rocks have double strength characteristics, which in turn determines that rock materials have the characteristics of hardening and softening. The difference in material characteristics determines that rock-like materials have many different mechanical characteristics from metallic materials, mainly in the compression hardness, isobaric yield characteristic, dilatancy, hardening and softening characteristics of rock-like materials, etc.
It can be seen that the mechanical characteristics of rock-like materials are very different from those of metallike materials. The characteristics of rock-like materials determine that the shear yield and failure of rocks must take into account the friction angle of the particles inside the rock, so the yield conditions and failure criteria of rock materials are different from those of metal materials.
The Mohr-Coulomb criterion is widely used for rock mass as the yield condition, which reflects the strength characteristics of rocks greatly and is suitable for building the constructive model of geotechnical materials, and is widely used in the geotechnical engineering industry now.
In the above equation, τ is the shear stress (shear strength) on the shear surface; σ is the normal stress on the shear surface; c is the bonding force; and φ is the internal friction angle.
In PFC, when the parallel bonding model is used, as shown in Fig. 1, the stiffness of the material with the parallel bonding model consists of both contact stiffness and bonding stiffness, and the parallel bonding component works together with the linear component to provide elastic interaction between the contact objects, once the contact stress between the particles exceeds the strength limit, the bonding breaks and the bonding is lost, and the model between particles degenerates into a linear contact model. The contact bond model can only transfer force, while the parallel bond model can transfer force and moment, so the parallel bond model can reflect the constructive relationship between the rock particles more accurately when simulating the contact between the rock particles.

Numerical test of particle flow discrete element triaxial loading and unloading 4.1 Principles of mechanics and computation
For the PFC software, the wall can only be given displacement and velocity, and there is no way to give pressure to the wall, so the stress from the triaxial test is converted to the displacement and velocity of the wall. In this model, there exists a cylindrical wall in contact with the edge of the spherical particles, since particles are spheres, then the contact mode between the particles and the wall is point contact, and the lateral displacement of the spherical particles occurs due to the application of axial pressure, and the extrusion of the surrounding wall occurs, and the lines between the contact points of the spherical particles and the wall and the centres of the spherical particles are always horizontal, so the force of the spherical particle on the wall is also horizontal, and let the force of each particle on the wall force is F i , then the sum of the force of each particle is the combined force on the particle. [7]   i F F (4)

Fig. 2. Schematic diagram of loading mechanism
Then the stresses generated at each boundary can be calculated by calculating the average stress of the particles on the wall according to Newton's third law, by dividing the combined force acting on the wall by the area where the interaction between the wall and the particles.
The strain can be calculated by the following equation.
In the above equation, L is the length of the specimen at a moment in the corresponding direction and L o is the original length of the specimen. For the cylindrical specimens in this paper, the following equations can be used to calculate the axial and normal strains.
The kinematic parameters of the wall are controlled to ensure that the axial and lateral stresses remain constant throughout the loading process of the specimen, which is performed by calling the get_gain and servo functions of the FISH language. A cycle is set and the sum of the stresses of each spherical particle on the lateral wall is collected once in each cycle to adjust the wall kinematic parameters, and minimize the difference between measured stress and demanded stress. According to the related literatures, the following algorithm can be used for the servo mechanism.
G in the above equation is the incremental function obtained by using the get_gain function, the maximum increment that causes the wall boundary force within a certain Δt is: In the above equation, N c is the number of particles in contact with the wall, and ΔF (ω) is the average stiffness of each particle contact surface, the average restraint stress is: In the above equation, A is the effective area where the cylindrical wall and the spherical particles are in contact and there is a force. At this point, the difference between the measured value of the restraint stress and the demand value needs to be checked. In practice, it is necessary to use the amplification factor method to derive the following expression.
Substituting equation (11) into equation (12), we can get: Introducing the incremental function G above, increment can be obtained by derivation.
According to the related literatures, we learn that in numerical simulations, ɑ=0.5.

Establishment of particle flow discrete element model
According to the principles introduced in the 4.1, the model can be built. First of all, the domain is created to determine the limits of the area where the model is to be built, and then three walls are set up, respectively the loading and download boards of infinite size and the side confining plate. After that, the particles are generated. First, a region is created using geometry, and then spherical particles are generated in the geometry region by using the ball command. At this point, it will be found that the locations of some particles are beyond the range of geometry and wall, in this case the run is likely to appear that the particles are extruded from the wall, so we can reduce the size of the sphere in the X, Y, and Z directions to make the particles shrink, so as to avoid the bad effect of this situation on the simulation. The model diagram is shown in Fig. 3.   Fig. 3. Schematic diagram of the establishment of PFC model In this model, the height of the model is 40 mm, and the width is 20 mm. In generating the spherical particles, the radius has a size between 0.075 mm and 0.1 mm. 10001 random numbers are used. Fig. 3 shows that 2819 particles were generated, and the line contact model was used between the particles.
After the particles were established, the principles command flow introduced in the 4.1 will be executed and the constrained stress state of the wall can be read, and the loading of the specimen will be performed.

Numerical simulation of elastic material specimen under triaxial loading and unloading
First, multiple plots need to be established to record the state curve of the specimen and the loading plates under the hist. After that, we use the command flow to load the specimen. First of all, the axial pressure is applied, the loading and download boards need to gradually reach the required loading pressure, so it is necessary to use the principles described in the 4.1 to use FISH language to form a circular structure, set the displacement within each step, and feedback the force of the particles on the wall at this time. When the demand pressure is not reached in one step, increasing the average velocity, and gradually approaching the demanded axial pressure. When the program gets the feedback that the axial pressure reaches the demand pressure, exiting the cycle and loading at this speed, thus obtaining a uniform axial loading speed, as shown in Fig. 4, the loading principle of the lateral confining pressure is the same as that of the axial pressure, that is, the cyclic structure is used to achieve the required test pressure, and unloading after loading is completed.
Making the specimen an elastic material, setting the specimen material mechanical parameters are as follows: tensile strength is 1×10 5 Pa, shear strength is 1×10 5 Pa, friction coefficient is 0.5, axial pressure is 2.3 MPa, confining pressure is 1 MPa, then performing triaxial loading and unloading.

Fig. 4. Velocity image of the wall under loading
After the command flow is finished, running the program to load the specimen. During the loading process, the displacement and contact force cloud charts of the particles can be viewed, as shown in Fig. 5,   Fig. 5. Cloud chart of displacement profile at a certain moment during loading After the PFC program is calculated, we can view the process curves during loading and unloading, as shown in Fig. 6, Fig. 7. The full stress-strain curve of the triaxial compression loading and unloading test of the elastic material specimen is shown in Fig. 6. Starting loading, the axial stress increases to the point D, and the axial displacement is 0, indicating that the rock is compact. Unloading is performed after loading to the point E, then the curve goes back to the point F. Unlike the conventional full stress-strain curve of elastic material specimen, the two curves of loading and unloading of the specimen shown in Fig. 6. do not overlap, because the microfractures and open structural planes between the constituent particles of the rock are gradually compacted during the loading process, the feature better reflects the compaction characteristics of the rock material.

Numerical simulation of elastoplastic material specimen under triaxial loading and unloading
Using the above principles, we can create elastoplastic specimen by changing material parameters, setting the specimen material mechanical parameters are as follows: tensile strength is 0.5×10 5 Pa, shear strength is 0.5×10 5 Pa, friction coefficient is 0.5, axial pressure is 3 MPa, confining pressure is 1 MPa, then performing triaxial loading and unloading.
Using the above principle to bring the force to the demand restraint stress, as shown in Fig. 8. After the loading is completed, velocity image of elastoplastic material under triaxial loading, full stressstrain curve of elastoplastic material under triaxial loading and unloading and axial strain curve of elastoplastic material under triaxial loading and unloading can be obtained, as shown in Fig. 9, Fig. 10, Fig. 11.
The full stress-strain curve often triaxial compression loading and unloading test of the elastic-plastic material specimen is shown in Fig. 10. In the early stage of loading, the stress-strain curve is convex upward. In the elastic phase, the stress and strain of rock is basically linear. In the post-peak phase, unloading starts at the point A, it can be seen that the unloading curve does not coincide with the loading curve, and the strain is not returns to 0, i.e., the point B. In this case, the recoverable deformation is called elastic deformation, expressed as ε e , which is the length indicated by BC. The unrecoverable deformation is called plastic deformation or residual deformation, expressed as ε p , i.e., the length from the origin to the point B. The plastic deformation of the specimen is generated and loading to the point A again, it can be seen that the loading and unloading curves of the rock cannot be completely overlapped, and the unloading curve always lags behind the loading curve, forming a closed plastic hysteresis loop with a sharp leaf shape. The area enclosed by the plastic hysteresis loop is the energy absorbed and consumed by the inelastic deformation part of the rock structural planes, which is called plastic work ΔW. The energy is used for adjustment of rock structure, structural planes compaction, or relative slip and movement of the structure. Continuing to load, it can be seen that the full stress-strain curve will have multiple steep falls after the peak, indicating that the rock produces multiple local brittle fractures.

Analysis of triaxial loading and unloading
test results Based on the above numerical simulation tests as well as the analysis results, the following conclusions can be drawn.
(l) The triaxial loading and unloading program can simulate the triaxial test well, and good conclusions can be drawn for elastic and elastoplastic materials, the loading and unloading principles in the 4.1 have good usability, and can be extended to obtain triaxial loading and unloading tests for multiple stress paths and multiple specimens.
(2) For the elastic material, it can be seen from the stress-strain curve that the residual stress after the loading and unloading process is 0. For the elastoplastic material, a large residual deformation is generated during the postpeak loading and unloading process, which is the same phenomenon as the triaxial test results, indicating that the test results are accurate.
(3) During the loading process, it can be found that the lateral displacements of particles at the periphery are larger and the particles closer to the axis have less lateral displacement. For the elastoplastic material, it is difficult to increase axial pressure when it is added to 3.56 MPa, because the material has already undergone plastic failure and the support capacity has been reduced.

Conclusion
The following conclusions can be drawn from the analysis of this paper.
(1) Numerical simulation is one of the important engineering aids under the increasingly complex engineering conditions, and its highly efficient computational capability provides an important reference for the prevention and control of engineering disasters.
(2) Particle Flow Code (PFC), as one of the important branches of discrete element method (DEM), is gradually being widely used in geotechnical engineering because of its ability to provide multiple contact relationships between particles and its powerful computational capability.
(3) Particle flow contact model is the core principle in PFC, with the development of PFC, the contact between particles has been diversified, such as rigid body contact, surface tension generated by liquid bridge and viscous force due to the fluid's own sticky etc.
(4) The triaxial test can be simulated by using PFC. By the algorithm of converting the force into the indenter velocity, the specimen consisting of particle flow is triaxial loaded, by setting the mechanical parameters of the material, the material is set as an elastic material and an elastic-plastic material respectively, and the stressstrain curves of the two materials are obtained by loading respectively. It can be seen from the stress-strain curves of these two materials that the elastic material residual strain is 0, while the elastoplastic material has a larger residual strain after loading and unloading repeatedly. It can be seen that the PFC can better simulate the triaxial test from the numerical simulation results. (5) Using this triaxial test model, it can be extended to realize numerical simulation of various tests by changing the shape of the specimen, such as hollow thin-walled cylinder, thick-walled cylinder, cube, etc., using these principles, some other loading and unloading processes can also be realized.
Regarding the calibration of some parameters, the author will further discuss in the future study. ←