Modeling and simulation of an acoustic well stimulation method

TL;DR

This paper develops a mathematical and numerical model to simulate acoustic well stimulation, aiming to optimize sound wave transmission for enhancing rock permeability around oil and gas wells.

Contribution

It introduces a novel Helmholtz equation-based model with impedance boundary conditions and non-reflecting boundaries, enabling efficient simulation and optimization of AWS.

Findings

Model accurately predicts acoustic energy transmission.

Optimal emission frequencies can be effectively identified.

Numerical examples demonstrate the method's effectiveness.

Abstract

This paper presents a mathematical model and a numerical procedure to simulate an acoustic well stimulation (AWS) method for enhancing the permeability of the rock formation surrounding oil and gas wells. The AWS method considered herein aims to exploit the well-known permeability-enhancing effect of mechanical vibrations in acoustically porous materials, by transmitting time-harmonic sound waves from a sound source device---placed inside the well---to the well perforations made into the formation. The efficiency of the AWS is assessed by quantifying the amount of acoustic energy transmitted from the source device to the rock formation in terms of the emission frequency and the well configuration. A simple methodology to find optimal emission frequencies for a given well configuration is presented. The proposed model is based on the Helmholtz equation and an impedance boundary condition that effectively accounts for the porous solid-fluid interaction at the interface between the rock formation and the well perforations. Exact non-reflecting boundary conditions derived from Dirichlet-to-Neumann maps are utilized to truncate the circular cylindrical waveguides considered in the model. The resulting boundary value problem is then numerically solved by means of the finite element method. A variety of numerical examples are presented in order to demonstrate the effectiveness of the proposed procedure for finding optimal emission frequencies.

Tables

Table 1. Table 1: Physical constants for crude oil and sandstone. The numerical values of the fluid constants were taken from  [ 5 ] . The numerical values of the porous solid constants, on the other hand, were taken from  [ 26 ] .

Constant	Value
Speed of sound in oil ($c$)	1524 m s-1
Oil density (${\rho }_0$)	1100 kg m-3
Oil shear viscosity ($\eta$)	1.2 Pa s
Oil bulk modulus ($B$)	3000 MPa
Rock formation permeability ($\kappa$)	$3\times {10}^{-13}$ m2
Rock formation porosity ($\varphi$)	$0.21$

Table 2. Table 2: Geometrical parameters of the well configurations considered. The dimensions of the transducer were selected according the device described in  [ 20 ] . The dimensions of a perforated well, on the other hand, were taken from  [ 11 ]

Parameter	Value
Well radius ($R$)	0.111 m
Perforated domain height ($2H$)	1.800 m
Transducer length	1.410 m
Transducer radius	0.054 m
Perforation radius ($r_0$)	0.020 m
Perforation depth	0.305 m
Perforation spacing (1st well)	0.257 m
Phasing angle (1st well)	$\pi /2$ rad
Perforation spacing (2nd well)	0.200 m
Phasing angle (2nd well)	$\pi /3$ rad
Perforation spacing (3rd well)	0.160 m
Phasing angle (3rd well)	$\pi /6$ rad

Equations

$$p$$

$$-\frac{i\omega }{c}p+{\rho }_0\mathrm{div}\mathbf{v}$$

$$-i\omega \mathbf{v}+\frac{1}{{\rho }_0}\nabla p$$

$$\Delta p+k^{2}p=0$$

$$\frac{\partial p}{\partial n}-\frac{ik}{\zeta }p=g$$

$$I_{\mathrm{abs}}=\frac{1}{2{\rho }_{0}c}{\int }_{{\Gamma }_w}\frac{|p{|}^2}{|\zeta {|}^2}\chi \mathrm{d}s.$$

$$I_{\mathrm{rad}}=\frac{1}{2{\rho }_0\omega }{\int }_{{\Gamma }_s}\mathrm{Im}\left\{p\bar{g}\right\}\mathrm{d}s.$$

$$Z=\frac{p}{\boldsymbol{v}\cdot\boldsymbol{n}}=\Bigg{(}\frac{\kappa\sqrt{i\omega m}}{\eta}\frac{H^{(1)}_{1}(\sqrt{i\omega m}r_{0})}{H^{(1)}_{0}(\sqrt{i\omega m}r_{0})}\Bigg{)}^{-1},$$

$$\mathbf{v}\cdot \mathbf{n}=-\frac{i}{\omega {\rho }_0}\frac{\partial p}{\partial n},$$

$$\frac{\partial p}{\partial n}-\frac{ikc{\rho }_0}{Z}p=0.$$

