About this Journal Submit a Manuscript Table of Contents
Mathematical Problems in Engineering
Volume 2011 (2011), Article ID 159389, 26 pages
http://dx.doi.org/10.1155/2011/159389
Research Article

Coupling the BEM/TBEM and the MFS for the Numerical Simulation of Wave Propagation in Heterogeneous Fluid-Solid Media

CICC, Department of Civil Engineering, University of Coimbra, Rua Luís Reis Santos, Pólo II da Universidade, 3030-788 Coimbra, Portugal

Received 31 March 2011; Accepted 24 July 2011

Academic Editor: Luis Godinho

Copyright © 2011 António Tadeu and Igor Castro. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

This paper simulates wave propagation in an elastic medium containing elastic, fluid, rigid, and empty heterogeneities, which may be thin. It uses a coupling formulation between the boundary element method (BEM)/the traction boundary element method (TBEM) and the method of fundamental solutions (MFS). The full domain is divided into subdomains, which are handled separately by the BEM/TBEM or the MFS, to overcome the specific limitations of each of these methods. The coupling is enforced by applying the prescribed boundary conditions at all medium interfaces. The accuracy, efficiency, and stability of the proposed algorithms are verified by comparing the results with reference solutions. The paper illustrates the computational efficiency of the proposed coupling formulation by computing the CPU time and the error. The transient analysis of wave propagation in the presence of a borehole driven in a cracked medium is used to illustrate the potential of the proposed coupling formulation.

1. Introduction

Various numerical methods have been proposed to simulate the propagation of waves in elastic and acoustic media, since analytical solutions are only known for simple and regular geometries (e.g., [16]). These techniques include the thin-layer method (TLM) [7], the boundary element method (BEM) [8], the finite element method (FEM) [9, 10], the finite difference method (FDM) [11], the ray tracing technique [12], and the method of fundamental solutions (MFS) [13].

Of these techniques, the FEM is the most widely used numerical method used by researchers and commercial software producers. It can be used to solve complex geometries, but it requires the full discretization of the media being modelled. This makes the FEM computationally unfeasible for very large scale models, such as those involving unbounded domains, unless substantial shortcuts are implemented. These may entail the use of coarse elements, low frequency simulations, or the introduction of boundary artefacts.

The BEM is one of the most suitable techniques for modelling wave propagation in homogeneous unbounded systems containing irregular interfaces and inclusions, because only the boundaries of the heterogeneities and interfaces need to be discretized and the far-field conditions are automatically satisfied [1416]. Despite this, the BEM still needs prior knowledge of fundamental solutions (Green’s functions) and also requires the correct integration of the resulting singular and hypersingular integrals to guarantee its efficiency. In addition, the number of boundary elements depends on the excitation frequency, and many boundary elements are needed to model high-frequency responses, a situation which leads to an undesirably high computational cost.

Furthermore, the simulation of wave propagation in the presence of very thin heterogeneities such as cracks leads to singular boundary element matrix systems, thus leading to the mathematical degeneration of the numerical formulation [17]. The dual boundary element method (DBEM) is one of the main boundary element formulations adopted to overcome this problem. Derivatives of the original BEM displacement formulation to produce a traction formulation first became necessary when fracture mechanics problems began to be addressed [18]. But these hybrid BEM formulations do not necessarily have to be used for solving such problems. Good results have been obtained in 2D examples of both elastodynamic and coupled-field problems involving stationary cracks when conventional, displacement-based BEM formulations were used in a transformed domain, with special treatment of the cracks [19, 20].

Using the DBEM, after the discretization of the inclusion’s surface, dipole loads are applied to the opposite surface, which is governed by the traction boundary integral equation [21], while monopole loads are applied to one part of the surface, which corresponds to applying the displacement boundary integral equation. In the case of a dimensionless empty crack, only a single line of boundary elements loaded with dipole loads is used to solve the problem, that is, by using only the traction boundary integral equation method [2224]. The appearance of hypersingular integrals is one of the difficulties posed by these formulations. In the particular case of 2D and 2.5D wave propagation in elastic and acoustic media, the resulting hypersingular kernels can be computed analytically [25].

Meshless techniques that require neither domain nor boundary discretization have recently become popular [26, 27]. The origin of the MFS has two sources and lies in the indirect BEM [28] and the general definition of a Green’s function [29]. The MFS copes with some of the mathematical complexity of the BEM and provides acceptable solutions for wave propagation problems at substantially lower computational cost [30, 31]. The MFS solution is based on a linear combination of fundamental solutions (Green’s functions), generated by a set of fictitious sources to simulate the scattered and refracted field produced by the heterogeneities. To avoid singularities, these virtual sources are placed at some distance from the inclusion’s boundary. The use of fundamental solutions allows the final solution to verify the unbounded boundary conditions automatically. Still the use of the MFS has its own limitations when thin inclusions such as cracks and inclusions with twisting (sinuous) boundaries are involved. The analysis would require the use of domain decompositions or/and the use of enriched functions, which increases computation costs [32]. The number of the virtual sources and their positions is another difficulty since the results are highly dependent on these parameters. Among the strategies that have been proposed to handle this problem is the verification of the solution’s accuracy by computing the solution at points other than the collocation ones, where the boundary and prescribed conditions are known a priori.

Researchers are currently trying to improve the results by coupling different methods so as to exploit the advantages of each technique and reduce their disadvantages, thereby speeding up analysis and ensuring efficiency, stability, accuracy, and flexibility.

BEM/FEM coupling has often been used, with each technique being applied to distinct subdomains [3335]. The two approaches most often used are a direct coupling and iterative coupling [3638]. Iterative coupling allows the subdomains to be analyzed separately, leading to smaller and better-conditioned systems of equations with independent discretizations being considered for each subdomain. Some authors have reported problems related to the convergence of ill-posed models, however. The coupling of meshless methods and the BEM is another approach. The coupling of the BEM/TBEM with the MFS to analyse acoustic wave propagation in the presence of multiple inclusions and thin heterogeneities is one example proposed by the authors. The full domain is first divided into subdomains which are modelled using the BEM/TBEM and the MFS. The subdomains are then coupled by imposing the required boundary conditions [39].

The paper extends that work with a formulation which couples the BEM/TBEM and the MFS to simulate the propagation of waves involving the fluid-solid interaction, as in the case of multielastic fluid layer systems, acoustic logging, and cross-hole surveying geophysical prospecting techniques [40, 41]. It is very often quite helpful to model the direct problem in order to better understand how waves propagate in the presence of such structures, particularly in cracked media and damaged zones, when it is not always easy to interpret the recorded results because of the unforeseen presence of those heterogeneities [42, 43].

The problem is formulated in the frequency domain. The waves generated by the virtual sources used by the MFS are seen as incident waves by the BEM/TBEM, while the BEM sees the collocation points used to impose the boundary conditions at the interfaces modelled by the MFS, as receivers. The approach is implemented for 2.5D problems in general. The accuracy of the proposed coupling algorithms, which use different combinations of BEM/TBEM and MFS formulations, is checked by means of a verification analysis using reference solutions.

The proposed coupling formulations for simulating wave propagation in the presence of fluid and elastic inclusions are described in the next section. The coupling formulations are first verified against solutions obtained using BEM/TBEM, taken as reference solutions. We then show the computational efficiency of the formulations by measuring the CPU time taken to compute the numerical responses provided by the different algorithms. Finally, the applicability of the proposed method is shown by means of a numerical example that simulates the propagation of waves generated by a line blast load in the vicinity of a fluid-filled borehole driven in a cracked elastic medium.

2. Boundary Integral Coupling Formulations

Consider two irregular two-dimensional cylindrical inclusions, 1 and 2, embedded in a homogeneous elastic medium (Medium 1) with density 𝜌1 (Figure 1) and allowing longitudinal (P-wave) and shear waves (S-wave) to travel at velocities 𝛼1 and 𝛽1, respectively. Medium 2, inside Inclusion 1, is fluid, has density 𝜌2, and permits pressure waves (P-wave) to travel at velocity 𝛼2. Inclusion 2 is elastic (Medium 3), has density 𝜌3, and allows longitudinal and shear waves to travel at velocities 𝛼3 and 𝛽3, respectively.

159389.fig.001
Figure 1: The geometry of the problem.

It is further assumed that this system is subjected to a dilatational line source placed at (𝑥𝑠,𝑦𝑠) whose amplitude varies sinusoidally in the third dimension (𝑧).

