Novel finite-element approach applied to borehole quadrupole dispersion analysis in stress-sensitive formationsFE analysis of quadrupole dispersion | Geophysics

ABSTRACT

Near a borehole, stress concentration effects may cause a complex spatial variation of elastic anisotropy. Stress-induced sonic anisotropy results when moduli and velocities are stress dependent and the state of stress is nonhydrostatic. For such cases, the appropriate material model is that of an elastic orthorhombic medium with material axes aligned with the principal axes of stress. We simulate the dispersion characteristics of quadrupole waves propagating along a borehole in a stress-sensitive formation. To include sensitivity to stress in the analysis, we derive a finite-element modeling approach which includes the stress-velocity relationship, the perturbed stress and velocity field near the borehole, and solutions to the full coupled anisotropic and 3D wave equation. We characterized the anisotropy around the borehole for different in situ stress conditions, and we derived the effects on quadrupole dispersion characteristics. Of particular interest was the phase velocity of the quadrupole wave, which at low frequencies provides an estimate of formation shear-wave velocity. Results indicated that the quadrupole mode dispersion was relatively insensitive to the hoop stress effects around a borehole giving greater confidence to shear-wave estimates made from these modes. The proposed ring-element formulation, which includes Fourier expansions of the field variables in the direction of the circumference, is well suited to this analysis and can be adapted to include other effects, such as a finite-length borehole.

acoustic, anisotropy, borehole geophysics, dispersion, finite element

INTRODUCTION

Quadrupole modes are guided waves that propagate along a borehole with the azimuthal acoustic profile $\cos (2 θ)$ in the propagating medium. In isotropic and homogeneous formations, the velocity of quadrupole waves moves asymptotically to the Scholte wave speed at high frequencies and to the formation shear speed at low frequencies. There is practically no energy directed into this wave at frequencies below the cut-off frequency velocity, where the notation of a cut-off frequency is used to refer to the frequency at which the phase velocity of the mode is equal to the formation shear velocity. Due to the asymptotic behavior, quadrupole acoustic logging can be used for estimating in situ shear-wave velocities in real formations. The use of this mode for shear-wave velocity estimation has been investigated for many years (e.g., White, 1967; Chen, 1989; Ellefsen, 1990; Wang and Tang, 2003; Zhang et al., 2010).

Because borehole acoustic waves only penetrate a few wavelengths or borehole radii into the formation, it may be suspected that near-borehole effects could dominate sonic measurements significantly. In rocks where sonic speed is coupled to stress, this phenomenon is sometimes referred to as acoustoelasticity, variations in stress concentration around the borehole (i.e., hoop stresses) affect sonic velocities near the borehole. This effect has been clearly shown in cross-dipole logging (e.g., Sinha and Kostek, 1996; Winkler et al., 1998; Huang et al., 1999; Sinha et al., 2000). The potential impact of this stress effect on quadrupole sonic logging in a wellbore is the subject of this study. Stress concentration effects cannot be removed from the reservoir and measurement environment. From a modeling point of view, we may include them or leave them out of our analyses. Hence, modeling can quantify their effect, assess their importance, and revise log interpretations where needed.

Acoustoelastic coupling is considered to be due to stiffening of intergranular contacts of the hertzian type, or closure of microcracks or cracklike components of the pore space. For an acoustoelastic material, a nonhydrostatic state of stress causes orthorhombic symmetry; i.e., the elastic properties are those of an orthorhombic medium, with material axes aligned with the principal axes of stress. The constitutive model that we use to capture this stress-moduli interdependency is originally proposed by Mavko et al. (1995).

In our numerical study, we consider a formation with properties similar to Berea sandstone, for which experiments show that velocities increase with increasing confining stress (Sayers et al., 1990; Fang et al., 2013). We consider a circular vertical borehole, for which we characterize the perturbed stress and velocity field and simulate the quadrupole wave under heterogeneous stress and velocity conditions. We establish the dispersion characteristics and demonstrate how effects related to the near-borehole stress field can be accounted for and the true formation shear velocity can be estimated.

FORMULATION

Propagation of acoustic waves in a borehole has been subject to extensive analytical works. A monograph around the topic has also been published (Tang and Cheng, 2004). Analytical methods work well in an unbounded linear elastic and isotropic homogeneous medium, but they have limitations when it comes to real formations. Real rocks are in many cases acoustoelastic, e.g., their velocities and stiffnesses increase with increasing confining stress (e.g., Nur and Simmons, 1969); hence, rocks are anisotropic and the formation is heterogeneous in the immediate proximity of the borehole. However, for reasons of modeling difficulty, the related impact of this effect on borehole acoustics is not fully implemented in most modeling methods.

We therefore derive a simulation model applicable for studying borehole acoustics in stress-sensitive formations that includes the elastic stress concentration around the borehole, and we assume orthorhombic material properties derived from a stress-dependent constitutive model (Mavko et al., 1995). We use a special purpose finite-element approach (Jørgensen, 1993), which has proven particularly useful in modeling stresses and strains in anisotropic elastic bodies of cylindrical geometry. The modeling approach was originally developed for materials science research, in which the study of composite materials has demanded numerical techniques for this particular class of problems (Jørgensen, 1994; Jørgensen et al., 1998).

Constitutive model for stress-dependent elasticity

The constitutive model (Mavko et al., 1995) regards a class of rocks in which elastic closure of a random distribution of compliant cracks or cracklike voids is stiffening the material under increasing confining stress. The derived elasticity (compliance) tensor $S_{i j k l}$ is isotropic (two independent parameters) under hydrostatic confining stress but orthorhombic (nine independent parameters), with material axes parallel to the directions of principal stresses, under a general triaxial loading.

Under the assumption that closure of microcracks, or cracklike components of the pore space, causes the elastic stiffening under increasing confining stress, the full stress-dependent compliance tensor $S_{i j k l} (σ)$ is established. Here, $σ$ refers to the full stress tensor and the compliance tensor is derived for an arbitrary state of stress.

Under very high confining stress, all cracks are closed and crack compliance no longer contributes to the effective medium compliance. This asymptotic and isotropic behavior plays a central role in the formulation. The compliance asymptote of the rock under very high hydrostatic confining stress

σ_{hyd}

⁠, i.e.,

S_{i j k l}^{0, ISO} = S_{i j k l}^{ISO} (σ_{hyd} \to \infty)

⁠, is measured in laboratory experiments. For lower levels of hydrostatic stress, additional compliance due to crack openings exists and the total compliance is the sum of the asymptotic compliance

