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

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.


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

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: where h 0 (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: px ,y ÀÁ > 0 and h x, y ðÞ ¼ 0 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: where p(x 0 ,y 0 ) is the normal pressure acting on location (x 0 ,y 0 ), K is the influence matrix, which represents the normal displacement at location (x, y) caused by unit normal pressure acting on location (x 0 ,y 0 ), and u e (x,y) is the elastic displacement at location (x, y). The influence matrix K can be expressed as: 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: 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: Love [12] first discretized Eqs. (4) and (5) as: where 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: where IFFT represents the inverse of Fourier transform. The FFT method reduces the number of operations from N 2 *M 2 to N*M*log(N*M) [1]. Therefore, when N and M are large, the FFT method can greatly reduce the calculation time.

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: 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: In Eq. (18), the creep compliance function J(t-τ) replaces the term 1/2G. Rearranging Eq. (18) gives: 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 0 ,y 0 , τ)dτ/ ∂τ can be substituted by a finite difference p(x 0 ,y 0 , τ +dτ) -p(x 0 ,y 0 , τ). Therefore, Eq. (19) becomes: 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: 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: Therefore, Eq. (21) can be discretized as: To implement FFT, Eq. (24) can be decoupled into two equations: and Eq. (26) can be solved by the FFT method, similar to Eqs. (13) and (14): 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.

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.
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: Eq. (28) implies that under constant stress, the strain rate does not change with time. The creep compliance can be expressed as: 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.   Another parameter, the relaxation time T, is defined as: 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: and 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 1 and G 2 , respectively. Under constant stress σ 0 , the strain can be obtained: The creep compliance J(t) is expressed as: The relaxation time T is defined as: 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: 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.  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. 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.

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: 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 2 ) 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.
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π/λ 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. [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.
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: Schematic of a power spectrum (Kang et al. [13]).

11
Using the Boundary Element Method to Simulate Visco-Elastic Deformations of Rough… DOI: http://dx.doi.org /10.5772/intechopen.96229 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 : 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].

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. 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.
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 1.
The parameters of the seven synthetic surface pairs.    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.

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.

Figure 11.
Aperture fields for different RMS values (Kang et al. [13] 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). 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.

Figures 12 and
Based on Figures 12 and 13, several conclusions can be drawn:   Table 4. Effect of surface parameters on the mean aperture and contact ratio. 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. 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.

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 Table 6.
Effect of surface parameters on the mean aperture and contact ratio.

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