Study on the Propagation Process Characteristics of Anisotropic Acoustic Waves in Shale Gas Well with the Reflection Rule of Lateral Fractures

Abstract

Based on the acoustic wave finite element (AWFE) method, one can establish an AWFE method and study the influence of mechanical parameters on the shale reservoir acoustic wave propagation characteristics. The different crack characteristics and different crack multi-physical coupling phenomena are studied by using the AWFE method on shale gas reservoir cracks. To calculate the shape and position along the crack near a side borehole, the model parameters are compared with the simulation results. The reflection waveform characteristics of adjacent cracks are studied by using the AWFE method. By considering the borehole axis of symmetry, for an acoustic impedance discontinuous interface on one side of the two-dimensional axisymmetric AWFE, one can establish a borehole cross-crack and an arc cross-crack reflection interface model with the AWFE method. By processing the waveform data received by different receiving points under the same source distance, the parameters, such as the reflection wave time and the distribution laws of the crack in the shale reservoir, are obtained. To verify the validity of the research method, the propagation of a reflected acoustic wave from the reservoir fracture by the filling with different media was also studied. The results show that a reflection wave arrival time changing with the source ordinate and present law, by side borehole fracture morphology, showed a suitable consistency. The well cross-crack angle range is 10~20°, according to the wave arrival time calculated by the side borehole fracture dip. For the acoustic signal propagation in the shale formation anisotropy, they found that an acoustic signal is always in the direction of the elastic modulus, with a further larger spread, a location of maximum amplitude, and a 45-degree direction to the axis of symmetry. In the lateral and longitudinal distance from an acoustic source of the same two receiver signal waves, the receiver vibration amplitude is bigger, and there is less attenuation. With the increase in the anisotropic index, the inside ovality amplitude distribution of the signal amplitude in this model is higher. When a side borehole has an arc crack and a reflected wave to time to obtain the coordinates of a reflecting interface and to compare with the results of the calculation model, the crack in the center position, and the reflection point coordinates, the relative error is less than 5%. Finally, the AWFE method could provide a new idea for the identification of the crack properties and also could be an inverse calculation of the position and morphological characteristics of fractures near the side borehole.

1. Introduction

With the large consumption of conventional energy, such as coal, oil, and other resources, the unconventional oil and gas reservoirs represented by shale gas reservoirs have attracted high attention around the world. The fracture is a direct flow channel of shale gas resources, whether it is hydraulic fracturing, fractures, or natural fractures. These fractures will greatly increase the oil and gas production of shale reservoirs. Therefore, it is of great significance to study the shape and location characteristics of shale reservoir fractures [1,2,3,4,5,6].

A lot of research has been conducted on sonic logging worldwide. On the one hand, the propagation characteristics of acoustic waves in the structure are studied [7,8,9,10,11,12,13,14,15,16]. Modave et al. [7], for example, and Miller et al. [8] demonstrate efficient finite element solutions of exterior acoustic problems with truncated computational domains surrounded by perfectly matched layers, combining projection-based model order reduction with an efficient time domain-impedance boundary condition formulation. Zeng et al. [9] proposed a physical model of a phased-array laser source for detecting surface defects. The parametric performance of the finite element simulation was optimized by comparing the amplitude and spectral characteristics of the ultrasonic waves generated by the phased-array laser source at different element spacings. Lin et al. [10] explored the feasibility and reliability of combining back-propagation neural networks (BPNN) and laser-generated surface acoustic waves (SAW) to detect the depth and length of subsurface cracks. The characteristic parameters of the emitted and reflected surface waves in the time and frequency domains were extracted by the Fourier transform and wavelet transform methods. Mohammadgholiha et al. [11] proposed a new method for frequency controllable acoustic transducers (FSAT) for guided wave detection and greatly simplified the hardware and reduced the cost of based systems. Tang et al. [12] classified shale microfracture processes into two types based on the characteristics of acoustic emission activity cycles and used digital image processing techniques to characterize the geometric features and inhomogeneous distribution of calcite minerals in shale at the microscale, providing an important theoretical guide for shale pneumatic fracture mining. Li et al. [13] combined acoustic emission (AE) monitoring and computed tomography to perform laboratory fracturing of shale samples containing multilayered facies and, based on the results, found that the pressure profile and acoustic emission response characteristics reflect well the growth behavior of hydraulic fractures in laminated shales. Wu et al. [15] analyzed the mechanical properties and damage modes of laminated shales under the coupling of seepage pressure and stress. It was found that the fractal dimension better reflected the fracture damage of shale. Guo et al. [16] considered the effects of shale anisotropy, shale stiffness degradation, and formation sliding on casing deformation and obtained the elastic modulus of shale rocks after obtaining water saturation at different times.

The above methods include the finite element method, the finite difference method, and physical simulation. On the other hand, the processing of well logging data is studied. Qiao et al. [17] have made a physical simulation model of aluminum and concrete, which can observe the refraction and reflection waves transmitted through the well wall and well side fracture. Che Xiaohua et al. [18,19] evaluated the geometrical characteristics of well side fractures by experiments and numerical simulations, respectively. Furthermore, according to the signal processing results, the dip angle and position of the fractures are effectively deduced. Wei Zhoutuo et al. [20] studied the well side fracture characteristics by numerical simulation. In addition, he put forward a method of extracting the reflected waves, which can effectively enhance the signal intensity of the reflected waves.

Based on the above research, this paper is based on the theory of wave propagation. The reflected wave characteristics were studied by the finite element method when there were different shapes and positions of fractures near the well. The reflected wave signals were extracted and processed, and the distribution of the fractures in the shale reservoir was determined.

2. Theoretical and Boundary Condition

2.1. The Equation of Acoustic Propagation in Shale Rock

According to the law of conservation of mass, the continuity equation in acoustic technic [21,22] can be expressed as:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho V_x\right)}{\partial x}+\frac{\partial \left(\rho V_y\right)}{\partial y}+\frac{\partial \left(\rho V_z\right)}{\partial z}=0$$

(1)