S_{i j k l}^{0, ISO}

and the compliance due to cracks

Δ S_{i j k l}^{ISO}

⁠; thus,

S_{i j k l}^{ISO} (σ_{hyd}) = Δ S_{i j k l}^{ISO} (σ_{hyd}) + S_{i j k l}^{0, ISO} .

(1)

The left side of equation can be established from a series of hydrostatic measurements and the unknown compliance under hydrostatic loading due to cracks can be isolated:

Δ S_{i j k l}^{ISO} (σ_{hyd}) = S_{i j k l}^{ISO} (σ_{hyd}) - S_{i j k l}^{0, ISO} .

(2)

The derivations thus far serve to quantify the nonlinear material behavior under strictly isotropic conditions. To reach beyond the isotropic conditions, Mavko et al. (1995) apply two physical assumptions: (1) normal and shear deformations of the crack decouple and (2) it is the crack normal component of stress alone that determines the crack closure. Under these assumptions, the normal and tangential crack compliance tensors

W^{N}

and

W^{T}

are (Mavko et al., 1995)

W^{N} (σ_{hyd}) = \frac{1}{2 π} Δ S_{j j k k}^{ISO} (σ_{hyd}),

(3)

and

W^{T} (σ_{hyd}) = \frac{Δ S_{j k j k}^{ISO} (σ_{hyd}) - Δ S_{j j k k}^{ISO} (σ_{hyd})}{4 Δ S_{j j k k}^{ISO} (σ_{hyd})} .

(4)

Finally, using the transformation properties of tensors, the stress-dependent compliance tensor under arbitrary stress

σ

reads

Δ S_{i j k l} (σ) = \int_{θ = 0}^{\frac{π}{2}} \int_{ϕ = 0}^{2 π} W^{N} (σ_{i j} m_{i} m_{j}) m_{i} m_{j} m_{k} m_{l} \sin (θ) d θ d ϕ + \int_{θ = 0}^{\frac{π}{2}} \int_{ϕ = 0}^{2 π} W^{T} (σ_{i j} m_{i} m_{j}) [δ_{i k} m_{j} m_{l} + δ_{i l} m_{j} m_{k} + δ_{j k} m_{i} m_{l} + δ_{j l} m_{i} m_{k} - 4 m_{i} m_{j} m_{k} m_{l}] \sin (θ) d θ d ϕ .

(5)

In equation ,

m_{i}

are the coordinates of the unit normal vector to the crack surface. In the finite-element formulation, we apply the stiffness tensor, which is found by inverting the total compliance; i.e.,

C_{i j k l} (σ) = {[Δ S_{i j k l} (σ) + S_{i j k l}^{0, ISO}]}^{- 1} .

(6)

Note that the latter inversion is very simple provided the compliance and stiffness matrices are evaluated in the material coordinate system.

3D modeling of stresses and strains

Anisotropy and spatial variations in anisotropy make the elastostatic problem challenging. Some observations of symmetry are, however, very helpful. First, because we consider reservoirs under considerable burial, we assume that the vertical direction (⁠ $z$ -direction) is one of the three orthogonal principal stress orientations. Second, we choose to align the global (x-, y-, z-) coordinates with the maximum and minimum horizontal stresses in the reservoir in all of our analyses, thus the global coordinate system (⁠ $x, y, z$ ⁠) has axes parallel with the in situ stresses; hence, $σ_{x x}, σ_{y y}, σ_{z z}$ are principal stresses.

We consider a vertical hole in the reservoir subject to in situ stresses (⁠ $σ_{x x}, σ_{y y}, σ_{z z}$ ⁠). We seek the stress distribution in the neighborhood of the hole, where a considerable stress concentration effect is present. The principal stresses and orthorhombic symmetry planes will have rotational symmetry about the vertical $z$ -axis only as long as the two horizontal in situ principal stresses are equal. For all other cases, the elastic anisotropy contradicts the basis of axial symmetry. Still, however, the vertical planes of mirror symmetry, indicated in the horizontal (⁠ $r, θ$ -plane in Figure 1) are noted. The two planes of symmetry go through the borehole center and are parallel with each of the horizontal in situ principal stresses.

Given those symmetries, we observe that the displacement field

d = (u, v, w)

⁠, with reference to polar coordinates (⁠

r, θ, z

⁠), can be arranged such that radial displacements

u

and axial displacements

w

are even functions of

2 θ

⁠, and circumferential displacements

v

are odd functions of

2 θ

⁠. So, without loss of generality we can write the displacement field as Fourier series (here truncated at

N

⁠):

u (r, θ, z) = \sum_{n = 0}^{N} u_{n} (r, θ, z) \cos (2 n θ), v (r, θ, z) = \sum_{n = 1}^{N} v_{n} (r, θ, z) \sin (2 n θ), w (r, θ, z) = \sum_{n = 0}^{N} w_{n} (r, θ, z) \cos (2 n θ) .

(7)

Very fast convergence of series 7 was demonstrated by Jørgensen (1993), and this makes the series expansions very useful for analysis concerning orthorhombic media and for our purpose in particular.

For the discretization we chose an eight-node isoparametric ring element, restricted to cylindrical domains (Jørgensen, 1993). The field variables, which are discretized by the finite elements, are the Fourier coefficient functions in equation , i.e., the generalized coordinates that uniquely specify the displacement field within element

e

⁠. We write the vector of generalized nodal coordinates in its short form as

d_{N_{e}} = {u_{0}^{i}, w_{0}^{i}, u_{1}^{i}, v_{1}^{i}, w_{1}^{i}, u_{2}^{i}, v_{2}^{i}, w_{2}^{i}, \dots ., u_{N}^{i}, v_{N}^{i}, w_{N}^{i}}_{e}^{T},

(8)

where

i = 1, 8

. The dimension of the nodal vector is

8 (2 + 3 N)

⁠. The element stiffness matrix

K_{s_{e}}

is established in the usual way by the volume integral:

K_{s_{e}} = \int_{V_{e}} = 2 \int_{0}^{π} \int_{A_{e}} B^{T} C B r d A d θ .

(9)

In equation ,

B

is the strain/displacement matrix;

C

is the constitutive matrix, i.e.,

C_{i j}

are the moduli using standard Voigt notation (see Appendix A);

V_{e}

is the volume of element

e

⁠; and

A_{e}

is the area of the vertical reference plane intersected by element

e

⁠. The dimension of the element stiffness matrix is

{[8 (2 + 3 N)]}^{2}

