Using the Boundary Element Method to Simulate Visco-Elastic Deformations of Rough Fractures

Abstract

In many engineering applications, such as tribology and rock mechanics, it is very important to understand the deformation of rough fractures to evaluate the safety and profitability of the project. Since a lot of materials can be characterized as visco-elastic materials, it is very significant to simulate the visco-elastic deformation of rough fractures. This chapter focuses on using the boundary element method to simulate visco-elastic deformations of rough fractures. First, the principles and procedures of the above-mentioned method will be introduced. Then, one example will be given in detail. This example investigates the effect of surface geometry on visco-elastic deformations of rough rock fractures under normal compressive stresses. The rock fracture surfaces are assumed to be self-affine, and synthetic rough surfaces are generated by systematically changing three surface roughness parameters: the Hurst exponent, root mean square roughness, and mismatch length. The results indicate that by decreasing the Hurst exponent or increasing the root mean square roughness or increasing the mismatch length, the fracture mean aperture increases, and the contact ratio (the number of contacting cells/total number of cells) increases slower with time. Finally, the limitations and possible future research directions will be briefly discussed.

1. Introduction

A lot of natural and engineering materials can be categorized as visco-elastic materials, such as rock, elastomers, and rubbers. In engineering applications, it is very important to understand and simulate the visco-elastic deformation of rough fractures. For example, in hydrocarbon extraction, we need to accurately simulate the visco-elastic deformation of rock fractures to predict production rates. In biomedical devices, we need to simulate the visco-elastic deformation of artificial joints to evaluate safety and effectiveness. Due to the geometrical complexity of rough fractures and the time-dependent properties of engineering materials, it is extremely difficult to obtain closed-form mathematical solutions. Thus, numerical models are required to simulate the time-dependent behavior of rough fractures.

The boundary element method (BEM) has been extensively used in solving rough surface contacting problems for distinct advantages compared with the traditional finite element method (FEM). First, it only requires discretization and calculation on the boundaries of the calculation domain, which is two-dimensional. On the contrary, FEM requires discretization and calculation for the whole calculation domain. As a result, to achieve the same stress calculation resolution, BEM requires much fewer numbers of elements and therefore, much less calculation time. In addition, since all the approximations are limited to the boundary, BEM has better stress calculation accuracy compared with FEM.

In recent years, researchers have been combining the BEM and fast numerical algorithms to achieve more efficient numerical simulations for contact problems. Stanley and Kato [1] published the first paper using the fast Fourier Transform (FFT) method to calculate the elastic deformation of rough surfaces under normal stresses. The FFT method makes the BEM simulation more efficient because FFT turns complicated convolution into simple matrix multiplication. Later, Polonsky and Keer [2] proposed the conjugate gradient (CG) method and combined it with the FFT method to further improve the efficiency. Liu et al. [3] improved the drawbacks of the FFT method proposed by Stanley and Kato [1]. Then, the CG and FFT methods have been applied to simulate plastic and visco-elastic deformations of rough fractures. Jacq et al. [4] and Sahlin et al. [5] considered perfect plasticity to simulate deformations of rough metal surfaces; and Wang et al. [6] considered strain-hardening plasticity.

For visco-elasticity, Chen et al. [7] first used the CG and FFT method to simulate visco-elastic deformations of rough fracture surfaces. They simulated three load-driven scenarios: rigid sphere indenting into PMMA surface, contact area evolution under constant load, and contact area evolution under harmonic cyclic load. Spinu and Cerlinca [8] applied different cut-off values for contact pressure to account for the plastic deformation of contacting asperities.

However, it appears that there is not much work that systematically simulates the visco-elastic deformation of rock fracture surfaces. Kang et al. [9] reported that for Musandam limestone fractures, the effect of mechanical compression on rock fracture time-dependent deformation is non-negligible, and should be systematically investigated. In addition, previous articles suggest that the fracture surface geometry has a significant effect on fracture time-dependent deformation. Therefore, we should systematically study the effect of surface geometry on rock fracture visco-elastic deformations.

Brown [10] proposed a simple probabilistic model to describe rock fracture surface geometry. In his model, the rock fracture surface geometry can be completely described by three key parameters: the Hurst exponent, the root mean square (RMS) roughness, and the mismatch length scale. In this research, his model will be used to generate synthetic fracture surface pairs, and the three key parameters will be changed systematically. The numerical method proposed by Chen et al. [7] will be used to simulate the visco-elastic deformation of synthetic fracture surfaces.

This chapter is organized as follows. Section 2 introduces and explains the principles and procedures of the numerical method. Section 3 provides a detailed example. The method for generating synthetic rough surfaces is introduced, and the effect of surface geometry parameters on the creep deformation is shown and discussed. Section 4 mentions the limitations of this method. Section 5 summarizes the findings.

2. BEM solution for visco-elastic deformations of rough fractures

2.1 Method for calculating fracture elastic deformation