where $V_x$, $V_y$, and $V_z$ are the components of the velocity vector in the $x,y,z$ Cartesian coordinate directions; $\rho$ denotes density; $x,y,z$ is the global coordinate system; and $t$ represents time.

To simplify the acoustic equation, one at first assumes that the fluid is compressible and non-flowing, and the following acoustic equation can be obtained, according to the Navier–Stokes and the continuity equations [21].

$$\nabla \cdot \left(\frac{1}{{\rho }_0}\nabla p\right)-\frac{1}{{\rho }_0{c}^2}\frac{{\partial }^{2}p}{\partial t^2}+\nabla \cdot \left[\frac{4\mu }{3{\rho }_0}\nabla \left(\frac{1}{{\rho }_0{c}^2}\frac{\partial p}{\partial t}\right)\right]=-\frac{\partial }{\partial t}\left(\frac{Q}{{\rho }_0}\right)+\nabla \cdot \left[\frac{4\mu }{3{\rho }_0}\nabla \left(\frac{Q}{{\rho }_0}\right)\right]$$

(2)

where ${\rho }_0$ is the fluid density, c is the acoustic velocity ($\sqrt{K/{\rho }_0}$) in the fluid media, $\mu$ is the dynamic viscosity, $K$ denotes the volume modulus, $t$ represents time, $p$ is the acoustic pressure, and $Q$ is the continuous equation quality source.

The acoustic pressure equation is as follows:

$$p\left(\overrightarrow{r},t\right)=\mathrm{Re}\left[p\left(\overrightarrow{r}\right)e^{j\omega t}\right]$$

(3)

where $p$ is the acoustic pressure, of which $j=\sqrt{-1}$, $\omega =2nf$, $f$ denote oscillation pressure frequency.

Then, Equation (2) can be approximately simplified to a non-homogeneous Helmholtz equation [21], as follows:

$$\begin{array}{l}\nabla \left(\frac{1}{{\rho }_0}\nabla p\right)+\frac{{\omega }^2}{{\rho }_0{c}^2}p+j\omega \nabla \cdot \left[\frac{4\mu }{3{\rho }_0}\nabla \left(\frac{1}{{\rho }_0{c}^2}p\right)\right]\\ =-j\omega \left(\frac{Q}{{\rho }_0}\right)+\nabla \left[\frac{4\mu }{3{\rho }_0}\nabla \left(\frac{Q}{{\rho }_0}\right)\right]\end{array}$$

(4)

2.2. Galerkin Method Acoustic Equation

By using the Galerkin method [21], from Equation (1) and Equation (4), a finite element equation of the acoustic wave can be obtained as follows:

$$\begin{array}{l}{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\frac{1}{{\rho }_0{c}^2}}w\frac{{\partial }^{2}p}{\partial t^2}dv+{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\nabla w\cdot \left(\frac{4\mu }{3{\rho }_0^2{c}^2}\nabla \frac{\partial p}{\partial t}\right)}dv+{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\nabla w\cdot \left(\frac{1}{{\rho }_0}\nabla p\right)}dv\\ -{\displaystyle \underset{{\mathsf{\Gamma }}_F}{∯}w\left(\frac{1}{{\rho }_0}+\frac{4\mu }{3{\rho }_0^2{c}^2}\frac{\partial }{\partial t}\right)}\widehat{n}\cdot \nabla pds+{\displaystyle \underset{{\mathsf{\Gamma }}_F}{∯}w}\frac{4\mu }{3{\rho }_0^2}\widehat{n}\cdot \nabla Qds\\ ={\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\frac{1}{{\rho }_0}}\frac{\partial Q}{\partial t}dv+{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\nabla w\left(\frac{4\mu }{3{\rho }_0^2}\nabla Q\right)}dv\end{array}$$

(5)

where $dv$ is the differential of the volume in the acoustic domain ${\mathsf{\Omega }}_F$; $ds$ is the differential of the acoustic domain boundary ${\mathsf{\Gamma }}_F$; and $\widehat{n}$ denotes a unit vector for the outer normal direction of the acoustic domain boundary ${\mathsf{\Gamma }}_F$.

According to the momentum conservation equation, the normal velocity at the boundary of the acoustic domain can be obtained as follows:

$$\frac{\partial v_{n,F}}{\partial t}=\widehat{n}\cdot \frac{\partial \overrightarrow{v}}{\partial t}=-\left(\frac{1}{{\rho }_0}+\frac{4\mu }{3{\rho }_0^2{c}^2}\frac{\partial }{\partial t}\right)\widehat{n}\cdot \nabla p+\frac{4\mu }{3{\rho }_0^2}\widehat{n}\cdot \nabla Q$$

(6)

Substituting Equations (1)–(6) into Equation (5), the following can be obtained:

$$\begin{array}{l}{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\frac{1}{{\rho }_0{c}^2}}w\frac{{\partial }^{2}p}{\partial t^2}dv+{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\nabla w\cdot \left(\frac{4\mu }{3{\rho }_0^2{c}^2}\nabla \frac{\partial p}{\partial t}\right)}dv+{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\nabla w\cdot \left(\frac{1}{{\rho }_0}\nabla p\right)}dv\\ +{\displaystyle \underset{{\mathsf{\Gamma }}_F}{∯}w}\frac{\partial v_{n,F}}{\partial t}ds={\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\frac{1}{{\rho }_0}}\frac{\partial Q}{\partial t}dv+{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\nabla w\cdot \left(\frac{4\mu }{3{\rho }_0^2}\nabla Q\right)}dv\end{array}$$

(7)

By using the method of normal displacement of the fluid particle, the acceleration can be expressed as follows:

$$\frac{\partial v_{n,F}}{\partial t}=\widehat{n}\cdot \frac{{\partial }^2{\overrightarrow{u}}_F}{\partial t^2}$$

(8)

where, ${\overrightarrow{u}}_F$ is the displacement of the flow constitution point.

Substituting Equations (1)–(7) into Equations (1)–(8), the following can be obtained:

$$\begin{array}{l}{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\frac{1}{{\rho }_0{c}^2}}w\frac{{\partial }^{2}p}{\partial t^2}dv+{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\nabla w\cdot \left(\frac{4\mu }{3{\rho }_0^2{c}^2}\nabla \frac{\partial p}{\partial t}\right)}dv+{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\nabla w\cdot \left(\frac{1}{{\rho }_0}\nabla p\right)}dv\\ +{\displaystyle \underset{{\mathsf{\Gamma }}_F}{∯}w\widehat{n}}\frac{{\partial }^2{\overrightarrow{u}}_F}{\partial t^2}ds={\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\frac{1}{{\rho }_0}}\frac{\partial Q}{\partial t}dv+{\displaystyle \underset{{\mathsf{\Omega }}_F}{\iiint }\nabla w\cdot \left(\frac{4\mu }{3{\rho }_0^2}\nabla Q\right)}dv\end{array}$$

(9)

2.3. Set Boundary Conditions and Build Models

The boundary of the acoustic field is a rigid boundary with zero pressure and a free boundary with zero velocity. Robin’s impedance boundary conditions can be expressed as follows [22]:

$$v_{n,F}\left(\overrightarrow{r}\right)-v_{n,s}\left(\stackrel{\rightharpoonup }{r}\right)=Y\left(\overline{r}\right)\rho \left(\overrightarrow{r}\right)$$

(10)

where $v_{n,F}$ is the normal velocity of the fluid particle on the boundary; $v_{n,s}$ is the normal velocity of the structure surface; $Y$ is the boundary admittance; and $Z$ is the boundary impedance and $Z=1/Y$.

Substituting Equations (1)–(10) into Equations (1)–(7), the following can be obtained:

$${\displaystyle \underset{{\mathsf{\Gamma }}_Z}{\iint }w}\frac{\partial v_{n,F}}{\partial t}ds={\displaystyle \underset{{\mathsf{\Gamma }}_Z}{\iint }w\left(Y\frac{\partial p}{\partial t}+\frac{\partial V_{n,S}}{\partial t}\right)ds}$$

(11)

The attenuation coefficient of the material absorption characteristics is defined as the ratio between the input acoustic power and the absorption acoustic power density, namely:

$$\mathsf{\alpha }=\frac{I_a}{I_{inc}}$$

(12)

where $I_a=\frac{1}{2}{P}_a\left(\widehat{n}\cdot {{\overrightarrow{v}}_a}^{\ast }\right)$; $I_{inc}=\frac{1}{2}{P}_{inc}\left(\widehat{n}\cdot {{\overrightarrow{v}}_{inc}}^{\ast }\right)$.

According to the reflection coefficient of the attenuated material, $\left|R\right|=\sqrt{I_{ref}/I_{inc}}=\sqrt{1-\partial }\left(0\le \left|R\right|\le 1\right)$.

Then, the equivalent impedance can be expressed as follows:

$$z=z_0\frac{1+\sqrt{1-\mathsf{\alpha }}}{1-\sqrt{1-\mathsf{\alpha }}}$$

(13)

where the acoustic impedance $z={\rho }_0{c}_0$

The above acoustic impedance boundary conditions are used in the modal analysis and the harmonic response analysis.

2.4. Acoustic Wave Absorption Boundary Conditions

In the study of acoustic problems, when the outside of the structural unit is an infinite uniformly non-viscous liquid, Equations (1) and (2) and Equations (3)–(5) can be expressed as follows [23]:

$$\nabla \cdot \left(\frac{1}{{\rho }_0}\nabla p\right)=\frac{1}{{\rho }_0{c}^2}\frac{{\partial }^{2}p}{\partial t^2}$$

(14)

It is necessary to satisfy the scattering conditions of the wave generated in the fluid, which can be expressed as follows:

$$\underset{r\to \infty }{\lim }{r}^{(\mathrm{d}-1)/2}\left(\frac{\partial p}{\partial r}+\frac{1}{c}\dot{p}\right)=0$$

(15)

where $r$ is the distance from the source, and $d$ is the number of dimensions.

When the finite element model is usually used to simulate the infinite medium, it is necessary to set up the absorption boundary in the structure unit. The two-dimensional analysis model can be expressed as follows:

$$P_n+\gamma P_n=-\frac{1}{c}\ddot{P}+\left(\frac{1}{2}\kappa -\frac{\gamma }{c}\right)P+\frac{1}{2}c{P}_{\lambda \lambda }+\left(\frac{1}{8}{\kappa }^{2}c+\frac{1}{2}\kappa \gamma \right)P$$

(16)

where $n$ is the outer normal direction of the boundary, $k$ is the boundary curvature, $P_n$ is the derivative of the upward pressure, $P_{\lambda \lambda }$ is the derivative of the pressure along a boundary, and $\gamma$ is the stability coefficient.

For 3D cases, this can be expressed as follows [21]:

$$\begin{array}{l}{\dot{P}}_n+\gamma P_n=-\frac{1}{c}\ddot{P}+\left(H-\frac{\gamma }{c}\right)\dot{P}+H\gamma P\\ +\frac{c}{2\sqrt{EG}}\left[{\left(\sqrt{\frac{G}{E}{P}_u}\right)}_u+{\left(\sqrt{\frac{G}{E}{P}_v}\right)}_v\right]+\frac{c}{2}\left(H^2-K\right)P\end{array}$$

(17)

where $n$ is the outer normal direction of the boundary; $u,v$ is the orthogonal curvilinear coordinate; $P_u,P_v$ is the derivative of pressure on the boundary surface; $H,K$ is the Gaussian curvature for the surface of two directions; and $E,G$ is the shale rock elastic modulus and shear modulus.

In combination with the related geometric relations in the analysis of the propagation path of the reflected wave in the above, the position coordinate $\left(x_p,y_p\right)$ of the reflection point is derived when the wave source is in the $\left(x_0,y_0\right)$ position:

$$\left(\sqrt{{\left(\frac{b^{\prime }{t}_p{v}_p}{2b+L^{\prime }}\right)}^2-{b^{\prime }}^2{\sin }^2\theta }+b^{\prime }\cos \theta \right)\sin \theta +0.1=x_p$$

(18)