$$\zeta=\Bigg{(}\frac{\rho_{0}c\kappa\sqrt{i\omega m}}{\eta}\frac{H^{(1)}_{1}(\sqrt{i\omega m}r_{0})}{H^{(1)}_{0}(\sqrt{i\omega m}r_{0})}\Bigg{)}^{-1}.$$

$$\Delta p+k^{2}p$$

$$\frac{\partial p}{\partial n}-\frac{ik}{\zeta }p$$

$$\frac{\partial p}{\partial n}$$

$$p(r,\theta ,z)=\sum _{n=-\infty }^{\infty }\sum _{m=1}^{\infty }{p}_{n,m}^{\pm }{v}_{n,m}(r,\theta )e^{\pm i(z\mp H)\sqrt{k^2-{\lambda }_{n,m}^2}}\text{mbox}in{\Omega }^{\pm },$$

$$v_{n,m}(r,\theta )=c_{n,m}{J}_n\left({\lambda }_{n,m}r\right)e^{in\theta }\text{mbox}and{\lambda }_{n,m}=\frac{j_{n,m}^{\prime }}{R},$$

$$c_{n,m}=\left\{\begin{array}[]{ccc}\displaystyle\frac{\lambda_{n,m}}{\sqrt{2\pi}\sqrt{\lambda^{2}_{n,m}R^{2}-n^{2}}J_{n}(\lambda_{n,m}R)}&\mbox{if}&\lambda_{n,m}>0,\vskip 6.0pt plus 2.0pt minus 2.0pt\\ \displaystyle\frac{1}{\sqrt{2\pi}R}&\mbox{if}&\lambda_{n,m}=0,\end{array}\right.$$

$$\Delta v_{n,m}+{\lambda }_{n,m}^2{v}_{n,m}=0\mathrm{in}\{r<R\},\frac{\partial v_{n,m}}{\partial n}=0\mathrm{on}\{r=R\},\text{mbox}and{\int }_{\{r<R\}}|v_{n,m}{|}^2=1.)$$

$$p_{n,m}^{\pm }={\int }_0^R{\int }_0^{2\pi }p(r,\theta ,\pm H)\overline{v_{n,m}(r,\theta )}r\mathrm{d}\theta \mathrm{d}r,-\infty <n<\infty ,m\ge 1,$$

$$\begin{array}{r}{\mathcal{T}}^{\pm }\left[p\right](\mathbf{x})=\sum _{n=-\infty }^{\infty }\sum _{m=1}^{\infty }i\sqrt{k^2-{\lambda }_{n,m}^2}{p}_{n,m}^{\pm }{v}_{n,m}(r,\theta ),\mathbf{x}=(r\cos \theta ,r\sin \theta ,\pm H)\in {\Gamma }^{\pm }.\end{array}$$

$$\Delta p+k^{2}p$$

$$\frac{\partial p}{\partial n}-\frac{ik}{\zeta }p$$

$$\frac{\partial p}{\partial n}$$

$$a(p,q)=f(q),\forall q\in H^1(\Omega ),$$

$$a(p,q)$$

$$f(q)$$

$$p_h(\mathbf{x})=\sum _{i=1}^N{p}_i{\varphi }_i(\mathbf{x}),\mathbf{x}\in \Omega ,$$

$$A\mathbf{p}=\mathbf{f},$$

$$g(\mathbf{x})=\frac{\partial G}{\partial n_{\mathbf{x}}}(\mathbf{x},\mathbf{y}),\mathbf{x}\in {\Gamma }_s,$$

$$G(\mathbf{x},\mathbf{y})=\sum _{n=-\infty }^{\infty }\sum _{m=1}^{\infty }\frac{v_{n,m}(r,\theta )\overline{v_{n,m}(\rho ,\vartheta )}}{\sqrt{{\lambda }_{n,m}^2-k^2}}{e}^{-\sqrt{{\lambda }_{n,m}^2-k^2}|z-\zeta |},$$

$$p(\mathbf{x})=G(\mathbf{x},\mathbf{y}),\mathbf{x}\in \tilde{\Omega }={\Omega }_w\setminus {\overline{\Omega }}_s,$$

$$E_h=\frac{\Vert p_h-{\Pi }_{h}p{\Vert }_{L^2\left(\Omega \right)}}{\Vert {\Pi }_{h}p{\Vert }_{L^2\left(\Omega \right)}},$$

$$Q_j$$

Full text

1 Introduction

