Abstract
In many rock engineering projects such as hydrocarbon extraction and geothermal energy utilization, the hydraulic and mechanical behavior of rock fractures significantly affects the safety and profitability of the project. In field conditions, the hydraulic and mechanical behavior of rock fractures changes with time (the rock fractures creep), and this creep is not negligible even under dry conditions. In addition, creep is strongly affected by the rock fracture surface geometry. However, there is not much literature systematically studying the effect of surface geometry on rock fracture creep under dry conditions. This paper presents the results of a numerical study considering the effect of surface geometry on rough fracture viscoelastic deformations. An in-house numerical code has been developed to simulate the viscoelastic deformation of rough fractures. In addition, another numerical code has been developed to generate synthetic rough fracture surfaces by systematically changing the surface roughness parameters: the Hurst exponent, mismatch length, and root mean square roughness. The results indicate that by increasing the Hurst exponent or decreasing the mismatch length or decreasing the root mean square roughness, the fracture mean aperture decreases, and the contact ratio (number of contacting cells/total number of cells) increases faster with time.
1. Introduction
In many rock engineering projects such as hydrocarbon extraction, geothermal energy utilization, or nuclear waste storage, the hydraulic and mechanical behavior of rock fractures has a strong effect on the safety and profitability of the project. In field conditions, under constant tectonic stresses, rock fractures may creep. Creep, the deformation with time under constant stress, occurs in many types of rocks. The creep of intact rock has been studied extensively, and some authors proposed sophisticated creep models by conceptualizing the rock creep as a combination of springs and dashpots [1–4]. Also, creep rate (i.e., the change of deformation with time) may accelerate, stay constant, or decelerate. Such creep mechanisms often occur in rock fractures, essentially involving a time-dependent change of aperture under constant load. Since creep, in time, affects hydraulic conductivity and thus affects the safety and profitability of rock engineering projects, there is significant interest in understanding the creep behavior of rock fractures.
Many papers have discussed numerical modeling and laboratory experiments on creep mechanisms of rock fractures. The creep of rock fractures containing fluids is a process of four interdependent mechanisms: mechanical compression, pressure solution, dissolution, and erosion. It appears that most experiments and numerical models focus on the mechanisms of pressure solution (at contacting asperities) and dissolution (on free surfaces). For dissolution, fracture permeability evolution under the influence of acidic solutions has been widely investigated [5–11]. A common conclusion is that for calcite rock, when the pH of the fluid in the fracture is less than 6, fracture aperture distribution becomes more heterogeneous, as a result of acidic dissolution. The effect of pressure solution on rock fracture creep has also been studied extensively [12–16]. A common finding is that pressure solution strongly affects the fracture mechanical and hydraulic aperture, especially for siliciclastic rock. In recent years, some numerical models, which simulate fracture-scale pressure solution, have been proposed [17–23].
Based on the literature review on rock fractures in general, there is not much work systematically studying the creep mechanism caused by mechanical compression. Kang et al. [24] reported that for Musandam limestone fractures, the creep deformation under dry conditions is not negligible. Under dry conditions, mechanical compression is the only creep mechanism suggesting that it should be systematically studied. In addition, previous experiments suggest that the fracture surface geometry has a dominant effect on the hydraulic and mechanical behavior of rock fractures [24–27]. Therefore, the effect of surface geometry on fracture viscoelastic deformation should be systematically investigated.
Brown [28] proposed a simple mathematical model to geometrically describe rock fracture surfaces. In this model, a rock fracture surface can be completely described by three parameters: the Hurst exponent, the mismatch length, and the root mean square roughness. In this research, Brown’s [28] model will be used to generate synthetic fracture surfaces: the three above-mentioned parameters will be varied systematically. In tribology, Chen et al. (2011) first numerically simulated the viscoelastic deformation of rough fracture surfaces. In this research, their method will be used to numerically determine the viscoelastic deformation for each surface geometry parameter combination.
The objective of this numerical study is to investigate the effect of surface geometry on rock fracture viscoelastic deformations. First, the method for surface geometry generation and viscoelastic deformation calculation will be explained in Sections 2.1–2.3. Then, the numerical results will be validated against analytical results for a given surface geometry in Section 2.4. Next, the viscoelastic deformation of different rough fracture surfaces (generated based on Brown’s model) will be compared and discussed in Sections 3 and 4. Finally, recommendations will be made in Section 5. The results of this numerical study will eventually provide a reference and basis for fracture viscoelastic deformation prediction.
2. Methodology
2.1. Methodology for Rough Fracture Surface Generation
Synthetic fracture surface pairs are generated based on Brown’s probabilistic model [28]. In this model, rough rock fracture surfaces can be approximated by a Gaussian height distribution and self-affine organization over a large range of length scales ([28–30]). The rock fracture surface can be described by three parameters: Hurst exponent (), root mean square roughness at a reference length scale (), and mismatch length scale () [28].
The surface roughness of rock fractures typically follows a self-affine fractal distribution ([29, 31]). A self-affine surface can be described as where is the surface elevation, is a constant for direction scaling, and is the Hurst exponent. The value is between 0 and 1; and a larger value corresponds to a smoother local surface profile.
can be well quantified from the power spectrum of the surface. The power spectrum is calculated by decomposing the surface into a sum of sinusoidal waves, each with its amplitude, phase, and wavelength. Figure 1 shows the process of wave decomposition. In the wave decomposition process, the Fourier transform is used to decompose the random rough surface into a series of sinusoidal wave components. The squared amplitude of each wave component () is defined as the power, and a plot of power versus the inverse of wavelength () is defined as the power spectrum (see Figure 2 for a schematic). For self-affine fracture surfaces, the power spectrum can be related to the wavenumber ():
The upper frequency end, , is defined as , where is the smallest representable or measurable size. The lower frequency end, , is defined as , where is the surface length.
The second parameter is the mismatch length . As shown in Figure 1, each waveform can be completely described by three parameters: wavelength (), amplitude (), and phase. In a natural rock joint, the two surfaces may not be perfectly matched; they may have relative shear displacements [28–30]. Rock fracture surface measurement results from Brown [28] and Matsuki et al. (1995) suggest that natural fractures are matched at long wavelengths, but the surfaces are not identical at short wavelengths. Based on these observations, Brown [28] proposed the parameter of critical wavelength , called the mismatch length scale. Above the mismatch length scale, the decomposed sinusoidal wavelengths are perfectly matched; they have the same phase, wavelength, and amplitude. Below the mismatch length scale, the decomposed sinusoidal wavelengths behave independently; they have identical relationship between and , but independent phases. Figure 3 illustrates the concept of the mismatch length scale.
(a)
(b)
The third parameter is the root mean square roughness (). This parameter is used to define the absolute scale of the surface height. can be expressed as
After generating the surface profile, its root mean square value, , is calculated. Then, the surface height is scaled linearly as where is the initial surface height, is the final surface height, and is the designated root mean square roughness. Figure 4 illustrates the concepts of (surface height) and (aperture).
In this research, rough fracture surface pairs will be generated based on Brown’s [28] model. The Hurst exponent, root mean square roughness, and mismatch length scale will be varied systematically to investigate the effect of surface roughness parameters on fracture viscoelastic deformations.
2.2. Methodology for Calculating Fracture Elastic Deformation
The method for calculating fracture elastic deformation is introduced in this section. An in-house contact simulation code has been developed, and the method is similar to that formulated by Polonsky and Keer [32]. Here, only the key mathematical aspects will be presented; the details can be found elsewhere [32]. Polonsky and Keer’s method has been widely used in tribology and rock mechanics.
First, the aperture distribution function is defined as where is the rigid-body movement between two surfaces under applied stress in the direction normal to the fracture surface (compression is defined as positive), is the initial aperture, and is the elastic displacement field between the two contacting surfaces inside the contacting area.
For two rough surfaces contacting each other, the boundary conditions below need to be satisfied: where is the contact pressure (in the normal direction) acting on location . This indicates that the contact pressure is zero at noncontacting areas and is larger than zero at contacting areas. Compressive stress is defined as positive.
The vertical elastic displacement can be calculated from the Boussinesq and Cerrutti (1885) equation: where is the influence matrix for normal displacement at surface location due to a normal unit stress applied at location and is the normal stress applied at location . Figure 5 shows the schematic of the Boussinesq and Cerrutti equation.
The influence matrix is defined as where and are the shear modulus and Poisson’s ratio, respectively.
The elastic displacement of two rough surfaces with partial contact usually cannot be solved analytically, but can be solved numerically. For numerical simulation, the analytical equations need to be discretized into numerical equations. First, as shown in Figure 6, the surface is discretized into a regular grid: where and are and coordinates, respectively; and are the grid spacings in the and directions, respectively; and and are the total number of grids in the and directions, respectively.
The aperture field can be calculated as
The boundary conditions become
Love [33] discretized the Boussinesq and Cerrutti equation as where
Figure 6 illustrates the indices.
Andrews [34] first used the fast Fourier transform (FFT) method to solve the Boussinesq and Cerrutti equation in rock fractures. Later, Stanley and Kato [35] used this method in tribology. Fast Fourier transform is an algorithm that calculates the discrete Fourier transform of a sequence. In engineering, people often use Fourier analysis to convert a signal from its original domain (space or time) to a representation in the frequency domain. Discrete Fourier transform is obtained by decomposing the sequence of values into components of different frequencies, and using the fast Fourier transform can save the computation time. The detailed explanation of the FFT method can be found in Stanley and Kato [35, 36]. Using the FFT method, can be expressed as where FFT-1 means inverse fast Fourier transform (IFFT). In mathematics, FFT changes the convolution (see Equation (13)) to a simple multiplication (see Equation (16)).
If the FFT method is not used, calculating Equation (10) requires operations, while by using the FFT method, calculating Equation (13) requires operations. Therefore, when and are large, using the FFT method can save computation time.
In addition, force balance needs to be achieved: where is the applied external force on the two fracture surfaces and represents the pressure acting on grid (see Equation (12) and Figure 5). The indices and are explained in Figure 5.
Equations (10), (11), (12), (16), and (17) are solved iteratively using the conjugate gradient (CG) algorithm proposed by Polonsky and Keer [32], where the details can be found.
2.3. Methodology for Calculating the Viscoelastic Deformation
Chen et al. (2011) [37] proposed a numerical method to solve the viscoelastic deformation of rough fracture surfaces. In this research, an in-house viscoelastic contact simulation code has been developed, and the method is similar to that formulated by Chen et al. (2011). Here, the key mathematical concepts will be presented. Other details can be found in Chen et al.’s work (2011).
Before describing the methodology, it is essential to review the basic concepts of viscoelasticity. Viscoelastic responses can be modeled by a group of dashpots and springs, in which the springs represent elasticity and dashpots represent viscosity. In this research, we assume that the material is linearly viscoelastic. This means that the stress/strain response scales linearly with the stress/strain input and follows the law of linear superposition. For any arbitrary series of stress or strain inputs, the outputs can be expressed as where and are the relaxation modulus function and creep compliance function, respectively. (unit: pressure) describes the time-dependent stress response to a step-change in strain, and (unit: 1/pressure) describes the time-dependent strain response to a step-change in stress. represents the time variable in the integration, and its range is (this means that is a constant and is the variable). Figure 7 explains the creep compliance and relaxation modulus.
(a)
(b)
The analytical solution for the Boussinesq and Cerrutti equation can be modified by adding time-dependent viscoelasticity. Adding the creep compliance function into the Boussinesq and Cerrutti equation (Equations (7) and (8)) yields
In Equation (19), the shear modulus (in Equation (8)) is replaced by the creep compliance .
Rearranging Equation (19) yields
Equation (20) cannot be solved analytically. However, it can be solved numerically, if it can be decomposed into a linear equation system. Compared with Equation (7), an additional integration over the time duration () is added in Equation (23). Therefore, the time duration also needs to be discretized into many small time steps. Here, the total time duration is discretized into time steps with a uniform interval . It is assumed that is sufficiently small that the pressure field at each time step is unchanged. In addition, according to the fundamental theorem of calculus, the term in Equation (20) can be expressed as . Based on the sufficiently small time step assumption and the fundamental theorem of calculus, Equation (20) can be rewritten as a summation: where .
Within each time duration, in a given surface, the pressure field does not change with time. Therefore, can be factored out from the integral. Equation (21) becomes
The integral part of Equation (22) can be discretized. From Section 2.2 (Equations (7), (8), and (13)), the following discretization relationship can be obtained:
By multiplying both sides of Equation (23) by , the left-hand side of Equation (26) becomes the integral part of Equation (22). Therefore, the integral part of Equation (22) can be discretized as where is the shear modulus.
Thus, Equation (22) can be rewritten as
Equation (25) can be decomposed into Equations (26) and (27):
It is worth noting that Equation (27) is analogous to Equation (13). Therefore, to save computation time, Equation (27) can also be solved using the FFT method, similar to Equation (16): where , , and in Equation (28) are equivalent to , , and in Equation (18), respectively.
Equations (26) and (28) will be used for numerical simulation. At each specific time step, Equations (10), (11), (12), (17), (26), and (28) will be solved iteratively using the CG method [32], and the pressure and displacement fields at each step can be obtained. Then, will be increased (a new time step will be added), and the new pressure and displacement fields will be solved.
Based on the above mathematical concepts, the numerical algorithm is developed. The main algorithm is summarized below: (1)Obtain the loading history , initial discretized aperture field , and the creep compliance function for the viscoelastic material(2)Calculate for each time step and the stiffness matrix (3)Use the method described in Section 2.2; calculate the initial contact state and the pressure distribution (). In this case, there is no creep effect, and the stress and displacement fields are calculated based on elasticity(4)Save the historical pressure distribution field(5)Now increase (apply a new time increment). Since the historical pressures are known, the pressure distribution field and displacement field at time step can be calculated by the CG method(6)If reaches (the total number of time steps), stop the program, and export the pressure and displacement field. Otherwise, go back to step 4
Figure 8 shows the flowchart of the main algorithm [38].
2.4. Validation of the Numerical Program
The results of the numerical program need to be validated against analytical solutions. In this research, the cases for validation will be the Maxwell and Standard Linear Solid (SLS) viscoelastic half-spaces indented by a spherical rigid body. Lee and Radok [39] developed the analytical solutions for the two above-mentioned cases, and the pressure distribution along the contacting region will be compared.
It is essential to explain the concept of the Maxwell and Standard Linear Solid (SLS) constitutive models. Figures 9(a) and 9(b) show the schematic for the Maxwell model and SLS model, respectively.
(a)
(b)
The Maxwell model consists of a spring and a dashpot. The spring represents elasticity, with a shear modulus of ; and the dashpot represents viscosity, with a viscosity of (unit: time pressure). Under constant stress, the strain can be expressed as
Equation (29) indicates that under constant stress, the strain rate (creep rate) is constant.
The creep compliance function can therefore be defined as where is the time. The relaxation time for the Maxwell model, , is defined as
The SLS model consists of two springs and one dashpot. The two springs represent elasticity, with an elastic modulus of and , respectively; and the dashpot represents viscosity, with a viscosity of (unit: time pressure). Under constant stress , the strain can be related to as
The creep compliance function can be expressed as where is the time. The relaxation time for the SLS model, , is defined as
In the numerical code for the Maxwell model, Equation (30) is implemented in Equation (25), and the pressure and displacement fields are solved as described in Methodology; for the SLS model, Equation (33) is implemented in Equation (25), and the pressure and displacement fields are solved as described in Methodology.
Lee and Radok [39] stated that if a spherical rigid body is indenting into a viscoelastic flat surface (see Figure 10) and the flat surface satisfies the Maxwell model, the pressure distribution field can be expressed as where is the pressure, is the distance from the center of the contacting region, is the total time duration, is the radius of the contacting region, is the spherical diameter, is the shear modulus, is the viscosity, and is the total load. Figure 10 illustrates the physical parameters.
Lee and Radok [39] also suggested that if the flat surface satisfies the SLS model, the pressure distribution field can be expressed as where is the pressure, is the distance from the center of the contacting region, is the total time duration, is the radius of the contacting region, is the spherical diameter, and are the shear moduli, is the viscosity, and is the total load. The physical parameters are also illustrated in Figure 10.
Johnson (1985) [40] solved Lee and Radok’s analytical equations. Figures 11 and 12 compare the analytical and numerical solutions for the Maxwell and SLS models, respectively. The analytical solutions were calculated by Johnson (1985), and the numerical solutions were calculated by the code developed in this research. In the comparison, the absolute values of and do not matter, because the pressure distribution and the radius of the contacting region are normalized.
The results in Figures 11 and 12 suggest that the difference between the numerical and analytical solutions is less than 10%. Therefore, the numerical results match well with the analytical results.
3. Results
3.1. Surface Generation Results
As described in Section 2.1, synthetic rough fracture surfaces are generated based on Brown’s model [28]. Brown [28] analyzed the surface profiles of 23 natural rock joints and calculated the Hurst exponent (), root mean square roughness (RMS), and mismatch length () of each rock joint. The rock types include tuff, granodiorite, granite, basalt, chalk, and siltstone. Table 1 summarizes his analysis results. RMS and are normalized by the profile length.
Based on the results in Table 1, it can be concluded that the Hurst exponent is normally between 0.50 and 0.90; the normalized RMS () is normally between 0.005 and 0.015; and the normalized () is normally between 0.02 and 0.2. The synthetic surface pairs described in Methodology (Section 2.1) are generated based on the above values. Specifically, seven pairs of the synthetic surfaces are generated with three different , RMS, and values. Table 2 summarizes the information of the seven pairs of the synthetic surfaces. For all the surfaces, the dimensions (profile lengths) in the and directions are 10 mm.
Table 2 shows that in surface pairs 1, 2, and 3, the value is varied; in surface pairs 2, 4, and 5, the value is varied; in surface pairs 2, 6, and 7, the RMS value is varied.
For each pair of the surfaces, the aperture field can be calculated. Figures 13–15 show the effect of different , , and RMS values on the aperture field, respectively. The following conclusions can be drawn: (1)From Figure 13, as decreases, the mean and standard deviation of the aperture increase(2)From Figure 14, as increases, the mean and standard deviation of the aperture increase(3)From Figure 15, as RMS increases, the mean and standard deviation of the aperture increase. The magnitude of the aperture field scales linearly with the RMS value (this can be explained from Equation (4) in Section 2.1)
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
Table 3 summarizes the mean and standard deviation of the aperture field for the seven surface pairs. Each of the calculated aperture field is implemented into the creep simulation code as the initial aperture field. The creep deformation of the aperture field is calculated as described in Sections 2.2 and 2.3.
3.2. Creep Simulation Results for the Maxwell Model
The creep deformations of the seven surface pairs are calculated using the Maxwell model. Table 4 summarizes the input parameters. The parameters are obtained from the testing results of the Vaca Muerta Shale published by Mighani et al. [41].
The following parameters are also defined: (1)Macroscopic stress , where is the total force applied on the fracture surface and is the area of the fracture surface ()(2). See the definition in Figure 6
Figure 16 summarizes the effect of different , , and RMS values on the creep deformation. In the simulation, is fixed at 10 MPa. The initial values of the average aperture and contact area correspond to the elastic deformation. The time duration is from 0 to .
(a) Mean aperture
(b) Contact ratio
The following conclusions can be drawn: (1)Based on the curves of surface pairs 1, 2, and 3, when increases, the average aperture decreases, and the contact ratio increases faster with time(2)Based on the curves of surface pairs 2, 4, and 5, when decreases, the average aperture decreases, and the contact ratio increases faster with time(3)Based on the curves of surface pairs 2, 6, and 7, when RMS decreases, the average aperture decreases, and the contact ratio increases faster with time(4)Under current surface parameters, macroscopic stresses, and the time durations, the contact ratio is generally less than 10%
Table 5 summarizes the above conclusions.
Figure 17 compares the contact area and local contact stress before and after the creep stage. The white areas correspond to noncontacting areas, and the colored areas correspond to contacting areas. The input parameters are , , , and , and the creep . The left figure corresponds to the stress distribution before the creep (the moment after the elastic deformation), and the right figure corresponds to the stress distribution after the time duration of . The scale of the color bar is 0 to 2000 MPa. After creep, the contact area increases, and the local contact stresses (the vertical stress in the contacting cells, see the definition of cell in Figure 6) decrease (the color becomes darker).
(a)
(b)
3.3. Creep Simulation Results for the Standard Linear Solid (SLS) Model
The creep deformations of the seven surface pairs are also calculated by the SLS model. Table 6 summarizes the input parameters. The parameters are obtained from the testing results of the Vaca Muerta Shale published by Mighani et al. [41].
Figure 18 summarizes the effect of different , , and RMS values on the creep deformation. In the simulation, is fixed at 10 MPa. The initial values of the average aperture and contact area correspond to the elastic deformation. The time duration is from 0 to . Compared with the Maxwell model, the time duration is extended from to to better illustrate the decreasing creep rate.
(a) Mean aperture
(b) Contact ratio
The following conclusions can be drawn: (1)Based on the curves of surface pairs 1, 2, and 3, when increases, the average aperture decreases, and the contact ratio increases faster with time(2)Based on the curves of surface pairs 2, 4, and 5, when decreases, the average aperture decreases, and the contact ratio increases faster with time(3)Based on the curves of surface pairs 2, 6, and 7, when RMS decreases, the average aperture decreases, and the contact ratio increases faster with time(4)Under current surface parameters, macroscopic stresses, and time durations, the contact ratio is generally less than 10%(5)Under current surface parameters, macroscopic stresses, and time durations, the creep rate (aperture decrease rate) decreases significantly with time. The main reason is that the SLS model assumes exponentially decreasing creep rate
Table 7 summarizes the above conclusions. The effects of , RMS, and on the creep rate under the SLS model are similar to those under the Maxwell model.
Figure 19 compares the contact area and local contact stress before and after the creep stage. The white areas correspond to noncontacting areas, and the colored areas correspond to contacting areas. The input parameters are , , , and , and the creep time . The left figure corresponds to the stress distribution before the creep (the moment after the elastic deformation), and the right figure corresponds to the stress distribution after the time duration of . The scale of the color bar is 0 to 2000 MPa. After creep, the contact area increases, and the local contact stresses (the vertical stress in the contacting cells, see the definition of cell in Figure 6) decrease (the color becomes darker).
(a)
(b)
4. Brief Discussion
4.1. Creep Deformation Normalization: The Maxwell Model
The creep curves (mean aperture and contact ratio versus normalized time) under different , , and RMS values can be normalized. For each surface pair (see the definition in Sections 2.1 and 3.1; each surface pair generates one aperture field), between time duration 0 and , the mean aperture, , and the contact ratio, , are normalized as where is the normalized mean aperture, is the mean aperture right after the elastic deformation (before creep), is the mean aperture right after the elastic deformation for surface pair 2 (, , and ), is the normalized contact ratio, is the contact ratio right after the elastic deformation (before creep), and is the contact ratio right after the elastic deformation for surface pair 2.
Figure 20 summarizes the normalized mean aperture and contact ratio for different , , and RMS values. is fixed at 10 MPa, and the time duration is from 0 to .
(a) Normalized mean aperture
(b) Normalized contact ratio
The results show that by normalizing the values as described in Equations (37) and (38), the curves for different surface parameters fall into a narrow region.
4.2. Creep Deformation Normalization: The SLS Model
The creep curves (mean aperture and contact ratio versus normalized time) under different , , and RMS values can be normalized similarly as described in Section 4.1. For each surface pair, between time duration 0 and , the mean aperture, , and the contact ratio, , are normalized as where is the normalized mean aperture, is the mean aperture right after the elastic deformation (before creep), is the mean aperture right after the elastic deformation for surface pair 2 (, , and ), is the normalized contact ratio, is the contact ratio right after the elastic deformation (before creep), and is the contact ratio right after the elastic deformation for surface pair 2.
Figure 21 summarizes the normalized mean aperture and contact ratio for different , , and RMS values. is fixed at 10 MPa, and the time duration is from 0 to .
(a) Normalized mean aperture
(b) Normalized contact ratio
The results show that by normalizing the values as described in Equations (39) and (40), the curves for different surface parameters fall into a narrow region.
4.3. Contact Pressure Evolution during the Creep Stage: One Example
Figure 22 compares the histogram of local contact stress (only for the contacting cells, see the definition of cell in Figure 6) before and after the creep stage. The horizontal axis represents the local contact stress in each contacting cell, and the vertical axis represents the normalized frequency (the number of contacting cells within the stress range/the total number of contacting cells). The SLS model is used as the viscoelastic model, and the time duration is . Other input parameters are , , , and . The blue histogram corresponds to the local stress histogram before the creep (the moment right after the elastic deformation), and the orange histogram corresponds to the local stress histogram after the time duration of . This plot implies that after creep, the local stresses in the contact region decrease (the histogram shifts leftwards).
From the histograms, it is worth noting that the local contact stress can be larger than 200 MPa. However, for both the Maxwell and SLS models, the input elastic modulus and viscosity are obtained from triaxial tests with confining and differential pressures less than 50 MPa. Therefore, if possible, creep data from triaxial tests with confining and differential pressures larger than 200 MPa should be used to obtain more accurate simulation results.
4.4. Limitations for the Numerical Method
In this simulation, the contacting asperities deform viscoelastically. This implies that there is no upper limit for the compressive stress in the contacting cells. For some synthetic surface pairs, the compressive stress in a few contacting cells reaches 1.5 GPa. In reality, under such high compressive stresses, the cells may deform plastically, and the contact ratio will further increase. Therefore, if plastic deformation is ignored, the compressive stress in the contacting cells may be overestimated, and the contact ratio may be underestimated. In addition, asperity breakage or cracking is ignored. In reality, under high compressive stress, the asperities may break, which will affect the contact ratio and contact stress distribution.
5. Conclusions and Recommendations
In this research, an in-house numerical code has been developed to generate synthetic rock fracture surfaces based on Brown’s [28] model, and another code has been developed to simulate viscoelastic deformations of rough fractures. Seven synthetic fracture surface pairs are generated using different values of the Hurst exponent (0.6, 0.8, and 1.0), root mean square roughness RMS (50 μm, 100 μm, and 150 μm), and mismatch length (1000 μm, 2000 μm, and 3000 μm). The viscoelastic deformations of the seven surface pairs are simulated based on the Maxwell and Standard Linear Solid (SLS) models. For both models, the following conclusions can be drawn: (1)When increases, the average aperture decreases, and the contact ratio increases faster with time(2)When decreases, the average aperture decreases, and the contact ratio increases faster with time(3)When RMS decreases, the average aperture decreases, and the contact ratio increases faster with time(4)For the surface parameters ( between 0.6 and 1.0, RMS between 50 and 150 μm, and between 1000 and 3000 μm), macroscopic stresses (between 5 and 20 MPa), and the time durations ( for the Maxwell model and for the SLS model) used in the simulation, the contact ratio is generally less than 10%(5)The creep curves for different surface parameters can be normalized so the curves fall into a very narrow region
While the results are very informative, future work would be beneficial. Currently, the numerical code only simulates the deformation using the Maxwell and SLS models. Other viscoelastic models, such as the power law model and the Burgers model, could be implemented. In addition, the code only considers the viscoelastic deformation, while in reality, the contacting asperities may deform plastically. The code could be modified to simulate the visco-elasto-plastic deformation of fracture surfaces. Last but not least, the code calculates the creep deformation under dry conditions. Pressure solution equations can be added in addition to the elastic deformation code to simulate the fracture creep deformation under wet conditions.
Data Availability
The data used to support the results of this research are available from the corresponding author upon request.
Conflicts of Interest
None of the authors have any conflicts of interest.
Acknowledgments
This work was partially funded by the Abu Dhabi National Oil Company (ADNOC) (project ID: 015-RCM-2015). The authors gratefully acknowledge Yves Bernabe, Brian Evans, and John Williams for their useful suggestions.