{'fraction_non_alphanumeric': 0.024967354400626796, 'fraction_numerical': 0.008461739357534604, 'mean_word_length': 5.166183574879227, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Closure behaviour of dislocated rock fractures subject to cyclic normal loadings

Understanding the closure behaviour of rock fractures is of great importance for various rock mechanics and rock engineering issues such as geothermal energy extraction, carbon dioxide capture and sequestration, and geological disposal of radioactive waste. The stress acting on a single bridging asperity results in the surface deformation of itself and the surrounding areas. The aggregate of such deformations exerted by all contacting asperities leads to the total deformation of a fracture. Therefore, the deformation behaviour of a stressed fracture is governed by the mechanical properties of the intact rock (elastic modulus and strength) and the geometrical characteristics of asperities (spatial distribution and local shapes). A contacting asperity may fail when the local stresses have reached some failure criterions, which is responsible to the hysteresis observed in cyclic loading experiments. In the present study, a high-precision contact model that used Boussinesq’s solution to relate contact stresses to surface displacements in two half-spaces was employed to calculate the deformation of rock fractures subject to normal loading-unloading cycles. The normal stress-displacement curves of cyclic loading experiment are in good agreement with the calculation results, and the measured surface damage area fits well the calculated plastic area, which efficiently verifies the validity of the numerical calculation model. An extended parametric study that considered more patterns of dislocation degree was implemented and the numerical calculation results show that the normal stress-displacement curves of different loading cycles show obvious hysterical behaviour. The hysteresis indicated that damages happened on partial contacting asperities, which tended to enhance the closure of a fracture. The hysteresis generated by a loading cycle is typically less than the previous one, showing a gradual weakening trend. More specifically, the damage of the contacted asperities accumulates successively with the increase of loading cycles, and the damage tends to cease after several cycles. The damage-induced permanent change of surface geometry varies the mechanical behaviour of fractures. As the number of loading cycle increases, the normal stiffness of the rock fracture also tends to increase, and eventually approaches a stable value after 4∼5 loading cycles. The numerical calculation provides quantitative estimations on the deformation and asperity damages, which are direct evidence for the hysteretic behaviour of rock fractures subject to cyclic loadings.

In the present study, a high-precision contact model that used Boussinesq's solution to relate contact stresses to surface displacements in two half-spaces was employed to calculate the deformation of rock fractures subject to normal loading-unloading cycles. The normal stressdisplacement curves of cyclic loading experiment are in good agreement with the calculation results, and the measured surface damage area fits well the calculated plastic area, which efficiently verifies the validity of the numerical calculation model. An extended parametric study that considered more patterns of dislocation degree was implemented and the numerical calculation results show that the normal stress-displacement curves of different loading cycles show obvious hysterical behaviour. The hysteresis indicated that damages happened on partial contacting asperities, which tended to enhance the closure of a fracture. The hysteresis generated by a loading cycle is typically less than the previous one, showing a gradual weakening trend. More specifically, the damage of the contacted asperities accumulates successively with the increase of loading cycles, and the damage tends to cease after several cycles. The damage-induced permanent change of surface geometry varies the mechanical behaviour of fractures. As the number of loading cycle increases, the normal stiffness of the rock fracture also tends to increase, and eventually approaches a stable value after 4~5 loading cycles. The numerical calculation provides quantitative estimations on the deformation and asperity damages, which are direct evidence for the hysteretic behaviour of rock fractures subject to cyclic loadings.

Introduction
It is an important topic in rock mechanics and rock engineering to study the closure behaviour and deformability of rock fractures subject to normal stresses [1]. The closure of fractures alters stress fields and aperture distributions, which significantly affect the mechanical, chemical, mass transport and their coupled processes [2][3][4]. In natural circumstances, the fracture surfaces between adjacent rocks are typically not completely mated caused by various geological processes. A lot of engineering practices show that the rock fractures are typically subject to varying stresses, such as the excavation and reinforcement of high slopes, the construction of underground caverns and the effect of blasting stress waves on tunnel support [5]. Therefore, rock fractures often undergo cyclic loading and unloading during the construction process and subsequent operation stages, and their deformation behaviour and mechanical properties need to be thoroughly studied.
Over the past decades, experimental and theoretical studies on the fracture deformation subject to normal stresses have been extensively conducted. In the early stage, Goodman et al. [6] and Bandis et al. [7] established the empirical relationship between normal stress and normal displacement based on experimental data. Goodman et al. [6] conducted the fracture normal deformation test and realized that with the increase of the normal stress, the deformation of rock fractures changes nonlinearly, which can be fitted by a hyperbolic model. Barton et al. [8] and Bandis et al. [7] further analysed the experimental data and then modified Goodman model to establish a new hyperbolic model (i.e., Barton-Bandis model). Subsequently, models based on the Hertzian elastic contact theory have been vigorously developed, of which the Greenwood-Williamson (GW) model is the most representative one [9]. The Greenwood-Williamson model can be used to quantitatively analyse the elastic contact problem between random rough surfaces. Brown and Scholz [10] extended the GW model to the contact of rough rock fracture surfaces and verified them through experiments. Hopkins [11] proposed a numerical simulation method based on the balance of contact force between two surface asperities, considering the impact of a single contact on the surrounding area and using cylinders with different heights to represent the geometric morphology of the asperities. The method using Boussinesq's solution to relate the contact stress to surface displacements in two half-spaces has been proved to be a physically effective algorithm. It considers the local plastic damage, and has been applied to the simulation of rock fracture compressive deformation [12][13].
Most of the previous models did not test the validity against fracture closure experimental data under the condition of cyclic loading. In this study, a high-precision contact model that used Boussinesq's solution to relate contact stresses to surface displacements in two half-spaces was employed to calculate the deformation of rock fractures subject to normal loading-unloading cycles. Then, a cyclic normal loading test was conducted to verify this model from both aspects of normal stress-displacement curve comparison and asperity damage analysis. The mechanism for hysteretic behaviour of rock fractures subject to normal loadings was investigated.

MethodologyA high-resolution contact model for fracture deformation analysis
Using the variational principle of contact mechanics, the fracture closure behavior subject to normal stress is simulated based on the minimum potential energy theory [9]. In this research, the contact model will be verified by a cyclic loading and unloading experiment. Here, the key points of this contact model are briefly presented.
Assuming that the asperity of the fracture is a kind of completely elastic-plastic material and the plastic deformation zone is limited to the local area, the elastic contact model is able to be applied to the elastic-plastic contact analysis. In this condition, the energy dissipation derived by plastic deformation is accumulated to the total complementary potential energy, which can be written as [14] * = 1 2 Ω n ( , ) where Ω is the surface contact area, n is the normal stress, ( , ) is the composite surface normal displacement at the contact area, * ( , ) is the total displacement of the two contacting asperities inside the assumed contact area based on geometrical interference, and ∆ is the composite > . The relationship between normal displacement and stress on the contact surface is given by the Boussinesq's solution of the normal load on the surface of two half-spaces, which can be expressed by Green's function [15].
In order to solve the problem of minimum complementary potential energy, the entire contact area is discretized into a structural mesh, and the size of each rectangular element is × . In this research, the element size is same with the measurement result of 3D optical profilometer. Then, the formula for calculating the total complementary potential energy can be discretized as where = [ 1 , 2 ,…, ,…, ] is the vector of stress, is the influence matrix and = [( * ) 1 , ( * ) 2 ,…,( * ) ,…,( * ) ] is the vector of normal displacement. When equation (3) is satisfied, the complementary potential energy reaches the minimum value.
(3) Taking into account the contact conditions and the applied load , a linear complementarity problem can be obtained. Afterwards, the stress and contact area can be found to satisfy the following condition where A is the element area. By using the constrained conjugate gradient algorithm, we can solve the contact pressure in an iterative way. Finally, the convergence conditions are reached, that is > 0 and the accumulation of equal to load . Note that, the parameter of was determined according to the empirical equation proposed by Momber that is related to the unconfined compressive strength c and [13].

Fracture Closure Experiment
A cubic granite sample was prepared with a side length of 80 mm, and the crystal grain size was between 0.1 and 0.2 mm. The intact sample was split into two halves of the same size along the central axis. Note that, the damages caused by the splitting test to the rock fracture surfaces were small that can be neglected. Thus, the two newly generated fracture surfaces have almost the same geometric characteristics. A MTS 815 testing system was used to conduct uniaxial compression tests on three standard cylinder samples to determine the values of Young's modulus , Poisson's ratio and uniaxial compressive strength . The mean values of testing results are = 41.2 GPa, = 0.29, = 135.9 MPa.
A high-precision 3D optical profilometer (Keyence VR3000) with a resolution of 1 μm was used to measure the fracture surfaces, and the obtained digital surfaces was converted into point cloud data with a point interval of 0.1 mm. The lower surface is presented in figure 1.
Before the normal cyclic loading experiment, the upper block of granite sample was shifted to the left by 4 mm offset. As shown in figure 2, the granite sample was placed in a universal material testing machine with its loading direction perpendicular to the fracture surface. A normal incremental load was applied to the sample at a rate of 0.1 mm/min until the maximum stress of 30 MPa was reached. After that, the stress was continuously decreased until 0 MPa at an unloading rate of 0.1 mm/min. Take the above loading and unloading process as a unit, four consecutive cycles were conducted. As for unmated fractures, the compressive stress can be easily accumulated on several major contacts, and the maximum stress the sample can sustain was much less than the compressive strength of intact rock blocks. Here, 30 MPa was selected as the maximum stress in order to maintain the integrity of the sample during the cyclic loadings. Two LVDTs with a precision of 1 μm were installed between the upper and lower compression plates to record the normal displacement on the left and right sides of the sample. An identical compressive test was conducted on the original intact sample before splitting. Subtracting the displacement of intact sample from the total displacement of the unmated rock fracture, the pure displacement of the fracture was obtained [15].

Comparison between experimental results and simulation results
The relationship between normal displacement and normal stress for the tested granite sample is shown in figure 3. Obviously, in each loading process, the displacement increases nonlinearly with the increasing normal stress, while the increasing rate gradually decreases. Comparing the four experimental loadings, the measured displacement at n = 30 MPa are 0.993 mm, 1.029 mm, 1.059 mm and 1.074 mm, respectively. In the corresponding four unloading stages, the residual value of displacement at n = 0 MPa are 0.643 mm, 0.756 mm, 0.798 mm and 0.813 mm, respectively. This phenomenon reflects the hysteresis. Taking two adjacent loading and unloading as a cycle, the end point of each cycle can not restore the starting point to form a closed loop, leaving a certain amount of irreversible deformation. More specifically, with the increase number of cycles, the maximum net displacement of each loading decreases (i.e., 0.993mm, 0.386mm, 0.303mm and 0.276mm), and the irreversible deformation of each loading also decreases (i.e., 0.643mm, 0.113mm, 0.042mm and 0.015mm). It can be observed from figure 3 that the slope of a latter loading curve is steeper than the former, which indicates that the stiffness of rock fracture is enhanced as the increasing number of cycles. Numerical simulations can reasonably reflect the nonlinear behaviour of displacements during four loading processes. The simulation result of the first stress-displacement loading curve matches well with the experiment data, while the displacement value generated during the subsequent three simulations are slightly smaller than the experimental results. The comparison of experimentally measured and numerically simulated fracture stiffness is shown in figure 4. For the first loading process, the measured stiffness curve matches well with the numerical simulation result, and the stiffness-stress curve is approximately a straight line. Obviously, the n from the second loading to the fourth loading increases nonlinearly with the increasing n , showing a bending downward trend. The experimentally measured n oscillates, which is caused by the insufficient accuracy of the experimental measurement device at large normal stresses. For both experimental and numerical cycles, n is improved as the increasing number of cycles, which is mainly caused by the variation of the contact geometry (contact area ratio).
In figure 5, (a)~(d) show the plastic deformation area generated by the simulated cyclic loadings respectively. The amount of plastic deformation produced by the first loading is the largest, and it occupies a dominant position in the total plastic deformation generated by the four cycles. As the number of cyclic simulation increases, the amount of plastic deformation produced by each subsequent loading gradually decreases. In figure 5, (e) and (f) are the plastic deformation accumulated by the four loading cycles (i.e., the summation of (a)~(d)) and the final damage produced by the experiment, respectively. From a macroscopic point of view, the total numerically calculated damage areas are in good agreement with the experimentally measured areas. Comparing the details, it can be found that local area of numerical result is smaller than the experimental measured results. For this unmated case, the fracture walls are mainly sustained by several contact areas, resulting in a significant concentration of stresses and therefore substantial damages. The failure of asperities of brittle materials (such as granite) will result in non-localized damage, which makes the measured damage area greater than the simulated ones.

Discussion
In this study, we verified the high-precision contact model for modelling rock fracture closure behaviour through comparison with a cyclic loading experiment. The results show that the normal stress-displacement curves of cyclic loading experiment are in good agreement with the simulation results. As the number of cyclic loading increases, the normal stiffness of the rock fracture calculated in each cycle is greater than the previous one.
In order to study the influence of different degree of matedness on cyclic loadings, we added another two numerical calculations with offset values of 2 mm and 6 mm while maintaining the upper and lower surfaces unchanged. As shown in figure 6 (a), with the increase of offset value, the maximum displacement corresponding to each loading also increases, and the amount of residual deformation increases as well. figure 6 (b) shows the evolution of the stiffness-stress curve during the four cyclic loadings with two different offset values. It can be observed that as the offset value increases the normal stiffness of rock fracture gradually decreases. Due to the increase of the offset value, the correlation between the upper and lower fracture surfaces is reduced, meanwhile the number of contact points also decreases, thus weakening the ability of the fracture surface to resist normal deformation and surface damage. The calculated plastic deformation area under different offset values are shown in figure 7. Obviously, as the mating degree between rock fracture surfaces decreases, the number of contact points gradually decreases, while the local deformation increases. This means that the fracture becomes more compliable as the mating degree decreases.

Conclusion
The main concluding remarks drawn from this study are summarized as follows.
(1) A high-precision contact model was verified by a cyclic experiment on a rock fracture subject to normal stress. From the comparison of the normal stress-displacement curves and the plastic deformation areas, the contact model shows a good agreement with the results of cyclic loading experiment.
(2) Both of the experimental results and numerical calculations show that the normal stressdisplacement curves of different loading cycles exhibit obvious hysterical behavior. The hysteresis indicated that damages happened on partial contacting asperities, which tended to enhance the closure of a fracture. The hysteresis generated by a loading cycle is typically less than the previous one, showing a gradual weakening trend.
(3) The permanent deformation caused by damage during cyclic loadings can change the mechanical properties of fractures. As the number of loading cycles increases, the normal stiffness of the rock fracture tends to increase.