The decrease of oil and gas recovery from a reservoir is clearly an important problem that affects the energy industry. One of the main causes of such problem is the local reduction of the reservoir permeability around producing wells due to the deposition of scales, precipitants and mud penetration during exploitation which, over time, give rise to an impermeable barrier to fluid flow [10]. Well stimulation methods play a prominent role in the exploitation of these essential natural resources as they are intended to increase the permeability of the reservoir, allowing the trapped fluid to flow toward the borehole and thus enhancing the productivity of the well. Various well stimulation methods are used in practice to cope with local deposits, including solvent and acid injection, treatment by mechanical scrapers and high pressure fracturing. Each one of these conventional methods have significant drawbacks and undesirable effects. Some of them, for instance, are expensive and produce damage to the well structure, while others are highly polluting, leading to harmful ecological effects associated with the contamination of underground water resources [10, 4]. The demonstrated effectiveness of mechanical vibrations on enhancing fluid flow through porous media [4, 12, 2], on the other hand, has led to the development of the so-called acoustic well stimulation (AWS) methods, which nowadays have broad acceptance by the hydrocarbon industry mainly due to the fact that they partially overcome the aforementioned issues.

This paper considers an AWS method based on the transmission of acoustic waves, emitted by a transducer submerged into the well, to the rock formation surrounding the well. The transducer is designed to trigger one of the physical processes known to enhance the permeability of the porous medium. Among such physical processes, we mention the reduction of the fluid viscosity by agitation and heating, stimulation of elastic waves on the well walls (to reduce the adherence forces in the layer between oil and rock formation), excitation of natural frequencies associated with the vibration of the fluid inside the porous medium, and the formation and collapse of cavitation bubbles near clogged pores of the rock formation. A variety of transducer designs have been proposed over the last three decades, which consider operation frequency and intensity ranges selected to target one (or several) of the aforementioned physical processes [21, 6, 19, 18].

This paper presents a mathematical model and a numerical procedure that allows us to find optimal emission frequencies for which the amount of energy transmitted from the transducer into the rock formation is maximized. The proposed methodology can potentially improve the performance of the whole class AWS methods considered, as the aforementioned physical processes take place within the porous medium. In detail, we develop a mathematical model based on the Helmholtz equation and an impedance boundary condition [7] that effectively accounts for the porous solid-fluid interaction at the interface between the rock formation and the well perforations [26]. Exact non-reflecting boundary conditions derived from Dirichlet-to-Neumann (DtN) maps are utilized to truncate the circular cylindrical waveguides considered in the model [9, 22, 23]. The resulting boundary value problem is numerically solved by means of the finite element (FE) method [25, 15, 13]. Optimal emission frequencies are then found by scanning the quotient of the emitted energy to the transmitted energy—toward the region of interest—over a range of frequencies. As expected, the optimal emission frequencies correspond to field distributions for which resonances occur inside the perforations.

The outline of this paper is as follows: The mathematical model is presented in Section 2. The DtN-FE method is then described and validated in Section 3. Section 4 provides numerical results for realistic well configurations. Section 5, finally, gives the concluding remarks of the present work.

2 Mathematical model

2.1 Geometry

A perforated well is created through two successive processes called drilling and completion. The former begins by drilling a borehole in the ground, which is covered by metal pipes that are attached to its walls by a layer of cement (cf. Figure 1). This part of the process, commonly referred to as casing, aims to stabilize the borehole structure. Once the well is cased, the completion process begins by shooting with explosives the portion of the casing that passes through the reservoir level—where the oil is trapped—forming small holes across the casing and the cement layer, and into the reservoir. These holes, referred to as perforations, are aimed at enabling the oil to flow from the reservoir into the well.

Upon completion, two different zones of the well can be identified; the zone containing the perforations, which we call the perforated domain, and the remaining part of the well, which we call the cylindrical domain. The perforated domain, denoted by $\Omega_{p}$, is assumed to be bounded. In addition, we assume that the cylindrical domain consists of two (semi-infinite) circular cylinders placed above and below the perforated domain, which we denote by $\Omega^{+}$ and $\Omega^{-}$, respectively. The model of a perforated well utilized in this paper then, corresponds to a locally perturbed circular cylinder defined as $\Omega_{w}=\Omega_{p}\cup\Omega^{+}\cup\Omega^{-}$. The interface between the perforated and the upper (resp. lower) cylindrical domains is denoted by $\Gamma^{+}$ (resp. $\Gamma^{-}$). Finally, the transducer (source) is assumed to occupy the bounded domain $\Omega_{s}\subset\Omega_{p}$ with boundary $\partial\Omega_{s}=\Gamma_{s}$. We refer to Figure 2 for the definition of all the relevant domains considered in the mathematical model.