The incident wave field generated by this source can be expressed in the frequency domain by means of the classic dilatational potential:𝜙inc𝑥,𝑦,𝜔,𝑘𝑧=i𝐴2𝐻0𝑘𝛼1𝑟,(2.1) in which 𝐻𝑛() are second kind Hankel functions of the order 𝑛, i=1, 𝑘𝛼1=𝜔2/𝛼21𝑘2𝑧 with Im(𝑘𝛼1)<0, 𝑟=(𝑥𝑥𝑠)2+(𝑦𝑦𝑠)2 and 𝑘𝑧 is the wavenumber along 𝑧.

Then, the displacement field can be expressed as𝑢inc𝑥𝑥,𝑦,𝑥𝑠,𝑦𝑠=,𝜔i𝐴2𝑘𝛼1𝐻1𝑘𝛼1𝑟𝜕𝑟,𝑢𝜕𝑥inc𝑦𝑥,𝑦,𝑥𝑠,𝑦𝑠=,𝜔i𝐴2𝑘𝛼1𝐻1𝑘𝛼1𝑟𝜕𝑟,𝑢𝜕𝑦inc𝑧𝑥,𝑦,𝑥𝑠,𝑦𝑠=𝐴,𝜔2𝑘𝑧𝐻0𝑘𝛼1𝑟.(2.2)

2.1. BEM/MFS Coupling Formulation

This section describes the coupling between the BEM and the MFS formulations used to obtain the wave field generated by the dilatational line source placed in the exterior medium, Medium 1. The first inclusion is modelled using the BEM while the other is solved with the MFS (see Figure 2).

159389.fig.002
Figure 2: Discretization of the system: position of virtual loads, collocation points and boundary elements.

Continuity of normal tractions and displacements and null tangential tractions are prescribed along the boundary of the fluid Inclusion 1. Three different boundary conditions may be ascribed to Inclusion 2′s surface: continuity of displacements and tractions (simulating an elastic inclusion); null tractions (an empty inclusion); null displacements (a rigid inclusion).

2.1.1. Fluid Inclusion 1 and Elastic Inclusion 2

Assuming Inclusion 1 to be bounded by a surface 𝑆1 and subjected to an incident field given by 𝑢inc, the boundary integral equation can be constructed by applying the reciprocity theorem (e.g., Manolis and Beskos [44]) leading to the following.

(a) Along the Exterior Domain of Inclusion 1 (Medium 1),
𝑐𝑖𝑗𝑢𝑖(1)𝑥0,𝑦0=,𝜔𝑆1𝑡1(1)𝑥,𝑦,𝐧𝑛1𝐺,𝜔(1)𝑖1𝑥,𝑦,𝑥0,𝑦0,𝜔d𝑠3𝑗=1𝑆1𝑢𝑗(1)(𝑥,𝑦,𝜔)𝐻(1)𝑖𝑗𝑥,𝑦,𝐧𝑛1,𝑥0,𝑦0,𝜔𝑑𝑠+𝑢inc𝑖𝑥0,𝑦0,𝑥𝑠,𝑦𝑠.,𝜔(2.3) In this equation, 𝑖,𝑗=1,2 correspond to the normal and tangential directions relative to the inclusion surface, while 𝑖,𝑗=3 correspond to the 𝑧 direction. In these equations, the superscript 1 represents the exterior domain; 𝐧𝑛1 is the unit outward normal along boundary 𝑆1, at (𝑥,𝑦), defined by the vector 𝐧𝑛1=(cos𝜃𝑛1,sin𝜃𝑛1). 𝐺(1)𝑖𝑗(𝑥,𝑦,𝑥0,𝑦0,𝜔) and 𝐻(1)𝑖𝑗(𝑥,𝑦,𝐧𝑛1,𝑥0,𝑦0,𝜔) define the fundamental solutions for displacements and tractions (Green’s functions), in direction 𝑗 on the boundary 𝑆1 at (𝑥,𝑦), caused by a unit point force in direction 𝑖 applied at the nodal point, (𝑥0,𝑦0) (see the appendix). 𝑢𝑗(1)(𝑥,𝑦,𝜔) corresponds to displacements in direction 𝑗 at (𝑥,𝑦), 𝑡𝑗(1)(𝑥,𝑦,𝐧𝑛1,𝜔) specifies the nodal tractions in direction 𝑗 on the boundary at (𝑥,𝑦) and 𝑢inc𝑖(𝑥0,𝑦0,𝑥𝑠,𝑦𝑠,𝜔) to the displacement incident field at (𝑥0,𝑦0) along direction 𝑖, when the source is located at (𝑥𝑠,𝑦𝑠). The coefficient 𝑐𝑖𝑗 is equal to 𝛿𝑖𝑗/2, with 𝛿𝑖𝑗 being the Kronecker delta, when the boundary is smooth.
Green’s functions for displacements along the 𝑥, 𝑦, and 𝑧 directions in the solid medium are listed in the appendix, and their derivation can be found in [45].
Equation (2.3) does not yet take into account the presence of the neighbouring Inclusion 2, which is modelled using the MFS. The MFS assumes that the response of this neighbouring inclusion is found as a linear combination of fundamental solutions simulating the displacement field generated by two sets of 𝑁𝑆 virtual sources. These virtual loads are distributed along the inclusion interface 𝑆2 at distances 𝛿 from that boundary towards the interior and exterior of the inclusion (lines 𝐶(1) and 𝐶(2) in Figure 2) in order to avoid singularities. Sources inside the inclusion have unknown amplitudes 𝑎(2)𝑛𝑗,𝑛_ext, while those placed outside the inclusion have unknown amplitudes 𝑎(2)𝑛𝑗,𝑛_int. In the exterior, and interior elastic media the scattered displacement fields are given by 𝑢𝑖(1)(𝑥,𝑦,𝜔)=𝑁𝑆3𝑛=1𝑗=1𝑎(2)𝑛𝑗,𝑛_ext𝐺(1)𝑗𝑖𝑥,𝑦,𝑥𝑛_ext,𝑦𝑛_ext,𝑢,𝜔𝑖(3)(𝑥,𝑦,𝜔)=𝑁𝑆3𝑛=1𝑗=1𝑎(2)𝑛𝑗,𝑛_int𝐺(3)𝑗𝑖𝑥,𝑦,𝑥𝑛_int,𝑦𝑛_int,,𝜔(2.4) where 𝐺(1)𝑗𝑖(𝑥,𝑦,𝑥𝑛_ext,𝑦𝑛_ext,𝜔) and 𝐺(3)𝑗𝑖(𝑥,𝑦,𝑥𝑛_int,𝑦𝑛_int,𝜔) are the fundamental solutions which represent the displacements at points (𝑥,𝑦) in Mediums 1 and 3, in direction 𝑖, caused by a unit point force in direction 𝑗 applied at the positions (𝑥𝑛_ext,𝑦𝑛_ext) and (𝑥𝑛_int,𝑦𝑛_int).  𝑛_ext and 𝑛_int are the subscripts that denote the load order number placed along lines 𝐶(1) and 𝐶(2).
The displacement field generated by this second inclusion can be viewed as an incident field that strikes the first inclusion. So (2.3) needs to be modified accordingly, 𝑐𝑖𝑗𝑢𝑖(1)𝑥0,𝑦0=,𝜔𝑆1𝑡1(1)𝑥,𝑦,𝐧𝑛1𝐺,𝜔(1)𝑖1𝑥,𝑦,𝑥0,𝑦0,𝜔d𝑠3𝑗=1𝑆1𝑢𝑗(1)(𝑥,𝑦,𝜔)𝐻(1)𝑖𝑗𝑥,𝑦,𝐧𝑛1,𝑥0,𝑦0,𝜔𝑑𝑠+𝑢inc𝑖𝑥0,𝑦0,𝑥𝑠,𝑦𝑠+,𝜔𝑁𝑆3𝑛=1𝑗=1𝑎(2)𝑛𝑗,𝑛_ext𝐺(1)𝑗𝑖𝑥0,𝑦0,𝑥𝑛_ext,𝑦𝑛_ext.,𝜔(2.5)

