Creep Behaviours of Argillaceous Sandstone: An Experimental and Modelling Study

: Understanding the creep behaviours of rocks is essential for the long-term stability of underground excavations in mining engineering. Creep behaviours are more important when the mining depth is greater, which leads to the emergence of weak rock masses and high in situ stresses. In this study, the creep behaviours of argillaceous sandstone (AS) were systematically investigated. For the experimental investigation, creep tests were conducted on AS with di ﬀ erent conﬁning pressures (3, 6, 9, 12, 15, and 18 MPa) using an MTS815.02 rock mechanics test system. The mechanical characteristics of AS were analysed. For the numerical study, a nonlinear creep model of AS under equal and di ﬀ erent conﬁning pressures was established based on rock creep theory and plastic theory. The results showed that conﬁning pressure could e ﬀ ectively improve the creep failure strength of AS, accelerating its creep deformation rate and process and reducing the ﬁnal expansion volume. The nonlinear creep model was embedded in the FLAC3D software, and the experimental and numerical results agreed well. The experimental investigation and proposed creep model can provide important guidance in underground mines for safe long-term stability of underground excavations.


Introduction
The long-term behaviors of rock masses around underground structures have received increasing attention from both academia and industry since these behaviors have significant influence on the structures' stability [1][2][3]. Creep behaviours refer to the continuous re-distribution of stress and displacement that occurs when the shear stress exceeds a threshold. For underground structures in weak rock masses or rock masses under high in situ stresses, creep behaviors are particularly common. Moreover, creep behaviors must be investigated for some underground structures that need to be operated for a long period of time, i.e., thousands to tens of thousands years for repositories for radioactive waste disposal [4].
into cylindrical specimens with a diameter of 50 mm and a height of 100 mm. The parallelism and flatness of the two ends of the processed specimens were controlled to be less than 0.02 mm and 0.5 mm, respectively. Figure 1 shows some AS samples used in this study. The basic physical and mechanical properties of the AS samples are listed in Table 1. The equipment used in the test was an MTS815.02 electro-hydraulic servo rock mechanics test system from the United States. Each index of the system was ensured to meet the requirements of uniaxial, triaxial, and creep tests. The test equipment and loading method of AS samples are shown in Figure 2. The circumferential deformation of the sample was measured by a circumferential extensometer. deform at a high deformation rate. The average deformation rates of the roof subsidence, the convergence of the two sides, and the floor heave were 14.4 mm/d, 12.9 mm/d, and 7.3 mm/d, respectively.
As introduced above, the AS samples were taken from the roof of the belt conveying roadway in Quandian Coal Mine 22. In the AS samples, clay materials account for 35.8%, while sand materials account for 64.2%. After the rock samples were collected in accordance with the International Society for Rock Mechanics (ISRM) recommendations for rock mechanics tests, the samples were processed into cylindrical specimens with a diameter of 50 mm and a height of 100 mm. The parallelism and flatness of the two ends of the processed specimens were controlled to be less than 0.02 mm and 0.5 mm, respectively. Figure 1 shows some AS samples used in this study. The basic physical and mechanical properties of the AS samples are listed in Table 1. The equipment used in the test was an MTS815.02 electro-hydraulic servo rock mechanics test system from the United States. Each index of the system was ensured to meet the requirements of uniaxial, triaxial, and creep tests. The test equipment and loading method of AS samples are shown in Figure 2. The circumferential deformation of the sample was measured by a circumferential extensometer.   Note: ρ represents density, E and µ represent Young's modulus and Poisson's ratio, respectively, σ c and σ t represents the uniaxial compressive strength and uniaxial tensile strength, c and ϕ represent cohesion and internal friction angle, respectively.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 3 of 23 for Rock Mechanics (ISRM) recommendations for rock mechanics tests, the samples were processed into cylindrical specimens with a diameter of 50 mm and a height of 100 mm. The parallelism and flatness of the two ends of the processed specimens were controlled to be less than 0.02 mm and 0.5 mm, respectively. Figure 1 shows some AS samples used in this study. The basic physical and mechanical properties of the AS samples are listed in Table 1. The equipment used in the test was an MTS815.02 electro-hydraulic servo rock mechanics test system from the United States. Each index of the system was ensured to meet the requirements of uniaxial, triaxial, and creep tests. The test equipment and loading method of AS samples are shown in Figure 2. The circumferential deformation of the sample was measured by a circumferential extensometer.