⁠. Note that the transformation properties of the orthorhombic material introduces factors of

\cos (2 (θ - θ_{1}))

⁠,

\cos (4 (θ - θ_{1}))

⁠,

\sin (2 (θ - θ_{1}))

⁠, and

\sin (4 (θ - θ_{1}))

C

⁠, where

(θ - θ_{1})

is the angle relative to the material axes and

θ

is the angle relative to the reference coordinate system (see Appendix A). Given the displacements in equation and the generalized coordinates in equation ,

B

contains all azimuthal modes

\sin (2 n θ)

and

\cos (2 n θ)

for

n = 1, N

⁠. One practical implication is that integration of the element stiffness matrix (equation ) involves integration of triple products of trigonometric functions. These integrals can be evaluated analytically following the simple scheme (where it is understood that

n, m

⁠, and

k

are positive even integers):

\int_{0}^{2 π} \cos (n θ) \cos (m θ) \cos (k θ) d θ = {\begin{matrix} \frac{π}{2} if | n \pm m | = k \\ 0 if | n \pm m | \neq k \end{matrix}}, \int_{0}^{2 π} \sin (n θ) \sin (m θ) \cos (k θ) d θ = {\begin{matrix} \mp \frac{π}{2} if | n \pm m | = k \\ 0 if | n \pm m | \neq k \end{matrix}}, \int_{0}^{2 π} \sin (n θ) \sin (m θ) \sin (k θ) d θ = 0, \int_{0}^{2 π} \cos (n θ) \cos (m θ) \sin (k θ) d θ = 0 .

(10)

Another implication is the coupling between modes, which inherently is related to anisotropy. To further assess the coupling, we consider the two element displacement vectors:

d_{n_{e}} = {0, \dots, 0, u_{n}^{i}, v_{n}^{i}, w_{n}^{i}, 0, . ., 0}_{e}^{T}, d_{m_{e}} = {0, \dots, 0, u_{m}^{i}, v_{m}^{i}, w_{m}^{i}, 0, . ., 0}_{e}^{T} .

(11)

In general,

d_{m_{e}}

and

d_{n_{e}}

may couple for

m \neq n

⁠; however, the coupling between modes is constrained by the following relation:

if | n - m | \geq 3 \Rightarrow d_{n_{e}}^{T} K_{s_{e}} d_{m_{e}} = 0 .

(12)

This

K_{s_{e}}

-orthogonality causes

K_{s_{e}}

to be strongly banded, and this applies also to the stiffness matrix of the full system. As a result, the full 3D field (equation ) can be solved with a computational effort similar to a 2D problem (Jørgensen, 1993). This statement applies to the numerical integration of the stiffness matrix as well, because the outer integral in equation can be evaluated analytically (see Appendix A; Jørgensen, 1993).

Another factor adding complexity is the stress sensitivity of the material causing an azimuthal variation in anisotropy; i.e., the principal axes of anisotropy vary along the circumference (see Figure 1 and the section “Stress distribution around a vertical borehole”). In Figure 1, four hatched angular segments are shown, each with the local material axes indicated with fine parallel lines. Within each segment, the stiffness tensor is evaluated following the method of Mavko et al. (1995), equations in this paper. Note that the local variation in the elasticity tensor of course complies with the overall mirror symmetry of the stress field. The outer integral of equation must take this anisotropy variation into account. This is overcome by dividing the outer integral into the sum

\int_{0}^{π} = (\int_{0}^{α_{1}} + \int_{π - α_{1}}^{π}) + (\int_{α_{1}}^{α_{2}} + \int_{π - α_{2}}^{π - α_{1}}) + \dots .

(13)

Each of the terms in parentheses has an analytical solution, as shown in Jørgensen (1994), but can also be evaluated using numerical integration. Importantly, the stiffness tensor is assumed constant within the subdivided integration bounds, but may vary between angular segments. In this way, the angular variation of the stiffness tensor within one element is included in the analysis.

The benefit of the applied method is the reduction in computational effort compared to conventional 3D finite-element analyses, which is the result of a reduced number of degrees of freedom, and a reduced bandwidth of the system stiffness matrix relative to conventional 3D finite-element modeling. This has to do with the fast convergence of the series in equation 7, which in turn has to do with the relatively weak coupling of modes (equation 12).

Sonic wave propagation along a borehole

All borehole-guided waves are eigenmodes and their shapes and associated frequencies can be derived as eigensolutions to a linear undamped eigenvalue problem. The eigenvalue problem is the wave equation formulated in the frequency domain excluding external forces. We now derive the eigenvalue problem and finally demonstrate a numerical formulation benchmark against known analytical solutions.

In the present finite-element formulation, we have chosen displacements

d = (u, v, w)

as fundamental degrees of freedom in the solid phase and we will choose acoustic pressure

p

in the fluid phase. We apply the same Fourier expansion for the pressure as for the radial displacements (see equation ); hence,

p (r, θ, z) = \sum_{n = 0}^{N} p_{n} (r, θ, z) \cos (2 n θ) .

(14)

We use the same eight-node finite element in the fluid phase, with the same shape functions

N^{i}

and the fluid element nodal pressure degrees of freedom therefore are

p_{N_{e}} = {p_{0}^{i}, p_{1}^{i}, p_{2}^{i}, \dots, p_{N}^{i}}_{e}^{T},

(15)

where

i = 1, 8

⁠. The vector defined by equation has the dimension

8 (1 + N)

⁠, whereas the element displacement vector (equation ) has the dimension

8 (2 + 3 N)

⁠. The dimension of the system displacement vector

d_{N}

is the total number of solid nodes times

(2 + 3 N)

⁠, whereas the dimension of the system nodal pressure vector

p_{N}

is the total number of pressure nodes times

(1 + N)

⁠.

N

still refers to the truncation number of the Fourier series (see equations and ).

The coupling between the fluid and solid means that the eigenvalue problem is a coupled one and takes the schematic form:

[[\begin{matrix} K_{s} & - Q \\ 0 & K_{f} \end{matrix}] - ω^{2} [\begin{matrix} M_{s} & 0 \\ Q^{T} & M_{f} \end{matrix}]] {\begin{matrix} d_{N} \\ p_{N} \end{matrix}} = {\begin{matrix} 0 \\ 0 \end{matrix}} .

(16)