(b) In the Interior Domain of Inclusion 1 (Medium 2),
𝑐𝑝(2)𝑥0,𝑦0=,𝜔𝑆1𝑞(2)𝑥,𝑦,𝐧𝑛1𝐺,𝜔𝑓(2)𝑥,𝑦,𝑥0,𝑦0,𝜔d𝑠𝑆1𝑝(2)(𝑥,𝑦,𝜔)𝐻𝑓(2)𝑥,𝑦,𝐧𝑛1,𝑥0,𝑦0,𝜔d𝑠.(2.6) In (2.6), the superscript 2 corresponds to the domain inside Inclusion 1, Medium 2. 𝐺𝑓(2)(𝑥,𝑦,𝑥0,𝑦0,𝜔) and 𝐻𝑓(2)(𝑥,𝑦,𝐧𝑛1,𝑥0,𝑦0,𝜔) are Green’s functions for pressure and the gradient of pressure on the boundary 𝑆1 at (𝑥,𝑦), caused by a unit point pressure at the nodal point, (𝑥0,𝑦0) (see the appendix). 𝑝(2)(𝑥,𝑦,𝜔) corresponds to the pressure at (𝑥,𝑦), 𝑞(2)(𝑥,𝑦,𝐧𝑛1,𝜔) specifies the nodal pressure gradients on the boundary at (𝑥,𝑦). The coefficient 𝑐 is equal to 1/2 when the boundary is smooth.

(c) In the Interior and Exterior Domains of Inclusion 2 (Media 1 and 3),
To determine the amplitudes of the unknown virtual point loads 𝑎(2)𝑛𝑗,𝑛_ext and 𝑎(2)𝑛𝑗,𝑛_int, it is also necessary to impose the continuity of displacements and tractions at interface 𝑆2, which is the boundary of Inclusion 2, along 𝑁𝑆 collocation points (𝑥col,𝑦col). This must be done so as to take into account the scattered field generated at Inclusion 1. The following two equations are thus defined: 𝑆1𝑡1(1)𝑥,𝑦,𝐧𝑛1𝐺,𝜔(1)𝑖1𝑥,𝑦,𝑥col,𝑦col,𝜔d𝑠3𝑗=1𝑆1𝑢𝑗(1)(𝑥,𝑦,𝜔)𝐻(1)𝑖𝑗𝑥,𝑦,𝐧𝑛1,𝑥col,𝑦col,𝜔d𝑠+𝑢inc𝑖xcol,𝑦col,𝑥𝑠,𝑦𝑠+,𝜔𝑁𝑆3𝑛=1𝑗=1𝑎(2)𝑛𝑗,𝑛_ext𝐺(1)𝑗𝑖𝑥col,𝑦col,𝑥𝑛_ext,𝑦𝑛_ext=,𝜔𝑁𝑆3𝑛=1𝑗=1𝑎(2)𝑛𝑗,𝑛_𝑖𝑛𝑡𝐺(3)𝑗𝑖𝑥col,𝑦col,𝑥𝑛_int,𝑦𝑛_int,,𝜔(2.7)𝑆1𝑡1(1)𝑥,𝑦,𝐧𝑛1,𝜔𝐺(1)𝑖1𝑥,𝑦,𝐧𝑛2,𝑥col,𝑦col,𝜔d𝑠3𝑗=1𝑆1𝑢𝑗(1)(𝑥,𝑦,𝜔)𝐻(1)𝑖𝑗𝑥,𝑦,𝐧𝑛1,𝐧𝑛2,𝑥col,𝑦col+,𝜔d𝑠𝑢inc𝑖xcol,𝑦col,𝐧𝑛2,𝑥𝑠,𝑦𝑠+,𝜔𝑁𝑆3𝑛=1𝑗=1𝑎(2)𝑛𝑗,n_ext𝐺(1)𝑗𝑖𝑥col,𝑦col,𝐧𝑛2,𝑥𝑛_ext,𝑦𝑛_ext=,𝜔𝑁𝑆3𝑛=1𝑗=1𝑎(2)𝑛𝑗,𝑛_int𝐺(3)𝑗𝑖𝑥col,𝑦col,𝐧𝑛2,𝑥𝑛_int,𝑦𝑛_int.,𝜔(2.8)
Green’s functions 𝐻(1)𝑖𝑗(𝑥,𝑦,𝐧𝑛1,𝐧𝑛2,𝑥col,𝑦col,𝜔) and 𝐺(1)𝑖𝑗(𝑥,𝑦,𝐧𝑛2,𝑥col,𝑦col,𝜔) are defined by applying the traction operator to 𝐻(1)𝑖𝑗(𝑥,𝑦,𝐧𝑛1,𝑥col,𝑦col,𝜔)and 𝐺(1)𝑖𝑗(𝑥,𝑦,𝑥col,𝑦col,𝜔), which can be obtained by combining the derivatives of the former Green’s functions, in order of 𝑥, 𝑦, and 𝑧, so as to obtain the stresses (see the appendix). In these equations, 𝐧𝑛2 is the unit outward normal to the boundary 𝑆2 at the collocation points (𝑥col,𝑦col).

(d) Final System of Equations.
The global solution is obtained by solving (2.5)–(2.8). This requires the discretization of the interface 𝑆1, which is the boundary of Inclusion 1. For the purposes of this work, this interface is discretized into 𝑁 straight boundary elements, with one nodal point in the middle of each element (see Figure 2).
The required integrations in (2.5)–(2.8) are evaluated in closed form when the element to be integrated is the loaded element [45, 46], while a numerical integration that uses a Gaussian quadrature scheme applies when the element to be integrated is not the loaded one.

The final integral equations are manipulated and combined so as to impose the continuity of normal tractions and displacements, and null tangential tractions along the boundary of Inclusion 1, and the continuity of displacements and tractions along the boundary of Inclusion 2, to establish a system of [(6𝑁𝑆+4𝑁)×(6𝑁𝑆+4𝑁)] equations. The relation 𝑢1(1)=(1/𝜌2)(𝜕𝑝(2)/𝜕𝐧𝑛1) is used to relate pressure gradients and displacements, while the normal pressure corresponds to normal tractions.

The solution of this system of equations gives the nodal tractions and displacements along the boundary 𝑆1 and the unknown virtual load amplitudes, 𝑎(2)𝑛𝑗,𝑛_ext and 𝑎(2)𝑛𝑗,𝑛_int, which allow the displacement field to be defined inside and outside the inclusions.

2.1.2. Empty Inclusion 2 (Null Tractions Along its Boundary)

In this case, null tractions are prescribed along the boundary 𝑆2. Thus, (2.5) and (2.6) are kept as before and (2.8) is simplified to 𝑆1𝑡1(1)𝑥,𝑦,𝐧𝑛1,𝜔𝐺(1)𝑖1𝑥,𝑦,𝐧𝑛2,𝑥col,𝑦col,𝜔d𝑠3𝑗=1𝑆1𝑢𝑗(1)(𝑥,𝑦,𝜔)𝐻(1)𝑖𝑗𝑥,𝑦,𝐧𝑛1,𝐧𝑛2,𝑥col,𝑦col+,𝜔d𝑠𝑢inc𝑖𝑥col,𝑦col,𝐧𝑛2,𝑥𝑠,𝑦𝑠+,𝜔𝑁𝑆3𝑛=1𝑗=1𝑎(2)𝑛𝑗,𝑛_ext𝐺(1)𝑗𝑖𝑥col,𝑦col,𝐧𝑛2,𝑥𝑛_ext,𝑦𝑛_ext,𝜔=0.(2.9)

The solution of the boundary integral along the surface (𝑆1) again requires its discretization into 𝑁 straight boundary elements, while the simulation of Inclusion 2 uses 𝑁𝑆 collocation points/virtual sources, following a procedure similar to the one described above. This leads to a system of [(3𝑁𝑆+4𝑁)×(3𝑁𝑆+4𝑁)] equations.

2.1.3. Rigid Inclusion 2 (Null Displacements Along its Boundary)

Null displacements on the surface of Inclusion 2 are now prescribed, which leads to (2.5) and (2.6) and to the following equation:𝑆1𝑡1(1)𝑥,𝑦,𝐧𝑛1𝐺,𝜔(1)𝑖1𝑥,𝑦,𝑥col,𝑦col,𝜔d𝑠3𝑗=1𝑆1𝑢𝑗(1)(𝑥,𝑦,𝜔)𝐻(1)𝑖𝑗𝑥,𝑦,𝐧𝑛1,𝑥col,𝑦col,𝜔d𝑠+𝑢inc𝑖𝑥col,𝑦col,𝑥𝑠,𝑦𝑠+,𝜔𝑁𝑆3𝑛=1𝑗=1𝑎(2)𝑛𝑗,𝑛_ext𝐺(1)𝑗𝑖𝑥col,𝑦col,𝑥𝑛_ext,𝑦𝑛_ext,𝜔=0.(2.10)