Experimental Scenarios
In this study, the creep characteristics of AS specimens under various confining pressures were investigated by the multi-stage loading method. The specific process was as follows: First, the compressive deviatoric strength of the specimens under different confining pressures was tested. Second, according to the compressive deviatoric strength, the specimens were subjected to multi-stage loading. Five-stage loading was employed for each confining pressure, and the loading interval for each stage was 48 h. Finally, the multi-stage loading curve was processed to obtain the stress-strain and mechanical parameters of specimens.
The selection of confining pressures was based on the engineering conditions. The depths of the underground roadway were 644.5-814.7 m. The average minimum in situ stress in the stability zone of the roadway was determined to be 18.23 MPa. Therefore, in order to analyse the creep behaviours in different positions along the roadway, the confining pressure C was set to 3 MPa, 6 MPa, 9 MPa, 12 MPa, 15 MPa, and 18 MPa.
The setup for the multi-stage loading was also determined carefully. The results of triaxial tests showed that the maximum deviatoric stress corresponding to the elastic deformation of specimens under different confining pressures was approximately 60% of the deviatoric strength σ c . Under different confining pressures, the first-stage loading σ 1 was determined to be 55-59% of the deviatoric strength.
The second-to fifth-stage loadings were increased progressively considering the first-stage loading, and each time, the loading was increased by 4 MPa. If the specimens under the fifth-stage load were not damaged, further loading was carried out with the same increment. If the tertiary creep stage occurred before the fifth-stage load was applied and the AS specimens were destroyed, no further loading was applied afterwards.
In the current study, the loading time per stage was determined to be 48 h. Several previous tests determined that the deformation rate of specimens after 24 h was less than 0.002 m/h at each stage when the tertiary creep stage was not reached. Thus, 48 h was considered long enough to reach a stable creep stage.
According to the results of triaxial compression tests, the creep test scenarios of AS samples under different confining pressures are shown in Table 2. In the table, σ i (i = I, II, . . . , V) represents creep loadings at different stages, and δ represents the ratio of creep loading σ i to deviatoric strength σ c .

Strain-Time Curves and Creep Characteristic Curves of AS
The strain-time curves of AS samples under different confining pressures were obtained by multi-stage loading tests. According to the principle of the multi-stage loading test, the strain-time curves for each loading stage could be obtained from the strain-time curves. The influence of confining pressures on the strain-time curves and strain-time curves for each loading stage of AS samples is shown in Figure 3. The following conclusions can be drawn from

1.
Regardless of the confining pressure, the deformation of AS samples at the initial loading stage increased with time, although at a decelerating rate. This stage is known as the primary creep stage. With continuous loading, the strain of the AS samples exhibited a linear response with time, and this stage is known as the secondary creep stage. At the final loading stage, which is known as the tertiary stage, the strain increased exponentially with time until the failure of samples. All AS samples exhibited three complete creep stages.

