Investigating the effect of unloading on artificial sandstone behaviour using the Discrete Element Method

. The Discrete Element Method (DEM) was used to simulate the mechanical behaviour of a reservoir sandstone. Triaxial tests were carried out using 3D-DEM to simulate the stress-strain behaviour of a sandstone with comparisons made between the numerical tests and the laboratory tests. The influence of isotropic unloading was investigated, which was found to have an impact on bond breakage and was successfully captured in the 3D shearing processes. It was found that bond breakages correlated strongly with the stress-strain behaviour of the sandstone affecting the peak strength. It was also found that unloading affected the bond breakages, which then changed the mechanical behaviour of the sandstone. The tangent stiffnesses of simulated virgin and cored samples under different confining stresses were compared. From the tangent stiffnesses, gross yield envelopes and the yielding surfaces for unloaded samples and virgin samples were plotted and analysed in detail.


Introduction
Many triaxial laboratory tests have been carried out to study the mechanical behaviour of sandstones. In particular, Alvarado et al. [1] developed a series of unloading tests and found that the cementation was damaged by isotropic unloading, representing coring. The difference between the cored and the virgin samples was revealed by the difference in the positions of yielding surfaces. Samples that have never been unloaded to zero confining stress are called 'virgin samples' while samples that have been unloaded to zero confining stress are called 'cored samples'. Using the Discrete Element Method (DEM), a parametric study on Castlegate and Saltwash sandstone was carried out by Cheung et al. [2] calibrating the stiffness and strength of the parallel bonds (Table 2), where in DEM parallel bonds act as a finite size bond, transmitting both a force and a moment over a circular cross-section.
In this study, parallel bonds were applied to the DEM sandstone model with cementation and the data were calibrated with the stress-strain behaviour of the Castlegate sandstone under triaxial shearing. Unloading tests were carried out based on the laboratory tests performed by Alvarado et al. [1]. The gross yield envelopes of both cored and virgin samples were investigated and compared to the envelopes obtained from the laboratory tests. The results of bond breakage caused by unloading are described quantitatively.

Modelling procedures
In the 3D-DEM model, the particles of sandstone were assumed to be spherical. The preparation of the cylindrical artificial sandstone samples was based on the laboratory triaxial tests with a height to diameter ratio of 2. During the shear test, vertical forces were generated by applying a velocity of 0.1m/s to the top and bottom platens, simulating the movement of the platens in a real triaxial test. The average unbalanced force over average contact force ratio was 1.22 x 10 -4 , which was below 10 -3 . Therefore, the system was in "quasi-static" equilibrium during shearing. Parallel bonds were applied at 29MPa for both the cored and virgin samples; this was then followed by the loading and unloading tests (Table 1). Virgin samples were all then unloaded or loaded to 10, 29 and 50 MPa. Parallel bond shear strength, S pb,S,shear, (Pa) 1.5 x10 9 1.5 x10 8 The particle Young's modulus, E c , was set to 20GPa and the damping coefficient was set to zero. A friction parameter μ of 0.7 was used; the input parameters were different from those of Cheung et al. [2] because the sample size and particle sizes were different in this study (Table 2).

Friction angles
The shear stresses and the confining stresses at the peak states have been obtained from the 3D simulations and are plotted in Figure 2. The figure shows that the peak friction angle is about 53.7° which is comparable to the laboratory result. Regardless of this good match of peak strength, real sandstone particles are far more irregular in shape. Due to the fact that all the uncemented DEM particles were spherical and that they will interact only through forces without moment, the internal friction angle of the uncemented DEM samples could be far lower than that of a real soil unless a much larger value ofμ or elongated-shaped particles are used.  Figure 3 describes the stress-strain behaviour of both cored and virgin samples. The DEM samples shearing at higher confining stresses gave higher peak stresses. Samples shearing at the same confining stress showed similar behaviour, such as V10 and C10, regardless of the different unloading. As the confining stress became larger, C50 reached a slightly higher peak than V50.

