Simulation of acoustic wave propagation in a borehole surrounded by cracked media using a finite difference method based on Hudson's approach

Logging and drilling data show that there are many aligned cracks in crustal rock (Crampin 1985). The cracks in formations are caused by directional tectonic stress. The characteristics of acoustic wave propagation are transversely isotropic in the cracked rock (Crampin 1978). With the development of oil and gas exploration, cracked media have become an important topic in reservoir research. Continental clastic and carbonate reservoirs have widely developed cracks. Cracks not only provide space for fluid but also obviously increase the permeability of the reservoir fluid. It is an important geophysical factor in oil exploration and exploitation.

Detection and identification of cracked reservoirs has become a difficult problem to solve. In recent years, researchers have paid a great deal of attention to cracked media.

There are two main effective-medium theories for cracked media. One was proposed by Schoenberg et al (1995) based on the normal and tangential compliance from linear slide compliance theory. The other one was proposed by Hudson (1990, 2000, 2010) based on wave scattering theory for high frequency. The former is not limited to pore shape and microstructure. It has better applicability to different cracks, but does not consider the factors of rock porosity and permeability. The latter can be applied to calculate the acoustic field of a cracked medium that contains uniform cracks, and the aspect ratio is small.

For cracked effective-medium theories, many researchers have carried out a lot of work. Andrey et al (2000a, 2000b, 2000c) gave a method for calculating crack parameters from seismic data. Cheng (1993) proposed an effective-medium theory for transversely isotropic (TI) cracked media. Kazatchenko et al (2004) simulated the full waveforms generated by an impulse source placed on the borehole axis with the method integral at the real axis. Grechka (2005) validated Hudson's approach by numerical simulation with the finite element method (FEM).

Zhang et al (2014) simulated seismic wave propagation in two-phase anisotropic media by velocity–stress wave equations. Rao and Wang (2015) analyzed seismic attenuation in fractured media based on numerical experiments for frequency-domain wave equations with three different parameters.

To study the influence of cracks, this paper simulated the borehole acoustic field propagation in vertical transverse isotropic (VTI) cracked formations. The paper discusses the influence of crack parameters on borehole acoustic field propagation.

Hudson's approach was derived from theory of the average wave-field scattering for high frequency. It is used to calculate the elastic modulus of an elastic medium containing thin coin-shaped cracks. This approach is applied to simulate acoustic field propagation in a cracked formation with the following conditions: crack density in the formation is less than 0.1, coin-shaped, the crack radius is small, and the crack interval is much smaller relative to the wavelength.

Hudson gave a method to calculate the elastic modulus. In isotropic media, if there are aligned coin-shaped cracks that are independently distributed, Hudson considered the medium as a transversely isotropic (TI) medium. If the isotropic medium is regarded as a background medium, the elastic modulus of the TI medium can be calculated from the elastic modulus of the background medium and crack parameters, and can be presented in the following form

$$\begin{eqnarray}&&C={C}^{0}+{C}^{1}+{C}^{2}+O(v{a}^{3})\end{eqnarray} \tag{ 1 }$$

where ${C}^{0}$ is the matrix formed by the elastic parameters of the background medium, ${C}^{1}$ is the first-order corrected matrix and regards each crack as having independent behavior, ${C}^{2}$ is the second-order corrected matrix and regards crack coupling as pairwise mutual, $O(v{a}^{3})$ is the higher-order perturbation produced by $v{a}^{3},$ $v$ is the crack density, and $a$ is the average radius of the cracks. In general, the higher-order perturbations can be ignored.

For the normal direction of cracks aligned along the Z-axis, cracked media exhibit vertical transversely isotropic symmetry. The background medium can be presented in the following form

$$\begin{eqnarray}{C}^{0}=\left[\begin{array}{ccccccc}\lambda +2\mu & & \lambda & \lambda & 0 & 0 & 0\\ \lambda & & \lambda +2\mu & \lambda & 0 & 0 & 0\\ \lambda & & \lambda & \lambda +2\mu & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \mu & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & \mu & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & \mu \end{array}\right],\end{eqnarray} \tag{ 2 }$$

where $\lambda$ is the Lame coefficient of the background medium and $\mu$ is the shear modulus. As ${c}_{ij}^{0},$ ${c}_{ij}^{1},$ ${c}_{ij}^{2}$ denote elements of ${C}^{0},$ ${C}^{1},$ ${C}^{2},$ the matrix ${C}^{0}$ can be represented as