2.
For creep tests with the same confining pressure, a certain amount of instantaneous deformation was observed when the stress loading at each stage was applied. Moreover, the deformation magnitude increased with increasing stress level. Instantaneous strain was the main component of the total strain. 3.
With increasing confining pressure, the failure load at the tertiary creep stage increased gradually. However, the failure load from the creep tests was less than the corresponding failure load from the conventional triaxial compressive tests.
In summary, Figure 3 shows that all AS samples exhibited four complete creep stages, namely, instantaneous deformation stage, primary creep stage, secondary creep stage and tertiary creep stage. Figure 4 shows a comparison of the strain-time curves at the final loading stage for all confining pressures. The confining pressure had significant influence on the strain-time curves at the tertiary creep stage. With increasing confining pressure, the duration of the tertiary creep stage decreased. In other words, the AS samples were destroyed after a shorter time when the confining pressure was increased.   is provided as an illustration, and similar results could be observed at other loading stages. As shown, the creep deformation rates of AS samples under various confining pressures followed a decreasing-stable-increasing pattern. Detailed conclusions are drawn from Figure 5 as follows: 1.
Based on different confining pressures, the deformation rates of AS samples showed a fast decrease at the beginning. The deformation rate became stable within 7 to 9 h. 2.
With increasing confining pressure, the average deformation rates in the stable stage were similar. However, evident differences were observed in the duration of the stable stage, which decreased gradually with increasing confining pressure. For example, when C = 6 MPa, the stable stage lasted approximately 30 h, and the average creep rate was 0.112 × 10 −3 mm/h. When C = 12 MPa, the stable stage lasted only approximately 5 h, and the average creep rate was 0.08 × 10 −3 mm/h. An interesting phenomenon was that the stable stage disappeared when C = 18 MPa. 3.
The influence of the confining pressure on the deformation rate at the acceleration stage was also evident. With increasing confining pressure, the duration of the acceleration stage decreased gradually. Meanwhile, the acceleration rate increased sharply. The nonlinear acceleration stage was due to the expansion of the internal cracks. At high confining pressure, the axial load of the AS samples was high, leading to a high speed of crack expansion. The crack expansion in turn accelerated the creep failure. Thus, the nonlinear acceleration stage of the deformation rate was observed. Figure 5 shows the creep deformation rates of AS samples under different confining pressures. The authors note that only the creep deformation rate at the final loading stage is provided as an illustration, and similar results could be observed at other loading stages. As shown, the creep deformation rates of AS samples under various confining pressures followed a decreasing-stableincreasing pattern. Detailed conclusions are drawn from Figure 5 as follows:

Creep Deformation Rate
1. Based on different confining pressures, the deformation rates of AS samples showed a fast decrease at the beginning. The deformation rate became stable within 7 to 9 h. 2. With increasing confining pressure, the average deformation rates in the stable stage were similar. However, evident differences were observed in the duration of the stable stage, which decreased gradually with increasing confining pressure. For example, when C = 6 MPa, the stable stage lasted approximately 30 h, and the average creep rate was 0.112 × 10 −3 mm/h. When C = 12 MPa, the stable stage lasted only approximately 5 h, and the average creep rate was 0.08 × 10 −3 mm/h. An interesting phenomenon was that the stable stage disappeared when C = 18 MPa. 3. The influence of the confining pressure on the deformation rate at the acceleration stage was also evident. With increasing confining pressure, the duration of the acceleration stage decreased gradually. Meanwhile, the acceleration rate increased sharply. The nonlinear acceleration stage was due to the expansion of the internal cracks. At high confining pressure, the axial load of the AS samples was high, leading to a high speed of crack expansion. The crack expansion in turn accelerated the creep failure. Thus, the nonlinear acceleration stage of the deformation rate was observed.

Volume Expansion of AS Samples
During the tests, the circumferential strains of AS samples under different confining pressures were obtained by a circumferential extensometer. Combined with the axial compressive strain data, the variation in the volumetric strain with time was obtained, as shown in Figure 6. Figure 6 illustrates that the volumetric strains of the AS samples under different confining pressures were basically similar. For each volumetric strain curve, three stages were observed: accelerated volumetric strain stage, stable volumetric strain stage, and dilatant volumetric strain stage. When the stress level was lower than 70% of the deviatoric strength, the volumetric

