`Mathematical Problems in EngineeringVolume 2011 (2011), Article ID 159389, 26 pageshttp://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

#### 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 (Figure 1) and allowing longitudinal (P-wave) and shear waves (S-wave) to travel at velocities and , respectively. Medium 2, inside Inclusion 1, is fluid, has density , and permits pressure waves (P-wave) to travel at velocity . Inclusion 2 is elastic (Medium 3), has density , and allows longitudinal and shear waves to travel at velocities and , respectively.

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: in which are second kind Hankel functions of the order , , with , and is the wavenumber along .

Then, the displacement field can be expressed as

##### 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).

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 and subjected to an incident field given by , 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),
In this equation, correspond to the normal and tangential directions relative to the inclusion surface, while correspond to the direction. In these equations, the superscript 1 represents the exterior domain; is the unit outward normal along boundary , at , defined by the vector . and define the fundamental solutions for displacements and tractions (Green’s functions), in direction on the boundary at , caused by a unit point force in direction applied at the nodal point, (see the appendix). corresponds to displacements in direction at , specifies the nodal tractions in direction on the boundary at and to the displacement incident field at along direction , when the source is located at . The coefficient is equal to , 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 at distances from that boundary towards the interior and exterior of the inclusion (lines and in Figure 2) in order to avoid singularities. Sources inside the inclusion have unknown amplitudes , while those placed outside the inclusion have unknown amplitudes . In the exterior, and interior elastic media the scattered displacement fields are given by where and 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 and .   and are the subscripts that denote the load order number placed along lines and .
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,

(b) In the Interior Domain of Inclusion 1 (Medium 2),
In (2.6), the superscript 2 corresponds to the domain inside Inclusion 1, Medium 2. and are Green’s functions for pressure and the gradient of pressure on the boundary at , caused by a unit point pressure at the nodal point, (see the appendix). corresponds to the pressure at , specifies the nodal pressure gradients on the boundary at . The coefficient is equal to 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 and , it is also necessary to impose the continuity of displacements and tractions at interface , which is the boundary of Inclusion 2, along collocation points . This must be done so as to take into account the scattered field generated at Inclusion 1. The following two equations are thus defined:
Green’s functions and are defined by applying the traction operator to and , 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, is the unit outward normal to the boundary at the collocation points .

(d) Final System of Equations.
The global solution is obtained by solving (2.5)–(2.8). This requires the discretization of the interface , 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 equations. The relation 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 and the unknown virtual load amplitudes, and , 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 . Thus, (2.5) and (2.6) are kept as before and (2.8) is simplified to

The solution of the boundary integral along the surface () 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 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:

The solution of these equations is once again obtained as described above, with a system of 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: 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 () 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 and , with radii of and (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.

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

The host elastic medium ( kg/m3) is homogeneous, permitting a -wave velocity of  m/s and an S-wave velocity of  m/s. The fluid in Inclusion 1 exhibits a density equal to  kg/m3 and allows P-wave velocity of  m/s, while the elastic Inclusion 2, of density  kg/m3, exhibits a dilatational and an -wave of velocities  m/s,  m/s, respectively. The excitation source is assumed to be inside the host medium at point . It is a harmonic dilatational line load whose amplitude varies sinusoidally in the third dimension according to  rad/m.

The responses are computed at receiver and , placed at and , 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 and 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 ) and pressure response (receiver ) 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.

Figure 4: Case  1: BEM, TBEM, MFS, and coupling formulations’ results: displacements at receiver () and pressure at receiver () when the system is excited by a blast load.
Figure 5: Case 2: BEM, TBEM, MFS, and coupling formulations’ results: displacements at receiver () and pressure at receiver () when the system is excited by a blast load.
Figure 6: Case 3: BEM, TBEM, MFS, and coupling formulations’ results: displacements at receiver () and pressure at receiver () 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)).

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 -wave and -wave velocities of 2630 m/s and 1416 m/s, respectively. The fluid-filled borehole is centred at with a radius of . Its fluid has a mass density of 1000 kg/m3 and permits a -wave speed of 1500 m/s. A null-thickness arc-shaped crack centred at has a radius of and a length of .

This system is illuminated by a wave field generated by a dilatational line load placed ( rad/m) from the crack at . The resulting displacement is obtained over a grid of 10140 receivers arranged along the and directions at equal intervals and placed from to and from to .

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: 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.

Figure 8: Displacements and pressures solutions for the excitation frequency 140 Hz.
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 to 20, as follows: the BEM/TBEM solutions were computed by discretizing the borehole and the crack interfaces with and boundary elements, respectively; the MFS solutions were obtained by simulating the borehole and the crack interface with and virtual sources/collocation points, respectively; the coupling MFS/TBEM solutions were obtained using and virtual sources/collocation points.

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 , 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 , parallel to the borehole axis (two-dimensional application), with a characteristic frequency of  Hz and which starts acting at  ms. A set of snapshots taken from computer animations is presented in Figure 12 to illustrate the resulting wave field at different time instants.

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.
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 where represents the amplitude; , corresponds to the time, is the time when the wavelet takes its maximum value, and is the characteristic (dominant) period of the Ricker wavelet. In the frequency domain, this pulse is written as where .

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 (with ). This procedure is later taken into account by rescaling the responses in the time domain with an exponential factor .

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

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 -wave and -wave velocities allowed in the host medium and its density remain constant at , , and , respectively. The fluid-filled borehole is centered at with a radius . Its medium has a mass density of , a -wave velocity of . 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 to and from to .

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 , 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 , the incident pulses are partly reflected back as -waves and -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 , but it is not easy to distinguish the two types of waves since they are almost coincident at this early stage. At , the reflected - and -waves are very well developed as they spread away from the crack. At , part of incident pulse that has just hit the borehole and been transmitted as -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 , 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 -waves and denote a delay in relation to the direct incident field, because of lower -wave speed inside the fluid.

The last snapshots ( and ) 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.

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

Definitions
-wave velocity,-wave velocity,, with ,, with ,, Hankel functions, functions, unit outward normal at , unit outward normal at , (the collocation point).

Solid Media Green’s Functions
One has
The derivatives of the above Green’s functions give the following tractions along the , , and directions, in the solid medium, with , , , and , . These expressions can be combined to obtain in the normal and tangential directions. In these equations, .

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 and . Along the boundary element, at , where the unit outward normal is defined by , and after the equilibrium of stresses, the following equations are expressed for , and generated by loads also applied along , and directions: with defining the unit outward normal at (the collocation point), , , , and .

Fluid Media Green’s Functions
One has where and .

Fluid Media Traction Green’s Functions
One has where and .

#### 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.
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.
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.
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.
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.
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.
12. A. Kulowski, “Algorithmic representation of the ray tracing technique,” Applied Acoustics, vol. 18, no. 6, pp. 449–469, 1985.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.
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.