The solution of these equations is once again obtained as described above, with a system of [(3𝑁𝑆+4𝑁)×(3𝑁𝑆+4𝑁)] equations. Other coupling combinations can be solved in the same way.

2.2. TBEM/MFS Coupling Formulation

The traction boundary element method (TBEM) can be proposed to simulate the scattered wave field in the vicinity of thin inclusions, for which the BEM formulation described above fails [47, 48]. This technique can be formulated by applying dipoles (dynamic doublets) instead of monopole loads. Replace the former (2.5) and (2.6), to give the following (2.11), while modelling the first inclusion:𝑎𝑖𝑗𝑢𝑖(1)𝑥0,𝑦0,𝜔+𝑐𝑖1𝑡1(1)𝑥0,𝑦0,𝐧𝑛1=,𝜔𝑆1𝑡1(1)𝑥,𝑦,𝐧𝑛1,𝜔𝐺(1)𝑖1𝑥,𝑦,𝐧𝑛2,𝑥0,𝑦0,𝜔d𝑠3𝑗=1𝑆1𝑢𝑗(1)(𝑥,𝑦,𝜔)𝐻(1)𝑖𝑗𝑥,𝑦,𝐧𝑛1,𝐧𝑛2,𝑥0,𝑦0+,𝜔d𝑠𝑢inc𝑖𝑥0,𝑦0,𝐧𝑛2,𝑥𝑠,𝑦𝑠+,𝜔𝑁𝑆3𝑛=1𝑗=1𝑎(2)𝑛𝑗,𝑛_ext𝐺(1)𝑗𝑖𝑥0,𝑦0,𝐧𝑛2,𝑥𝑛_ext,𝑦𝑛_ext,,𝜔𝑎𝑝(2)𝑥0,𝑦0,𝜔+𝑐𝑞(2)𝑥0,𝑦0,𝐧𝑛1=,𝜔𝑆1𝑞(2)𝑥,𝑦,𝐧𝑛1,𝜔𝐺𝑓(2)𝑥,𝑦,𝐧𝑛2,𝑥0,𝑦0,𝜔d𝑠𝑆1𝐻𝑓(2)𝑥,𝑦,𝐧𝑛1,𝐧𝑛2,𝑥0,𝑦0𝑝,𝜔(1)(𝑥,𝑦,𝜔)d𝑠.(2.11) As noted by Guiggiani [49], the coefficients 𝑎𝑖𝑗 and 𝑎 are zero for piecewise straight boundary elements. The factors 𝑐𝑖𝑗 and 𝑐 are constants, defined as above.

Equations (2.7)) and (2.8) can be kept the same for modelling the second inclusion.

The solutions of these equations are defined as before by discretizing the boundary surface (𝑆1) into 𝑁 straight boundary elements, with one nodal point in the middle of each element. The integrations in (2.11) are performed through a Gaussian quadrature scheme when the element being integrated is not the loaded one. When the element being integrated is the loaded one, the integrals become hypersingular. An indirect approach is used for the analytical solution of those hypersingular integrals. This consists of defining the dynamic equilibrium of an isolated semicylinder, above each boundary element (see [47, 48]).

Manipulating (2.7), (2.8) and (2.11) as described above, the cavity and the rigid inclusions and their combinations, can also be modelled.

2.3. Combined (TBEM+BEM)/MFS Coupling Formulation

The displacement and traction formulations can be combined so as to solve the problems described above. This has the advantage of allowing the solution to be defined when Inclusions 1 or/and 2 are thin. In these cases, part of the boundary surface of the inclusion is loaded with monopole loads (formulation in displacements), while the remaining part is loaded with dipoles (formulation in tractions).

3. Verification of the Coupling Algorithms

The proposed coupling algorithms (MFS/BEM and MFS/TBEM) described are verified against BEM and MFS solutions by solving the elastic field produced by two circular inclusions embedded in an unbounded elastic medium, centred at (0.0m,20.0m) and (22.0m,5.0m), with radii of 5.0m and 6.0m (see Figure 3). Three separate problems are solved by combining different types of Inclusion 2, namely, an elastic inclusion (Case  1), a cavity (Case  2), and a rigid inclusion (Case  3). Inclusion 1 is always filled with fluid.

159389.fig.003
Figure 3: Two circular inclusions (fluid and elastic) embedded in an unbounded elastic medium.

The host elastic medium (𝜌1=2140 kg/m3) is homogeneous, permitting a P-wave velocity of 𝛼1=4208 m/s and an S-wave velocity of 𝛽1=2656 m/s. The fluid in Inclusion 1 exhibits a density equal to 𝜌2=1000 kg/m3 and allows P-wave velocity of 𝛼2=1500 m/s, while the elastic Inclusion 2, of density 𝜌3=2250 kg/m3, exhibits a dilatational and an S-wave of velocities 𝛼3=2630 m/s, 𝛽3=1416 m/s, respectively. The excitation source is assumed to be inside the host medium at point (10.0m,17.0m). It is a harmonic dilatational line load whose amplitude varies sinusoidally in the third dimension according to 𝑘𝑧=0.2 rad/m.

The responses are computed at receiver 𝑅1 and 𝑅2, placed at (15.0m,10.0m) and (1.0m,19.0m), respectively. The computations are performed in the frequency domain from 1 Hz to 200 Hz.

All the illustrated simulations used interior and exterior virtual sources respectively placed at distances 0.9×𝑟 and 1.1×𝑟 from the centre of the inclusion, with 𝑟 being radii.

Figures 4, 5, and 6 present the real (solid line) and imaginary (dashed line) parts of the displacements 𝑢𝑥, 𝑢𝑦, and 𝑢𝑧 (receiver 𝑅1) and pressure response (receiver 𝑅2) for the three cases. The lines correspond to the BEM responses, that is, when the inclusions are each modelled with 200 boundary elements. Different BEM/TBEM, MFS, and coupling solutions are indicated by the marked points and labelled “BEM/TBEM”, “MFS”, “MFS/BEM”, and “MFS/TBEM”. 200 boundary elements and virtual sources are used in the MFS and coupling solutions for each inclusion. An analysis of the results shows very good agreement between the proposed coupling solutions and both the BEM and MFS models’ solutions.

fig4
Figure 4: Case  1: BEM, TBEM, MFS, and coupling formulations’ results: displacements at receiver (𝑅1) and pressure at receiver (𝑅2) when the system is excited by a blast load.
fig5
Figure 5: Case 2: BEM, TBEM, MFS, and coupling formulations’ results: displacements at receiver (𝑅1) and pressure at receiver (𝑅2) when the system is excited by a blast load.
fig6
Figure 6: Case 3: BEM, TBEM, MFS, and coupling formulations’ results: displacements at receiver (𝑅1) and pressure at receiver (𝑅2) when the system is excited by a blast load: BEM.

4. Computational Efficiency of the Coupling Algorithms

The computational efficiency of the proposed coupling formulations is illustrated by calculating at a grid of receivers the displacements caused by an empty crack of null thickness, placed in the vicinity of a fluid-filled borehole (Figure 7(a)).

fig7
Figure 7: Numerical application used to illustrate the computational efficiency of the proposed algorithm: (a) geometry of a fluid-filled borehole with a null-thickness empty crack in its vicinity and position of the blast load; (b) boundary elements used by the BEM/TBEM model; (c) position of virtual loads and collocation points used by the MFS model; (d) position of virtual loads, collocation points (MFS), and boundary elements (TBEM) used by the proposed MFS/TBEM coupling formulation.

The host medium, with a density of 2250 kg/m3, allows P-wave and S-wave velocities of 2630 m/s and 1416 m/s, respectively. The fluid-filled borehole is centred at (0.0m,0.0m) with a radius of 0.05m. Its fluid has a mass density of 1000 kg/m3 and permits a P-wave speed of 1500 m/s. A null-thickness arc-shaped crack centred at (0.0m,0.0m) has a radius of 0.10m and a length of (1.2𝜋/32)m.

This system is illuminated by a wave field generated by a dilatational line load placed 0.05m (𝑘𝑧=0 rad/m) from the crack at (0.15m,0.0m). The resulting displacement is obtained over a grid of 10140 receivers arranged along the 𝑥 and 𝑦 directions at equal intervals and placed from 𝑥=0.10m to 𝑥=0.25m and from 𝑦=0.15m to 𝑦=0.15m.