Volume Expansion of AS Samples
During the tests, the circumferential strains of AS samples under different confining pressures were obtained by a circumferential extensometer. Combined with the axial compressive strain data, the variation in the volumetric strain with time was obtained, as shown in Figure 6. basically maintained at a stable value. When the deviatoric stress level reached 80-85% of the deviatoric strength, the microcracks in AS samples began to expand and gradually penetrate. At this point, the specimen entered the dilatant volumetric strain stage. Moreover, the expansion rate increased gradually with increasing time until AS samples were destroyed. Due to the confining pressure constraint, the final volumetric strain of the AS samples gradually decreased with increasing confining pressure.   Figure 6 illustrates that the volumetric strains of the AS samples under different confining pressures were basically similar. For each volumetric strain curve, three stages were observed: accelerated volumetric strain stage, stable volumetric strain stage, and dilatant volumetric strain stage. When the stress level was lower than 70% of the deviatoric strength, the volumetric deformation of AS samples was mainly compressive. With the gradual increase of stress level, when it reached approximately 75% of the deviatoric strength, the volume deformation of AS samples was basically maintained at a stable value. When the deviatoric stress level reached 80-85% of the deviatoric strength, the microcracks in AS samples began to expand and gradually penetrate. At this point, the specimen entered the dilatant volumetric strain stage. Moreover, the expansion rate increased gradually with increasing time until AS samples were destroyed. Due to the confining pressure constraint, the final volumetric strain of the AS samples gradually decreased with increasing confining pressure. Figure 7 shows the variations in creep failure strength with confining pressure. The figure displays that with increasing confining pressure, the creep failure strength increased linearly. When the confining pressure was increased from 3 MPa to 18 MPa, the creep failure strength increased from 33.7 MPa to 56.8 MPa, which was a 68.55% increase. Thus, the confining pressure could effectively improve the creep failure strength of AS.

Creep Failure Strength
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 23 33.7 MPa to 56.8 MPa, which was a 68.55% increase. Thus, the confining pressure could effectively improve the creep failure strength of AS.

Establishment of Nonlinear Creep Model
In this paper, the establishment of the nonlinear creep model was conducted by a combination of elemental components. As discussed above, the creep results revealed four stages during the creep tests: instantaneous elastic deformation, primary creep stage, secondary creep stage, and tertiary creep stage. Among them, the instantaneous elastic deformation and primary creep stages could be described by a Hook elastic element ( Figure 8I) and Kelvin viscoelastic model ( Figure 8II), respectively. When the stress level reached a certain threshold, the surrounding rock began to enter the secondary creep stage, which could be described by the Bingham model ( Figure 8III) (Liu et al. 2017). When the stress level was higher than a certain threshold, the creep curves showed the characteristics of accelerated deformation. For conventional rheological elements to accurately describe the accelerated creep stage of AS was difficult. In this paper, the nonlinear volume plasticity Bingham (NVPB) model ( Figure 8IV) was introduced to describe the tertiary creep stage. Based on the above analysis, this paper proposes a nonlinear viscoelastic-plastic creep model, as shown in Figure 8, which would be suitable for the description of AS.