Before explaining the method for visco-elastic deformation calculation, it is essential to introduce the method for elastic deformation calculation. The author has developed an in-house numerical code, which is similar to the algorithm proposed by Polonsky and Keer [2]. In this section, only the key mathematical concepts will be shown; the details can be found in their work [2]. It is worth noting that only the compressive stress (stress normal to the fracture surface) is considered; the shear stress (stress parallel to the fracture surface) is not considered.

First, the aperture (surface gap between two rough surfaces) distribution h (x,y) needs to be defined:

$$\mathrm{h}\left(\mathrm{x},\mathrm{y}\right)=h_0\left(x,y\right)-u_e\left(x,y\right)-\delta$$E1

where h₀(x,y) is the initial aperture distribution, u_e(x,y) is the elastic deformation of fracture surfaces, and δ is the rigid body displacement between two surfaces under applied normal stress. Here, compressive stress and fracture closure are defined as positive.

The boundary conditions are expressed as:

$$\mathrm{p}\left(\mathrm{x},\mathrm{y}\right)=0\mathit{\text{and}}h\left(x,y\right)>0$$E2

$$\mathrm{p}\left(\mathrm{x},\mathrm{y}\right)>0\mathit{\text{and}}h\left(x,y\right)=0$$E3

where p(x,y) is the contacting stress (normal to the surface) acting on location (x,y). Eqs. (2) and (3) indicate that the contacting stress is larger than zero at contacting regions, and is zero at non-contacting regions.

Boussinesq and Cerrutti [11] stated that the vertical displacement u_e (x,y) can be calculated as:

$$u_e\left(x,y\right)=\underset{-\infty }{\overset{+\infty }{\int }}\underset{-\infty }{\overset{+\infty }{\int }}K\left(x,y,x^{\prime },y^{\prime }\right)p\left({\mathrm{x}}^{\prime },y^{\prime }\right)d{x}^{\prime }d{y}^{\prime }$$E4

where p(x′, y′) is the normal pressure acting on location (x′, y′), K is the influence matrix, which represents the normal displacement at location (x, y) caused by unit normal pressure acting on location (x′, y′), and u_e (x,y) is the elastic displacement at location (x, y). The influence matrix K can be expressed as:

$$K\left(x,y,x^{\prime },y^{\prime }\right)=\frac{1-\nu }{2\mathit{\pi G}}\frac{1}{\sqrt{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2}}$$E5

where G is the shear modulus, and v is the Poisson’s ratio.

As mentioned in the introduction section, it is difficult to obtain the analytical solution for rough surface deformation under normal stress. However, the numerical solution can be obtained. To obtain the numerical solution, the fracture surface area needs to be discretized into rectangular grids:

$$x_i=i\mathrm{\Delta }x,i=1,2,\dots ,N$$E6

$$y_j=j\mathrm{\Delta }y,j=1,2,\dots ,M$$E7

where x_i, y_j are x and y coordinates, respectively; N and M are total number of grids in x- and y-direction, respectively; and Δx and Δy are the grid dimensions in x- and y-direction, respectively. After discretization, the aperture distribution function and boundary conditions can be expressed as:

$$h_{i,j}={\left(h_0\right)}_{i,j}+{\left(u_e\right)}_{i,j}-\delta$$E8

$$p_{i,j}=0\mathit{\text{if}}{h}_{i,j}>0$$E9

$$p_{i,j}>0\mathit{\text{if}}{h}_{i,j}=0$$E10

Love [12] first discretized Eqs. (4) and (5) as:

$${\left(u_e\right)}_{i,j}=\sum _{l=1}^M\sum _{k=1}^N{K}_{i,k,j,l}\times p_{k,l}$$E11

$$K_{i,k,j,l}=\frac{1-\nu }{2\mathit{\pi G}}\left(aln\frac{c+\sqrt{a^2+c^2}}{d+\sqrt{a^2+d^2}}+bln\frac{d+\sqrt{b^2+d^2}}{c+\sqrt{b^2+c^2}}+cln\frac{a+\sqrt{a^2+c^2}}{b+\sqrt{b^2+c^2}}+dln\frac{b+\sqrt{b^2+d^2}}{a+\sqrt{a^2+d^2}}\right)$$E12

where