Computational efficiency was evaluated by determining the CPU time taken to compute the solution for the full grid of receivers by the BEM/TBEM, the MFS and the MFS/TBEM, at two specific frequencies: 140Hz and 9000 Hz.

As there are no known analytical solutions, the BEM/TBEM solution for 770 boundary elements is used as reference solution. The crack is discretized as an open line and loaded with dipole loads (210 TBEM boundary elements), while the fluid-filled borehole boundary is discretized using a classical closed surface and loaded with monopole loads (560 BEM boundary elements) (see Figure 7(b)).

Figures 8 and 9 illustrate the real and imaginary part of the reference solutions for both excitation frequencies.

fig8
Figure 8: Displacements and pressures solutions for the excitation frequency 140 Hz.
fig9
Figure 9: Displacements and pressures solutions for the excitation frequency 9000 Hz.

The MFS is less efficient at modelling thin inclusions such as cracks when good accuracy is required. The approach used here to model the displacement around the crack is based on the decomposition of the inner domain into two different subdomains, as illustrated in Figure 7(c). The interface between these two subdomains will be circular and contain the crack, 𝑇, and a fictitious interface, 𝐹. In order to correctly describe the behaviour of the null-thickness crack, null tractions are ascribed to both sides of interface 𝑇 and continuity of displacements and tractions is imposed along the interface 𝐹. The distances between the virtual sources and the boundary have been defined by computing the errors along the boundary outside the collocation points, where prescribed conditions are known. The error along the boundary is computed as the integral of the error surface, which is defined by the difference between the responses and the prescribed conditions along the boundary. The final positions of the virtual sources were those that led to the slowest boundary errors. This is the same as the procedure proposed in [31], where a stability analysis is presented.

The MFS/TBEM coupling model discretizes the crack with boundary elements loaded with dipole loads (TBEM), while the fluid-filled borehole is modelled using a set of virtual point sources (MFS), whose positions are defined as explained above. The collocation points are evenly distributed along the wall surfaces.

The errors yielded by the methods within the domain are assessed by comparing the responses obtained with those provided by the reference solution, the BEM/TBEM solution, found using 770 boundary elements. A global domain error is defined by computing the integration of the volume generated by the absolute value of the difference between the reference and the different model responses at the grid of receivers. To evaluate the computational efficiency, the CPU time taken by the three computational models to compute the solutions at the grid of receivers placed in the exterior medium (displacements) and at the grid of receivers placed within the borehole (pressures) is registered. All solutions were computed on a laptop computer with an Intel Core Duo CPU E6750.

Figure 10 illustrates the global domain error registered versus CPU time required by each formulation, for the two frequencies computed above and varying the number of degrees of freedom, that is, changing the number of boundary elements and virtual sources/collocation points. For each formulation, the number of degrees of freedom varies according to the value of 𝑚=1 to 20, as follows: the BEM/TBEM solutions were computed by discretizing the borehole and the crack interfaces with 10𝑚 and 4𝑚 boundary elements, respectively; the MFS solutions were obtained by simulating the borehole and the crack interface with 10𝑚 and 20𝑚 virtual sources/collocation points, respectively; the coupling MFS/TBEM solutions were obtained using 10𝑚 and 4𝑚 virtual sources/collocation points.

fig10
Figure 10: Global domain error versus CPU time for solving the system composed of a fluid-filled borehole placed in the vicinity of a crack.

The global domain errors shown in Figure 10 are displayed in a logarithmic scale to allow an easier interpretation of the results. An analysis of the responses shows that the BEM/TBEM and the MFS/TBEM register smaller errors as the number of degrees of freedom increases. The MFS does not exhibit a permanent trend and its behaviour fluctuates, particularly for larger numbers of virtual sources, since the global equation system may become ill-conditioned. The results show that the coupled MFS/TBEM formulation is the algorithm that requires the least CPU time for the same accuracy. In both cases, for the same CPU time, the coupled MFS/TBEM solution has the smallest global domain error, except when a very small number of degrees of freedom are used.

5. Numerical Application

The applicability of the proposed coupling formulations for solving more complex systems is illustrated by calculating the wave field in the vicinity of a fluid-filled circular borehole, with a radius of 0.1016m, driven in a cracked medium, as illustrated in Figure 11. The system is subjected to a dilatational line source pulse, modelled as a Ricker wavelet placed at (0.1m,0.3m), parallel to the borehole axis (two-dimensional application), with a characteristic frequency of 30000 Hz and which starts acting at 𝑡=0 ms. A set of snapshots taken from computer animations is presented in Figure 12 to illustrate the resulting wave field at different time instants.

159389.fig.0011
Figure 11: Geometry of a fluid-filled borehole driven in cracked medium: position of virtual loads, collocation points (MFS), and boundary elements (TBEM) used by the proposed MFS/TBEM coupling formulation.
fig12
Figure 12: Snapshots illustrating the displacement, 𝑢𝑥 and 𝑢𝑦, and the pressure generated by a line blast load, modelled as a Ricker pulse with a characteristic frequency of 30000 Hz.

The responses in the time domain are computed by applying an inverse (fast) Fourier transform to the responses in the frequency domain 𝜔. The Ricker pulse modelled is expressed in the time domain by𝑢(𝜏)=𝐴12𝜏2e𝜏2,(5.1) where 𝐴 represents the amplitude; 𝜏=(𝑡𝑡𝑠)/𝑡0, 𝑡 corresponds to the time, 𝑡𝑠 is the time when the wavelet takes its maximum value, and 𝜋𝑡0 is the characteristic (dominant) period of the Ricker wavelet. In the frequency domain, this pulse is written as 𝑈(𝜔)=𝐴2𝑡0𝜋ei𝜔𝑡𝑠Ω2eΩ2,(5.2) where Ω=𝜔𝑡0/2.

The Fourier transformation is computed by adding together a finite number of terms. The frequency increment, Δ𝜔, needs to be small enough to avoid the aliasing phenomena. These are almost completely eliminated by the introduction of complex frequencies with a small imaginary part of the form 𝜔𝑐=𝜔i𝜂 (with 𝜂=0.7Δ𝜔). This procedure is later taken into account by rescaling the responses in the time domain with an exponential factor e𝜂𝑡.

The computations are performed in the frequency domain for frequencies ranging from 140Hz to 71680Hz, with a frequency increment of 140Hz, which determines a total time window of 7.14ms.

The results were computed using the MFS/TBEM coupling model. The empty crack is discretized using a number of boundary elements defined by the relation between the wavelength and the length of the boundary elements, which was set at 10. A minimum of 10 boundary elements were used. The inclusion is simulated by the MFS, using a minimum of 40 virtual loads/collocation points. The number of virtual sources/collocation points increases with the frequency, according to the relation, defined above, between the wavelength and the distance between collocation points. Figure 11 illustrates the position of the virtual sources, collocation points, and boundary elements.

The P-wave and S-wave velocities allowed in the host medium and its density remain constant at 4208m/s, 2656m/s, and 2140kg/m3, respectively. The fluid-filled borehole is centered at (0.0m,0.0m) with a radius 0.1016m. Its medium has a mass density of 1000kg/m3, a P-wave velocity of 1500m/s. A null-thickness crack is embedded in the vicinity of this elastic inclusion.

The resulting displacement (in elastic medium) and pressure (into the fluid-filled borehole) are obtained over a two-dimensional grid of 10120 receivers arranged along the 𝑥 and 𝑦 directions at equal intervals and placed in the vicinity of the inclusion and crack from 𝑥=0.4m to 𝑥=0.4m and from 𝑦=0.4m to 𝑦=0.4m.

A set of snapshots taken from computer animations is presented in Figure 12 to illustrate the resulting wave field in both the fluid inside the borehole and in the vicinity of the crack, at different time instants. In this figure, the left and the centre columns present the horizontal displacements (𝑢𝑥) and the vertical displacements (𝑢𝑦) in the elastic medium, while the right column exhibits the pressure inside the borehole. In these plots, a colour gradient between the red and the blue represents responses from positive to negative values.

In the first plots, at 𝑡=0.01ms, the pulse excited by the dilatational source can be seen travelling in the elastic medium without perturbations as it has not yet reached the fluid-filled borehole. The fluid inside the borehole has not yet suffered any pressure variation. The differences of the component displacements in the horizontal and vertical direction can be seen clearly.