2.2 Acoustic waves

The transducer is herein modeled as a time-harmonic vibrating surface $\Gamma_{s}$ that operates at a fixed frequency $f=\omega/2\pi$, where $\omega>0$ denotes the angular frequency in radians. Being excited by a single time-harmonic source, the pressure $P$, the density $\varrho$, and the velocity $\boldsymbol{V}$ fields eventually reach a stationary (time-harmonic) regime for which $P(\boldsymbol{x},t)={\mathrm{Re}\,}\left\{p(\boldsymbol{x})\operatorname{e}^{-i\omega t}\right\}$, $\varrho(\boldsymbol{x},t)={\mathrm{Re}\,}\left\{\rho(\boldsymbol{x})\operatorname{e}^{-i\omega t}\right\}$, and $\boldsymbol{V}(\boldsymbol{x},t)={\mathrm{Re}\,}\left\{\boldsymbol{v}(\boldsymbol{x})\operatorname{e}^{-i\omega t}\right\}$, where $t>0$ denotes the time variable and $p$, $\rho$ and $\boldsymbol{v}$ denote the amplitudes of the pressure, the density and the velocity, respectively, which only depend on the position $\boldsymbol{x}$. The linearized equations of state and conservation of mass and momentum in this case, read as [7, 17]

[TABLE]

where $c>0$ and $\rho_{0}>0$ denote the speed of sound and the equilibrium density of the fluid that fills the well, respectively. Suitably combining equations (1a), (1b) and (1c) we then obtain that $p$ satisfies the Helmholtz equation

[TABLE]

in the domain $\widetilde{\Omega}=\Omega_{w}\setminus\overline{\Omega_{s}}$ occupied by the fluid, where $k=\omega/c$ denotes the wavenumber.

Note that dissipation effects can be easily taken into account by considering a complex wavenumber with spatial absorption depending on the equilibrium density and the shear and bulk viscosities [17]. For presentation simplicity, however, we only consider real wavenumbers.

2.3 Boundary conditions

Throughout this paper we consider boundary conditions of the form

[TABLE]

on the surfaces of the well ($\Gamma_{w}=\partial\Omega_{w}$) and on the transducer ($\Gamma_{s}$), where $\zeta\in\mathbb{C}$ denotes the dimensionless surface impedance and the function $g$ corresponds to the excitation prescribed on the surface $\Gamma_{s}$ of the transducer. The dimensionless impedance takes the form $\zeta=\chi+\operatorname{i}\xi$, where $\chi$ and $\xi$ ($\chi,\xi:\Gamma_{w}\cup\Gamma_{s}\to\mathbb{R}$) are known as the resistive (real) and reactive (imaginary) parts of the impedance, respectively. The dimensionless impedance $\zeta$ and the pressure field $p$ are related to the time-averaged energy flux through $\Gamma_{w}$ by the formula [7]

[TABLE]

The time-averaged acoustic energy radiated by transducer, on the other hand, is given by

[TABLE]

The spatial dependence of the dimensionless impedance $\zeta$ in (3) is determined by the mechanical properties of the various materials that are in direct contact with the fluid. Being the casing made of metal (see Section 2.1)—which is usually modeled as a sound hard (Neumann) boundary condition—the admittance $1/\zeta$ is taken equal to zero over the cylindrical domain and the cased portion of the perforated domain. The sound hard boundary condition ($1/\zeta=0$) is also used on the transducer $\Gamma_{s}$. In order to determine suitable impedance values to be used over the boundary of the perforations, in turn, we follow the analytical calculations presented by J. E. White in [26] for the wall impedance at the interface between a liquid and a porous material. According to these calculations, the wall impedance $Z$—defined as the quotient of the pressure amplitude to the normal velocity amplitude on the boundary of the perforation— is given by

[TABLE]

where $H^{(1)}_{0}$ and $H^{(1)}_{1}$ denote the Hankel functions of the first kind and order zero and one, respectively [1], $r_{0}$ is the radius of the perforation, $\kappa$ is the permeability of the porous medium, $\eta$ is the shear viscosity of the fluid, and $m=\phi\eta/(\kappa B)$, being $\phi$ the porosity and $B$ the bulk modulus of the fluid in the pore space. It is important to highlight that the impedance model (6) is valid under the assumption that $r_{0}$ is smaller than the wavelength $\lambda=2\pi/k$. For the sake of completeness, the analytical derivations leading to (6) are reproduced in A. On the other hand, in order to link $Z$ with the dimensionless surface impedance $\zeta$ we get, from the momentum conservation equation (1c), the relation