$$\begin{eqnarray}&&\left\{\begin{array}{c}{c}_{11}^{0}=\lambda +2\mu \\ {c}_{13}^{1}=\lambda \\ {c}_{33}^{1}=\lambda +2\mu \\ {c}_{44}^{1}=\mu \end{array}\right.,\end{eqnarray} \tag{ 3 }$$

and the first-order corrected matrix ${C}^{1}$ can be presented as

$$\begin{eqnarray}&&\left\{\begin{array}{c}{c}_{11}^{1}=-\displaystyle \frac{{\lambda }^{2}}{\mu }\varepsilon {U}_{3}\\ {c}_{13}^{1}=-\displaystyle \frac{\lambda (\lambda +2\mu )}{\mu }\varepsilon {U}_{3}\\ {c}_{33}^{1}=-\displaystyle \frac{{(\lambda +2\mu )}^{2}}{\mu }\varepsilon {U}_{3}\\ {c}_{44}^{1}=-\mu \varepsilon {U}_{1}\end{array}\right.,\end{eqnarray} \tag{ 4 }$$

where $\varepsilon$ is the crack density and ${U}_{1},$ ${U}_{3}$ are determined by the elastic parameters of the background medium. For dry cracked media, the forms of ${U}_{1},$ ${U}_{3}$ are given by the functions

$$\begin{eqnarray}&&\left\{\begin{array}{c}{U}_{1}=\displaystyle \frac{16(\lambda +2\mu )}{3(3\lambda +4\mu )}\\ {U}_{3}=\displaystyle \frac{4(\lambda +2\mu )}{3(\lambda +\mu )}\end{array}\right..\end{eqnarray} \tag{ 5 }$$

For saturated fluids, ${U}_{1},$ ${U}_{3}$ are given by the functions

$$\begin{eqnarray}&&\left\{\begin{array}{c}{U}_{1}=\displaystyle \frac{16(\lambda +2\mu )}{3(3\lambda +4\mu )}\displaystyle \frac{1}{1+M}\\ {U}_{3}=\displaystyle \frac{4(\lambda +2\mu )}{3(\lambda +\mu )}\displaystyle \frac{1}{\kappa +1}\end{array}\right.\end{eqnarray} \tag{ 6 }$$

where $\kappa$ and $M$ are parameters related to the nature of the crack inclusions and background, and they can be given as the following functions

$$\begin{eqnarray}&&\left\{\begin{array}{c}\kappa =\displaystyle \frac{\left[K^{\prime} +\displaystyle \frac{4}{3}\mu ^{\prime} \right](\lambda +2\mu )}{\pi \alpha \mu (\lambda +\mu )}\\ M=\displaystyle \frac{\mu ^{\prime} }{\pi \alpha \mu }\displaystyle \frac{\lambda +2\mu }{3\lambda +4\mu }\end{array}\right.,\end{eqnarray} \tag{ 7 }$$

where $K^{\prime}$ is the bulk modulus of the crack inclusions, $\mu ^{\prime}$ is the shear modulus of the crack inclusions, $\mu ^{\prime} =0$ as the crack inclusion is fluid and $\alpha$ is the average aspect ratio of the cracks. ${C}^{2}$ can be presented as the following expression

$$\begin{eqnarray}&&\left\{\begin{array}{c}{c}_{11}^{2}=\displaystyle \frac{q}{15}\displaystyle \frac{{\lambda }^{2}}{\lambda +2\mu }{(\varepsilon {U}_{3})}^{2}\\ {c}_{13}^{2}=\displaystyle \frac{q}{15}\lambda {(\varepsilon {U}_{3})}^{2}\\ {c}_{33}^{2}=\displaystyle \frac{q}{15}(\lambda +2\mu ){(\varepsilon {U}_{3})}^{2}\\ {c}_{44}^{2}=\displaystyle \frac{2}{15}\displaystyle \frac{\mu (3\lambda +8\mu )}{\lambda +2\mu }{(\varepsilon {U}_{1})}^{2}\end{array}\right.,\end{eqnarray} \tag{ 8 }$$

where $q$ was given by the function

$$\begin{eqnarray}&&q=15\displaystyle \frac{{\lambda }^{2}}{{\mu }^{2}}+28\displaystyle \frac{\lambda }{\mu }+28.\end{eqnarray} \tag{ 9 }$$

The crack density can be calculated by the following expression

$$\begin{eqnarray}&&\varepsilon =\displaystyle \frac{3\phi }{4\pi \alpha },\end{eqnarray} \tag{ 10 }$$

where $\phi$ is the crack porosity. The elastic parameters of the cracked medium ${c}_{ij}$ can be obtained as

$$\begin{eqnarray}&&{c}_{ij}={c}_{ij}^{0}+{c}_{ij}^{1}+{c}_{ij}^{2}.\end{eqnarray} \tag{ 11 }$$

Using the relationship between the acoustic velocity and the elastic parameters, we can give the P-wave velocity in the axial direction as ${v}_{{pz}}=\sqrt{\tfrac{{c}_{33}}{\rho }}\,,$ the P-wave velocity in radial direction as ${v}_{{pr}}=\sqrt{\tfrac{{c}_{11}}{\rho }}\,,$ the S-wave velocity as ${v}_{{sz}}=\sqrt{\tfrac{{c}_{55}}{\rho }}$ and the anisotropy coefficient of the P-wave as $\gamma =\tfrac{{c}_{11}-{c}_{33}}{2{c}_{33}}$.

In cylindrical coordinates, the three-dimensional wave equation can be transformed into two-dimensions for VTI media with cracks about the axial symmetry distribution. Ignoring external forces, the stress-displacement elastic wave equation can be presented as the following expression

$$\begin{eqnarray}&&\left\{\begin{array}{c}\displaystyle \frac{{\partial }^{2}u}{\partial {t}^{2}}=\displaystyle \frac{1}{\rho }\left(\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}(r{\tau }_{rr})+\displaystyle \frac{\partial {\tau }_{rz}}{\partial z}\right),\\ \displaystyle \frac{{\partial }^{2}w}{\partial {t}^{2}}=\displaystyle \frac{1}{\rho }\left(\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}(r{\tau }_{rz})+\displaystyle \frac{\partial {\tau }_{zz}}{\partial z}\right),\\ {\tau }_{rr}={c}_{11}\displaystyle \frac{\partial u}{\partial r}+{c}_{13}\left(\displaystyle \frac{u}{r}+\displaystyle \frac{\partial w}{\partial z}\right),\\ {\tau }_{zz}={c}_{33}\displaystyle \frac{\partial w}{\partial z}+{c}_{13}\left(\displaystyle \frac{u}{r}+\displaystyle \frac{\partial u}{\partial r}\right),\\ {\tau }_{rz}={c}_{44}\left(\displaystyle \frac{\partial u}{\partial z}+\displaystyle \frac{\partial w}{\partial r}\right),\end{array}\right.\end{eqnarray} \tag{ 12 }$$

where $u$ and $w$ respectively indicate displacement in the $r$ and $z$ directions, ${\tau }_{rr}$ and ${\tau }_{zz}$ are normal stress in the $r$ and $z$ directions, $\rho$ is the bulk density and ${c}_{11},$ ${c}_{13},$ ${c}_{33},$ ${c}_{44}$ are elastic constants of cracked media. If $\displaystyle \frac{\partial u}{\partial t}={v}_{r},$ $\displaystyle \frac{\partial w}{\partial t}={v}_{z},$ then we can transform equation (12) into the following expression

$$\begin{eqnarray}&&\left\{\begin{array}{c}\displaystyle \frac{\partial {v}_{r}}{\partial t}=\displaystyle \frac{1}{\rho }\left(\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}(r{\tau }_{rr})+\displaystyle \frac{\partial {\tau }_{rz}}{\partial z}\right),\\ \displaystyle \frac{\partial {v}_{z}}{\partial t}=\displaystyle \frac{1}{\rho }\left(\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}(r{\tau }_{rz})+\displaystyle \frac{\partial {\tau }_{zz}}{\partial z}\right),\\ \displaystyle \frac{\partial {\tau }_{rr}}{\partial t}={c}_{11}\displaystyle \frac{\partial {v}_{r}}{\partial r}+{c}_{13}\left(\displaystyle \frac{{v}_{r}}{r}+\displaystyle \frac{\partial {v}_{z}}{\partial z}\right),\\ \displaystyle \frac{\partial {\tau }_{zz}}{\partial t}={c}_{33}\displaystyle \frac{\partial {v}_{z}}{\partial z}+{c}_{13}\left(\displaystyle \frac{{v}_{r}}{r}+\displaystyle \frac{\partial {v}_{r}}{\partial r}\right),\\ \displaystyle \frac{\partial {\tau }_{rz}}{\partial t}={c}_{44}\left(\displaystyle \frac{\partial {v}_{r}}{\partial z}+\displaystyle \frac{\partial {v}_{z}}{\partial r}\right),\end{array}\right.\end{eqnarray} \tag{ 13 }$$

so, the function (13) is the speed-stress equation in cylindrical coordinates.

In cylindrical coordinates, we set grids for the computing space. ${\rm{\Delta }}r$ and ${\rm{\Delta }}z$ are computing steps in the radial and axial directions respectively. $M$ and $N$ are the numbers of meshes in the radial and axial directions respectively. So, $(M\times {\rm{\Delta }}r)\times (N\times {\rm{\Delta }}Z)$ is the calculation area. We transform the wave equation (13) into a high-order staggered grid difference format as the following functions

$$\begin{eqnarray}\left\{\begin{array}{rcl}{U}_{i+\displaystyle \frac{1}{2},j}^{k+\displaystyle \frac{1}{2}} & = & {U}_{i+\displaystyle \frac{1}{2},j}^{k-\displaystyle \frac{1}{2}}+\displaystyle \frac{{\rm{\Delta }}t}{\rho }\left[\displaystyle \frac{{R}_{i,j}^{k}}{r}+{D}_{r}^{-}{R}_{i,j}^{k}+{D}_{z}^{-}{H}_{i+\displaystyle \frac{1}{2},j+\displaystyle \frac{1}{2}}^{k}\right]\\ {V}_{i,j+\displaystyle \frac{1}{2}}^{k+\displaystyle \frac{1}{2}} & = & {V}_{i,j+\displaystyle \frac{1}{2}}^{k-\displaystyle \frac{1}{2}}+\displaystyle \frac{{\rm{\Delta }}t}{\rho }\left[\displaystyle \frac{{H}_{i+\displaystyle \frac{1}{2},j+\displaystyle \frac{1}{2}}^{k}}{r}+{D}_{r}^{+}{H}_{i,j}^{k}+{D}_{z}^{+}{T}_{i,j}^{k}\right]\\ {T}_{i,j}^{k+1} & = & {T}_{i,j}^{k}+{c}_{11}(i,j){\rm{\Delta }}t{D}_{r}^{+}{U}_{i+\displaystyle \frac{1}{2},j}^{k+\displaystyle \frac{1}{2}}\\ & & +{c}_{13}(i,j){\rm{\Delta }}t\left[\left(\displaystyle \frac{{U}_{i+\displaystyle \frac{1}{2},j}^{k+\displaystyle \frac{1}{2}}}{r}\right)+{D}_{z}^{-}{V}_{i,j+\displaystyle \frac{1}{2}}^{k+\displaystyle \frac{1}{2}}\right],\\ {R}_{i,j}^{k+1} & = & {R}_{i,j}^{k}+{c}_{33}(i,j){\rm{\Delta }}t{D}_{r}^{-}{U}_{i+\displaystyle \frac{1}{2},j}^{k+\displaystyle \frac{1}{2}}\\ & & +{c}_{13}(i,j){\rm{\Delta }}t\left[\left(\displaystyle \frac{{U}_{i+\displaystyle \frac{1}{2},j}^{k+\displaystyle \frac{1}{2}}}{r}\right)+{D}_{r}^{-}{U}_{i,j+\displaystyle \frac{1}{2}}^{k+\displaystyle \frac{1}{2}}\right]\\ {H}_{i+\displaystyle \frac{1}{2},j+\displaystyle \frac{1}{2}}^{k+1} & = & {H}_{i+\displaystyle \frac{1}{2},j+\displaystyle \frac{1}{2}}^{k}+{c}_{44}(i,j){\rm{\Delta }}t\\ & & \times \left[{D}_{z}^{+}{U}_{i+\displaystyle \frac{1}{2},j}^{k+\displaystyle \frac{1}{2}}+{D}_{r}^{-}{V}_{i,j+\displaystyle \frac{1}{2}}^{k+\displaystyle \frac{1}{2}}\right]\end{array}\right.\end{eqnarray} \tag{ 14 }$$

where ${U}_{i+\displaystyle \frac{1}{2},j}^{k+\displaystyle \frac{1}{2}}$ and ${V}_{i+\displaystyle \frac{1}{2},j}^{k+\displaystyle \frac{1}{2}}$ denote ${v}_{r}$ and ${v}_{z}$ respectively, ${R}_{i,j}^{k},$ ${T}_{i,j}^{k}$ and ${H}_{i+\displaystyle \frac{1}{2},j+\displaystyle \frac{1}{2}}^{k}$ denote ${\tau }_{rr},$ ${\tau }_{zz}$ and ${\tau }_{zr}$ respectively, ${D}_{r}^{-},$ ${D}_{r}^{+}$ are high order difference operators with backward and forward direction respectively, which are about the first order derivation in the radial direction. Similarly, ${D}_{z}^{-}$ and ${D}_{z}^{+}$ denote the operators in the axial direction. They can be presented as the following expression

$$\begin{eqnarray}\left\{\begin{array}{rcl}{D}_{r}^{-}f(i{\rm{\Delta }}r) & = & \displaystyle \frac{1}{{\rm{\Delta }}r}\displaystyle \sum _{n=1}^{n=L}{c}_{n}^{(L)}\{f[i{\rm{\Delta }}r+(n-1){\rm{\Delta }}r]\\ & & -f[(i-n){\rm{\Delta }}r]\}\\ {D}_{r}^{+}f(j{\rm{\Delta }}r) & = & \displaystyle \frac{1}{{\rm{\Delta }}r}\displaystyle \sum _{n=1}^{n=L}{c}_{n}^{(L)}\{f[(j+n){\rm{\Delta }}r]\\ & & -f[j{\rm{\Delta }}r-(n-1){\rm{\Delta }}r]\}\end{array}\right.,\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}\left\{\begin{array}{rcl}{D}_{z}^{-}f(i{\rm{\Delta }}z) & = & \displaystyle \frac{1}{{\rm{\Delta }}z}\displaystyle \sum _{n=1}^{n=L}{c}_{n}^{(L)}\{f[i{\rm{\Delta }}z+(n-1){\rm{\Delta }}z]\\ & & -f[(i-n){\rm{\Delta }}z]\}\\ {D}_{z}^{+}f(j{\rm{\Delta }}z) & = & \displaystyle \frac{1}{{\rm{\Delta }}z}\displaystyle \sum _{n=1}^{n=L}{c}_{n}^{(L)}\{f[(j+n){\rm{\Delta }}z]\\ & & -f[j{\rm{\Delta }}z-(n-1){\rm{\Delta }}z]\}\end{array}\right.,\end{eqnarray} \tag{ 16 }$$

where ${c}_{n}^{L}$ is the differential coefficient, $i=1,$ 2,..., $M$ is the number of discrete points in the radial direction, j = 1, 2,..., $N$ is the number of discrete points in the axial direction.

We simulated the acoustic field in a borehole in cylindrical coordinates based on the functions in (14). The borehole radius r = 0.12 m, step ${\rm{\Delta }}r={\rm{\Delta }}z=0.005\,{\rm{m}}$ and the point source was placed on the Z axis at the center of the borehole, as shown in figure 1.

Zoom In Zoom Out Reset image size

5.1. Acoustic source function

The first order derivative of Gaussian function was modeled as the sound source, and can be presented as the following expression

$$\begin{eqnarray}&&g(t)=-2\xi (t-{t}_{s})\exp [-\xi {(t-{t}_{s})}^{2}]\end{eqnarray} \tag{ 17 }$$

where $\xi$ is a parameter related to the bandwidth of the pulse, ${{t}}_{{s}}$ is a time translation parameter, ${{f}}_{0}$ is the center frequency of the sound source, ${{f}}_{{\rm{\max }}}$ is maximum frequency, and their values are presented in table 1.

5.2. Dispersion with the size of the grid and time step

If the mesh size is too large, the sound waves will cause dispersion. Conversely, if the mesh size is too small, the computational efficiency will be decreased. Alford et al (1974) gave an empirical formula to decide the grid size

$$\begin{eqnarray}&&{\rm{\Delta }}r\leqslant \displaystyle \frac{{{V}}_{{\rm{\min }}}}{{{\rm{1}}{\rm{0}}{f}}_{{\rm{\max }}}},\end{eqnarray} \tag{ 18 }$$

where ${\rm{\Delta }}r$ is the maximum size of the grid, ${{V}}_{{\rm{\min }}}$ is the minimum P-wave velocity and ${{f}}_{{\rm{\max }}}$ is the maximum frequency. The time step can be determined according to the following formula

$$\begin{eqnarray}&&{\rm{\Delta }}t\leqslant \displaystyle \frac{{\rm{\Delta }}r}{\sqrt{2}{{V}}_{{\rm{\max }}}}\end{eqnarray} \tag{ 19 }$$

where ${{V}}_{{\rm{\max }}}$ is the maximum speed of the wave. In this paper, we calculated using the grid spacing interval ${\rm{\Delta }}r={\rm{\Delta }}z\,=0.005\,{\rm{m}}$ and the time interval ${\rm{\Delta }}t=0.4\,\mu {\rm{s}},$ in line with formulas (18) and (19).

Table 1.  Parameters of the sound source function.

${{f}}_{0}$(kHz)	${{f}}_{{\rm{\max }}}$ (kHz)	$\xi$	${{t}}_{{s}}$ (ms)
10.6	20.3	0.22	0.056

5.3. Boundary conditions

We used an exponential function (Alford et al 1974) in an expanding attenuation band for the upper and lower boundaries, as well as the outer boundary in the radial direction. We considered adding an attenuation band out of the boundary, the physical quantity of which attenuates exponentially in this field. In addition, the axis of the borehole and several grids near the axis can be assumed as r < 0. Then, it can be calculated according to the symmetry of the sound field as the following conditions

$$\begin{eqnarray}&&\left\{\begin{array}{c}U(-r)=-U(r)\\ V(-r)=V(r)\\ T(-r)=T(r)\\ R(-r)=R(r)\\ H(-r)=H(r)\end{array}\right..\end{eqnarray} \tag{ 20 }$$

The elastic properties influenced by cracks are connected mainly to the shape and density of the cracks. This paper pays attention to these two factors to study the acoustic wave propagation characteristics in a borehole surrounded by vertical transversely isotropic cracked media.

6.1. Influence of crack density

According to Hudson's approach, the elastic modulus of dry rock is connected to crack density, so we simulated the borehole acoustic field with different cracked densities in order to study the influence of crack density on dry rock. It is numerically simulated using a finite-difference technique with second order in time and tenth order in space. The petrophysical parameters of the background medium are shown in table 2.

Table 2.  Petrophysical parameters of background media.

Parameters	${c}_{11}^{0}$(GPa)	${c}_{13}^{0}$(GPa)	${c}_{33}^{0}$(GPa)	${c}_{44}^{0}$(GPa)	$\rho$(kg m−3)	$\phi$
Value	32.218	1.5175	32.218	8.5212	2630	0.01

As shown in figure 2, the curves illustrate the relationship between the velocity of the P-wave and S-wave, as well as the anisotropy factor and crack density. In figure 2, the abscissa is the crack density and the ordinate is the velocity of the wave. In the figure, ${v}_{{pr}}$ is the velocity of the P-wave in the radial direction, and ${v}_{{pz}}$ is the velocity of the P-wave in the axial direction. As can be seen from figure 2, the velocity of the P-wave in both axial and radial directions decrease as the crack density varies from 0 to 0.1, but more rapidly in the axial direction while the crack density increases. The velocity of the S-wave decreases slowly with the increase in crack density.

Zoom In Zoom Out Reset image size

As shown in figure 3, the curve illustrates the relationship between the anisotropy factor of the P-wave and crack density. The abscissa is the crack density and the ordinate is the anisotropy coefficient of the P-wave. As can be seen from figure 3, the anisotropy of the medium increases linearly while the crack density increases.

Zoom In Zoom Out Reset image size

As shown in figure 4, panels (a), (b) and (c) are snapshots of the borehole acoustic wave field when the cracked density $\varepsilon =0,$ 0.03, 0.05 respectively at time t = 280 μs (take the grid units as coordinate units). From the three panels we can see there are spherical waves reflected by the wall of the borehole. Furthermore, there are two arcs linking the slashes which are symmetrically distributed in the wellbore outside. The outer arc is the P-wave of formation and the inside arc is the S-wave. Comparing these three figures, the velocities of the P-wave and the S-wave decrease with decreasing crack density, but more rapidly in the axial direction. The S-wave velocity decreases slowly with increasing crack density. So, the velocity of the P-wave is more sensitive than the S-wave.

Zoom In Zoom Out Reset image size

As shown in figure 5, panels (a), (b) and (c) are waveforms received at the axis from a distance between the sound source 1 m, when the cracked density $\varepsilon =0,$ 0.03, 0.05 respectively. From the three panels we can see there are P-waves, S-waves, pseudo-Rayleigh waves and Stoneley-waveforms in the time axis. Comparing these three figures, the velocities of the P-wave and S-wave decrease with increasing crack density, while the wave amplitude of the P-wave, S-wave, pseudo-Rayleigh waves and Stoneley-wave decrease with increasing crack density.

Zoom In Zoom Out Reset image size

6.2. Influence of the aspect ratio of cracks

According to Hudson's approach, the elastic modulus of fluid-saturated rock is connected with the aspect ratio of the cracks. So we simulated the borehole acoustic field with cracks of different aspect ratios in order to study the influence of the aspect ratio of cracks on fluid-saturated rock. It is numerically simulated using the finite-difference technique with second order in time and tenth order in space. The petrophysical parameters of the background medium are shown in table 2.

As shown in figure 6, the curves illustrate the relationship between the anisotropy factor of the P-wave, S-wave and aspect ratio of the cracks. From figure 6, it can be seen that the velocity of the P-wave in the radial direction is greater than in the axial direction, the velocity in both the radial and axial directions linearly increases with increasing aspect ratio in the range of 0–0.1. Furthermore, from the curve of the relationship with the S-wave and the aspect ratio, it can be seen that the velocity of the S-wave increases with increasing aspect ratio, more quickly in the range of 0–0.02 and after that almost linearly increasing.

Zoom In Zoom Out Reset image size

Figure 7 shows the relationship curve between the anisotropy coefficient and the aspect ratio of the cracks. In figure 7, the abscissa is the aspect ratio and the ordinate is the anisotropy coefficient of the P-wave. It can be seen that the anisotropy coefficient decreases with increasing aspect ratio.

Zoom In Zoom Out Reset image size

As shown in figure 8, panels (a), (b) and (c) are snapshots of the borehole acoustic wave field where the cracked density $\alpha =0,$ 0.03, 0.05 respectively at time t = 280 μs (take the grid units as coordinate units). Comparing these three panels, the distance of the P-wave and S-wave wave field spread increases with increasing aspect ratio. From the shape of the P-wave field in the snapshots, with an increase in the aspect ratio it varys from oblate to round. Obviously, the anisotropy factor of the P-wave decreases with the aspect ratio.

Zoom In Zoom Out Reset image size

As shown in figure 9, panels (a), (b) and (c) are the are waveforms received at the axis from a distance between the sound source 1 m and aspect ratio α = 0, 0.03, 0.05 respectively. Comparing these three panels, the velocity of the P-wave and S-wave increases with increasing aspect ratio, while the wave amplitude of the P-wave, S-wave, pseudo-Rayleigh waves and Stoneley-wave does not obviously change obviously with increasing aspect ratio.

Zoom In Zoom Out Reset image size

This paper simulated a borehole acoustic field using the finite difference method based on Hudson's approach. It studied the propagation characteristics of a borehole acoustic field surrounded by vertical transversely isotropic cracked media and also the effects of crack density and crack aspect ratio on the propagation of acoustic waves.

Numerical simulation shows that, in dry rock, the velocity of the P-wave and S-wave decreases with increasing crack density, while it decreases more rapidly for the S-wave. The attenuation of the wave amplitude increases with increasing crack density. In the saturated fluid medium, the velocity of the P-wave and S-wave increases with the aspect ratio of cracks. The anisotropy factor of the P-wave decreases with the aspect ratio. The influence of the aspect ratio of the cracks on attenuation is not obvious.

It should be noted that the simulation satisfied the presupposition of effective-medium theories in Hudson's approach.

We are grateful to the School of Geology Engineering and Geomatics, Chang'an University, where this study was carried out in the framework of the research program 'Study on the characteristics of acoustic wave in a borehole surrounded by unconsolidated sandstone' subsidized by the Natural Science Foundation of China.