 Open access peer-reviewed chapter

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

Written By

Hao Kang

Submitted: December 5th, 2020 Reviewed: January 26th, 2021 Published: February 25th, 2021

DOI: 10.5772/intechopen.96229

From the Edited Volume

## Recent Developments in the Solution of Nonlinear Differential Equations

Edited by Bruno Carpentieri

Chapter metrics overview

View Full Metrics

## 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.

### Keywords

• visco-elastic deformation
• fast Fourier transform
• Boussinesq’s solution
• linear viscoelasticity
• rough fracture
• self-affine

## 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  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  proposed the conjugate gradient (CG) method and combined it with the FFT method to further improve the efficiency. Liu et al.  improved the drawbacks of the FFT method proposed by Stanley and Kato . Then, the CG and FFT methods have been applied to simulate plastic and visco-elastic deformations of rough fractures. Jacq et al.  and Sahlin et al.  considered perfect plasticity to simulate deformations of rough metal surfaces; and Wang et al.  considered strain-hardening plasticity.

For visco-elasticity, Chen et al.  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  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.  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  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.  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 . In this section, only the key mathematical concepts will be shown; the details can be found in their work . 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:

hxy=h0xyuexyδE1

where h0(x,y) is the initial aperture distribution, ue(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:

pxy=0andhxy>0E2
pxy>0andhxy=0E3

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  stated that the vertical displacement ue (x,y) can be calculated as:

uexy=++KxyxypxydxdyE4

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 ue (x,y) is the elastic displacement at location (x, y). The influence matrix K can be expressed as:

Kxyxy=1ν2πG1xx2+yy2E5

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:

xi=iΔx,i=1,2,,NE6
yj=jΔy,j=1,2,,ME7

where xi, yj 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:

hi,j=h0i,j+uei,jδE8
pi,j=0ifhi,j>0E9
pi,j>0ifhi,j=0E10

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

uei,j=l=1Mk=1NKi,k,j,l×pk,lE11
Ki,k,j,l=1ν2πGalnc+a2+c2d+a2+d2+blnd+b2+d2c+b2+c2+clna+a2+c2b+b2+c2+dlnb+b2+d2a+a2+d2E12

where

a=xixk+x2,b=xixkx2,c=yjyl+y2,d=yjyly2E13

As mentioned before, Stanley and Kato  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:

uei,j=IFFTFFTKi,k,j,l×FFTpk,lE14

where IFFT represents the inverse of Fourier transform. The FFT method reduces the number of operations from N2 * M2 to N*M*log(N*M) . 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:

Ftotal=k=1Nl=1Mpk,lE15

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

### 2.2 Method for calculating fracture visco-elastic deformation

As described before, Chen et al.  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. . In this section, only the key mathematical aspects will be introduced; the rest can be found in their work .

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:

σt=0tEtττdtE16
εt=0tJtττdtE17

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:

uexyt=0tJtττ++pxyτ1νπxx2+yy2dxdyE18

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

uexyt=0t++Jtτ1νπxx2+yy2pxyττdxdyE19

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 Nt 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:

uexyαt=α=1α++Jαtαt1νπxx2+yy2pxyαpxyα1dxdyE20

where α = 1, 2, …, Nt.

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:

uexyαt=α=1αpxyαpxyα1++Jαtαt1νπxx2+yy2dxdyE21

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:

++1ν2πGxx2+yy2dxdyDiscretizel=1Mk=1NKi,k,j,lE22

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

++Jαtαt1νπxx2+yy2dxdyDiscretizel=1Mk=1N2GJααtKi,k,j,lE23

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

ueijαt=α=1αl=1Mk=1N2GJααtKi,k,j,lpk,l,αpk,l,α1E24

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

ueijαt=α=1αueαE25

and

ueα=l=1Mk=1N2GJααtKi,k,j,lpk,l,αpk,l,α1E26

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

ueα=IFFTFFT2GJααtKi,k,j,l×FFTpk,l,αpk,l,α1E27

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.

### 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  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. ). 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. ). (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 σ0, the strain can be obtained:

εt=σ01G+tηE28

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

Jt=1G+tηE29

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

T=η/GE30

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:

ptr=2πR1υ0tGettG/ηddta2tr21/2dtE31

and

at=31νRP41G+tη3E32

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 G1 and G2, respectively. Under constant stress σ0, the strain can be obtained:

σt=G1G2G1+G2εt+G1ηG1+G2tdtE33

The creep compliance J(t) is expressed as:

Jt=1G1+1etG2/ηG2E34

The relaxation time T is defined as:

T=η/G2E35

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:

ptr=2πR1υ0tG1G1+G2G2+G1ettG1+G2/ηddta2tr21/2dtE36

and

at=31νRP41G1+1G21etG2/η3E37

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  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 . In Figures 4 and 5, rh is the contacting region at time zero, ph 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. ). Figure 5.Numerical and analytical solutions for the Maxwell model (Kang et al. ).

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 . Brown’s probabilistic model assumes that the surface is self-affine, and the surface height distribution follows Gaussian distribution . 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:

zxεHzεxE38

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 (A2) 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. ). Figure 7.Schematic of a power spectrum (Kang et al. ).

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

Cqq21+HE39

In Figure 7, the q has an upper bound and a lower bound. For the lower bound, qmin = 2π / λL, where λL is the surface dimension; for the upper bound, qmax = 2π/λ1, where λ1 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.  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  proposed a parameter: critical wavelength λc, which is also called the mismatch length scale. Brown  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.

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:

σ2=qminqmaxCqdqE40

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:

zdes=ziniσdesσiniE41

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

### 3.2 Generated synthetic surface pairs

Brown  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λcRMS
λc/LAbsolute value (μm)RMS/LAbsolute value (μm)
1100.60.110000.00550
2100.80.110000.00550
3101.00.110000.00550
4100.80.220000.00550
5100.80.330000.00550
6100.80.110000.010100
7100.80.110000.015150

### 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. ). (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. ). (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. ). (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 911, 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 911) is considered as the initial aperture field.

Surface pair No.Hλc (μm)RMS (μm)Average aperture (μm)Standard deviation of aperture (μm)
10.610005063.4114.29
20.810005037.308.57
31.010005021.895.14
40.820005055.9415.01
50.830005066.1020.12
60.8100010074.5917.15
70.81000150111.8925.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.  are used as the input parameters, and those properties are summarized in Table 3.

ParametersValue
Shear modulus, G (GPa)7.0
Poisson’s ratio, υ0.25
Viscosity, η (GPa*s)2.0 × 107
Relaxation time, τ = η / G (s)2.857 × 106

### 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. ). The time duration is normalized by τ. Figure 13.Contact ratio changing with time (Kang et al. ). 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.

ParametersAverage apertureContact ratio
Initial valueDecrease rateInitial valueIncrease 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. ). (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.  are used as the input parameters, and those properties are summarized in Table 5.

ParametersValue
Shear modulus, G1 (GPa)7.0
Shear modulus, G2 (GPa)7.0
Poisson’s ratio, υ0.25
Viscosity, η (GPa*s)2.0 × 107
Relaxation time, τ = η / G2 (s)2.857 × 106

### 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. ). The time duration is normalized by τ. Figure 16.Contact ratio changing with time (Kang et al. ). 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.

ParametersAverage apertureContact ratio
Initial valueDecrease rateInitial valueIncrease 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. ). (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 . 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 , 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.