$$\left(\sqrt{{\left(\frac{b^{\prime }{t}_p{v}_p}{2b+L^{\prime }}\right)}^2-{b^{\prime }}^2{\sin }^2\theta }+b^{\prime }\cos \theta \right)\cos \theta -b^{\prime }+a\tan {\theta }^{\prime }+y_0=y_p$$

(19)

When the inclination angle is negative, the coordinate $\left(x_p,y_p\right)$ of the reflection point is calculated as follows:

$$\left(b^{\prime }\cos \theta -\sqrt{{\left(\frac{b^{\prime }{t}_p{v}_p}{2b-L^{\prime }}\right)}^2-{b^{\prime }}^2{\sin }^2\theta }\right)\sin \theta +0.1=x_p$$

(20)

$$b^{\prime }-\left(b^{\prime }\cos \theta -\sqrt{{\left(\frac{b^{\prime }{t}_p{v}_p}{2b+L^{\prime }}\right)}^2-{b^{\prime }}^2{\sin }^2\theta }\right)\cos \theta +a\tan {\theta }^{\prime }+y_0=y_p$$

(21)

2.5. Perfect Match Layer Boundary Condition

Perfect match layer boundary conditions can absorb the incident wave in any direction other than the parallel direction of the PML layer. The properties of the PML material can be expressed as follows [21]:

$$\nabla p=-j\omega {\rho }_0\overrightarrow{v}$$

(22)

$$\nabla \cdot \overrightarrow{v}=-\frac{1}{{\rho }_0{c}^2}j\omega p$$

(23)

$$\overline{\overline{\rho }}={\rho }_0\left[\mathsf{\Lambda }\right]$$

(24)

$$\overline{\overline{\mathsf{\zeta }}}=\frac{1}{{\rho }_0{c}^2}\left[\mathsf{\Lambda }\right]$$

(25)

where $\left[\mathsf{\Lambda }\right]$ represents a diagonal matrix of constants.

In the vector calculation, the control equation can be expressed as follows:

$$\nabla p=-j\omega \overline{\overline{\rho }}\cdot \overrightarrow{v}$$

(26)

$$\nabla \cdot \left({\overline{\overline{\mathsf{\zeta }}}}^{-1}\cdot \overrightarrow{v}\right)=-j\omega p$$

(27)

If the correction operator is expressed as

$$\nabla S=\frac{1}{s_x}\frac{\partial }{\partial x}\widehat{x}+\frac{1}{s_y}\frac{\partial }{\partial y}\widehat{y}+\frac{1}{s_z}\frac{\partial }{\partial z}\widehat{z}$$

(28)

then the acoustic propagation of the PML layer can be expressed as follows:

$${\nabla }_s\left(\frac{1}{{\rho }_0}{\nabla }_{s}p\right)+\frac{{\omega }^2}{{\rho }_0{c}^2}p=0$$

(29)

In order to realize the effect of the complete matching layer boundary condition, $\left[\mathsf{\Lambda }\right]$ can be expressed as follows:

$$s_i=1+\frac{{\mathsf{\sigma }}_i}{j\omega }\hspace{1em}\left(i=1,2,3\right)$$

(30)

where ${\mathsf{\sigma }}_i$ is an attenuation coefficient.

Then, $\left[\mathsf{\Lambda }\right]$ multiplied by Equations (3)–(28) can be expressed as follows:

$$\nabla \left(\frac{1}{{\rho }_0}\left[\mathsf{\Psi }\right]\cdot \nabla p\right)+\frac{{\omega }^2}{{\rho }_0{c}^2}{s}_x{s}_y{s}_{z}p=0$$

(31)

In the finite element calculation, if the boundary condition of the complete matching layer is used, the excitation source cannot be set at the boundary, in order to prevent the calculation error.

3. Case Study: Propagation Process Characteristics of Anisotropic Acoustic Waves in Shale Gas Well

When there is a medium in the well, the propagation path of the reflected wave in reflected sonic logging is shown in Figure 1 [24]. The sonic propagation velocities in the well and the formation medium are $v_1$ and $v_{\mathrm{p}}$, respectively. The distance between the acoustic transmitter and the acoustic receiver is L, according to the law of acoustic wave propagation:

$$\frac{v_1}{\sin {\theta }_1}=\frac{v_p}{\sin {\theta }_2}$$

(32)

$$\frac{v_1}{\sin {\theta }_5}=\frac{v_p}{\sin {\theta }_4}$$

(33)

When the well fracture pattern is parallel to the well, the well radius is a, the distance between the well centerline and the fracture is X, and the reflected wave travel time is determined by the geometric relation:

$$t_{pp}=\frac{a}{v_1\cos {\theta }_5}+\frac{a}{v_1\cos {\theta }_1}+\frac{\sqrt{{L^{\prime }}^2+4{\left(X-a\right)}^2}}{v_p}$$

(34)

$$L^{\prime }=L-a\tan {\theta }_1-a\tan {\theta }_5$$

(35)

When the reflection interface is not parallel with the well direction, assume that the crack angle is θ, the distance between the intersection point of the extension line of the reflecting interface and the centerline of the well and the acoustic signal transmission point is b; then:

$${\theta }_2={\theta }_3-\theta$$

(36)

$${\theta }_4={\theta }_3+\theta$$

(37)

$$\cot {\theta }_3=\frac{L^{\prime }+2b^{\prime }}{L}\tan \theta$$

(38)

The receiving time of the reflected wave is deduced from the geometric relation of the triangle:

$$t_{pp}=\frac{a}{v_1\cos {\theta }_5}+\frac{a}{v_1\cos {\theta }_1}+\frac{\sqrt{{\left(L^{\prime }+b^{\prime }\right)}^2+{b^{\prime }}^2-2b^{\prime }\left(b^{\prime }+L^{\prime }\right)\cos 2\theta }}{v_p}$$

(39)

$$b^{\prime }=b-a\cot \theta +a\tan {\theta }_1$$

(40)

Then:

$$t_{pp}=\frac{a}{v_1}\left[\frac{1}{\cos {\theta }_5}+\frac{1}{\cos {\theta }_1}\right]+\frac{\sqrt{{L^{\prime }}^2+4b^{\prime }\left(L^{\prime }+b^{\prime }\right){\sin }^2\theta }}{v_p}$$