At 𝑡=0.04ms, the incident pulses are partly reflected back as P-waves and S-waves after hitting the crack, propagating away from the crack on the right in the elastic medium and creating a shadow zone behind it. This is already perceptible at 𝑡=0.07ms, but it is not easy to distinguish the two types of waves since they are almost coincident at this early stage. At 𝑡=0.10ms, the reflected P- and S-waves are very well developed as they spread away from the crack. At 𝑡=0.04ms, part of incident pulse that has just hit the borehole and been transmitted as P-waves into the fluid can also be seen in the pressure field generated within the borehole. Also note the diffracted wave field moving around the crack, once the incident pulses reach its ends. These waves generate refracted waves that travel along both sides of the crack as guided waves.

By 𝑡=0.07ms, the waves have hit the first crack, which is on the left side of the borehole. The waves that pass through the borehole fluid (pressure) are in their initial development stages as P-waves and denote a delay in relation to the direct incident field, because of lower P-wave speed inside the fluid.

The last snapshots (𝑡=0.10ms and 𝑡=0.15ms) show the first reflected waves continuing to propagate in the unbounded medium. Multiple reflections of waves are visible as they impinge upon the crack surfaces. The wave energy trapped between cracks and within the fluid borehole generates a complex wave field due to the multiple reflections and refractions. It can be seen that the pulses that have travelled around the exterior of the inclusion appear before those that have propagated through the fluid borehole, since waves may travel more slowly through this heterogeneity. The multiple reflected and diffracted pulses on the crack surfaces and within the fluid borehole will continue until the total energy had dissipated.

6. Conclusions

Coupled formulations between the boundary element method (BEM)/traction boundary element method (TBEM) and the method of fundamental solutions (MFS) have been developed and proposed for simulating wave propagation involving solid-fluid interaction in media containing multiple inclusions. The proposed coupling formulations overcome the limitations posed by each method individually and require less computational effort, while maintaining reasonable accuracy.

The formulations have been verified against referenced solutions. The wave field generated by fluid, rigid, free, and elastic heterogeneities embedded in an unbounded homogeneous elastic medium and subjected to waves originated by dilatational loads (blast loads) has been simulated. The results were found to closely match the behaviour of the conventional direct BEM or TBEM solutions.

The coupling formulation between the MFS and the TBEM was proposed to overcome the problems posed by thin inclusion, such as cracks. The simulation of the wave propagation in the vicinity of a fluid-filled borehole driven in a cracked medium has been presented to illustrate the stability and efficiency of the proposed coupling formulations.

Appendix

2.5D Green’s Functions for Unbounded Elastic and Fluid Media

Definitions
𝛼P-wave velocity,𝛽S-wave velocity,𝑘𝛼=𝜔2/𝛼2𝑘2𝑧, with Im(𝑘𝛼)<0,𝑘𝛽=(𝜔2/𝛽2𝑘2𝑧), with Im(𝑘𝛽)<0,𝑟=(𝑥𝑥0)2+(𝑦𝑦0)2, 𝐻𝑛𝛼=𝐻𝑛(𝑘𝛼𝑟),𝐻𝑛𝛽=𝐻𝑛(𝑘𝛽𝑟) Hankel functions,𝐵𝑛=𝑘𝑛𝛽𝐻𝑛𝛽𝑘𝑛𝛼𝐻𝑛𝛼B𝑛 functions,𝐧𝑛1=(cos𝜃𝑛1,sin𝜃𝑛1) unit outward normal at (𝑥,𝑦),𝐧𝑛2=(cos𝜃𝑛2,sin𝜃𝑛2) unit outward normal at (𝑥0,𝑦0), (the collocation point).

Solid Media Green’s Functions
One has𝐺𝑥𝑥=14i𝜌𝜔2𝜔2𝛽2𝐻0𝛽1𝑟𝐵1+𝑥𝑥0𝑟2𝐵2,𝐺𝑦𝑦=14i𝜌𝜔2𝜔2𝛽2𝐻0𝛽1𝑟𝐵1+𝑦𝑦0𝑟2𝐵2,𝐺𝑧𝑧=14i𝜌𝜔2𝜔2𝛽2𝐻0𝛽𝑘2𝑧𝐵0,𝐺𝑥𝑦=𝐺𝑦𝑥=14i𝜌𝜔2𝑥𝑥0𝑟𝑦𝑦0𝑟𝐵2,𝐺𝑥𝑧=𝐺𝑧𝑥=i𝑘𝑧14i𝜌𝜔2𝑥𝑥0𝑟𝐵1,𝐺𝑦𝑧=𝐺𝑧𝑦=i𝑘𝑧14i𝜌𝜔2𝑦𝑦0𝑟𝐵1.(A.1)
The derivatives of the above Green’s functions give the following tractions along the 𝑥, 𝑦, and 𝑧 directions, in the solid medium, 𝐻𝑟𝑥𝛼=2𝜇22𝛽2𝜕𝐺𝑟𝑥+𝛼𝜕𝑥22𝛽21𝜕𝐺𝑟𝑦+𝜕𝑦𝜕𝐺𝑟𝑧𝜕𝑧cos𝜃𝑛1+𝜇𝜕𝐺𝑟𝑦+𝜕𝑥𝜕𝐺𝑟𝑥𝜕𝑦sin𝜃𝑛1,𝐻𝑟𝑦𝛼=2𝜇22𝛽21𝜕𝐺𝑟𝑥+𝜕𝑥𝜕𝐺𝑟𝑧+𝛼𝜕𝑧22𝛽2𝜕𝐺𝑟𝑦𝜕𝑦sin𝜃𝑛1+𝜇𝜕𝐺𝑟𝑦+𝜕𝑥𝜕𝐺𝑟𝑥𝜕𝑦cos𝜃𝑛1,𝐻𝑟𝑧=𝜇𝜕𝐺𝑟𝑥+𝜕𝑧𝜕𝐺𝑟𝑧𝜕𝑥cos𝜃𝑛1+𝜇𝜕𝐺𝑟𝑦+𝜕𝑧𝜕𝐺𝑟𝑧𝜕𝑦sin𝜃𝑛1,(A.2) with 𝐧𝑛1=(cos𝜃𝑛1,sin𝜃𝑛1), 𝐻𝑟𝑡=𝐻𝑟𝑡(𝑥,𝑦,𝐧𝑛1,𝑥0,𝑦0,𝜔), 𝐺𝑟𝑡=𝐺𝑟𝑡(𝑥,𝑦,𝑥0,𝑦0,𝜔), and 𝑟, 𝑡=𝑥,𝑦,𝑧. These expressions can be combined to obtain 𝐻𝑖𝑗(𝑥,𝑦,𝐧𝑛1,𝑥0,𝑦0,𝜔) in the normal and tangential directions. In these equations, 𝜇=𝜌𝛽2.