In equation 16, $K_{s}$ is the system solid stiffness matrix, $K_{f}$ is the system fluid stiffness matrix, $M_{s}$ is the system solid mass matrix, and $M_{f}$ is the system fluid mass matrix, $Q$ is the system coefficient of coupling, and $ω$ is the angular frequency. The total number of degrees of freedom is $(1 + N) L_{f} + (2 + 3 N) L_{s}$ ⁠, where $L_{s}$ and $L_{f}$ are the number of solid and fluid element nodes, respectively, and $N$ is the Fourier series truncation number.

We have already discussed the derivation of the element solid stiffness matrix (Jørgensen, 1993). We now continue and define, at the element level, the solid mass matrix

M_{s}

the fluid-to-solid coupling matrix,

Q

⁠, and finally the fluid stiffness and mass matrices

K_{f}

and

M_{f}

⁠. We apply the principle of virtual work, which states that the internal virtual work (IVW) is equal to the external virtual work (EVW) when equilibrated forces undergo virtual displacements. The IVW done on a solid element is the work done by the elastic forces going through the virtual displacements

δ d_{N_{e}}^{T}

IVW = δ d_{N_{e}}^{T} K_{s_{e}} d_{N_{e}} .

(17)

The EVW is the work done by external forces going through the virtual displacements

δ d_{N_{e}}^{T}

⁠. In the present case, external forces include inertia forces and forces acting on the surface of the element. Let us first consider the EVW done by the inertia forces

{EVW}_{inertia}

⁠. The solid element mass matrix

M_{s}

follows from expressing the virtual work done by the inertia forces undergoing virtual displacement

δ d_{N_{e}}^{T}

{EVW}_{inertia} = - ω^{2} δ d_{N_{e}}^{T} M_{s_{e}} d_{N_{e}} .

(18)

We arrive at the expression in equation by expanding the displacement field variables

(u, v, w)

in terms of the nodal degrees of freedom (equation ) and the shape functions

N^{i}

⁠. Then, for each of the radial, circumferential, and axial modes we define the submatrices:

M_{s_{u_{n}}} = \int_{V_{e}} ρ_{s} N^{i} N^{j} \cos^{2} (2 n θ) d V, M_{s_{v_{n}}} = \int_{V_{e}} ρ_{s} N^{i} N^{j} \sin^{2} (2 n θ) d V, M_{s_{w_{n}}} = \int_{V_{e}} ρ_{s} N^{i} N^{j} \cos^{2} (2 n θ) d V .

(19)

In equation ,

ρ_{s}

is the specific mass of the solid. The solid element mass matrix corresponding to the form given in equation reads

M_{s_{e}} = [\begin{matrix} M_{s_{u_{0}}} & 0 & \dots & \dots & \dots & 0 \\ 0 & M_{s_{w_{0}}} & 0 & \dots & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & \dots & \dots & 0 & M_{s_{v_{N}}} & 0 \\ 0 & \dots & \dots & \dots & 0 & M_{s_{w_{N}}} \end{matrix}] .

(20)

The dimension of the solid element mass matrix is

{[8 (2 + 3 N)]}^{2}

⁠. For solid elements with one of the element sides coinciding with the fluid-to-solid interface, the fluid-to-solid coupling needs to be derived. Again, we express the virtual work done by the acoustic pressure due to virtual surface displacements

δ d_{N_{e}}^{T}

in the form

{EVW}_{fluid} = δ d_{N_{e}}^{T} Q_{e} p_{N_{e}} .

(21)

The acoustic pressure acts on the element side(s) that are coinciding with the interface between the fluid and solid phases

Γ_{I}

in the direction normal to the interface. Hence, the acoustic pressure performs virtual work on the solid element when the fluid-to-solid interface undergoes virtual displacements in the direction normal to the interface. Let the unit normal vector to the fluid surface be

n_{f}

⁠. The fluid-solid surface displacements follow from the expansion

N_{Γ_{I}, d}^{i} δ d_{N_{e}}

where

N_{Γ_{I}, d}^{i}

are the displacement shape functions restricted to the fluid-solid interface

Γ_{I}

⁠. Note that if a solid element has one side coinciding with the fluid-solid interface, then only the three shape functions belonging to the three element nodes on the interface are nonzero along

Γ_{I}

⁠. Similarly, for the adjacent fluid element, only the three shape functions belonging to the three element nodes on the interface are nonzero along

Γ_{I}

⁠, and similarly we expand the pressure degree of freedom as

N_{Γ_{I}, p}^{i} δ p_{N_{e}}

⁠:

{EVW}_{fluid} = δ d_{N_{e}}^{T} \int_{Γ_{I}} N_{Γ_{I}, d}^{i} n_{f} N_{Γ_{I}, p}^{j} d Γ p_{N_{e}} .

(22)

The coupling coefficient matrix

Q_{e}

follows from combining equations and :

Q_{e} = \int_{Γ_{I}} N_{Γ_{I}, d}^{i} n_{f} N_{Γ_{I}, p}^{j} d Γ .

(23)

The dimension of the matrix coupling one solid element to one fluid element through three shared nodes is

{[3 (1 + N)]}^{2}

⁠.

The fluid stiffness and mass matrices are derived by writing the wave equation in its weak form and applying the finite-element discretization to it. The weak form is acquired by multiplying the acoustic wave equation with a virtual pressure variation followed by integration over the domain of one element:

0 = \int_{V_{e}} δ p \frac{1}{ρ_{f}} (\frac{1}{v_{f}^{2}} \frac{\partial^{2} p}{\partial t^{2}} - \nabla^{2} p) d V .

(24)

Integration by parts leads to

0 = \int_{V_{e}} \frac{1}{ρ_{f}} {(\nabla δ p)}^{T} \nabla p d V + \int_{V_{e}} \frac{1}{ρ_{f}} (δ p) \ddot{p} d V - \int_{Γ_{e}} \frac{1}{ρ_{f}} δ p \nabla p \cdot n_{f} d Γ .

(25)

In equation ,

Γ_{e}

is the contour of element

e

⁠. The first integral is associated with the IVW due to a variation in pressure, and we seek to write it in the form

\int_{V_{f}} \frac{1}{ρ_{f}} {(\nabla δ p)}^{T} \nabla p d V = δ p_{N_{e}}^{T} K_{f_{e}} p_{N_{e}} .

(26)

Associated with the azimuthal modes

p_{n}

⁠, we write

K_{f_{p_{n}}} = \int_{V_{e}} \frac{1}{ρ_{f}} [(\frac{\partial N^{i}}{\partial r} \frac{\partial N^{j}}{\partial r} + \frac{\partial N^{i}}{\partial z} \frac{\partial N^{j}}{\partial z}) \cos^{2} (2 n θ) + \frac{N^{i} N^{j}}{r^{2}} 4 n^{2} \sin^{2} (2 n θ)] d V .