(41)

When the reflection dip angle is opposite, the intersection of the reflection interface extends, and the well is at the top of the acoustic source. The traveling time of the reflected wave is:

$$b^{\prime }=b-a\cot \theta -a\tan {\theta }_1$$

(42)

$$t_{pp}=\frac{a}{v_1}\left[\frac{1}{\cos {\theta }_5}+\frac{1}{\cos {\theta }_1}\right]+\frac{\sqrt{{L^{\prime }}^2+4b^{\prime }\left(b^{\prime }-L^{\prime }\right){\sin }^2\theta }}{v_p}$$

(43)

3.1. Calculation of Propagation Path of Reflected Acoustic Wave

In the actual application of sonic logging, there are mud and well logging tools in the well. The acoustic signals transmitted by the phased array acoustic sources are reflected by the borehole and the formation and then received by the receiving array in the borehole. The received signals include slide waves, direct waves, and reflected waves.

The main content of this paper is the simulation of acoustic reflection logging by the finite element method. Taking the borehole as the boundary, a two-dimensional finite element simulation of the fractured side of the reservoir is carried out. Cross-shaped fracture and arc-shaped fracture finite element models are established and are shown in Figure 2 and Figure 3. The locations of each coordinate point are shown in

Table 1. Key point coordinates of cross-shaped fracture.

	a	b	c	d	e
x	0.1	6.38	0.1	5.04	3.1
y	0	0	20	20	9

and

Table 2. Key point coordinates of arc-shaped fracture.

	a	b	c	d	e
x	0.1	8	0.1	8	4
y	0	0	20	20	9

. The sizes of ${\theta }_1$ and ${\theta }_2$ are 10 and 20 degrees, respectively. Several acoustic transmitters and receivers are located along the center of the borehole. The free boundary on the right is used to simulate the reflecting interface. The Y dotted line represents the centerline of the borehole, and the well depth is 20 m. The distance between the borehole centerline and the borehole wall is 0.1 m. Both sides of the wall of the well are the plane elements of the mud in the well and the plane element of the shale reservoir. The material parameters are shown in

Table 3. Material parameter table.

Material Attributes	Sonic Velocity/(m/s)	Density/(kg/m3)
Shale reservoir	3675	2500
Well	2900	2000

.

3.2. Excitation Source

In the aspect of the excitation source selection, the acoustic source signal is a Gauss sine wave with a center frequency of 10 kHz [25], and the expression of the acoustic source signal is:

$$S(t)=A{e}^{-k{(t-t_0)}^2}\cos \left(2\pi f_0\left(t-t_0\right)\right)$$

(44)

$$k=\frac{5{\pi }^2{b}^2{f}_0^2}{q\times \ln \left(10\right)}$$

(45)

In the upper case, the amplitude A is 1, B is the normalized bandwidth (0.8), q is the attenuation (6 dB), $f_0$ is the center frequency of the acoustic source (10 kHz), and $t_0$ is the delay time (150). The waveform of the source signal is shown in Figure 4.

3.3. Shale Gas Well AWFE Parameter Setting

The physical model of the shale formation structure is established based on ANSYS. A parameter of the model is defined, and a physical model of the structure is established. At the same time, to reduce the computer calculation burden, a certain proportion can be passed to reduce the size of the model and to set up some special boundary conditions to simulate the infinite formation. The build radius for the circular area of 5 m is the simulated formation; the setting of the physical performance parameters of the elastic modulus is 30 GPa. Poisson’s ratio is 0.22. The density is 2500 kg/m³. In the middle of a circular area, the radius R is set as 20 mm round holes to simulate the shale gas well eye. Then, one can launch acoustic signals in the shale gas wells, and different receiving points are set up in the stratum to study the shale reservoir acoustic wave propagation (SRAWP) in the formation.

3.3.1. Boundary Conditions and Loading of AWFE

In the process of the acoustic field calculation, the finite element method is used to simulate the acoustic wave propagation in the infinite stratum, and the absorption effect of the signal wave can be realized by using the absorbing boundary condition. In the simulation of the finite element acoustic field, the Fluid29 and Fluid30 elements are usually used to simulate the fluid part, and then, the Fluid129 and Fluid130 elements are used together. The solid model can only be contacted with Fluid29 and Fluid30 units. The Fluid129 and Fluid130 elements surround the internal structural units and fluid units.

Using AWFE, the Fluid29 and Fluid129 elements are used to simulate and calculate, respectively. The AWFE mesh boundary needs to be two-dimensional and axisymmetric. Some studies have shown that the radius of the lateral absorbing boundary surrounded by this unit should be specified as the real constant rad.

This needs to be achieved: $rad\ge (\frac{D}{2}+0.2\frac{c}{f})$. This is condition [26], where D is the model radius, C is the wave velocity, and f is the frequency of the acoustic wave propagation.

3.3.2. Mesh Generation and Material Parameter Setting

The model is a complete axisymmetric structure. The mapping mesh method is used. The model diagram is shown in Figure 5 after the mesh division. At the source of the acoustic wave, the mesh is finer, the length of the unit is about 0.2 mm, the grid at the boundary of the model is coarse, and the number of units is 59,287. The model represents the stratigraphic unit and the absorption unit of the two layers successively, from the inside to the outside. The parameter KEYOPT2 of the two-layer absorption unit is set to 0 and 1, respectively.

In terms of material properties, in addition to the structural elements, there are also contact units, fluid units, and boundary absorption units in the model. The fluid unit needs density and acoustic velocity as the material attribute FLUID129 unit only needs acoustic velocity. The boundary admittance β (absorption coefficient) of the fluid-structure interface is used to express the absorption of the acoustic velocity, and the value β is usually measured by experiments. The cell parameter settings are shown in

Table 4. Unit parameter setting.

	Density (kg/m3)	Sound Velocity (m/s)	Admittance Coefficient β
Contact unit FLUID29KEYOPT(2) = 0	1.21	340	/
Fluid unit FLUID29KEYOPT(2) = 1	1.21	340	/
Absorption boundary FLUID129	1350	1500	1