Establishment of Nonlinear Creep Model
In this paper, the establishment of the nonlinear creep model was conducted by a combination of elemental components. As discussed above, the creep results revealed four stages during the creep tests: instantaneous elastic deformation, primary creep stage, secondary creep stage, and tertiary creep stage. Among them, the instantaneous elastic deformation and primary creep stages could be described by a Hook elastic element ( Figure 8I) and Kelvin viscoelastic model ( Figure 8II), respectively. When the stress level reached a certain threshold, the surrounding rock began to enter the secondary creep stage, which could be described by the Bingham model ( Figure 8III) (Liu et al. 2017). When the stress level was higher than a certain threshold, the creep curves showed the characteristics of accelerated deformation. For conventional rheological elements to accurately describe the accelerated creep stage of AS was difficult. In this paper, the nonlinear volume plasticity Bingham (NVPB) model ( Figure 8IV) was introduced to describe the tertiary creep stage. Based on the above analysis, this paper proposes a nonlinear viscoelastic-plastic creep model, as shown in Figure 8, which would be suitable for the description of AS. 2017). When the stress level was higher than a certain threshold, the creep curves showed the characteristics of accelerated deformation. For conventional rheological elements to accurately describe the accelerated creep stage of AS was difficult. In this paper, the nonlinear volume plasticity Bingham (NVPB) model ( Figure 8IV) was introduced to describe the tertiary creep stage. Based on the above analysis, this paper proposes a nonlinear viscoelastic-plastic creep model, as shown in Figure 8, which would be suitable for the description of AS.  In the above models, E K 1 and E K 2 represent the elastic moduli of elastic elements in Hook and Kelvin-Voigt bodies, respectively; η 1 , η 2 , and η 3 represent the viscous coefficients of viscous elements in Kelvin, Bingham, and NVPB bodies, respectively; N is the rheological index in the NVPB body, indicating the accelerated rheological rate of rock, n > 0; σ s 1 and σ s 2 represent the crack initiation stress and the expansion stress of the rock, respectively.
In Figure 8, the physical meaning and constitutive relation of the first three parts are clear and are not discussed in detail in this study. Interested readers can referred to published papers for more detailed information [26,27]. For the NVPB model, the viscosity coefficient η 3 is assumed to satisfy the following equation.
where t 0 is the relative reference time with a value of 1 and η c is the equivalent viscosity coefficient of the NVPB model, which is closely related to the properties of rock materials. Then, the constitutive equation of the NVPB model is The corresponding constitutive equation is where H(σ 0 − σ s 2 ) is a plastic switching function, which can be expressed as follows: Equation (4) shows that when σ 0 ≤ σ s 2 , the model can reflect the instantaneous deformation, primary creep stage and secondary creep stage of AS; when σ 0 > σ s 2 , the model can fully represent the tertiary creep stage. Figure 9 shows the strain-time curves of the NVPB model under different rheological indices n. As the rheological index n increased, the deformation rate of the AS samples increased gradually. This property could be used to describe the accelerated creep of AS samples under different confining pressures. According to the series-parallel relationship between components, the following equations can be obtained.  When σ < σ s 1 , the model degenerates into a Hook body and Kelvin elements in series, and its constitutive relationship is as follows: When σ s 1 ≤ σ < σ s 2 , the model is mainly composed of Hook, Kelvin, and Bingham elements. The constitutive relationship is as follows: When σ ≥ σ s 2 , the model contains all four elements, and its constitutive relationship is as follows: The simultaneous Equations (5)-(7) are used to obtain the constitutive equation of the nonlinear creep model of AS, as follows: According to Equation (8), the nonlinear creep equation of AS can be obtained as follows:

Nonlinear Creep Model in Three Dimensions
Equations (8) and (9) are actually one-dimensional creep equations for AS. For underground roadways, a one-dimensional constitutive model cannot accurately describe the stress-strain distribution of the surrounding rock masses. Therefore, the model needs to be expanded to three dimensions so that it can be used in 3D stress/displacement analysis.
When the rock is loaded in three dimensions, the stress state of any point in the rock sample can be divided into a stress spherical tensor and a stress deviatoric tensor. The stress spherical tensor can influence only the volume of the cell, and the stress deviatoric tensor can influence the shape of the cell. The stress tensor and strain tensor can be expressed as follows: where σ m and S ij represent the stress spherical tensor and stress deviatoric tensor, respectively, and ε m and e ij represent the strain spherical tensor and strain deviatoric tensor, respectively. σ m , S ij and ε m can be calculated as follows: For the elastomer, the Hook element expression in the three-dimensional force state can be derived from the Hook body expression in the one-dimensional force state as follows: where K and G are the bulk modulus and shear modulus of the rock mass, respectively, as follows: The constitutive relation of the elastomer can be obtained by introducing Equation (11) into Equation (12): Assume that rock is a homogeneous material, and that creep deformation is caused only by the stress deviatoric tensor. In this case, the creep equation of a three-dimensional viscoelastic body can be obtained by directly replacing σ 0 in a one-dimensional creep equation with (S ij ) 0 , in which (S ij ) 0 is the constant deviatoric stress imposed on the rock sample during creep tests. Therefore, the three-dimensional constitutive relationship of the viscoelastic body is as follows: For a viscoplastic body, the creep behaviors are mainly reflected in the shear deformation. Thus, the constitutive relation under the three-dimensional stress state can be expressed as (Brotóns et al., 2014): where g and g 0 are the initial values of the yield function and yield function, respectively; F is the plastic potential function; and φ( g g 0 ) is the model switching function under three-dimensional conditions, which satisfies the conditions: Based on the above analysis, substitute Equations (13)-(15) into Equation (9); thus, the nonlinear viscoelastic-plastic creep equation of AS under a three-dimensional stress state is obtained: In the current study, the intermediate and minimum principal stresses were equal, i.e., σ 2 = σ 3 . Therefore, the stress spherical tensor and the stress deviatoric tensor can be expressed as follows: Introducing Equation (19) into Equation (18), the nonlinear creep model for AS with σ 2 = σ 3 is obtained as follows:

Cracking Stress σ s 1 and Creep Dilatancy Stress σ s 2
According to the creep test results, when the loading stress reached approximately 80% of the deviatoric stress before failure, the volume of AS samples began to change from compressive to expansive. Therefore, 80% of the loading stress was taken as the value of σ s 1 . When the loading stress reached approximately 85% of the final loading stress in the expansion stage, the sample began to expand. Therefore, 85% of the final loading stress in the expansion stage was taken as σ s 2 value.
The values of creep parameters σ s 1 and σ s 2 under different confining pressure conditions are shown in Table 3. Note: σ ve represents the stress when the volume changed from compressive to expansive, σ ds represents the stress level in the volume expansion stage.

Determination of the Remaining Creep Parameters
When σ C < σ s 1 , the parameters to be fitted include K, G 0 , G 1 and η 1 ; when σ s 1 ≤ σ C < σ s 2 , the parameters to be fitted are K, G 0 , G 1 , η 1 and η 2 ; when σ C ≥ σ s 2 , the creep parameters to be fitted are K, G 0 , G 1 , η 1 , η 2 and the rheological index n. In this study, the least squares method was used to fit the parameters. First, the fitting function was determined.
Use b 1 to b 7 to represent the approximate values of parameters K, G 0 , G 1 , η 1 , η 2 , η c , and n and define: By expanding the Taylor series and extracting only the linear part, we can obtain: Then, the following equations can be obtained: According to Equation (24), the objective function of the least squares fitting problem is obtained: If ∂Q ∂δ k = 0, (k = 1, 2, 3, 4, 5, 6, 7) (26) Then 7 j=1 a k j δ j = C k , (k = 1, 2, 3, 4, 5, 6, 7) where a k j and C k in Equation (27) are defined as follows: By solving Equation (28), the δ j values can be obtained, which yield the approximate values of the creep parameters. Equation (28) can be solved as follows: • Arbitrarily select initial values of a set of coefficients b 1 to b 7 ; • Calculate a k j , C k ; • Substitute a k j and C k into (2 − 46) (28) to obtain δ 1 to δ 7 ; • Take b 1 + δ 1 , b 2 + δ 2 , . . . , b 7 + δ 7 as the approximate values of the regression coefficients; The nonlinear creep model parameters of AS samples under different confining pressures and loading levels were obtained by the above methods. As shown in Table 4, the model could reflect the creep behaviors of AS. A comparison between the experimental data and the fitting value is shown in Figure 10.