[TABLE]

which combined with the definition of $Z$ yields

[TABLE]

From (3) with $g=0$ and (7), we obtain that $Z=\rho_{0}c\zeta$. Therefore, the dimensionless surface impedance to be utilized in (3) on the surface of the perforations is given by

[TABLE]

2.4 Boundary value problem

We are now in position to put together the boundary value problem to be solved in what follows of this paper. The time-harmonic pressure field $p:\widetilde{\Omega}\to\mathbb{C}$, which is driven by the transducer submerged into the well, satisfies

[TABLE]

where the dimensionless impedance $\zeta$ is given by (8) on the boundary of perforations and it equals infinity (i.e., $1/\zeta=0$) everywhere else on $\Gamma_{w}$ (see Section 2.3). In order for the boundary value problem (9) to be well-posed, $p$ has to satisfy a certain radiation condition—which differs from the classical Sommerfeld condition—that is expressed in terms of the propagative modes associated with the upper and lower unbounded cylindrical domains $\Omega^{+}$ and $\Omega^{-}$ [9, 22, 23].

3 Dirichlet-to-Neumann Finite Element Method

3.1 The DtN map

In what follows we present a DtN-FE method for the numerical solution of (9). Notice that standard finite element (FE) methods do not directly apply to this problem due to the unboundedness of the domain $\widetilde{\Omega}$. The DtN-FE method is based on the DtN operators $\mathcal{T}^{\pm}$ that map the boundary values $p|_{\Gamma^{\pm}}$ on $\Gamma^{\pm}$ into the corresponding normal derivatives $\partial p/\partial n|_{\Gamma^{\pm}}$ on $\Gamma^{\pm}$ [9, 22, 3]. As these DtN maps provide exact non-reflecting boundary conditions on $\Gamma^{\pm}$ they allow us to write a boundary value problem posed on the bounded domain $\Omega=\widetilde{\Omega}\setminus\overline{(\Omega^{+}\cup\Omega^{-})}=\Omega_{p}\setminus\overline{\Omega}_{s}$ that is equivalent to (9) and is suitable to be solved by FE methods (or any other standard numerical method for solving PDEs).

In order to provide explicit expressions for the DtN maps, we first introduce a cylindrical coordinate system $(r,\theta,z)$, with $r\geq 0,0\leq\theta\leq 2\pi$ and $z\in\mathbb{R}$, upon which the upper and lower cylindrical domains can be expressed as $\Omega^{\pm}=\{r<R,\pm z>H\}\subset\mathbb{R}^{3}$, where $H>0$ denotes the truncation height and $R>0$ denotes the radius of the well. The series representation of the desired DtN maps are then obtained by applying the method of separation of variables to solve the Helmholtz equation in the domains $\Omega^{\pm}$ with Neumann boundary condition on the surface $\{r=R\}$. Enforcing the radiation condition—by eliminating both down-going (resp. up-going) and exponentially growing solutions in $\Omega^{+}$ (resp. $\Omega^{-}$)—we obtain the following Fourier-Bessel series for the pressure field [22]

[TABLE]

where, letting $j^{\prime}_{n,m}\geq 0$ denote the $m$-th non-negative zero of the derivative of the Bessel function of first kind $J_{n}$, we have that

[TABLE]

with

[TABLE]

correspond to the normalized Neumann-Laplace eigenfunctions and eigenvalues of the circle $\{r<R\}\subset\mathbb{R}^{2}$, respectively (i.e., they satisfy

[TABLE]

The Fourier coefficients $p^{\pm}_{n,m}$ in (10), in turn, are given by

[TABLE]

where $p(r,\theta,\pm H)=p|_{\Gamma^{\pm}}$. Taking normal derivative of (10) on $\Gamma^{\pm}$ (with unit normal vectors pointing toward $\Omega^{\pm}$) we finally arrive at the following expression for the DtN maps

[TABLE]

3.2 Equivalent boundary value problem

Using the continuity of the pressure field and its normal derivative across $\Gamma^{\pm}$ we thus obtain the following equivalent boundary value problem

[TABLE]

for the pressure field in the bounded domain $\Omega$.

Multiplying the Helmholtz equation (12a) across by a test function $q\in H^{1}(\Omega)$ and integrating by parts, we arrive at the variational (or weak) formulation of (12), which is expressed as follows: Find $p\in H^{1}(\Omega)$ such that

[TABLE]

where

[TABLE]

The well-posedness of the variational problem (13) can be easily established following the analysis presented in [9].

3.3 Finite element discretization

The discretization of the variational formulation (14) by finite elements is straightforward. We consider a family of regular tetrahedral meshes $\mathcal{T}_{h}$ of the domain $\Omega$, such that $\overline{\Omega}=\bigcup_{T\in\mathcal{T}_{h}}T$ ($\Omega$ is assumed to be a tetrahedral domain) where $h=\max\{\mathrm{diam}\,T:T\in\mathcal{T}_{h}\}$, with $\mathrm{diam}\,T=\max\{|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}|:\boldsymbol{x}_{1},\boldsymbol{x}_{2}\in T\}$. Using standard linear Lagrange elements, the approximate solution $p_{h}$ of (14) is expressed as