(27)

The fluid element stiffness matrix reads

K_{f_{e}} = [\begin{matrix} K_{f_{p_{0}}} & 0 & \dots & \dots & \dots & 0 \\ 0 & K_{f_{p_{1}}} & 0 & \dots & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & \dots & \dots & 0 & K_{f_{p_{N - 1}}} & 0 \\ 0 & \dots & \dots & \dots & 0 & K_{f_{p_{N}}} \end{matrix}] .

(28)

The second integral in equation is associated with the virtual work done by the inertia forces due to a variation in pressure. Again, we seek to write it in the form

\int_{V_{e}} \frac{1}{ρ_{f}} (δ p) \ddot{p} d V = - ω^{2} δ p_{N_{e}}^{T} M_{f_{e}} p_{N_{e}} .

(29)

To this end, we define the particular fluid mass matrix corresponding to azimuthal mode

p_{n}

⁠:

M_{f_{p_{n}}} = \int_{V_{e}} \frac{1}{ρ_{f} {v_{f}}^{2}} N^{i} N^{j} \cos^{2} (2 n θ) d V .

(30)

Hereafter, the fluid element mass matrix reads

K_{f_{e}} = [\begin{matrix} M_{f_{p_{0}}} & 0 & \dots & \dots & \dots & 0 \\ 0 & M_{f_{p_{1}}} & 0 & \dots & \dots & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots \\ 0 & \dots & \dots & 0 & M_{f_{p_{N - 1}}} & 0 \\ 0 & \dots & \dots & \dots & 0 & M_{f_{p_{N}}} \end{matrix}] .

(31)

The dimension of the fluid element mass and stiffness matrices is

{[8 (1 + N)]}^{2}

⁠. In equation , the coupling coefficient matrix [

Q

] expresses the boundary conditions at the solid/fluid interface, i.e., that pressure and normal velocities are continuous. In the derivation of the coupling coefficient matrix given by equation , we applied the first condition. Now, for the second displacement boundary condition, we write

\dot{d_{f}} \cdot n_{f} = \dot{d_{s}} \cdot n_{f},

(32)

where

n_{f}

is the unit normal vectors to the fluid-solid interfaces and

\dot{()}

indicates the time derivative. We rewrite the latter equation by making use of the equation of motion

ρ_{f} \ddot{d_{f}} = - \nabla p .

(33)

Taking the time derivative of equation and substituting equation into it gives

\frac{1}{ρ_{f}} \nabla p \cdot n_{f} = - {\ddot{d}}_{s} \cdot n_{f} .

(34)

We now focus on the part of element boundary

Γ_{e}

⁠, which coincides with the fluid-solid interface

Γ_{I}

⁠. Considering this part of the third integral in equation , we write

- \int_{Γ_{I}} \frac{1}{ρ_{f}} δ p \nabla p \cdot n_{f} d Γ = \int_{Γ_{I}} δ p {\ddot{d}}_{s} \cdot n_{f} d Γ = - ω^{2} δ p_{N_{e}}^{T} \int_{Γ_{I}} N_{Γ_{I}, p}^{j} n_{f} N_{Γ_{I}, d}^{i} d Γ d_{N_{e}} .

(35)

The right-most expression is the complementary virtual work done by the inertia forces going through a virtual variation in pressure, and we identify the transpose of the coupling matrix by writing

- ω^{2} δ p_{N_{e}}^{T} \int_{Γ_{I}} N_{Γ_{I}, p}^{j} n_{f} N_{Γ_{I}, d}^{i} d Γ d_{N_{e}} = - ω^{2} δ p_{N_{e}}^{T} Q_{e}^{T} d_{N_{e}} \Rightarrow Q_{e}^{T} = \int_{Γ_{I}} N_{Γ_{I}, p}^{j} n_{f} N_{Γ_{I}, d}^{i} d Γ .

(36)

Be aware that in the equation above,

Γ_{I}

refers to the part of the fluid element boundary that coincides with the boundary between the fluid and solid phases, i.e., the side of the fluid element that connects to a solid element.

Test of numerical approach

To derive the dispersion characteristics numerically, we evaluate the eigenvalue problem for different wavenumbers

k_{z}

⁠. The eigenvalue associated with the quadrupole eigenmode; i.e., the squared angular eigenfrequency

ω^{2}

can be associated with a point on the dispersion curve because

(f, V_{phase}) = (\frac{ω}{2 π}, \frac{ω}{k_{z}})

⁠. Using the present finite-element formulation, the vertical dimension of the model determines the wavelength

λ

⁠; i.e., boundary conditions are imposed such that the height of the model

H

spans one half-wavelength

H = λ / 2

⁠. Hence, the wavenumber is

k_{z} = 2 π / λ = π / H

(see Figure 2). The relevant boundary conditions depend on the specific mode under consideration. For a quadrupole half-period, they read

for n = 0, 1, . ., N : u_{n} (r, θ z = H) = 0; v_{n} (r, θ z = H) = 0; p_{n} (r, θ z = H) = 0 u_{n} (r, θ z = H - λ / 2) = 0; v_{n} (r, θ z = H - λ / 2) = 0; p_{n} (r, θ z = H - λ / 2) = 0 u_{n} (r = R, θ, z) = 0; v_{n} (r = R, θ, z) = 0; for n = 1, 2, . ., N : p_{n} (r = 0, z) = 0 .

(37)

In addition to the boundary conditions, we impose the symmetry condition:

w (r, θ z = H - λ / 4) = 0 .

(38)

The latter condition serves to fix the model in space. The pressure condition at the centerline

(r = 0)

⁠, which is not a real boundary condition but rather a condition reflecting continuity, makes coupling to the monopole mode possible. The quadrupole mode shape, the finite-element grid, and associated axial boundary conditions are shown in Figure 2. To benchmark the numerical model against an analytical solution, we solve a specific problem studied by Tang and Cheng (2004) with the following model parameters:

r_{w} = 0.1 m, V_{p} = 4000 m / s; V_{s} = 2300 m / s; V_{f} = 1500 m / s, ρ_{s} = 2500 kg / m^{3}; ρ_{f} = 1000 kg / m^{3} .

(39)

For this isotropic case, the quadrupole mode does not couple to any other modes. We can therefore truncate the series at

N = 1

and also eliminate all degrees of freedom with the azimuthal profile

n = 0

