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

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.


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., 1-6 .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 [14][15][16] .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 22-24 .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 33-35 .The two approaches most often used are a direct coupling and iterative coupling 36-38 .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 fluidfilled borehole driven in a cracked elastic medium.

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.It is further assumed that this system is subjected to a dilatational line source placed at x s , y s whose amplitude varies sinusoidally in the third dimension z .
The incident wave field generated by this source can be expressed in the frequency domain by means of the classic dilatational potential: 2 and k z is the wavenumber along z.
Then, the displacement field can be expressed as

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 .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 .
• O (x s , y s ) Nodal points Virtual loads Figure 2: Discretization of the system: position of virtual loads, collocation points and boundary elements.

Fluid Inclusion 1 and Elastic Inclusion 2
Assuming Inclusion 1 to be bounded by a surface S 1 and subjected to an incident field given by u inc , the boundary integral equation can be constructed by applying the reciprocity theorem e.g., Manolis and Beskos 44 leading to the following.

2.3
In this equation, i, j 1, 2 correspond to the normal and tangential directions relative to the inclusion surface, while i, j 3 correspond to the z direction.In these equations, the superscript 1 represents the exterior domain; n n1 is the unit outward normal along boundary S 1 , at x, y , defined by the vector n n1 cos θ n1 , sin θ n1 .G 1 ij x, y, x 0 , y 0 , ω and H 1 ij x, y, n n1 , x 0 , y 0 , ω define the fundamental solutions for displacements and tractions Green's functions , in direction j on the boundary S 1 at x, y , caused by a unit point force in direction i applied at the nodal point, x 0 , y 0 see the appendix .u 1 j x, y, ω corresponds to displacements in direction j at x, y , t x, y, n n1 , ω specifies the nodal tractions in direction j on the boundary at x, y and u inc i x 0 , y 0 , x s , y s , ω to the displacement incident field at x 0 , y 0 along direction i, when the source is located at x s , y s .The coefficient c ij is equal to δ ij /2, with δ ij being the Kronecker delta, when the boundary is smooth.
Green's functions for displacements along the x, y, and z 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 NS virtual sources.These virtual loads are distributed along the inclusion interface S 2 at distances δ from that boundary towards the interior and exterior of the inclusion lines C 1 and C 2 in Figure 2 in order to avoid singularities.Sources inside the inclusion have unknown amplitudes a 2 nj,n ext , while those placed outside the inclusion have unknown amplitudes a 2 nj,n int .In the exterior, and interior elastic media the scattered displacement fields are given by where G 1 ji x, y, x n ext , y n ext , ω and G 3 ji x, y, x n int , y n int , ω are the fundamental solutions which represent the displacements at points x, y in Mediums 1 and 3, in direction i, caused by a unit point force in direction j applied at the positions x n ext , y n ext and x n int , y n int .n ext and n int are the subscripts that denote the load order number placed along lines C 1 and C 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,

2.6
In 2.6 , the superscript 2 corresponds to the domain inside Inclusion 1, Medium 2.
x, y, −n n1 , x 0 , y 0 , ω are Green's functions for pressure and the gradient of pressure on the boundary S 1 at x, y , caused by a unit point pressure at the nodal point, x 0 , y 0 see the appendix .p 2 x, y, ω corresponds to the pressure at x, y , q 2 x, y, n n1 , ω specifies the nodal pressure gradients on the boundary at x, y .The coefficient c 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 a 2 nj,n ext and a 2 nj,n int , it is also necessary to impose the continuity of displacements and tractions at interface S 2 , which is the boundary of Inclusion 2, along NS collocation points x col , y 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: x, y, n n1 , x col , y col , ω ds

2.8
Green's functions x col , y col , ω , which can be obtained by combining the derivatives of the former Green's functions, in order of x, y, and z, so as to obtain the stresses see the appendix .In these equations, n n2 is the unit outward normal to the boundary S 2 at the collocation points x col , y col .

(d) Final System of Equations.
The global solution is obtained by solving 2.5 -2.8 .This requires the discretization of the interface S 1 , which is the boundary of Inclusion 1.For the purposes of this work, this interface is discretized into N 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 6NS 4N × 6NS 4N equations.The relation u ∂n n1 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 S 1 and the unknown virtual load amplitudes, a 2 nj,n ext and a 2 nj,n int , which allow the displacement field to be defined inside and outside the inclusions.

Empty Inclusion 2 (Null Tractions Along its Boundary)
In this case, null tractions are prescribed along the boundary S 2 .Thus, 2.5 and 2.6 are kept as before and 2.8 is simplified to

2.9
The solution of the boundary integral along the surface S 1 again requires its discretization into N straight boundary elements, while the simulation of Inclusion 2 uses NS collocation points/virtual sources, following a procedure similar to the one described above.This leads to a system of 3NS 4N × 3NS 4N equations.

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: x, y, n n1 , x col , y col , ω ds

2.10
Mathematical Problems in Engineering 9 The solution of these equations is once again obtained as described above, with a system of 3NS 4N × 3NS 4N equations.Other coupling combinations can be solved in the same way.

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: x 0 , y 0 , n n1 , ω x, y, n n1 , ω G 1 i1 x, y, n n2 , x 0 , y 0 , ω ds x, y, n n2 , x 0 , y 0 , ω ds x, y, −n n1 , n n2 , x 0 , y 0 , ω p 1 x, y, ω ds.

2.11
As noted by Guiggiani 49 , the coefficients a ij and a are zero for piecewise straight boundary elements.The factors c ij and c 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 S 1 into N 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.

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 .