[TABLE]

where $N$ is the number of nodes of the mesh and $\{\phi_{1},\phi_{2},\ldots,\phi_{N}\}$ is the nodal basis of the finite dimensional function space V_{h}=\big{\{}q\in H^{1}(\Omega):\ q\in\mathcal{C}^{0}(\Omega),\ q\mid_{T}\in\mathcal{P}_{1}(T),\ \forall\,T\in\mathcal{T}_{h}\big{\}}\subset H^{1}(\Omega) where $\mathcal{C}^{0}(\Omega)$ denotes the set of continuous functions in $\Omega$, and $\mathcal{P}_{1}(T)$ denotes the set of polynomials of degree at most one defined in $T$. A system of equations for the node values $p_{i}$, $i=1,\ldots,N$ in (15) is obtained by substituting $p$ by $p_{h}$ in (14) and taking test functions $q_{h}$ from the nodal basis of $V_{h}$. Doing so, and further replacing the bilinear form $a$ by an approximate bilinear form $\tilde{a}$, given by (14a) but with the DtN maps $\mathcal{T}^{\pm}$ in the last two integrals expressed in terms of truncated series representations, we obtain the linear system

[TABLE]

where $A_{ij}=\tilde{a}(\phi_{i},\phi_{j})$, $1\leq i,j\leq N$, $\boldsymbol{p}=[p_{1},\ldots,p_{N}]^{T}$ and $\boldsymbol{f}=[f(\phi_{1}),\ldots,f(\phi_{N})]^{T}$. In order to ensure the uniqueness of the solution of the linear system, it suffices to consider truncated series representations of the DtN maps that include all the modes satisfying $|\lambda_{n,m}|\leq k$ [14].

Remark 3.1**.**

It is worth mentioning that one of the main advantages of the proposed absorbing boundary conditions over perfectly matched layers (PMLs) lies in the fact that the absorbing boundaries $\Gamma^{\pm}$ can be placed arbitrarily close to the region of interest (near the perforations and the transducer) provided that a sufficiently large number of modes are considered in the truncated series representations of the DtN maps. Off-the-shelf PMLs that absorb only propagative modes, on the other hand, would have to be placed far away enough from the region of interest so that all the evanescent modes are sufficiently attenuated, leading to larger computational domains and larger linear systems. Alternatively, PMLs that absorb both propagative and evanescent modes can also be used, provided that the mesh is properly refined to account for the frequency increment within the absorbing layers [16].

3.4 Validation

In this section, we present a numerical experiment devised to validate the proposed DtN-FE method. We thus consider a test geometry consisting of a non-perforated well and a spherical transducer, given by $\Omega_{w}=\left\{\boldsymbol{x}=(r\cos\theta,r\sin\theta,z)\in\mathbb{R}^{3}:r<R\right\}\subset\mathbb{R}^{3}$ and $\Omega_{s}=\left\{\boldsymbol{x}\in\mathbb{R}^{3}:|\boldsymbol{x}-\boldsymbol{y}|<\delta\right\}$, respectively, where $\Omega_{s}$ is centered at a point $\boldsymbol{y}\in\Omega_{w}$ and $\delta>0$ is small enough so that $\overline{\Omega_{s}}\subset\Omega_{w}$. On the spherical surface of the transducer, we prescribe the excitation

[TABLE]

where $G$ is the Green’s function of the infinite cylinder with homogeneous Neumann boundary conditions, which can be expressed in terms of the Neumann-Laplace eigenfunctions [24] as

[TABLE]