$$a=x_i-x_k+\raisebox{1ex}{$∆x$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,b=x_i-x_k-\raisebox{1ex}{$∆x$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,c=y_j-y_l+\raisebox{1ex}{$∆y$}\!\left/ \!\raisebox{-1ex}{$2$}\right.,d=y_j-y_l-\raisebox{1ex}{$∆y$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$$E13

As mentioned before, Stanley and Kato [1] first the FFT method to solve Eq. (11) to make the calculation more efficient. The FFT method turns complicated convolution into simple matrix multiplication. By using the FFT method, Eq. (11) becomes:

$${\left(u_e\right)}_{i,j}=\mathit{\text{IFFT}}\left[\mathit{FFT}\left(K_{i,k,j,l}\right)\times \mathit{FFT}\left(p_{k,l}\right)\right]$$E14

where IFFT represents the inverse of Fourier transform. The FFT method reduces the number of operations from N² * M² to N*M*log(N*M) [1]. Therefore, when N and M are large, the FFT method can greatly reduce the calculation time.

The force balance over the entire fracture surface needs to be satisfied:

$$F_{\mathit{\text{total}}}=\sum _{k=1}^N\sum _{l=1}^M{p}_{k,l}$$E15

Eqs. (8)–(10), (14), and (15) are solved iteratively using the CG method proposed by Polonsky and Keer [2].

2.2 Method for calculating fracture visco-elastic deformation

As described before, Chen et al. [7] first combined the FFT and CG method to simulate visco-elastic deformations of rough fractures. The author has developed an in-house numerical code, which is similar to the algorithm described by Chen et al. [7]. In this section, only the key mathematical aspects will be introduced; the rest can be found in their work [7].

In this simulation, the rock materials are assumed to be linear viscoelastic. Therefore, is it essential to introduce the concept of linear viscoelasticity first. For linear viscoelastic materials, the stress/strain response scales linearly with the strain/stress input, and the behavior follows the rule of linear superposition. Mathematically, the stress/strain at time t can be expressed as:

$$\sigma \left(t\right)={\int }_0^{t}E\left(t-\tau \right)\frac{\mathit{d\epsilon }\left(\tau \right)}{\mathit{dt}}\mathit{d\tau }$$E16

$$\epsilon \left(t\right)={\int }_0^{t}J\left(t-\tau \right)\frac{\mathit{d\sigma }\left(\tau \right)}{\mathit{dt}}\mathit{d\tau }$$E17

where J(t) and E(t) are the creep compliance function and the relaxation modulus function, respectively. J(t) represents the time-dependent strain change with a step change in stress, and E(t) represents the time-dependent stress change with a step change in strain. Based on Eq. (17), the Boussinesq and Cerrutti equation can be modified to represent linear viscoelasticity by adding the creep compliance function:

$$u_e\left(x,y,t\right)={\int }_0^{t}J\left(t-\tau \right)\frac{\partial }{\partial \tau }\left[{\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }\frac{p\left(x^{\prime },y^{\prime },\tau \right)\left(1-\nu \right)}{\pi \sqrt{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2}}d{x}^{\prime }d{y}^{\prime }\right]\mathit{d\tau }$$E18

In Eq. (18), the creep compliance function J(t-τ) replaces the term 1/2G. Rearranging Eq. (18) gives:

$$u_e\left(x,y,t\right)={\int }_0^t{\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }\frac{J\left(t-\tau \right)\left(1-\nu \right)}{\pi \sqrt{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2}}\frac{\partial p\left(x^{\prime },y^{\prime },\tau \right)}{\partial \tau }d{x}^{\prime }d{y}^{\prime }\mathit{d\tau }$$E19

Eq. (19) cannot be solved analytically for rough fracture surfaces. However, if the time integration term can be de-coupled with the pressure integration term, Eq. (19) will become a linear equation system, and can therefore be solved numerically. To de-couple the time integration term, the time duration t is discretized into N_t time steps. The time interval is uniform, and is termed as Δt. The time interval is assumed to be sufficiently small that the pressure distribution field within each time interval does not change. In addition, based on the fundamental theorem of calculus, the term ∂p(x′, y′, τ) dτ/ ∂τ can be substituted by a finite difference p(x′, y′, τ + dτ) – p(x′, y′, τ). Therefore, Eq. (19) becomes:

$$u_e\left(x,y,\alpha ∆t\right)=\sum _{{\alpha }^{\prime }=1}^{\alpha }\left\{{\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }\frac{J\left(\alpha ∆t-{\alpha }^{\prime }∆t\right)\left(1-\nu \right)}{\pi \sqrt{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2}}\left[p\left(x^{\prime },y^{\prime },{\alpha }^{\prime }\right)-p\left(x^{\prime },y^{\prime },{\alpha }^{\prime }-1\right)\right]d{x}^{\prime }d{y}^{\prime }\right\}$$E20

where α = 1, 2, …, N_t.

In addition, within each time interval, the pressure distribution field does not change. Therefore, the pressure distribution field can be removed from the integration term:

$$u_e\left(x,y,\alpha ∆t\right)=\sum _{{\alpha }^{\prime }=1}^{\alpha }\left[p\left(x^{\prime },y^{\prime },{\alpha }^{\prime }\right)-p\left(x^{\prime },y^{\prime },{\alpha }^{\prime }-1\right)\right]\left\{{\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }\frac{J\left(\alpha ∆t-{\alpha }^{\prime }∆\mathrm{t}\right)\left(1-\nu \right)}{\pi \sqrt{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2}}d{x}^{\prime }d{y}^{\prime }\right\}$$E21

Eq. (21) indicates that the time integration term is de-coupled with the pressure integration term. The pressure integration term can then be discretized, similar to Eq. (11). From Eqs. (4), (5), and (11), the Boussinesq equation can be discretized as:

$${\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }\frac{\left(1-\nu \right)}{2\mathit{\pi G}\sqrt{{\left(x-x^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2}}d{x}^{\prime }d{y}^{\prime }\stackrel{\mathit{\text{Discretize}}}{\to }\sum _{l=1}^M\sum _{k=1}^N{K}_{i,k,j,l}$$E22

Based on Eq. (22), the integration part of Eq. (21) can then be discretized:

$${\int }_{-\infty }^{+\infty }{\int }_{-\infty }^{+\infty }\frac{J\left(\alpha ∆t-{\alpha }^{\prime }∆t\right)\left(1-\nu \right)}{\pi \sqrt{{\left(x-{\mathrm{x}}^{\prime }\right)}^2+{\left(y-y^{\prime }\right)}^2}}d{x}^{\prime }d{y}^{\prime }\stackrel{\mathit{\text{Discretize}}}{\to }\sum _{l=1}^M\sum _{k=1}^{N}2\mathit{GJ}\left(\left(\alpha -{\alpha }^{\prime }\right)∆t\right)K_{i,k,j,l}$$E23

Therefore, Eq. (21) can be discretized as:

$$u_e\left(i,j,\alpha ∆t\right)=\sum _{{\alpha }^{\prime }=1}^{\alpha }\sum _{l=1}^M\sum _{k=1}^{N}2\mathit{GJ}\left(\left(\alpha -{\alpha }^{\prime }\right)∆t\right)K_{i,k,j,l}\left(p_{k,l,{\alpha }^{\prime }}-p_{k,l,{\alpha }^{\prime }-1}\right)$$E24

To implement FFT, Eq. (24) can be decoupled into two equations:

$$u_e\left(i,j,\alpha ∆t\right)=\sum _{{\alpha }^{\prime }=1}^{\alpha }{\left(u_e\right)}_{{\alpha }^{\prime }}$$E25

and

$${\left(u_e\right)}_{{\alpha }^{\prime }}=\sum _{l=1}^M\sum _{k=1}^{N}2\mathit{GJ}\left(\left(\alpha -{\alpha }^{\prime }\right)∆t\right)K_{i,k,j,l}\left(p_{k,l,{\alpha }^{\prime }}-p_{k,l,{\alpha }^{\prime }-1}\right)$$E26

Eq. (26) can be solved by the FFT method, similar to Eqs. (13) and (14):

$${\left(u_e\right)}_{{\alpha }^{\prime }}=\mathit{\text{IFFT}}\left[\mathit{FFT}\left(2\mathit{GJ}\left(\left(\alpha -{\alpha }^{\prime }\right)∆t\right)K_{i,k,j,l}\right)\times \mathit{FFT}\left(p_{k,l,{\alpha }^{\prime }}-p_{k,l,{\alpha }^{\prime }-1}\right)\right]$$E27

Within each time step, Eqs. (8)–(10), (15), (25), and (27) are solved using the CG method. The pressure distribution field is obtained and stored. Then, a new time step will be added (α will be increased by one), and the new deformation and pressure fields will be solved based on the historical pressure fields. Figure 1 summarizes the main calculation algorithm based on the above mathematical concepts.

Figure 1.

Summary of the calculation algorithm (Kang et al. [13]).

2.3 Model validation

Before simulating visco-elastic deformations of rough rock fractures, it is essential to validate the numerical code against analytical solutions. In this research, the analytical solutions provided by Radok and Lee [14] will be used for validation. In their solutions, a rigid spherical indenter is indented into a flat visco-elastic surface; and the visco-elastic models for the flat surface are the Maxwell and Standard Linear Solid (SLS) model. Figure 2 illustrates the geometry setup for the analytical solution, and Figure 3 shows the concepts of the Maxwell and SLS model.

Figure 2.

Geometry setup for the analytical solution (Kang et al. [13]). R is the radius of the spherical rigid indenter, P is the total load, δ is the indentation depth, t is the time duration, and a(t) is the radius of the contacting region.

Figure 3.

Concepts of the Maxwell and SLS model (Kang et al. [13]). (a): Schematic of the Maxwell model; (b): Schematic of the SLS model.

The Maxwell model consists of a dashpot and a spring. The dashpot represents viscosity, with a viscosity of η; the spring represents elasticity, with a shear modulus of G. Under constant stress σ₀, the strain can be obtained:

$$\epsilon \left(t\right)={\sigma }_0\left(\frac{1}{G}+\frac{t}{\eta }\right)$$E28

Eq. (28) implies that under constant stress, the strain rate does not change with time. The creep compliance can be expressed as:

$$J\left(t\right)=\frac{1}{G}+\frac{t}{\eta }$$E29

Another parameter, the relaxation time T, is defined as:

$$T=\eta /G$$E30

In the numerical simulation, Eq. (29) will be implemented into Eq. (27), and the displacement and pressure field will be solved as described in Sections 2.1 and 2.2. For the geometry setup shown in Figure 2, the analytical solution for the contacting region radius and pressure field can be obtained:

$$p\left(t,r\right)=\frac{2}{\mathit{\pi R}\left(1-\upsilon \right)}{\int }_0^{t}G{e}^{-\left(t-t^{\prime }\right)G/\eta }\frac{d}{d{t}^{\prime }}{\left[a^2\left(t^{\prime }\right)-r^2\right]}^{1/2}d{t}^{\prime }$$E31

and

$$a\left(t\right)=\sqrt[3]{\frac{3\left(1-\nu \right)\mathit{RP}}{4}\left(\frac{1}{G}+\frac{t}{\eta }\right)}$$E32

where p is the pressure field, t is the total time duration, υ is the Poisson’s ratio, and r is the distance from the center of the contacting region.

The SLS model consists of one dashpot and two springs. The dashpot represents viscosity, with a viscosity η; the two springs represent elasticity, with a shear modulus of G₁ and G₂, respectively. Under constant stress σ₀, the strain can be obtained:

$$\mathrm{\sigma }\left(\mathrm{t}\right)=\frac{G_1{G}_2}{G_1+G_2}\epsilon \left(t\right)+\frac{G_1\eta }{G_1+G_2}\frac{\mathit{d\epsilon }\left(t\right)}{\mathit{dt}}$$E33

The creep compliance J(t) is expressed as:

$$\mathrm{J}\left(\mathrm{t}\right)=\frac{1}{G_1}+\frac{1-e^{-\mathrm{t}{G}_2/\eta }}{G_2}$$E34

The relaxation time T is defined as:

$$T=\eta /G_2$$E35

In the numerical simulation, Eq. (34) will be implemented into Eq. (27), and the displacement and pressure field will be solved as described in Sections 2.1 and 2.2. For the geometry setup shown in Figure 2, the analytical solution for the contacting region radius and pressure field can be obtained:

$$p\left(t,r\right)=\frac{2}{\mathit{\pi R}\left(1-\upsilon \right)}{\int }_0^t\frac{G_1}{G_1+G_2}\left[G_2+G_1{e}^{-\left(t-t^{\prime }\right)\left(G_1+G_2\right)/\eta }\right]\frac{d}{d{t}^{\prime }}{\left[a^2\left(t^{\prime }\right)-r^2\right]}^{1/2}d{t}^{\prime }$$E36

and

$$a\left(t\right)=\sqrt[3]{\frac{3\left(1-\nu \right)\mathit{RP}}{4}\left(\frac{1}{G_1}+\frac{1}{G_2}\left(1-e^{-t{G}_2/\eta }\right)\right)}$$E37

where p is the pressure field, t is the total time duration, υ is the Poisson’s ratio, and r is the distance from the center of the contacting region.

Johnson [15] solved Eqs. (31), (32), (36), and (37), Figures 4 and 5 compare the numerical and analytical solutions for the SLS and Maxwell models, respectively. The solid lines are the numerical solutions obtained by the author, and the dashed lines are the analytical solutions solved by Johnson [15]. In Figures 4 and 5, r_h is the contacting region at time zero, p_h is the maximum contacting pressure at time zero, and T is the relaxation time.

Figure 4.

Numerical and analytical solutions for the SLS model (Kang et al. [13]).

Figure 5.

Numerical and analytical solutions for the Maxwell model (Kang et al. [13]).

Figures 4 and 5 indicate the deviation between the numerical and analytical results is less than 10%. Therefore, the numerical code can be used to simulate the visco-elastic deformations of rough fractures. For the two validation cases, the numerical simulation accuracy is not strongly dependent on the total number of elements, but on the time interval Δt. The deviation between numerical and analytical solutions will be smaller if the time interval Δt is reduced.

3. One example: visco-elastic deformations of rough rock fractures

3.1 Brief introduction of Brown’s (1995) model

In this chapter, synthetic fracture surface pairs are generated based on Brown’s model [10]. Brown’s probabilistic model assumes that the surface is self-affine, and the surface height distribution follows Gaussian distribution [10]. The surface geometry can be completely described by three parameters: the Hurst exponent H, the mismatch length λ_c, and the root mean square roughness RMS.

Mathematically, a self-affine surface is defined as:

$$z\left(x\right)\sim {\epsilon }^{-H}z\left(\mathit{\epsilon x}\right)$$E38

where H is the Hurst exponent, z is the surface height, and ε is a constant for scaling at the x-direction. The H value is between 0 and 1, and it describes local roughness. A smaller H value corresponds to a rougher local surface profile.

The H value can be obtained from the power spectrum of surface height. The power spectrum of a surface can be obtained by decomposing the surface profile into a series of sinusoidal waves via Fourier transform, and each sinusoidal wave has its own amplitude A, wavelength λ, and phase. Figure 6 shows the schematic of the decomposition process. The power (A²) is defined as the square of the amplitude A; and the plot of power versus the wavelength number (the inverse of wavelength, which is 2π/λ) is defined as the power spectrum. Figure 7 shows the schematic of power spectrum.

Figure 6.

Schematic of wave decomposition via Fourier transform (Kang et al. [13]).

Figure 7.

Schematic of a power spectrum (Kang et al. [13]).

For a self-affine fracture surface, the power C (=A²) can be related to the wavelength number q (=2π/λ) as:

$$\mathrm{C}\left(\mathrm{q}\right)\sim q^{-2\left(1+H\right)}$$E39

In Figure 7, the q has an upper bound and a lower bound. For the lower bound, q_min = 2π / λ_L, where λ_L is the surface dimension; for the upper bound, q_max = 2π/λ₁, where λ₁ is the surface measurement resolution.

The second parameter is the mismatch length, λ_c. As illustrated in Figure 6, each wave component has its own wavelength λ. Glover et al. [16] and Brown [10, 17, 18] stated that for most natural rock joints, the two surfaces have relative shear displacements. At long wavelengths, the wave components match well; at short wavelengths, the wave components are not identical. Based on the above observation, Brown [10] proposed a parameter: critical wavelength λ_c, which is also called the mismatch length scale. Brown [10] assumed that above the mismatch wavelength, the decomposed wave components of two surfaces match perfectly; they have the same amplitudes, wavelengths, and phases. On the contrary, below the mismatch wavelength, the decomposed wave components of two surfaces do not match; they have the same amplitudes and wavelengths, but the phases are independent. Figure 8 illustrates the concept of the mismatch wavelength.

Figure 8.

Illustration of the mismatch wavelength (Kang et al. [13]).

The third parameter is the root mean square roughness, RMS. It represents the absolute scale of surface asperity elevation. Mathematically, the RMS is defined as:

$${\sigma }^2={\int }_{q_{\mathit{min}}}^{q_{\mathit{max}}}C\left(q\right)\mathit{dq}$$E40

where C is the power, q is the wavelength number, and σ is the RMS value. When generating the synthetic surface, the surface heights are normalized by its own RMS value, σ_ini, and then multiplied by the designated RMS value, σ_des:

$$z_{\mathit{des}}=z_{\mathit{ini}}\frac{{\sigma }_{\mathit{des}}}{{\sigma }_{\mathit{ini}}}$$E41

where z_ini is the initial surface height and z_des is the surface height after linear scaling. In this chapter, only the key mathematical concepts of Brown’s [10] model is introduced; other details can be found in [10].

3.2 Generated synthetic surface pairs

Brown [10] measured the Hurst exponent H, mismatch length λ_c, and RMS for 23 natural rock joints. His measurement results imply that the H value is normally between 0.5 and 1.0; the normalized λ_c value (λ_c/fracture profile length) is normally between 0.02 and 0.2, and the normalized RMS value (RMS/fracture profile length) is normally between 0.005 and 0.015. Based on the above conclusion, seven synthetic fracture surface pairs are generated, with different H, λ_c, and RMS values. Table 1 summarizes the parameters of the seven synthetic surface pairs. It is worth noting that surface pair No. 2 is the reference surface pair.

Surface Pair No.	Profile length L (mm)	H	λc		RMS
			λc/L	Absolute value (μm)	RMS/L	Absolute value (μm)
1	10	0.6	0.1	1000	0.005	50
2	10	0.8	0.1	1000	0.005	50
3	10	1.0	0.1	1000	0.005	50
4	10	0.8	0.2	2000	0.005	50
5	10	0.8	0.3	3000	0.005	50
6	10	0.8	0.1	1000	0.010	100
7	10	0.8	0.1	1000	0.015	150

Table 1.

The parameters of the seven synthetic surface pairs.

Table 1 shows that between surface pairs 1, 2, and 3, the H value is varied; between surface pairs 2, 4, and 5, the λ_c value is varied; between surface pairs 2, 6, and 7, the RMS value is varied. For each surface pair, the aperture distribution field can be plotted. Figure 9 plots the aperture fields for surface pairs 1, 2, and 3; Figure 10 plots the aperture fields for surface pairs 2, 4, and 5, and Figure 11 plots the aperture fields for surface pairs 2, 6, and 7.

Figure 9.

Aperture fields for different H values (Kang et al. [13]). (a): aperture field for surface pair 1; (b): aperture field for surface pair 2; (c): aperture field for surface pair 3. The color bar scales are identical.

Figure 10.

Aperture fields for different λ_c values (Kang et al. [13]). (a): aperture field for surface pair 2; (b): aperture field for surface pair 4; (c): aperture field for surface pair 5. The color bar scales are identical.

Figure 11.

Aperture fields for different RMS values (Kang et al. [13]). (a): aperture field for surface pair 2; (b): aperture field for surface pair 6; (c): aperture field for surface pair 7. The color bar scales scale linearly with the RMS value.

Based on Figures 9–11, we have the following observations:

1. According to Figure 9, when H increases, the average and standard deviation of the aperture decreases;

2. According to Figure 10, when λ_c deceases, the average and standard deviation of aperture decreases;

3. According to Figure 11, the average and standard deviation of aperture scales linearly with the RMS value.

Table 2 summarizes the mean and standard deviation of aperture for each surface pair. In the numerical code, each calculated aperture field (shown in Figures 9–11) is considered as the initial aperture field.

Surface pair No.	H	λc (μm)	RMS (μm)	Average aperture (μm)	Standard deviation of aperture (μm)
1	0.6	1000	50	63.41	14.29
2	0.8	1000	50	37.30	8.57
3	1.0	1000	50	21.89	5.14
4	0.8	2000	50	55.94	15.01
5	0.8	3000	50	66.10	20.12
6	0.8	1000	100	74.59	17.15
7	0.8	1000	150	111.89	25.72

Table 2.

The average and standard deviation of seven synthetic surface pairs.

3.3 Creep simulation results for the Maxwell model

The author uses the Maxwell model to calculate the visco-elastic deformation of seven synthetic surface pairs. The mechanical properties of Vaca Muerta Shale measured by Mighani et al. [19] are used as the input parameters, and those properties are summarized in Table 3.

Parameters	Value
Shear modulus, G (GPa)	7.0
Poisson’s ratio, υ	0.25
Viscosity, η (GPa*s)	2.0 $\times$ 10^7
Relaxation time, τ = η / G (s)	2.857 $\times$ 10^6

Table 3.

Input parameters for the Maxwell model.

Before showing the results, two parameters are introduced: macroscopic stress σ and contact ratio:

1. The macroscopic stress σ = total force applied to the fracture/fracture surface area;

2. Contact ratio = 100 * (the number of grids in contact/total number of grids).

Figures 12 and 13 show the mean aperture and contact ratio evolving with time for seven synthetic surface pairs, respectively. The total time duration is 2τ, and the macroscopic stress σ = 10 MPa. The initial changes of the mean aperture and contact ratio correspond to fracture elastic deformation.

Figure 12.

Mean aperture changing with time (Kang et al. [13]). The time duration is normalized by τ.

Figure 13.

Contact ratio changing with time (Kang et al. [13]). The time duration is normalized by τ.

Based on Figures 12 and 13, several conclusions can be drawn:

1. As H decreases, the mean aperture increases, and the contact ratio increases slower with time;

2. As RMS increases, the mean aperture increases, and the contact ratio increases slower with time;

3. As λ_c increases, the mean aperture increases, and the contact ratio increases slower with time.

4. Under current macroscopic stress, time duration, and surface parameters, the contact ratio is generally less than 9.5%.

Table 4 summarizes the effect of surface parameters on the mean aperture and contact ratio.

Parameters	Average aperture		Contact ratio
	Initial value	Decrease rate	Initial value	Increase rate
H↓	↑	↑	↓	↓
λc↑	↑	↑	↓	↓
RMS↑	↑	↑	↓	↓

Table 4.

Effect of surface parameters on the mean aperture and contact ratio.

Figure 14 shows the contact region and local contacting stress evolution of surface pair 3 before and after the creep stage. The macroscopic stress is 10 MPa and the creep time duration is 2τ. The colored regions and white regions correspond to the contacting regions and non-contacting regions, respectively. The color bar scale is 2000 MPa. After the creep stage, the area of contacting regions becomes larger, and the local contacting stress reduces. However, even after the creep stage, the contact ratio is still less than 9.5%. Under the same time duration, if η is reduced, the contact area will increase more rapidly.

Figure 14.

Contact region and local contacting stress evolution before and after the creep stage (Kang et al. [13]). (a): before the creep stage; (b): after the creep stage. In both x- and y-directions, the number of grids is 512. The contact area increase is qualitatively shown.

3.4 Creep simulation results for the SLS model

The author also uses the SLS model to calculate the visco-elastic deformation of seven synthetic surface pairs. The mechanical properties of Vaca Muerta Shale measured by Mighani et al. [19] are used as the input parameters, and those properties are summarized in Table 5.

Parameters	Value
Shear modulus, G1 (GPa)	7.0
Shear modulus, G2 (GPa)	7.0
Poisson’s ratio, υ	0.25
Viscosity, η (GPa*s)	2.0 $\times$ 10^7
Relaxation time, τ = η / G2 (s)	2.857 $\times$ 10^6

Table 5.

Input parameters for the SLS model.

Figures 15 and 16 show the mean aperture and contact ratio evolving with time for seven synthetic surface pairs, respectively. The total time duration is 5τ, and the macroscopic stress σ = 10 MPa. The total time duration is extended from 2τ to 5τ to show the time-decaying creep rate. The initial changes of the mean aperture and contact ratio correspond to fracture elastic deformation.

Figure 15.

Mean aperture changing with time (Kang et al. [13]). The time duration is normalized by τ.

Figure 16.

Contact ratio changing with time (Kang et al. [13]). The time duration is normalized by τ.

Based on Figures 15 and 16, several conclusions can be drawn:

1. As H decreases, the mean aperture increases, and the contact ratio increases slower with time;

2. As RMS increases, the mean aperture increases, and the contact ratio increases slower with time;

3. As λ_c decreases, the mean aperture increases, and the contact ratio increases slower with time.

4. Under current macroscopic stress, time duration, and surface parameters, the contact ratio is generally less than 7.0%.

5. Under current macroscopic stress, time duration, and surface parameters, the creep rate decreases significantly with time. This is mainly because the SLS model assumes an exponentially decaying creep rate.

Table 6 summarizes the effect of surface parameters on the mean aperture and contact ratio.

Parameters	Average aperture		Contact ratio
	Initial value	Decrease rate	Initial value	Increase rate
H↓	↑	↑	↓	↓
λc↑	↑	↑	↓	↓
RMS↑	↑	↑	↓	↓

Table 6.

Effect of surface parameters on the mean aperture and contact ratio.

Figure 17 shows the contact region and local contacting stress evolution of surface pair 3 before and after the creep stage. The macroscopic stress is 10 MPa and the creep time duration is 5τ. The colored regions and white regions correspond to the contacting regions and non-contacting regions, respectively. The color bar scale is 2000 MPa. After the creep stage, the area of contacting regions becomes larger, and the local contacting stress reduces. However, even after the creep stage, the contact ratio is still less than 7.0%. Under the same time duration, if η is reduced, the contact area increase will increase more rapidly.

Figure 17.

Contact region and local contacting stress evolution before and after the creep stage (Kang et al. [13]). (a): before the creep stage; (b): after the creep stage. In both x- and y-directions, the number of grids is 512. The contact area increase is qualitatively shown.

4. Limitations of the method

In this numerical method, the contacting asperities deform visco-elastically, and there is no upper limit on the local contacting stress. For some synthetic surfaces, the contacting stress in a few cells exceed 1.3 GPa. In reality, under such high contacting stresses, the asperities may deform plastically. Ignoring the plastic deformation will underestimate the contact ratio and overestimate the local contacting stress. In addition, asperity breakage is ignored in this numerical method. Under high contacting stresses, asperities may break, which will further change the contacting regions and the contacting stress distribution [20]. Furthermore, the effect of shear stress on fracture visco-elastic deformations is also not considered. In engineering applications (especially in oil and gas production), fractures may subject to shear stress, which may significantly change the visco-elastic deformations.

5. Conclusions

This chapter explains how to use the boundary element method to calculate visco-elastic deformations of rough fractures. Fast numerical algorithms (CG and FFT) are implemented to further improve the efficiency. In addition, one example, which investigates the effect of surface geometry on visco-elastic deformations of rough rock fractures, is given. In this example, the author builds two in-house numerical codes: one code generates synthetic fracture surface pairs using Brown’s probabilistic model [10], and the other simulates the visco-elastic deformations of the synthetic surface pairs. Seven synthetic surface pairs are generated by systematically changing the values of the root mean square roughness RMS (50 μm, 100 μm, and 150 μm), mismatch length λ_c (1000 μm, 2000 μm, and 3000 μm), and Hurst exponent H (0.6, 0.8, and 1.0). Then, the author simulates the visco-elastic deformation of the seven surface pairs by using the Standard Linear Solid (SLS) and the Maxwell model. The following key conclusions can be drawn:

1. As RMS increases, the average aperture increases, and the contact ratio increases slower with time;

2. As λ_c increases, the average aperture increases, and the contact ratio increases slower with time;

3. As H decreases, the average aperture increases, and the contact ratio increases slower with time;

4. For the macroscopic stress (10 MPa), time durations (5τ for the SLS model and 2τ for the Maxwell model), and the surface roughness parameters (RMS between 50 and 150 μm, λ_c between 1000 and 3000 μm, H between 0.6 and 1.0) used in the examples, the contact ratio is less than 9.5%.

While the results are useful, future work would be helpful. First, more surface roughness parameter values can be used so a quantitative relationship between surface parameters and contact ratio or average aperture can be obtained. In addition, other visco-elastic models, such as the Burgers model and the Power Law model, can be implemented. Furthermore, in this simulation, the plastic deformation of contacting asperities is not considered. As a result, the local contacting stress may be overestimated. The plastic deformation of contacting asperities can be considered so the results can be more realistic. Last but not least, the effect of shear stress can be simulated to make the results more applicable.

Acknowledgments

The authors gratefully acknowledge Farrokh Sheibani, Stephen Brown, Herbert Einstein, and John Germaine for their useful suggestions.

Conflict of interest

The authors declare no conflict of interest.