from the problem without loss of generality. The finite-element mesh comprises 200 fluid elements and 4000 solid elements. Regular discretization (20 elements) applies to the axial direction (⁠

z

⁠). For the radial direction, 10 elements discretize the radial direction in the fluid phase (regular discretization) while 200 elements discretize the solid phase (element size increases with radial distance). The outer radius of the model is

R = 25 m

⁠. The quadrupole becomes the dominant eigenmode (lowest eigenvalue) when the azimuthal profile corresponding to

n = 1

is the only one included in the analysis. The method known as “inverse iteration” has been used throughout, and because this method always converges to the dominant eigensolution, it converges to the solution we are after.

The phase velocity is frequency dependent, and a plot of pairs (⁠ $f, V_{phase}$ ⁠) reveals the dispersion characteristics (Figure 3). The simulated dispersion curve corresponds very well to the analytical counterpart; thus, the error related to discretization can be considered minimal.

In the present analysis, which is purely based on real number arithmetic, the cut-off frequency is identified as the frequency where transition to leaky modes occurs, i.e., where the pure shear deformation of horizontal polarization is excited, and the phase velocity approaches the true formation shear-wave velocity (⁠ $V_{S} = \sqrt{G / ρ_{s}}$ ⁠), which in this case is $2300 m / s$ ⁠. Mode shapes are shown for different frequencies, and the gradual change in mode shape toward pure shear near the cut-off frequency is displayed in Figure 4. Note that mode shapes are shown as dimensionless.

Stress distribution around a vertical borehole

Simulation of the effect of the direction and magnitude of principal stresses on elastic moduli and dispersion curves provides the framework for understanding how velocities are perturbed near the borehole.

We use Berea sandstone for our study case. For this rock, a nonlinear stress-strain relationship of the type previously described in the section “Constitutive model for stress-dependent elasticity” applies, and we can use the $V_{P}$ and $V_{S}$ versus the $σ_{hyd}$ relations published by Fang et al. (2013) (Figure 5). Because velocity and compliance are uniquely related, we can calculate the isotropic elastic compliances as a function of hydrostatic stress. Now, for any state of stress, which includes three principal stresses of different magnitude, we can derive the resulting full orthorhombic stiffness tensor as described in equations 1–6. The asymptotic compliance is approximated by the compliance at the highest hydrostatic stress point. This means the $V_{P}$ and $V_{S}$ velocities at 50 MPa (see Figure 5) are used to calculate the linear elastic compliances $S_{i j k l}^{0, ISO}$ and the trend curves (red and blue curves) are used to calculate compliances at lower confining stress $S_{i j k l}^{ISO}$ ⁠.

We consider a vertical hole subject to in situ stresses (⁠ $σ_{x x}, σ_{y y}, σ_{z z}$ ⁠). The stress distribution in the solid is established by imposing the in situ stresses on the outer circumference of the model $r = R$ and then solving the static problem. Stress-free conditions at the fluid-solid interface are applied, and the vertical displacements are fixed at $z = H$ and $z = 0$ ⁠. We truncate the Fourier series at $N = 3$ ⁠. The static boundary value problem is a nonlinear one due to the nonlinear stress-strain relationship, and the solution is derived using an incremental approach; i.e., we impose the in situ stresses in increments and solve the linearized problem at each loading step. At each loading step, the stiffness matrix is updated in accordance with the concurrent stress state. Importantly, the angular variation in stress, and therefore in the stiffness tensor, is taken into account (see Figure 1 and equation 13). At the last loading step, the distributed stiffness tensor in the model reflects the stress variation in the model as defined by the constitutive model for stress-dependent elasticity (equations 1–6).

To demonstrate the significance of stress perturbations on the near-borehole stress and velocity field, we consider two cases: (1) a case where horizontal stresses are isotropic and (2) a case where strong horizontal stress anisotropy exists:

case 1 : σ_{z z} = 14.5 MPa; σ_{x x} = 14.5 MPa; σ_{y y} = 14.5 MPa, case 2 : σ_{z z} = 14.5 MPa; σ_{x x} = 14.5 MPa; σ_{y y} = 5.0 MPa .

(40)

In an unstressed or in a hydrostatic state of stress, the stiffness tensor has isotropic symmetry properties. Hence, away from the borehole, the elastic medium corresponding to case 1, may be considered isotropic. However, near the borehole, both cases require full orthorhombic material representation, i.e., nine independent elastic constants

C_{11}

⁠,

C_{22}

⁠,

C_{33}

⁠,

C_{13}

⁠,

C_{23}

⁠,

C_{12}

⁠,

C_{12}

⁠,

C_{55}

⁠, and

C_{66}

⁠. Case 1 preserves axial symmetry, and field variables and elastic parameters are solely functions of radius. However, the variations in principal stresses in the radial, axial, and circumferential directions near the borehole make orthorhombic symmetry the highest order of symmetry applicable for describing the elasticity tensor for this case. For both cases, the elastic parameters vary spatially and so do the material axes. The framework for the anisotropy is the principal stresses

σ_{1}

and

σ_{2}

(Figure 6). Clearly, symmetry about the axis applies when the two horizontal stresses are equal in magnitude

σ_{x x} = σ_{y y}

⁠; e.g., the stresses are functions of radius but not of

θ

(Figure 6a and 6b). For the other case, where

σ_{x x} = σ_{z z} > σ_{y y}

⁠, the anisotropy is more complex as shown in Figure 7a and 7b.

Velocity distribution and quadrupole wave dispersion

Because the orthorhombic elasticity tensor is uniquely related to the stress tensor and velocities are analytical functions of moduli, we may easily derive the perturbed velocity field around the borehole. In Figure 6, we compare the derived $V_{P 1} = \sqrt{C_{11} / ρ_{s}}$ and $V_{P 2} = \sqrt{C_{22} / ρ_{s}}$ to the principal stresses (Figure 6a and 6b).

Of particular interest is the shear velocity that comprises the upper limit of the quadrupole phase velocity. Unlike the homogeneous case, this limiting velocity is not known a priori when stresses and moduli vary in space. In the stress-sensitive medium, the slowest shear wave is one that is polarized in the direction of minimum stress and propagating along one of the two other orthogonal axes. The speed of the wave that is polarized in the direction of least horizontal compression (⁠ $i = 2$ ⁠) and propagating in the vertical direction (⁠ $i = 3$ ⁠) is $V_{S 32} = \sqrt{C_{44} / ρ_{s}}$ ⁠, whereas the speed of the wave propagating along the alternate horizontal material axis is $V_{S 12} = \sqrt{C_{66} / ρ_{s}}$ ⁠.