Finite Difference Form of the Nonlinear Creep Model
According to the components of the nonlinear model, the total deviatoric stress and the total deviatoric strain tensors can be expressed as follows: where e H ij , e K ij , e B ij and e N ij represent the partial strains of Hook, Kelvin, Bingham, and NVPB, respectively; S K ij , S H ij , S B ij and S N ij represent the deviatoric (partial) stress of Hook, Kelvin, Bingham, and NVPB, respectively.
According to the superposition principle, the total deviatoric strain rate of the nonlinear creep model is as follows: The Hook element deviatoric stress is expressed as The stress-strain relationship of the combined model of the Hook element and Kelvin element is as follows: The strain rate of the Bingham viscoplastic element in three-dimensional form is as follows: Furthermore, the deviatoric strain rate of the Bingham viscoplastic element is as follows: Based on the assumption in plastic mechanics, the plastic deformation of rock would not change the volume, and the volumetric stress results only from elastic volumetric strain; thus, the following equation can be obtained: where K B and K H both represent the bulk modulus of the rock.
Based on the principle of centre difference, Equation (33) can be transformed into where ∆t represents a time increment step and S ij and e ij represent the average deviating stress and bias strain in a time increment, respectively. Among these parameters, where the superscript symbols N and O represent the new and old values in a computing time step, respectively.
Substituting Equation (38) into Equation (37), the following equation can be obtained: The first part of the spherical strain expression of Equation (39) can be calculated as follows: By substituting Equations (37) and (40) into Equation (31), the following results can be obtained: When σ < σ s 1 , When σ s 1 ≤ σ < σ s 2 ,

Numerical Implementation in FLAC3D
FLAC3D numerical simulation of the nonlinear creep model of AS was achieved using the C++ programming language to connect it to the secondary development interface of FLAC3D. Its process included main three steps: modifying the header file (.h file), modifying the program source file (.cpp file) and generating the dynamic connection or linking library file (.dll format). The primary purpose of the above modifications was to ensure that the embedded files and parameters were correctly identified by the FLAC3D software. In the current secondary development process, the model header file was named HKBN.h, the file was named HKBN, and the model ID was changed to 336.
(2) Modifying the source file (.cpp file) At this stage, the Properties function, the Strata function, the Getproperty function, the Initialize function, and the Run function were modified so that essential functions, such as identification and initialization, could be used. Its main contents included the following steps: Modify the model structure: The constitutive model is actually an empty function. The main purpose is to assign the initial value to the derived new member in the custom HKBN.h header file. In general, the initial value is set to 0.
Customize the Properties function: This function converts creep parameters in the nonlinear creep model into a series of strings. During FLAC3D operation, new parameters can be assigned to the relevant definitions through the string.
Modify the Strata function, Getproperty function and SetProperty function: The Strata function is the unit status indicator, which realizes the precise positioning of the unit state. The main function of the Getproperty function and SetProperty function is to realize the assignment of creep parameters to the model parameters in the derived class.
Modify the Initialize function: The main function of this function is to initialize the unit status indicator and related parameter assignments and assign values to the private variables of the custom nonlinear creep model defined in the derived class.
Modify the Run function: The main purpose of this function is to realize the cumulative calculation process of defining a new model; that is, to calculate the strain increment in each cycle and convert it into a stress increment. Then, the new stress value of the calculation unit body is superimposed.
Modify the Save.Restore function: The calculation results can be saved after this modification. (3) Generating dynamic link library (.dll format) file of nonlinear creep model The dynamic link library is a necessary part of the normal call of the nonlinear creep model. The main steps are as follows. First, the user defined map (UDM) folder is named HKBN and copied to the FLAC3D installation location. Then, the Visual Studio software is run, and the udm.veproj file is opened. Second, the UDM command is executed to modify the output file named Debug/HKBN.dll. The file is a dynamic connection library of the nonlinear creep model. Finally, the solution command is run to check whether the dynamic library file has been successfully exported. If the export is successful, it is copied to the FLAC3D installation location and invoked directly during simulation. Otherwise, it is checked and modified, and the above steps are repeated after modification. The entire process is shown in Figure 11.