Stress-strain behaviour and bond breakage
The stress-strain behaviour shown in Figure 3 has a strong relationship with bond breakage. In Figure 4, the initial bond number of the unloaded (cored) samples was found to be lower than that of the virgin samples, showing that unloading damages the bonds between particles. Moreover, the bond numbers of all the cored samples remained the same until an axial strain of 1% was reached, indicating that the initial increasing deviator stress was controlled by the strength of the cementing materials. Sample V50 had the highest bond numbers which might due to the higher confining pressure giving more bonded surfaces around the particles. C10 had the most remaining bonds when compared to the other samples at the strain-softening stage. The bond number of V29 reduced faster than C29 at the second stage and eventually reached the same state after 5% axial strain. With the bonds gradually breaking during shearing, the remaining bond numbers were similar for both the cored and virgin samples under the same confining pressure.

Tangent stiffness
The tangent stiffnesses of both cored and virgin samples are shown in Figures 5 and 6. Unlike Alvarado et al. [1] who used a Y2 method to identify the yield point, the gross yield points in this paper were obtained from the tangent stiffness plots at the end of the plateaus, where there was a sudden drop at around 0.6% axial strain, indicating the start of the major plastic deformation. Figure 4 indicates that these gross yield points happen when significant bond breakage starts. In the laboratory, the gross yield points were also obtained from the stiffness degradation curves. The gross yield point was different from the true yield point as pointed by Alvarado et al. [1], since the true yield points occurred at very small strains below which the behaviour was both linear and elastic. It can be seen from Figures 5 and 6 that below 0.10% axial strain, the tangent stiffnesses of the virgin samples were slightly lower than those of the cored samples and the plateaus were more obvious in the virgin samples than in the cored samples.

Yielding surface
In the laboratory test, Alvarado et al. [1] used the Y2 method, which was defined by Jardine [5], and the yield points were identified when the direction of the strain increment changed. According to these yield points, it was found that the gross yielding envelopes and the yielding surfaces were similar to the laboratory results under relatively higher confining stresses. The stress paths of the three cored and virgin samples tested under 10, 29 and 50MPa confining pressures are shown in Figures 7 and 8 respectively. The critical state points were obtained from the end points of the stress-strain curves in Figure 3. Peak points at the top of the stress paths, form a peak envelope (dotted line) that curves at medium stress levels. The shape of the gross yield envelope for the cored samples is close to elliptical and is similar to the results of the real sandstone [3]. The gross yield envelope for the virgin samples, however, is much larger than that of the cored samples, which is similar to the results in the laboratory. The difference between the two yield envelopes was due to the unloading process, indicating that there were cementation bond breakages during the unloading of the cored samples.   It should be noted that bond breakage is essential to reach the CSL, which was also studied by Cheng et al. [4] who pointed out that particle breakage also had a significant impact on the shearing behaviour of soil. The small amount of bond breakages revealed in the DEM simulations were difficult to identify in the laboratory tests. The laboratory tests also had difficulty in determining the yield points at high confining pressures due to the limitations of the apparatus.

Conclusions
This paper studied the mechanical behaviour of a reservoir sandstone based on the discrete element method. Qualitative comparisons were made between the numerical results and laboratory results mainly referring to the studies by Alvarado et al. and Cheung et al. [1,2]. The results could potentially help engineers to have a better understanding of the shearing behaviour of the reservoir sandstone from a micromechanical perspective.
The main conclusions of this study are: 1) The stress-strain behaviour of sandstones under triaxial shearing was confining pressure dependent. It was found that the microstructural parameters, in this case, the cementation bond breakage during unloading, had a great impact on the shearing behaviour of a sandstone.
2) There was a relationship between bond breakage and the stress-strain behaviour of a sandstone. The parallel bonds started to break when the deviator stresses reached the peak states. Isotropic unloading has been found to affect the behaviour of the sandstone, especially bond breakage. The breakage of bonds also reduced the tangent stiffnesses. However for real sandstone, there might be particle breakages (on top of cementation bond breakage) that further disturb the original structure of the sandstone, leading to a different tangent stiffness, when compared to the DEM data, due to increased locked stresses and frictional forces between the particles.
3) The gross yield envelopes between cored and virgin samples were compared. It was found that the gross yield envelope of the virgin samples was larger than that of the cored samples, both of them enclosed the curing stress (29MPa). It is difficult to tell the difference between the Y2 surface and the gross yield surface of the virgin samples.