.

3.4. Exciting Load

The AWFE analysis is carried out in ANSYS, and the transient dynamic analysis method is used to simulate the excitation and propagation of the guided wave. The time period of the calculation is 0–0.5 ms. The 10 kHz sinusoidal signal with five cycles in a single cycle is chosen as the excitation signal of the wave source. A cyclic signal is shown in Figure 6.

4. Results and Discussion

4.1. Effect of Rock Mechanics Parameters on Acoustic Wave Velocity

In order to facilitate the discussion of the results, this paper mainly analyzes the variation of the acoustic wave propagation velocity in shale strata when the elastic modulus changes at 4 GPa intervals. For the acoustic source settings, as shown in Figure 7, the exact rate of simulation is applied at the central borehole of the shale formation model to simulate the propagation of the borehole acoustic signals in the formation. A reception point is set at the boundary of the formation unit, and the position of the reception point is shown in Figure 8.

When the elastic modulus is 30 GPa and Poisson’s ratio is 0.22, the waveform diagram of the receiving point is shown in Figure 9. It can be seen that the acoustic wave signal is received in the shale formation at this position at about 1.38 ms, and the received waveform is approximately the same as that of the source signal. Compared with the maximum amplitude of the source signal, the signal frequency is basically unchanged, and the maximum amplitude is reduced by 96.9 kyr. According to the position relation, the velocity of the acoustic wave propagation is about 3616 m/s.

Comparing the waveform of the receiving point with the different elastic moduli, the amplitude attenuation is roughly between 96% and 97%. As shown in Figure 10, according to the time of the receiving point, the propagation velocity of the acoustic wave signal in the shale formation is obtained under different elastic moduli. The results are shown in

Table 5. Acoustic velocity in shale formations with different elastic moduli.

E/GPa	18	22	26	30	34
Propagation time/ms	1.75	1.61	1.48	1.38	1.28
Propagation speed/(m/s)	2851	3099	3371	3616	3898

.

It can be found that when the modulus of elasticity changes in 18 GPa–34 GPa, the velocity of the acoustic wave in the shale strata is in the range of 2800 m/s~3900 m/s, and the relationship between the velocity of the acoustic wave propagation and the modulus of elasticity is shown in Figure 11.

According to the fitting results, the velocity of the acoustic wave in the shale strata increases with the increase of the elastic modulus of the formation, and the expression of the relationship between them can be obtained as approximately:

$$V=65.275E+1669.85$$

(46)

In addition to the elastic modulus, Poisson’s ratio is another important rock mechanics parameter in shale formation. The Poisson’s ratio of shale strata is generally about 0.25. The variation law of acoustic velocity is mainly analyzed when the Poisson’s ratio varies from 0.2–0.28. The calculation and setting, including the selection of the acoustic source and receiving point, the excitation signal, and so on, are the same as the elastic modulus analysis above.

As shown in Figure 12, comparing the waveform of the receiving point with different Poisson’s ratios in shale formation, it can be seen that when the Poisson’s ratio changes in the range of 0.2–0.28, the amplitude attenuation of the receiving point remains the same; the amplitude is about the same, compared with the amplitude of the signal wave source. When the Poisson’s ratio changes in a certain range, the receiving time of the signal wave reception point is basically the same. According to the calculation of the reception point time, the propagation velocity of the acoustic wave signal in shale formation is about 3616 m/s at different Poisson’s ratios. It is worth noting that when the Poisson’s ratio varies from 0.2 to 0.28, the velocity of the acoustic wave in the shale strata remains basically unchanged and is always around 3616 m/s. That is, the Poisson’s ratio has little effect on the velocity of the acoustic wave propagation in shale formation.

4.2. Analysis of Acoustic Wave Propagation Characteristics in Isotropic Shale Formation

Using the ANSYS finite element method to simulate the acoustic wave propagation characteristics in shale formation is an important means of understanding the specific propagation of the acoustic wave in acoustic logging.

The elastic modulus is 30 GPa and Poisson’s ratio is 0.22. The propagation characteristics of acoustic waves in shale strata are analyzed by calculation. After calculating that, the source excites at 0.5 ms; that is, after the wave source excites a cycle, the waveform diagram of the acoustic wave propagation in shale strata is shown in Figure 13. At this point, the maximum magnitude of the formation is 4.74 × 10⁻⁷. Compared with the maximum value of 1 × 10⁻⁵ of the source signal, it attenuates at about 95.26, which is consistent with the phenomenon of energy attenuation caused by damping during the acoustic wave propagation.

The waveforms of the acoustic waves propagating in the shale strata at different times are compared. Taking the time nodes as 0.5 ms, 1.5 ms, 2.5 ms, and 3.5 ms, respectively, we can find that the largest values in the strata at different times are always near the wave source, which is related to the attenuation degree of the amplitude of the acoustic wave propagation process, and at about three periods; that is, the wave source signal is relayed to the boundary point of the model stratum for 1.5 ms, and the maximum value is that when the signal wave is affected by the absorbing boundary condition, when the signal wave propagates to the formation boundary point, compared with the wave source signal attenuation of nearly 95.28%.

4.3. Analysis of Acoustic Waves Propagation Characteristics in Anisotropic Shale Formation

In order to describe the anisotropy degree of shale formation [27], the anisotropy index e can be defined as:

$$e=\left|1-E_x/E_y\right|$$

(47)

where $E_x$ is the elastic modulus or the Poisson’s ratio is parallel to the isotropic surface, and $E_y$ is the elastic modulus or Poisson’s ratio perpendicular to the isotropic surface. The influence of the anisotropy of the shale formation on acoustic wave propagation is considered, and the influence of the anisotropic exponent e on the acoustic wave propagation characteristics is mainly considered. The setting of the physical performance parameters is shown in

Table 6. Parameter setting of anisotropic media.

Parameter Type	Value
$\rho$$/(kg/m^3$)	2500
$E_x$/GPa	24.91
$E_y$/GPa	14.093
${\mu }_x$	0.324
${\mu }_y$	0.367

.