with $\boldsymbol{x}=(r\cos\theta,r\sin\theta,z)$ and $\boldsymbol{y}=(\rho\cos\vartheta,\rho\sin\vartheta,\zeta)$. It is easy to verify, from the definition of the Green’s function, that

[TABLE]

is in fact the exact solution of (9) for the test geometry considered. This exact solution (17) is then compared with approximate solutions obtained by means of the DtN-FE method described in Section 3 for various mesh sizes $h>0$. In order to compare both the exact and the approximate solution, we define the relative error

[TABLE]

where $\Pi_{h}p$ denotes the Lagrange interpolation of the exact solution using the tetrahedral mesh $\mathcal{T}_{h}$.

The results of this numerical experiment are presented in Figure 3, which displays the relative numerical errors (18) for the test problem with $R=0.5$, $\boldsymbol{y}=(0,0.25,0)$ and $\delta=0.2$. The unbounded computational domain $\widetilde{\Omega}$ was truncated by introducing artificial boundaries $\Gamma^{\pm}$ placed at $z=\pm H$, with $H=1.5$. Clearly, the numerical solution converges to the exact solution as the grid size tends to zero at a rate that is slightly faster than the expected second-order rate.

4 Numerical simulations

This section presents numerical simulations of the AWS method modeled in this paper. The values of the relevant physical constants of the fluid and the porous material—needed to evaluate the wavenumber $k=\omega/c$ and the surface impedance $\zeta$ in (8)—are displayed in Table 1. In detail, the fluid is assumed to be crude oil, with physical constants taken from [2], and the porous rock formation is assumed to be sandstone, with permeability and porosity values obtained from [27].

In order to properly simulate the operation of the AWS method, the excitation $g$ on the surface of the transducer has to be suitably prescribed. For that purpose, the transducer is modeled as a constant-amplitude time-harmonic vibrating surface $\Gamma_{s}$ with $g=1\ \mathrm{N\,m}$. A more sophisticated transducer model can be easily incorporated into the simulations by considering more general functions $g\in H^{-1/2}(\Gamma_{s})$.

The generic well configuration to be considered in the simulations is depicted in Figure 4, which includes the definition of the relevant geometrical parameters. Three particular well configurations are initially considered, with specific geometrical parameters provided in Table 2. The 1st, 2nd and 3rd well configurations include $N_{p}=6,8$ and $10$ perforations, respectively. The resulting computational domains, which were meshed using Gmsh [8], are shown in Figure 5.

Next, we compute the energy transmission through the surface of the perforations and the energy emitted by the transducer using formulae (4) and (5), respectively, for a certain range of frequencies $f=\omega/(2\pi)$. In order to find (local) optimal emission frequencies, we look for local maxima of the individual and consolidated energy transmission factors, that is,

[TABLE]

respectively—which are dimensionless quantities—as functions of the excitation frequency $f=\omega/(2\pi)$. The pressure field $p$ in (19a) and (19b) corresponds to the solution of (12) and $\Gamma_{p}^{j}\subset\Gamma_{p}$, $j=1,\cdots,N_{p}$, denotes the surface of the $j$-th perforation (the perforations are sorted from top to bottom). Note that in virtue of the conservation of energy principle and the fact that $0\leq Q\leq 1$, we have that the quantity $100\times Q$ corresponds to the percent of energy effectively transmitted to the porous reservoir rock through the perforations.

Figure 6 displays the consolidated energy transmission factor $Q$ as a function of the excitation frequency $f=\omega/(2\pi)$ for the three well configurations laid out in Table 2 and Figure 5. In these numerical simulations, the transducer is placed exactly at the center of the perforated domain. Sharp peaks of the energy transmission factor $Q$—many of them reaching values close to the upper bound $Q=1$—are observed at various frequencies for the well configurations considered (e.g., the peaks values around $f=0.895,1.585,2.79,3.695$ and 5.525 kHz). The existence of these peaks is explained by resonance phenomena taking place inside the perforations. As illustrated by the pressure field at a peak frequency displayed in Figure 7, the factor $Q$ attains its local maxima at “resonance” frequencies, for which the associated pressure field exhibits inordinate large amplitudes inside the perforations.