Numerical Model Verification
To verify the accuracy of the secondary development of the nonlinear creep model, a numerical calculation model that was consistent with the creep test of AS was established. The size of the simulation model was a 50 mm × 100 mm cylinder with a total of 14,280 units and 14,686 nodes. A vertical uniform stress was applied at the top of the model, and a fixed circumferential stress was applied on both sides to simulate the confining pressure. A displacement constraint was applied at the bottom in the vertical direction. Displacement constraints, specific numerical models and boundary conditions are shown in Figure 12a. The parameters of the model were the material parameters of AS. In this simulation, the axial load was applied by means of fixed confining pressure and stepwise loading. The confining pressure was set to 18 MPa, and the creep time of each stage was 48 h. The axial load deviatoric stresses of successive stages were 39.6 MPa, 43.9 MPa, 48.2 MPa, 52.5 MPa, and 56.8 MPa.
In the numerical simulation, the rock parameters were taken as the average values from the conventional triaxial compression tests. The variations in the specimen axial displacement with time under different loading stresses are shown in Figure 12b.

Conclusions
In this study, the creep behaviors of AS were investigated both experimentally and numerically. Multi-stage creep tests were conducted on AS samples under different confining pressures (3,6,9,12,15,and 18 MPa). The creep characteristics of AS were analyzed. Based on the experimental results, a new nonlinear creep model was proposed. The nonlinear creep model was implemented in the FLAC3D software, and the procedure for the model parameter determination was presented. A comparison between the numerical and experimental results was also conducted for numerical modelling validation. Based on the above results, the following conclusions can be drawn: not only related to the confining pressure, but also closely affiliated with factors such as loading rate and initial damage degree etc. In this paper, as the dominant factor in the engineering case, the confining pressure is used as a variable to carry out creep test research, which has certain limitations and will be supplemented in the future.

Conclusions
In this study, the creep behaviors of AS were investigated both experimentally and numerically. Multi-stage creep tests were conducted on AS samples under different confining pressures (3,6,9,12,15,and 18 MPa). The creep characteristics of AS were analyzed. Based on the experimental results, a new nonlinear creep model was proposed. The nonlinear creep model was implemented in the FLAC3D software, and the procedure for the model parameter determination was presented. A comparison between the numerical and experimental results was also conducted for numerical modelling validation. Based on the above results, the following conclusions can be drawn:

1.
During the creep tests, the AS samples under various confining pressures exhibited four stages: instantaneous deformation, primary creep, secondary creep, and tertiary creep. With increasing confining pressure, the acceleration rate of the AS samples increased, and the creep failure strength gradually increased. However, due to the confining effect of confining pressure, the final volume gradually decreased.

2.
The nonlinear creep model of AS was constructed by adding another element, the NVPB element, to the traditional Nishihara model. The three-dimensional form of the creep model was deduced to simulate the rock creep problem under the equal confining pressure of surrounding rock masses. The parameters of the creep model could be obtained by fitting the least squares method.

3.
Based on the principle of finite difference, the finite difference form of the nonlinear creep model of AS was derived, and the model was embedded in the FLAC3D software. The numerical calculation model of AS was established. The creep behaviors of specimens under a confining pressure of 18 MPa were calculated and verified with experimental data. A good agreement was observed between the numerical and experimental results, indicating the feasibility of the creep model. 4.
The creep and rupture process of rock is a process of damage accumulation essentially, which is not only related to the confining pressure, but also closely affiliated with factors such as loading rate and initial damage degree etc. In this paper, as the dominant factor in the engineering case, the confining pressure is used as a variable to carry out creep test research, which has certain limitations and will be supplemented in the future.