In order to study the application of the numerical simulation method in anisotropic media and to simplify the calculation, the section of the shale stratigraphic plane in the vertical direction is taken to be simulated and analyzed. The size of the model is 5 × 5 m, and the number of grid units is 137,920. The parameters of each unit are set up as shown in Table 4. The excitation source is set up in the middle of the model, and the excitation signal wave of the wave source frequency is a 10 kHz sinusoidal signal with five cycles.

After the above parameters are set, the acoustic wave finite element simulations are performed on the anisotropic shale formation to analyze the acoustic propagation characteristics of the anisotropic formation. After the excitation of one cycle in the wellbore of the model well, the sonic propagation displacement cloud at 0.5 ms is shown in Figure 14. The acoustic wave diagram shows that the acoustic signal is diffused from the wave source position to the periphery. Due to the periodic variation of the amplitude of the wave source signal, at 0.5 ms, the strength distribution of the model stratigraphic amplitude is consistent with the periodic variation of the wave source signal; with the increasing of the signal source distance, the signal intensity tends to decrease.

In terms of propagation distance, it is parallel to the shale bedding plane; that is, the isotropic plane parallel to the x-axis direction travels farther because of the existence of interlaminar gaps and micro-cracks perpendicular to the bedding plane in the actual shale formation. The reason for this is also consistent with the results of the relationship between acoustic velocity and the rock mechanics parameters analyzed in the previous section. The maximum amplitude is in the direction of 45° with the axis, and the maximum amplitude is 7.08 × 10⁻⁷, compared to the wave source. The signal is attenuated by 92.92%, which is due to the amplitude superposition of the acoustic signal in the different directions of the anisotropic material.

As shown in Figure 14, the waveforms propagating in the anisotropic media at the times of 0.5 ms, 1.5 ms, 2.5 ms and 3.5 ms are compared with the waveforms at the different periods of time when the wave source is excited, and it can be seen that the wave pattern of the acoustic wave propagation in the anisotropic medium is obviously different from that in the isotropic medium. The velocity of the signal wave propagating along the transverse direction is obviously faster than that in the longitudinal direction. The position of the maximum amplitude is not circular and uniform, and the position is 45° with the symmetry axis.

By comparing different periodic waveforms with the passage of time, the wave source signal propagates gradually along the radial direction around; the amplitude of the acoustic wave is distributed alternately, which is consistent with the wave source signal. Because of the attenuation of the amplitude in the acoustic wave propagation process, the maximum amplitude is near the wave source, and the maximum amplitude is 7.08 × 10⁻⁷ at the end of the different excitation periods, which is caused by the combined cause of the periodic excitation signal wave and the amplitude attenuation in the position of the source.

In order to analyze the specific acoustic vibration of each point in the anisotropic media, the radial displacement–time history curves of the points parallel to the stratigraphic plane and (0, 4) are taken, respectively, as shown in Figure 15. The time point of the receiving acoustic wave signal is about 1.13 ms, and the time of the receiving acoustic wave signal is about 1.60 ms. According to the source distance, the average velocity of the acoustic wave propagation is 3530 m/s and 2493 m/s, respectively, combined with each anisotropic medium. The mechanical parameters correspond to the direction. The results are in agreement with the previous analysis of the relationship between acoustic velocity and the rock mechanics parameters. In addition, comparing the amplitude of the two receiving points, the amplitude value of the coordinate (4, 0) is obviously larger than that of coordinate (0, 4), and the acoustic wave signal perpendicular to the angle of the stratigraphic plane in this model attenuates more quickly.

4.4. Analysis of Anisotropy Index on Acoustic Wave Propagation in Shale Formation

In order to study the influence of the anisotropy degree of the shale reservoir on the acoustic wave propagation characteristics, different anisotropic indices are taken, and the corresponding parameters are set as shown in

Table 7. Anisotropic exponential parameter setting.

$\mathit{e}$	0.3	0.4	0.5	0.6
$E_x/GPa$	40	40	40	40
$E_y/GPa$	28	24	20	16

.

The acoustic source setting is the same as above. The acoustic propagation cloud images of different anisotropic exponential models at the 1 ms moment of acoustic source excitation are taken, as shown in Figure 16. It can be seen that the maximum amplitude of the model is higher, and the maximum amplitude of the model increases with the increase of the anisotropic exponent, and the position of the maximum amplitude is 45° with the axis of symmetry. It can be seen from the cloud map that the intensity distribution of the acoustic waves is ellipsoidal due to the heterogeneity, and the ellipticity increases with the increase of the heterogeneity indices. This is due to the difference of the mechanical parameters in different directions, which leads to the difference of the velocity of the acoustic wave signal and then presents different ellipticity.

In the deep well direction, the acoustic wave transmitter is set at a depth of 3–14 m, and the acoustic receiver is set at a depth of 6–17 m. Figure 17 shows the waveform of the received signal wave when the source distance is 3 m. The waveform wave after superimposing the signal wave data is shown in Figure 18.

In Figure 17, with the source distance of 3 m, the received signal wave is a slide wave signal with a very small amplitude, which then receives the direct and total reflected waves transmitted from the well. By Figure 18, the weak signal wave is caused by the mesh error in the finite element calculation at about 1–1.5 ms, which does not affect the recognition of the reflected wave signal.

In addition, when the source longitudinal coordinates are near 7 m, the intensity of the reflected signal wave is weaker, which is due to the reflected signal path passing through the fracture at the side of the well, causing the further diffusion of the signal wave. As a whole, the receiving time of the reflected wave signal can be clearly observed. When the emitter position is below 7 m, the receiving time of the reflected wave signal becomes smaller with the increase in the longitudinal coordinates of the source. Furthermore, when the emitter position is above 8 m, the receiving time of the reflected wave increases with the increase in the longitudinal coordinates of the source. As the dotted line in the picture shows, the variation of the arrival time of the reflected wave is consistent with the shape of the formation fracture. The receiving time at different acoustic source locations is shown in

Table 8. Wave arrival time at different source locations.