Far away from the wellbore, the stress conditions are homogeneous; e.g., $C_{44} = C_{66} < C_{55}$ for the anisotropic case and $C_{44} = C_{55} = C_{66}$ for the isotropic case. Near the wellbore, shear velocities are influenced by the stress concentration (see Figure 8).

The spatial distributions of $V_{S 12}$ and $V_{S 32}$ for the anisotropic (a and b) and isotropic (c and d) in situ stress conditions are shown in Figure 8. For the isotropically stressed formation, the S-wave velocity attains its minimum at the borehole wall at $2300 m / s$ ⁠, which is to be compared to $2355 m / s$ farther away from the borehole.

As expected, the variation is much greater for the anisotropic case, where the range is as wide as $2100 - 2400 m / s$ for either of $V_{S 12}$ and $V_{S 32}$ ⁠. For comparison, the horizontal P-wave velocities computed for isotropic in situ stress conditions vary between 3540 and 3750 m/s (radial and circumferential directions, respectively, see Figure 7c and 7d) and the horizontal P-wave velocities computed for anisotropic in situ stress conditions vary between 3500 and 3900 m/s. We now direct our attention to how the quadrupole phase velocity and dispersion characteristics are affected by the perturbed velocity field. To this end, we again focus on the eigenvalue problem (equation 16), and again we solve it repeatedly for different wavenumbers $k_{z}$ ⁠. That is, we alter the vertical dimension of the model step by step and solve the boundary value problem (equation 16) using the boundary conditions (equation 37) and the symmetry condition (equation 38). Clearly, the solid stiffness matrix is evaluated by solving exactly the static problem described above, and the remaining matrices in the formulation of equation 16 are evaluated as described in the section “Sonic wave propagation along a borehole.” Each eigenpair comprises a mode and a frequency, and we focus consistently on the quadrupole eigenmode (see Figure 2a). In a strict sense, the quadrupole eigenmode is not characterized by the acoustic profile $\cos (2 θ)$ due to coupling related to anisotropy. The term “quadrupole eigenmode” refers, in the anisotropic case, to the eigenmode dominated by the azimuthal acoustic profile $\cos (2 θ)$ ⁠. Numerically, the quadrupole eigenpair (eigenmode and frequency) is again found by the iterative algorithm known as inverse iteration, this time with a shift. One advantage provided with this method is the convergence to the eigenvalue nearest to a predefined estimate (the shift value). An approximation to the quadrupole eigenvalue is easily obtained by truncating the series, given in equation 7, at $N = 1$ ⁠, and eliminating the degrees of freedom associated to the axial symmetric profile $n = 0$ from the eigenvalue problem (equation 16). In this way, a reduced eigenvalue problem is obtained, i.e., one where coupling is eliminated; however, the lowest eigenvalue for each wavenumber $k_{z}$ is an approximation to the quadrupole squared eigenfrequency. To each quadrupole eigenpair, we associate the mode’s phase velocity $V_{phase} = ω / k_{z}$ and frequency $f = ω / 2 π$ and plot the dispersion curve (⁠ $f, V_{phase}$ ⁠). These curves are plotted in Figure 9 for the isotropic as well as for the anisotropic cases. Also in Figure 9, the results of simulations are shown in which the stress concentration effect on the borehole hoop stress has been removed from the analysis. For the latter cases (dotted lines), the velocity fields are homogeneous and reflect conditions far away from the borehole. The curves converge toward the true formation shear-wave speed $V_{S 32} = V_{S 12}$ as expected. The curves influenced by near-borehole stresses are consistently below their undisturbed counterparts, except at the cut-off frequency where they coincide. Intuitively, this is to be expected because in either cases, the shear mode is excited in the entire solid phase and the speed of propagation is $V_{S 32} = V_{S 12}$ ⁠. However, it should be noted that shear velocity estimates in field measurements are based on extrapolations of quadrupole dispersion curves to the cut-off frequency because the mode is generally excited above the cut-off frequency only. To this end, we note that for the anisotropic case, the influence of the borehole hoop stress concentration causes the phase velocity trend to be around 1% less than without the stress influence. This would imply that the quadrupole dispersion is essentially insensitive to near-wellbore effects for the cases considered, even when the spatial variation of shear velocities is quite significant near the borehole. This is consistent with approximations derived from perturbation analyses (Sinha et al., 2003) and recent experimental results (Hsu et al., 2011).

Information about anisotropy can be extracted from analysis of quadrupole modes. The quadrupole mode, which is associated with azimuthal variation $\cos (2 n θ)$ ⁠, couples to modes with azimuthal variation $\cos (4 n θ)$ ⁠, $\cos (6 n θ)$ …, as well as to axisymmetrical modes. The strength of coupling reflects the degree of anisotropy. Also, the azimuthal directions (⁠ $θ, θ + π / 2$ ⁠) of maximum amplitude coincide with the directions of maximum and minimum in situ horizontal stress.

CONCLUSION

We have derived an integrated numerical finite-element approach, which includes couplings between nonlinear elastic deformation, near-borehole stress redistribution, and borehole acoustics. The coupling between stress and elasticity means that spatial variations in elastic anisotropy apply near a borehole, where stress concentration effects are present. We demonstrate how the method is constructed to include heterogeneous anisotropy effects for the case of a vertical wellbore in a stress regime where two in situ principal stress directions are horizontal and one is vertical. We have used the method to study the dispersion characteristics of quadrupole waves traveling along a borehole. Our analysis suggests that the dispersion characteristics of quadrupole waves are essentially insensitive to the spatial variation of elastic anisotropy immediately around the borehole. Hence, extrapolation of a quadrupole dispersion curve may provide estimates of the true formation shear velocity, independently of the spatial variation of this parameter in the close proximity of the borehole. This finding, which is in line with experimental results and approximations derived from perturbation analysis, is here confirmed in a full 3D analysis.

The presented method is based on a ring-type finite element, thus a radial/axial plane is taken as a plane of reference. The field variation out of this plane, i.e., the displacement field in the rock and the pressure field in the fluid, is restricted to a subspace formed by a truncated Fourier series. Full 3D analysis, at a computational effort that is remarkably reduced compared to conventional 3D finite-element modeling, can thus be accomplished for a small but very important class of problems by the proposed method. It is noteworthy that this class of problems includes the possibility of taking the effect of the bottom of the hole into account in the modeling. Hence, the method has potential for investigating applications for imaging ahead of the bit during drilling.