Solid Media Traction Green’s Functions
These Green’s functions can be seen as the combination of the derivatives of the equations (A.1) and (A.2), in order along 𝑥, 𝑦, and 𝑧, so as to obtain stresses 𝐺𝑖𝑗(𝑥,𝑦,𝐧𝑛2,𝑥col,𝑦col,𝜔) and 𝐻𝑖𝑗(𝑥,𝑦,𝐧𝑛1,𝐧𝑛2,𝑥col,𝑦col,𝜔). Along the boundary element, at (𝑥,𝑦), where the unit outward normal is defined by 𝐧𝑛1=(cos𝜃𝑛1,sin𝜃𝑛1), and after the equilibrium of stresses, the following equations are expressed for 𝑥, 𝑦 and 𝑧 generated by loads also applied along 𝑥, 𝑦 and 𝑧 directions: 𝐺𝑥𝑟𝛼=2𝜇22𝛽2𝜕𝐺𝑥𝑟+𝛼𝜕𝑥22𝛽21𝜕𝐺𝑦𝑟+𝜕𝑦𝜕𝐺𝑧𝑟𝜕𝑧cos𝜃𝑛2+𝜇𝜕𝐺𝑦𝑟+𝜕𝑥𝜕𝐺𝑥𝑟𝜕𝑦sin𝜃𝑛2,𝐺𝑦𝑟𝛼=2𝜇22𝛽21𝜕𝐺𝑥𝑟+𝜕𝑥𝜕𝐺𝑧𝑟+𝛼𝜕𝑧22𝛽2𝜕𝐺𝑦𝑟𝜕𝑦sin𝜃𝑛2+𝜇𝜕𝐺𝑦𝑟+𝜕𝑥𝜕𝐺𝑥𝑟𝜕𝑦cos𝜃𝑛2,𝐺𝑧𝑟=𝜇𝜕𝐺𝑥𝑟+𝜕𝑧𝜕𝐺𝑧𝑟𝜕𝑥cos𝜃𝑛2+𝜇𝜕𝐺𝑦𝑟+𝜕𝑧𝜕𝐺𝑧𝑟𝜕𝑦sin𝜃𝑛2,𝐻𝑥𝑟𝛼=2𝜇22𝛽2𝜕𝐻𝑥𝑟+𝛼𝜕𝑥22𝛽21𝜕𝐻𝑦𝑟+𝜕𝑦𝜕𝐻𝑧𝑟𝜕𝑧cos𝜃𝑛2+𝜇𝜕𝐻𝑦𝑟+𝜕𝑥𝜕𝐻𝑥𝑟𝜕𝑦sin𝜃𝑛2,𝐻𝑦𝑟𝛼=2𝜇22𝛽21𝜕𝐻𝑥𝑟+𝜕𝑥𝜕𝐻𝑧𝑟+𝛼𝜕𝑧22𝛽2𝜕𝐻𝑦𝑟𝜕𝑦sin𝜃𝑛2+𝜇𝜕𝐻𝑦𝑟+𝜕𝑥𝜕𝐻𝑥𝑟𝜕𝑦cos𝜃𝑛2,𝐻𝑧𝑟=𝜇𝜕𝐻𝑥𝑟+𝜕𝑧𝜕𝐻𝑧𝑟𝜕𝑥cos𝜃𝑛2+𝜇𝜕𝐻𝑦𝑟+𝜕𝑧𝜕𝐻𝑧𝑟𝜕𝑦sin𝜃𝑛2(A.3) with 𝐧𝑛2=(cos𝜃𝑛2,sin𝜃𝑛2) defining the unit outward normal at (𝑥0,𝑦0) (the collocation point), 𝐺𝑡𝑟=𝐺𝑡𝑟(𝑥,𝑦,𝐧𝑛2,𝑥col,𝑦col,𝜔), 𝐺𝑡𝑟=𝐺𝑡𝑟(𝑥,𝑦,𝑥0,𝑦0,𝜔), 𝐻𝑡𝑟=𝐻𝑡𝑟(𝑥,𝑦,𝐧𝑛1,𝐧𝑛2,𝑥col,𝑦col,𝜔), 𝐻tr=𝐻tr(𝑥,𝑦,𝐧𝑛1,𝑥0,𝑦0,𝜔) and 𝑟,𝑡=𝑥,𝑦,𝑧.

Fluid Media Green’s Functions
One has𝐺𝑓i=4𝐻0𝑘𝛼𝑟,𝐻𝑓=i4𝑘𝛼𝐻0𝑘𝛼𝑟𝜕𝑟𝜕𝐧𝑛1,(A.4) where 𝐺𝑓=𝐺𝑓(𝑥,𝑦,𝑥0,𝑦0,𝜔) and 𝐻𝑓=𝐻𝑓(𝑥,𝑦,𝐧𝑛1,𝑥0,𝑦0,𝜔).

Fluid Media Traction Green’s Functions
One has𝐺𝑓=i4𝑘𝛼𝐻1𝑘𝛼𝑟𝜕𝑟𝜕𝐧𝑛2,𝐻𝑓=i4𝑘𝛼𝑘𝛼𝐻2𝑘𝛼𝑟𝜕𝑟𝜕𝑥2𝜕𝑥𝜕𝐧𝑛1+𝜕𝑟𝜕𝑥𝜕𝑟𝜕𝑦𝜕𝑦𝜕𝐧𝑛1+𝐻1𝑘𝛼𝑟𝑟𝜕𝑥𝜕𝐧𝑛1𝜕𝑥𝜕𝐧𝑛2+i4𝑘𝛼𝑘𝛼𝐻2𝑘𝛼𝑟𝜕𝑟𝜕𝑥𝜕𝑟𝜕𝑦𝜕𝑥𝜕𝐧𝑛1+𝜕𝑟𝜕𝑦2𝜕𝑦𝜕𝐧𝑛1+𝐻1𝑘𝛼𝑟𝑟𝜕𝑦𝜕𝐧𝑛1𝜕𝑦𝜕𝐧𝑛2,(A.5) where 𝐺𝑓=𝐺𝑓(𝑥,𝑦,𝐧𝑛2,𝑥0,𝑦0,𝜔) and 𝐻𝑓=𝐻𝑓(𝑥,𝑦,𝐧𝑛1,𝐧𝑛2,𝑥0,𝑦0,𝜔).

Acknowledgment

The research work presented herein was supported by the Portuguese Foundation for Science and Technology (FCT), under Grant SFRH/BD/37425/2007.