Verification of the Coupling Algorithms
The All the illustrated simulations used interior and exterior virtual sources respectively placed at distances 0.9 × r and 1.1 × r from the centre of the inclusion, with r being radii.
Figures 4, 5, and 6 present the real solid line and imaginary dashed line parts of the displacements u x , u y , and u z receiver R 1 and pressure response receiver R 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.

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 .The host medium, with a density of 2250 kg/m 3 , allows P-wave and S-wave velocities of 2630 m/s and 1416 m/s, respectively.The fluid-filled borehole is centred at 0.0 m, 0.0 m with a radius of 0.05 m.Its fluid has a mass density of 1000 kg/m 3 and permits a P-wave speed of 1500 m/s.A null-thickness arc-shaped crack centred at 0.0 m, 0.0 m has a radius of 0.10 m and a length of 1.2π/32 m.
This system is illuminated by a wave field generated by a dilatational line load placed 0.05 m k z 0 rad/m from the crack at 0.15 m, 0.0 m .The resulting displacement is obtained over a grid of 10140 receivers arranged along the x and y directions at equal intervals and placed from x −0.10 m to x 0.25 m and from y −0.15 m to y 0.15 m.
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: 140 Hz 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.
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, T , and a fictitious interface, F. In order to correctly describe the behaviour of the null-thickness crack, null tractions are ascribed to both sides of interface T and continuity of displacements and tractions is imposed along the interface F. 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 m 1 to 20, as follows: the BEM/TBEM solutions were computed by discretizing the borehole and the crack interfaces with 10m and 4m boundary elements, respectively; the MFS solutions were obtained by simulating the borehole and the crack interface with 10m and 20m virtual sources/collocation points, respectively; the coupling MFS/TBEM solutions were obtained using 10m and 4m virtual sources/collocation points.
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.

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.1016 m, 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.1 m, −0.3 m , parallel to the borehole axis two-dimensional application , with a characteristic frequency of 30000 Hz and which starts acting at t 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.
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 A represents the amplitude; τ t − t s /t 0 , t corresponds to the time, t s is the time when the wavelet takes its maximum value, and πt 0 is the characteristic dominant period of the Ricker wavelet.In the frequency domain, this pulse is written as where Ω ωt 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 ω c ω − iη with η 0.7Δω .This procedure is later taken into account by rescaling the responses in the time domain with an exponential factor e ηt .
The computations are performed in the frequency domain for frequencies ranging from 140 Hz to 71680 Hz, with a frequency increment of 140 Hz, which determines a total time window of 7.14 ms.
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 4208 m/s, 2656 m/s, and 2140 kg/m 3 , respectively.The fluid-filled borehole is centered at 0.0 m, 0.0 m with a radius 0.1016 m.Its medium has a mass density of 1000 kg/m 3 , a P-wave velocity of 1500 m/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 x and y directions at equal intervals and placed in the vicinity of the inclusion and crack from x −0.4 m to x 0.4 m and from y −0.4 m to y 0.4 m.
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 u x and the vertical displacements u y 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 t 0.01 ms, the pulse excited by the dilatational source can be seen travelling in the elastic medium without perturbations as it has not yet reached the fluidfilled 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 t 0.04 ms, 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 t 0.07 ms, but it is not easy to distinguish the two types of waves since they are almost coincident at this early stage.At t 0.10 ms, the reflected P-and S-waves are very well developed as they spread away from the crack.At t 0.04 ms, 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 t 0.07 ms, 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 t 0.10 ms and t 0.15 ms 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.

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.
The derivatives of the above Green's functions give the following tractions along the x, y, and z directions, in the solid medium, with n n1 cos θ n1 , sin θ n1 , H rt H rt x, y, n n1 , x 0 , y 0 , ω , G rt G rt x, y, x 0 , y 0 , ω , and r, t x, y, z.These expressions can be combined to obtain H ij x, y, n n1 , x 0 , y 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 x, y, and z, so as to obtain stresses G ij x, y, n n2 , x col , y col , ω and H ij x, y, n n1 , n n2 , x col , y col , ω .Along the boundary element, at x, y , where the unit outward normal is defined by n n1 cos θ n1 , sin θ n1 , and after the equilibrium of stresses, the following equations are expressed for x, y and z generated by loads also applied along x, y and z directions: where G f G f x, y, x 0 , y 0 , ω and H f H f x, y, n n1 , x 0 , y 0 , ω .

Fluid Media Traction Green's Functions
One has where G f G f x, y, n n2 , x 0 , y 0 , ω and H f H f x, y, n n1 , n n2 , x 0 , y 0 , ω .

Figure 1 :
Figure 1: The geometry of the problem.

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

Figure 4 :
Figure 4: Case 1: BEM, TBEM, MFS, and coupling formulations' results: displacements at receiver R 1 and pressure at receiver R 2 when the system is excited by a blast load.

FrequencyFigure 5 :
Figure 5: Case 2: BEM, TBEM, MFS, and coupling formulations' results: displacements at receiver R 1 and pressure at receiver R 2 when the system is excited by a blast load.

FrequencyFigure 6 :
Figure 6: Case 3: BEM, TBEM, MFS, and coupling formulations' results: displacements at receiver R 1 and pressure at receiver R 2 when the system is excited by a blast load: BEM.

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

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

Figure 11 :
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 :
Figure 12: Snapshots illustrating the displacement, u x and u y , and the pressure generated by a line blast load, modelled as a Ricker pulse with a characteristic frequency of 30000 Hz.