ACKNOWLEDGMENTS

The authors wish to thank Maersk Oil and Gas AS for sponsoring this work and the Earth Resources Laboratory (ERL) at MIT for hosting it. The authors also wish to thank the students, staff, and faculty at MIT-ERL for many fruitful discussions.

TRANSFORMATION OF ORTHORHOMBIC STIFFNESS PROPERTIES

In the following, we will refer to the Cartesian coordinate system

(x, y, z)

using index notation

x_{i} = (x_{1}, x_{2}, x_{3})

⁠. In a Cartesian reference coordinate system

x_{i}

aligned with the axes of material symmetry

x_{i}^{'}

⁠, we may write orthorhombic stiffness tensor in the form

{\begin{matrix} σ_{11} \\ σ_{22} \\ σ_{33} \\ σ_{23} \\ σ_{13} \\ σ_{12} \end{matrix}} = [\begin{matrix} C_{11} & C_{12} & C_{13} & 0 & 0 & 0 \\ C_{12} & C_{22} & C_{23} & 0 & 0 & 0 \\ C_{13} & C_{23} & C_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & C_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & C_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & C_{66} \end{matrix}] {\begin{matrix} ϵ_{11} \\ ϵ_{22} \\ ϵ_{33} \\ 2 ϵ_{23} \\ 2 ϵ_{13} \\ 2 ϵ_{12} \end{matrix}} .

(A-1)

We can transform equation into cylindrical coordinates and let the angle

θ

denote the rotation about axis

x_{3}

relative to axis

x_{1}

⁠:

{\begin{matrix} σ_{r r} \\ σ_{θ θ} \\ σ_{z z} \\ σ_{θ z} \\ σ_{r z} \\ σ_{r θ} \end{matrix}} = [\begin{matrix} {\bar{C}}_{11} & {\bar{C}}_{12} & {\bar{C}}_{13} & 0 & 0 & {\bar{C}}_{16} \\ {\bar{C}}_{12} & {\bar{C}}_{22} & {\bar{C}}_{23} & 0 & 0 & {\bar{C}}_{26} \\ {\bar{C}}_{13} & {\bar{C}}_{23} & {\bar{C}}_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & {\bar{C}}_{44} & {\bar{C}}_{45} & 0 \\ 0 & 0 & 0 & {\bar{C}}_{45} & {\bar{C}}_{55} & 0 \\ {\bar{C}}_{16} & {\bar{C}}_{26} & 0 & 0 & 0 & {\bar{C}}_{66} \end{matrix}] {\begin{matrix} ϵ_{r r} \\ ϵ_{θ θ} \\ ϵ_{z z} \\ 2 ϵ_{θ z} \\ 2 ϵ_{r z} \\ 2 ϵ_{r θ} \end{matrix}} .

(A-2)

In equation ,

{\bar{C}}_{11} = U_{1} + U_{2} \cos (2 θ) + U_{3} \cos (4 θ); {\bar{C}}_{12} = U_{4} - U_{3} \cos (4 θ); {\bar{C}}_{13} = U_{6} - U_{7} \cos (2 θ); {\bar{C}}_{16} = - \frac{1}{2} U_{2} \sin (2 θ) - U_{3} \sin (4 θ); {\bar{C}}_{22} = U_{1} - U_{2} \cos (2 θ) + U_{3} \cos (4 θ); {\bar{C}}_{23} = U_{6} + U_{7} \cos (2 θ); {\bar{C}}_{26} = - \frac{1}{2} U_{2} \sin (2 θ) + U_{3} \sin (4 θ); {\bar{C}}_{33} = C_{33}; {\bar{C}}_{44} = U_{8} - U_{9} \cos (2 θ); {\bar{C}}_{45} = - U_{9} \sin (2 θ); {\bar{C}}_{55} = U_{8} + U_{9} \cos (2 θ); {\bar{C}}_{66} = U_{5} - U_{3} \cos (4 θ); U_{1} = \frac{1}{8} (3 C_{11} + 3 C_{22} + 2 C_{12} + 4 C_{66}); U_{2} = \frac{1}{2} (C_{11} - C_{22}); U_{3} = \frac{1}{8} (C_{11} + C_{22} - 2 C_{12} - 4 C_{66}); U_{4} = \frac{1}{8} (C_{11} + C_{22} + 6 C_{12} - 4 C_{66}); U_{5} = \frac{1}{8} (C_{11} + C_{22} - 2 C_{12} + 4 C_{66}); U_{6} = \frac{1}{2} (C_{23} + C_{13}); U_{7} = \frac{1}{2} (C_{23} - C_{13}); U_{8} = \frac{1}{2} (C_{55} + C_{44}); U_{9} = \frac{1}{2} (C_{55} - C_{44}) .

(A-3)

In the analysis presented in this paper, the material axes

x_{1}^{'}

and

x_{2}^{'}

(or equivalently,

x^{'}

and

y^{'}

⁠) rotate about the vertical axis of reference,

x_{3}

(or

z

⁠), due to the perturbed stress field near the borehole. If

θ_{1}

denotes the local angle of rotation (see Figure 1), then

θ - θ_{1}

replaces

θ

in the equations above.

Novel finite-element approach applied to borehole quadrupole dispersion analysis in stress-sensitive formations

ABSTRACT

INTRODUCTION

FORMULATION

Constitutive model for stress-dependent elasticity

3D modeling of stresses and strains

Sonic wave propagation along a borehole

Test of numerical approach

Stress distribution around a vertical borehole

Velocity distribution and quadrupole wave dispersion

CONCLUSION

ACKNOWLEDGMENTS

TRANSFORMATION OF ORTHORHOMBIC STIFFNESS PROPERTIES

Data & Figures

About GSW

Resources

Explore

Connect

Novel finite-element approach applied to borehole quadrupole dispersion analysis in stress-sensitive formations Open Access

ABSTRACT

INTRODUCTION

FORMULATION

Constitutive model for stress-dependent elasticity

3D modeling of stresses and strains

Sonic wave propagation along a borehole

Test of numerical approach

Stress distribution around a vertical borehole

Velocity distribution and quadrupole wave dispersion

CONCLUSION

ACKNOWLEDGMENTS

TRANSFORMATION OF ORTHORHOMBIC STIFFNESS PROPERTIES

Data & Figures

About GSW

Resources

Explore

Connect

This Feature Is Available To Subscribers Only

Novel finite-element approach applied to borehole quadrupole dispersion analysis in stress-sensitive formations