The large correlation between the location of the peaks for the various well configurations observed in Figure 6, on the other hand, can be explained by the resonant frequencies of an individual perforation. In fact, large values of $Q_{j}$ are expected to occur at the resonant frequencies of the $j$-th perforation. Since the same perforation radius, the same perforation length, and the same location of the transducer are utilized in the three configurations considered, all the perforations are expected to resonate collectively at approximately the same frequency. Therefore, the factors $Q_{j}$, $j=1,\dots,N_{p}$ attain simultaneously local maxima at these “resonance” frequencies. To look into that in more detail, we present Figure 8—which displays the individual factors $Q_{j}$, $j=1,\ldots,N_{p}$—where it can be clearly observed that the factors $Q_{j}$ attain collectively local maxima at certain frequencies that indeed correspond to the largest peak values of $Q$ observed in Figure 6.

Although the simulation results presented above seemingly indicate the existence of optimal frequencies for which nearly 100% ($Q\approx 1$) of energy transmission is achieved, in practice, uncertain variations in the shape of the perforations might result in an overall reduction of the peak values of $Q$. To briefly study the effect of small shape variations on the location of the local maxima of $Q$, we consider perturbations of the three aforementioned configurations, which are generated by introducing random changes in the perforation radius, the perforation length, and the location of the transducer. Figure 9 displays the $Q$ factors obtained for the new well configurations, where it can be observed a much weaker correlation between the location of the peak values, as compared to the results presented in Figure 6. This weaker correlation is further explained by the results displayed in Figure 10, which show that, as expected, the factors $Q_{j}$, $j=1,\dots,N_{p}$ do not attain their local maxima at the same frequencies. Despite this fact, remarkably large peak values of $Q$ ($Q\approx 1$) are still observed. Nearly perfect transmission is achieved in this case by excitation of “resonant” frequencies associated with just a few perforations, for which the local energy transmission factors $Q_{j}$ lies well above 50% (e.g., the plot at the top of Figure 10—corresponding to the first well configuration—around $f=3.24$ kHz, where $Q_{1}=0.81$). Figure 11 displays the pressure field at one of the peak values of $Q$, where large pressure amplitude values inside some of the perforations are again observed. We thus finally conclude that, in principle, it would be possible to achieve nearly perfect transmission for realistic well configurations, provided the model assumptions are satisfied.

5 Concluding remarks

A mathematical model—based upon the Helmholtz equation and the use of a suitable impedance boundary condition—and a DtN-FE procedure are presented for the numerical simulation of an AWS method. The existence of optimal emission frequencies, associated with acoustic resonance phenomena, is demonstrated by means of numerical simulations for a variety of realistic well configurations. We believe that the proposed methodology and the numerical results presented in this work provide valuable information for design and optimization of the AWS method as its performance can be significantly improved by properly selecting the operating frequencies of the AWS device (transducer).

Appendix A White’s wall impedance model

Let us consider a circular cylinder of radius $r_{0}>0$ which is assumed to be filled with a liquid and surrounded everywhere by an unbounded porous material. The pressure $P$ and the average flow velocity $V$ in the radial direction are related by Darcy’s law

[TABLE]

where $\eta$ denotes the shear viscosity of the fluid, and $\kappa$ denotes the permeability of the porous material. Note that it is assumed in (20) that both the elastic expansion of the tube and the average direction of the flow through the pore space, are radial. The equation of conservation of mass

[TABLE]

together with the compressibility relation

[TABLE]

where $B$ denotes the bulk modulus of the fluid in the pore space and $\phi$ denotes the porosity, lead to

[TABLE]

Combining equations (20) and (21) we arrive at

[TABLE]

where $m=\phi\eta/(\kappa B)$. Further assuming that the velocity and pressure fields in the porous material are time-harmonic, i.e., $P(r,t)={\mathrm{Re}\,}\{p(r)\operatorname{e}^{-i\omega t}\}$ and $V(r,t)={\mathrm{Re}\,}\{v(r)\operatorname{e}^{-i\omega t}\}$, we obtain that the pressure amplitude $p$ satisfies the Bessel differential equation

[TABLE]

Looking for bounded outgoing-wave solutions at infinity fulfilling the boundary condition $p(r_{0})=p_{0}$ at the interface between the fluid and the porous material ($r=r_{0}$), we arrive at

[TABLE]

where $H^{(1)}_{0}$ denotes the Hankel function of the first kind and order zero [1]. From Darcy’s law (20), on the other hand, we obtain that the velocity amplitude $v$ is given by

[TABLE]

Combining (22) and (23), it is straightforward to evaluate the wall impedance $Z$, which is defined as the quotient of the pressure amplitude $p$ to the radial velocity amplitude $v$ at the surface of the cylinder, that is,

[TABLE]

This wall impedance, given in terms of the frequency $\omega$, accounts for the effect that the porous medium has on the fluid dynamics inside the cylinder.