Acoustic source ordinate/m	3	4	5	6	7	8
$t_{pp}$/ms	2.61	2.43	2.25	2.10	1.92	1.93
Acoustic source ordinate/m	9	10	11	12	13	14
$t_{pp}$/ms	2.02	2.11	2.20	2.30	2.38	2.47

.

In this model, when ${\theta }_2={90}^{°}$, according to the speed of the acoustic, ${\theta }_1={55.36}^{°}$ is obtained, and this angle is the critical angle ${\theta }^{\prime }$. In the critical case, the signal wave travels farthest in the borehole, and when ${\theta }_1>{\theta }^{\prime }$, the signal wave propagates in the borehole only. Therefore, when the reflection wave propagation path is $0<{\theta }_1<{\theta }^{\prime }$ and a is the borehole radius, the propagation time of the reflected wave in the well is satisfied:

$$\frac{2a}{v_1}<t_0<\frac{2a}{\cos {\theta }^{\prime }{v}_1}$$

(48)

As the model size and the selection of the source distance changes, the ${\theta }_1$ of the reflected wave signal changes as well. When ${\theta }_1={\theta }_5=\raisebox{1ex}{${\theta }^{\prime }$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, the reflection wave receiving the time data is divided into two parts, taking part of the results of the two points (0, 10) and (0, 11) into the Equations (4), (8), and (9) and taking part of the results of the two points (0, 4) and (0, 5) into the Equations (4), (11) and (12), respectively; when calculated, the well fracture angles are 10.13 degrees and 20.19 degrees, and the b is 18.34 and 13.46.

The error rates of the dip angles are 1.3% and 0.95%; the error rates of b are 1.3% and 0.51%. The source of error is mainly from the reading of the receiving time and the approximate calculation of the angle of incidence in the well. In a word, the calculation results reflect the shape and location of the fracture near the wellbore and verify the reliability of the simulation results.

As before, the source position varies from 3–14 m, and the signal receiving point varies from 6–17 m, keeping the source distance at 3 m. The waveform after superimposing the signal wave data is shown in Figure 19. From Figure 8, when the well fracture shape is an arc, with the change of the acoustic source location, the morphological change of the reflected wave receiving time is presented in an arc. In addition, when the acoustic source is near 8–9 m, the reception time of the reflected signal wave is almost the same. When the source is near 8–9 m, the propagation path of the reflected signal is approximately 0 degrees through the fracture cutting angle. The wave receiving time of the reflected wave is shown in

Table 9. Arrival time of reflected wave.

Acoustic source ordinate/m	3	4	5	6	7	8
$t_{pp}$/ms	2.80	2.72	2.60	2.51	2.45	2.41
Acoustic source ordinate/m	9	10	11	12	13	14
$t_{pp}$/ms	2.41	2.47	2.50	2.58	2.63	2.70

. According to the variation law of the receiving time of the reflected wave, the location and shape feature of the well side fracture can be basically determined.

According to the different source points, the coordinates of the reflection points of the fractures beside the well are calculated separately. The reflected interface coordinates are calculated from the receiving time of the reflected wave, as shown in

Table 10. Reflection interface coordinates.

$x_p$/m	4.748	4.525	4.357	4.225	4.127
$y_p$/m	5.278	6.638	7.324	8.035	8.550
$x_p$/m	4.034	4.172	4.183	4.368	4.442
$y_p$/m	9.983	11.239	11.791	13.045	13.854

. In a word, the location and shape of the fracture near the well can be effectively deduced according to the receiving time of the reflected wave. The reliability of the finite element method is further verified.

5. Conclusions

The fracture is a direct flow channel of shale gas resources; so, it is very urgent to study the characteristics of reservoir fractures. In this paper, based on the theory of acoustic wave propagation, the reflected wave characteristics are studied by AWFE in the presence of side fractures. The receiving time of a reflected wave is analyzed according to the propagation path of the reflected wave. In addition, the finite element method is used to study the reflected wave signals under different shapes of fractures near the well. The conclusions are as follows:

(1)

The amplitude of the acoustic signal in the well is higher than that in the formation, and the amplitude decreases gradually with the increase in the propagation time. The velocity of a slide wave propagating along the wall of a well is the fastest, and the signal intensity of the slide wave is obviously weaker than that of the signal transmitted in the well.

(2)

When there is a cross fracture near the well, the receiving time of the reflected wave is obtained by processing the waveform data received by different receiving points under the same source distance condition. The variation of the arrival time of the reflected wave is consistent with the shape of the formation fracture.

(3)

When there is an arc-shaped fracture near the well, the coordinates of the reflecting interface can effectively deduce the location and shape feature of the side borehole fracture. According to the basic principle of reflection wave logging, the coordinate calculation formula for the crack reflection point is derived. By comparing these coordinates with the shape and position of the reflecting interface in the model, the method can effectively reverse the location and morphological features of the fracture near the well.

(4)

When the modulus of elasticity changes from 18 GPa to 34 GPa, the velocity of the acoustic waves in the shale strata is in the range of 2800 m/s–3900 m/s, and the velocity of the acoustic waves in the shale strata increases with the increase of the elastic modulus of the shale formation. When the Poisson’s ratio changes from 0.2–0.28, the velocity of the acoustic waves in the shale changes. The propagation velocity of an acoustic wave in a shale formation remains basically unchanged and is always around 3616 m/s. That is, Poisson’s ratio has little effect on the velocity of acoustic wave propagation in a shale formation.

(5)

Analyzing the acoustic wave propagation in shale strata at different times, it is found that, with the passage of time, when the intensity of acoustic signal is decaying continuously, the maximum value of the model is 4.74 × 10⁻⁷, the amplitude of the signal attenuation is about 95.26, and the absorbing boundary condition can effectively absorb the acoustic wave signal.

(6)

It is found that when the acoustic wave signal propagates further in the direction of the high elastic modulus and the position of the maximum amplitude is 45° with the symmetry axis, the acoustic wave propagation velocities of the two points are 3530 m/s and 2493 m/s, respectively, and the amplitude of vibration of the receiving point parallel to the direction of the stratification plane is larger; that is, the attenuation is less. With the increase of the heterogeneity index, the ellipticity of the signal amplitude distribution in the study model is higher.