References

  1. M. D. Trifunac, “Surface motion of a semi-cylindrical alluvial valley for incident plane SH waves,” Bulletin of the Seismological Society of America, vol. 61, pp. 1755–1770, 1971.
  2. Y. H. Pao and C. C. Mow, Diffraction of Elastic Waves and Dynamic Stress Concentrations, Crane and Russak, 1973.
  3. H. L. Wong and M. D. Trifunac, “Surface motion of semi-elliptical alluvial valley for incident plane SH-waves,” Bulletin of the Seismological Society of America, vol. 64, pp. 1389–1403, 1974.
  4. V. W. Lee, “Three-dimensional diffraction of elastic waves by a spherical cavity in an elastic half-space, I: closed-form solutions,” Soil Dynamics and Earthquake Engineering, vol. 7, no. 3, pp. 149–161, 1988. View at Publisher · View at Google Scholar
  5. V. W. Lee and J. Karl, “Diffraction of SV waves by underground, circular, cylindrical cavities,” Soil Dynamics and Earthquake Engineering, vol. 11, no. 8, pp. 445–456, 1992. View at Publisher · View at Google Scholar
  6. F. J. Sánchez-Sesma and U. Iturrarán-Viveros, “Scattering and diffraction of SH waves by a finite crack: an analytical solution,” Geophysical Journal International, vol. 145, no. 3, pp. 749–758, 2001. View at Publisher · View at Google Scholar
  7. G. Kausel, “Thin-layer method: formulation in the time domain,” International Journal for Numerical Methods in Engineering, vol. 37, no. 6, pp. 927–941, 1994.
  8. M. H. Aliabadi, Ed., The Boundary Element Method: Appl. in Solids and Structures, John Wiley & Sons, New York, NY, USA, 2002.
  9. F. Ihlenburg, Finite Element Analysis of Acoustic Scattering, vol. 132 of Applied Mathematical Sciences, Springer-Verlag, New York, NY, USA, 1998. View at Publisher · View at Google Scholar
  10. L. L. Thompson, “A review of finite-element methods for time-harmonic acoustics,” Journal of the Acoustical Society of America, vol. 119, no. 3, pp. 1315–1330, 2006. View at Publisher · View at Google Scholar
  11. L. Savioja, T. Rinne, and T. Takala, “Simulation of room acoustics with a 3-D finite difference mesh,” in Proceedings of the International Computer Music Conference (ICMC'94), pp. 463–466, Aarhus, Denmark, 1994. View at Publisher · View at Google Scholar
  12. A. Kulowski, “Algorithmic representation of the ray tracing technique,” Applied Acoustics, vol. 18, no. 6, pp. 449–469, 1985. View at Publisher · View at Google Scholar
  13. C. S. Chen, A. Karageorghis, and Y. S. Smyrlis, Eds., The Method of Fundamental Solutions: A Meshless Method, Dynamic Publishers, 2008.
  14. A. A. Stamos and D. E. Beskos, “3-D seismic response analysis of long lined tunnels in half-space,” Soil Dynamics and Earthquake Engineering, vol. 15, no. 2, pp. 111–118, 1996. View at Publisher · View at Google Scholar
  15. A. J. B. Tadeu, J. M. P. António, and E. Kausel, “3D scattering of waves by a cylindrical irregular cavity of infinite length in a homogeneous elastic medium,” Computer Methods in Applied Mechanics and Engineering, vol. 191, no. 27-28, pp. 3015–3033, 2002. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  16. J. António and A. Tadeu, “3D seismic response of a limited valley via BEM using 2.5D analytical Green's functions for an infinite free-rigid layer,” Soil Dynamics and Earthquake Engineering, vol. 22, no. 8, pp. 659–673, 2002. View at Publisher · View at Google Scholar
  17. D. N. Dell'erba, M. H. Aliabadi, and D. P. Rooke, “Dual boundary element method for three-dimensional thermoelastic crack problems,” International Journal of Fracture, vol. 94, no. 1, pp. 89–101, 1998. View at Publisher · View at Google Scholar
  18. T. A. Cruse, “Fracture mechanics,” in Boundary Element Methods in Mechanics, D. E. Beskos, Ed., pp. 333–365, North Holland, Amsterdam, The Netherlands, 1987.
  19. P. S. Dineva and G. D. Manolis, “Scattering of seismic waves by cracks in multi-layered geological regions I. Mechanical model,” Soil Dynamics and Earthquake Engineering, vol. 21, no. 7, pp. 615–625, 2001. View at Publisher · View at Google Scholar
  20. P. S. Dineva and G. D. Manolis, “Scattering of seismic waves by cracks in multi-layered geological regions II. Numerical results,” Soil Dynamics and Earthquake Engineering, vol. 21, no. 7, pp. 627–641, 2001. View at Publisher · View at Google Scholar
  21. M. H. Aliabadi, “A new generation of boundary element methods in fracture mechanics,” International Journal of Fracture, vol. 86, no. 1-2, pp. 91–125, 1997. View at Publisher · View at Google Scholar
  22. K. Takakuda, “Diffraction of plane harmonic waves by cracks,” Bulletin of the Japan Society of Mechanical Engineers, vol. 26, no. 214, pp. 487–493, 1983. View at Publisher · View at Google Scholar
  23. T. A. Cruse, Boundary Element Analysis in Computational Fracture Mechanics, vol. 1 of Mechanics: Computational Mechanics, Kluwer Academic, Dodrecht, The Netherlands, 1988.
  24. J. Sládek and V. Sládek, “A boundary integral equation method for dynamic crack problems,” Engineering Fracture Mechanics, vol. 27, no. 3, pp. 269–277, 1987. View at Publisher · View at Google Scholar
  25. A. Tadeu, L. Godinho, J. António, and P. Amado Mendes, “Wave propagation in cracked elastic slabs and half-space domains-TBEM and MFS approaches,” Engineering Analysis with Boundary Elements, vol. 31, no. 10, pp. 819–835, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  26. G. Fairweather, A. Karageorghis, and P. A. Martin, “The method of fundamental solutions for scattering and radiation problems,” Engineering Analysis with Boundary Elements, vol. 27, no. 7, pp. 759–769, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  27. L. Godino, A. P. Mendes, A. Tadeu et al., “Numerical simulation of ground rotations along 2D topographical profiles under the incidence of elastic plane waves,” Bulletin of the Seismological Society of America, vol. 99, no. 2, pp. 1147–1161, 2009. View at Publisher · View at Google Scholar
  28. M. A. Jawson and G. T. Symm, Integral Equation Methods in Potential theory and Elastostatics, Academic Press, London, UK, 1977.
  29. M. D. Greenberg, Application of Green’s Functions in Science and Engineering, Prentice Hall, Englewood Cliffs, NJ, USA, 1971.
  30. L. Godinho, A. Tadeu, and N. A. Simões, “Accuracy of the MFS and BEM on the analysis of acoustic wave propagation and heat conduction problems,” in Advances in Meshless Methods, S. Jan and S> Vladimir, Eds., Tech Science Press, 2006.
  31. A. Tadeu, J. António, and L. Godinho, “Defining an accurate MFS solution for 2.5D acoustic and elastic wave propagation,” Engineering Analysis with Boundary Elements, vol. 33, no. 12, pp. 1383–1395, 2009. View at Publisher · View at Google Scholar
  32. C. J. S. Alves and V. M. A. Leitão, “Crack analysis using an enriched MFS domain decomposition technique,” Engineering Analysis with Boundary Elements, vol. 30, no. 3, pp. 160–166, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  33. D. Soares Jr., W. J. Mansur, and O. Von Estorff, “An efficient time-domain FEM/BEM coupling approach based on FEM implicit Green's functions and truncation of BEM time convolution process,” Computer Methods in Applied Mechanics and Engineering, vol. 196, no. 9–12, pp. 1816–1826, 2007. View at Publisher · View at Google Scholar
  34. A. Warszawski, D. Soares Jr., and W. J. Mansur, “A FEM-BEM coupling procedure to model the propagation of interacting acoustic-acoustic/acoustic-elastic waves through axisymmetric media,” Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 45–48, pp. 3828–3835, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  35. Z. C. He, G. R. Liu, Z. H. Zhong, G. Y. Zhang, and A. G. Cheng, “A coupled ES-FEM/BEM method for fluid-structure interaction problems,” Engineering Analysis with Boundary Elements, vol. 35, no. 1, pp. 140–147, 2011. View at Publisher · View at Google Scholar
  36. D. Soares Jr., O. Von Estorff, and W. J. Mansur, “Iterative coupling of BEM and FEM for nonlinear dynamic analyses,” Computational Mechanics, vol. 34, no. 1, pp. 67–73, 2004.
  37. D. Soares Jr., J. A. M. Carrer, and W. J. Mansur, “Non-linear elastodynamic analysis by the BEM: an approach based on the iterative coupling of the D-BEM and TD-BEM formulations,” Engineering Analysis with Boundary Elements, vol. 29, no. 8, pp. 761–774, 2005. View at Publisher · View at Google Scholar
  38. D. Soares Jr., O. Von Estorff, and W. J. Mansur, “Efficient nonlinear solid-fluid interaction analysis by an iterative BEM/FEM coupling,” International Journal for Numerical Methods in Engineering, vol. 64, pp. 1416–1431, 2005. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  39. A. Tadeu, J. António, and I. Castro, “Coupling the BEM/TBEM and the MFS for the numerical simulation of acoustic wave propagation,” Engineering Analysis with Boundary Elements, vol. 34, no. 4, pp. 405–416, 2010. View at Publisher · View at Google Scholar
  40. J. António, A. Tadeu, and P. A. Mendes, “Simulation of wave propagation in a fluid-filled borehole embedded in a cracked medium using a coupled BEM/TBEM formulation,” Bulletin of the Seismological Society of America, vol. 99, no. 6, pp. 3326–3329, 2009. View at Publisher · View at Google Scholar
  41. A. Rodríguez-Castellanos, E. Flores, F. J. Sánchez-Sesma, C. Ortiz-Alemán, M. Nava-Flores, and R. Martin, “Indirect boundary element method applied to fluid—solid interfaces,” Engineering, vol. 31, pp. 470–477, 2002.
  42. G. Dresen, S. Stanchits, and E. Rybacki, “Borehole breakout evolution through acoustic emission location analysis,” International Journal of Rock Mechanics and Mining Sciences, vol. 47, no. 3, pp. 426–435, 2010. View at Publisher · View at Google Scholar
  43. P. Qu, R. Shen, L. Fu, and Z. Wang, “Time delay effect due to pore pressure changes and existence of cleats on borehole stability in coal seam,” International Journal of Coal Geology, vol. 85, no. 2, pp. 212–218, 2011. View at Publisher · View at Google Scholar
  44. G. D. Manolis and D. E. Beskos, Boundary Element Methods in Elastodynamics, Chapman & Hall, London, UK, 1988.
  45. A. J. B. Tadeu and Kausel, “Green's functions for two-and-a-half-dimensional elastodynamic problems,” Journal of Engineering Mechanics, vol. 126, no. 10, pp. 1093–1097, 2000. View at Publisher · View at Google Scholar
  46. A. J. B. Tadeu, P. F. A. Santos, and E. Kausel, “Closed-form integration of singular terms for constant, linear and quadratic boundary elements. Part 1. SH wave propagation,” Engineering Analysis with Boundary Elements, vol. 23, no. 8, pp. 671–681, 1999. View at Publisher · View at Google Scholar
  47. A. Tadeu, P. A. Mendes, and J. António, “The simulation of 3D elastic scattering produced by thin rigid inclusions using the traction boundary element method,” Computers and Structures, vol. 84, no. 31-32, pp. 2244–2253, 2006. View at Publisher · View at Google Scholar
  48. P. A. Mendes and A. Tadeu, “Wave propagation in the presence of empty cracks in an elastic medium,” Computational Mechanics, vol. 38, no. 3, pp. 183–199, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  49. M. Guiggiani, “Formulation and numerical treatment of boundary integral equations with hypersingular kernels,” in Singular Integrals in Boundary Element Methods, V. Sladek and J. Sladek, Eds., Computational Mechanics Publications, Boston, Mass, USA, 1998.