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

Authors: António Tadeu; Igor Castro
Journal: Mathematical Problems in Engineering
 (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.


---

## Body

## 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., [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–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 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.Figure 1
The geometry of the problem.It is further assumed that this system is subjected to a dilatational line source placed at(xs,ys) 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.1)ϕinc(x,y,ω,kz)=-iA2H0(kα1r),
in which Hn(⋯) are second kind Hankel functions of the order n, i=-1, kα1=ω2/α12-kz2 with Im⁡(kα1)<0, r=(x-xs)2+(y-ys)2 and kz is the wavenumber along z.Then, the displacement field can be expressed as(2.2)uxinc(x,y,xs,ys,ω)=iA2kα1H1(kα1r)∂r∂x,uyinc(x,y,xs,ys,ω)=iA2kα1H1(kα1r)∂r∂y,uzinc(x,y,xs,ys,ω)=A2kzH0(kα1r).


### 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 Figure2).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 surfaceS1 and subjected to an incident field given by uinc, 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),
(2.3)cijui(1)(x0,y0,ω)=∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,x0,y0,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,x0,y0,ω)ds+uiinc(x0,y0,xs,ys,ω).
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; nn1 is the unit outward normal along boundary S1, at (x,y), defined by the vector nn1=(cos⁡θn1,sin⁡θn1). Gij(1)(x,y,x0,y0,ω) and Hij(1)(x,y,nn1,x0,y0,ω) define the fundamental solutions for displacements and tractions (Green’s functions), in direction j on the boundary S1 at (x,y), caused by a unit point force in direction i applied at the nodal point, (x0,y0) (see the appendix). uj(1)(x,y,ω) corresponds to displacements in direction j at (x,y), tj(1)(x,y,nn1,ω) specifies the nodal tractions in direction j on the boundary at (x,y) and uiinc(x0,y0,xs,ys,ω) to the displacement incident field at (x0,y0) along direction i, when the source is located at (xs,ys). The coefficient cij is equal to δij/2, with δij being the Kronecker delta, when the boundary is smooth.
Green’s functions for displacements along thex, 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 S2 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 anj,n_ext(2), while those placed outside the inclusion have unknown amplitudes anj,n_int⁡(2). In the exterior, and interior elastic media the scattered displacement fields are given by
(2.4)ui(1)(x,y,ω)=∑n=1NS∑j=13anj,n_ext(2)Gji(1)(x,y,xn_ext,yn_ext,ω),ui(3)(x,y,ω)=∑n=1NS∑j=13anj,n_int⁡(2)Gji(3)(x,y,xn_int⁡,yn_int⁡,ω),
where Gji(1)(x,y,xn_ext,yn_ext,ω) and Gji(3)(x,y,xn_int⁡,yn_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 (xn_ext,yn_ext) and (xn_int⁡,yn_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.5)cijui(1)(x0,y0,ω)=∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,x0,y0,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,x0,y0,ω)ds+uiinc(x0,y0,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)Gji(1)(x0,y0,xn_ext,yn_ext,ω).(b) In the Interior Domain of Inclusion 1 (Medium 2),
(2.6)cp(2)(x0,y0,ω)=∫S1q(2)(x,y,-nn1,ω)Gf(2)(x,y,x0,y0,ω)ds-∫S1p(2)(x,y,ω)Hf(2)(x,y,-nn1,x0,y0,ω)ds.
In (2.6), the superscript 2 corresponds to the domain inside Inclusion 1, Medium 2. Gf(2)(x,y,x0,y0,ω) and Hf(2)(x,y,-nn1,x0,y0,ω) are Green’s functions for pressure and the gradient of pressure on the boundary S1 at (x,y), caused by a unit point pressure at the nodal point, (x0,y0) (see the appendix). p(2)(x,y,ω) corresponds to the pressure at (x,y), q(2)(x,y,nn1,ω) 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 loadsanj,n_ext(2) and anj,n_int⁡(2), it is also necessary to impose the continuity of displacements and tractions at interface S2, which is the boundary of Inclusion 2, along NS collocation points (xcol,ycol). This must be done so as to take into account the scattered field generated at Inclusion 1. The following two equations are thus defined:
(2.7)∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,xcol,ycol,ω)ds+uiinc(xcol,ycol,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)Gji(1)(xcol,ycol,xn_ext,yn_ext,ω)=∑n=1NS∑j=13anj,n_int(2)Gji(3)(xcol,ycol,xn_int⁡,yn_int⁡,ω),(2.8)∫S1t1(1)(x,y,nn1,ω)G̅i1(1)(x,y,nn2,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)H̅ij(1)(x,y,nn1,nn2,xcol,ycol,ω)ds+u̅iinc(xcol,ycol,nn2,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)G̅ji(1)(xcol,ycol,nn2,xn_ext,yn_ext,ω)=∑n=1NS∑j=13anj,n_int⁡(2)G̅ji(3)(xcol,ycol,nn2,xn_int⁡,yn_int⁡,ω).
Green’s functionsH̅ij(1)(x,y,nn1,nn2,xcol,ycol,ω) and G̅ij(1)(x,y,nn2,xcol,ycol,ω) are defined by applying the traction operator to Hij(1)(x,y,nn1,xcol,ycol,ω)and Gij(1)(x,y,xcol,ycol,ω), 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, nn2 is the unit outward normal to the boundary S2 at the collocation points (xcol,ycol).(d) Final System of Equations.
The global solution is obtained by solving (2.5)–(2.8). This requires the discretization of the interface S1, 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 u1(1)=-(1/ρ2)(∂p(2)/∂nn1) 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 boundaryS1 and the unknown virtual load amplitudes, anj,n_ext(2) and anj,n_int⁡(2), 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 boundaryS2. Thus, (2.5) and (2.6) are kept as before and (2.8) is simplified to
(2.9)∫S1t1(1)(x,y,nn1,ω)G̅i1(1)(x,y,nn2,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)H̅ij(1)(x,y,nn1,nn2,xcol,ycol,ω)ds+u̅iinc(xcol,ycol,nn2,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)G¯ji(1)(xcol,ycol,nn2,xn_ext,yn_ext,ω)=0.The solution of the boundary integral along the surface (S1) 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.


#### 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:(2.10)∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,xcol,ycol,ω)ds+uiinc(xcol,ycol,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)Gji(1)(xcol,ycol,xn_ext,yn_ext,ω)=0.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.


### 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:(2.11)aijui(1)(x0,y0,ω)+ci1t1(1)(x0,y0,nn1,ω)=∫S1t1(1)(x,y,nn1,ω)G¯i1(1)(x,y,nn2,x0,y0,ω)ds-∑j=13∫S1uj(1)(x,y,ω)H¯ij(1)(x,y,nn1,nn2,x0,y0,ω)ds+u¯iinc(x0,y0,nn2,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)G¯ji(1)(x0,y0,nn2,xn_ext,yn_ext,ω),ap(2)(x0,y0,ω)+cq(2)(x0,y0,-nn1,ω)=∫S1q(2)(x,y,-nn1,ω)G¯f(2)(x,y,nn2,x0,y0,ω)ds-∫S1H¯f(2)(x,y,-nn1,nn2,x0,y0,ω)p(1)(x,y,ω)ds.
As noted by Guiggiani [49], the coefficients aij and a are zero for piecewise straight boundary elements. The factors cij 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 (S1) 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.


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


## 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 Figure2).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 surfaceS1 and subjected to an incident field given by uinc, 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),
(2.3)cijui(1)(x0,y0,ω)=∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,x0,y0,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,x0,y0,ω)ds+uiinc(x0,y0,xs,ys,ω).
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; nn1 is the unit outward normal along boundary S1, at (x,y), defined by the vector nn1=(cos⁡θn1,sin⁡θn1). Gij(1)(x,y,x0,y0,ω) and Hij(1)(x,y,nn1,x0,y0,ω) define the fundamental solutions for displacements and tractions (Green’s functions), in direction j on the boundary S1 at (x,y), caused by a unit point force in direction i applied at the nodal point, (x0,y0) (see the appendix). uj(1)(x,y,ω) corresponds to displacements in direction j at (x,y), tj(1)(x,y,nn1,ω) specifies the nodal tractions in direction j on the boundary at (x,y) and uiinc(x0,y0,xs,ys,ω) to the displacement incident field at (x0,y0) along direction i, when the source is located at (xs,ys). The coefficient cij is equal to δij/2, with δij being the Kronecker delta, when the boundary is smooth.
Green’s functions for displacements along thex, 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 S2 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 anj,n_ext(2), while those placed outside the inclusion have unknown amplitudes anj,n_int⁡(2). In the exterior, and interior elastic media the scattered displacement fields are given by
(2.4)ui(1)(x,y,ω)=∑n=1NS∑j=13anj,n_ext(2)Gji(1)(x,y,xn_ext,yn_ext,ω),ui(3)(x,y,ω)=∑n=1NS∑j=13anj,n_int⁡(2)Gji(3)(x,y,xn_int⁡,yn_int⁡,ω),
where Gji(1)(x,y,xn_ext,yn_ext,ω) and Gji(3)(x,y,xn_int⁡,yn_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 (xn_ext,yn_ext) and (xn_int⁡,yn_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.5)cijui(1)(x0,y0,ω)=∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,x0,y0,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,x0,y0,ω)ds+uiinc(x0,y0,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)Gji(1)(x0,y0,xn_ext,yn_ext,ω).(b) In the Interior Domain of Inclusion 1 (Medium 2),
(2.6)cp(2)(x0,y0,ω)=∫S1q(2)(x,y,-nn1,ω)Gf(2)(x,y,x0,y0,ω)ds-∫S1p(2)(x,y,ω)Hf(2)(x,y,-nn1,x0,y0,ω)ds.
In (2.6), the superscript 2 corresponds to the domain inside Inclusion 1, Medium 2. Gf(2)(x,y,x0,y0,ω) and Hf(2)(x,y,-nn1,x0,y0,ω) are Green’s functions for pressure and the gradient of pressure on the boundary S1 at (x,y), caused by a unit point pressure at the nodal point, (x0,y0) (see the appendix). p(2)(x,y,ω) corresponds to the pressure at (x,y), q(2)(x,y,nn1,ω) 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 loadsanj,n_ext(2) and anj,n_int⁡(2), it is also necessary to impose the continuity of displacements and tractions at interface S2, which is the boundary of Inclusion 2, along NS collocation points (xcol,ycol). This must be done so as to take into account the scattered field generated at Inclusion 1. The following two equations are thus defined:
(2.7)∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,xcol,ycol,ω)ds+uiinc(xcol,ycol,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)Gji(1)(xcol,ycol,xn_ext,yn_ext,ω)=∑n=1NS∑j=13anj,n_int(2)Gji(3)(xcol,ycol,xn_int⁡,yn_int⁡,ω),(2.8)∫S1t1(1)(x,y,nn1,ω)G̅i1(1)(x,y,nn2,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)H̅ij(1)(x,y,nn1,nn2,xcol,ycol,ω)ds+u̅iinc(xcol,ycol,nn2,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)G̅ji(1)(xcol,ycol,nn2,xn_ext,yn_ext,ω)=∑n=1NS∑j=13anj,n_int⁡(2)G̅ji(3)(xcol,ycol,nn2,xn_int⁡,yn_int⁡,ω).
Green’s functionsH̅ij(1)(x,y,nn1,nn2,xcol,ycol,ω) and G̅ij(1)(x,y,nn2,xcol,ycol,ω) are defined by applying the traction operator to Hij(1)(x,y,nn1,xcol,ycol,ω)and Gij(1)(x,y,xcol,ycol,ω), 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, nn2 is the unit outward normal to the boundary S2 at the collocation points (xcol,ycol).(d) Final System of Equations.
The global solution is obtained by solving (2.5)–(2.8). This requires the discretization of the interface S1, 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 u1(1)=-(1/ρ2)(∂p(2)/∂nn1) 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 boundaryS1 and the unknown virtual load amplitudes, anj,n_ext(2) and anj,n_int⁡(2), 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 boundaryS2. Thus, (2.5) and (2.6) are kept as before and (2.8) is simplified to
(2.9)∫S1t1(1)(x,y,nn1,ω)G̅i1(1)(x,y,nn2,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)H̅ij(1)(x,y,nn1,nn2,xcol,ycol,ω)ds+u̅iinc(xcol,ycol,nn2,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)G¯ji(1)(xcol,ycol,nn2,xn_ext,yn_ext,ω)=0.The solution of the boundary integral along the surface (S1) 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.


### 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:(2.10)∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,xcol,ycol,ω)ds+uiinc(xcol,ycol,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)Gji(1)(xcol,ycol,xn_ext,yn_ext,ω)=0.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.


## 2.1.1. Fluid Inclusion 1 and Elastic Inclusion 2

Assuming Inclusion 1 to be bounded by a surfaceS1 and subjected to an incident field given by uinc, 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),
(2.3)cijui(1)(x0,y0,ω)=∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,x0,y0,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,x0,y0,ω)ds+uiinc(x0,y0,xs,ys,ω).
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; nn1 is the unit outward normal along boundary S1, at (x,y), defined by the vector nn1=(cos⁡θn1,sin⁡θn1). Gij(1)(x,y,x0,y0,ω) and Hij(1)(x,y,nn1,x0,y0,ω) define the fundamental solutions for displacements and tractions (Green’s functions), in direction j on the boundary S1 at (x,y), caused by a unit point force in direction i applied at the nodal point, (x0,y0) (see the appendix). uj(1)(x,y,ω) corresponds to displacements in direction j at (x,y), tj(1)(x,y,nn1,ω) specifies the nodal tractions in direction j on the boundary at (x,y) and uiinc(x0,y0,xs,ys,ω) to the displacement incident field at (x0,y0) along direction i, when the source is located at (xs,ys). The coefficient cij is equal to δij/2, with δij being the Kronecker delta, when the boundary is smooth.
Green’s functions for displacements along thex, 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 S2 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 anj,n_ext(2), while those placed outside the inclusion have unknown amplitudes anj,n_int⁡(2). In the exterior, and interior elastic media the scattered displacement fields are given by
(2.4)ui(1)(x,y,ω)=∑n=1NS∑j=13anj,n_ext(2)Gji(1)(x,y,xn_ext,yn_ext,ω),ui(3)(x,y,ω)=∑n=1NS∑j=13anj,n_int⁡(2)Gji(3)(x,y,xn_int⁡,yn_int⁡,ω),
where Gji(1)(x,y,xn_ext,yn_ext,ω) and Gji(3)(x,y,xn_int⁡,yn_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 (xn_ext,yn_ext) and (xn_int⁡,yn_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.5)cijui(1)(x0,y0,ω)=∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,x0,y0,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,x0,y0,ω)ds+uiinc(x0,y0,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)Gji(1)(x0,y0,xn_ext,yn_ext,ω).(b) In the Interior Domain of Inclusion 1 (Medium 2),
(2.6)cp(2)(x0,y0,ω)=∫S1q(2)(x,y,-nn1,ω)Gf(2)(x,y,x0,y0,ω)ds-∫S1p(2)(x,y,ω)Hf(2)(x,y,-nn1,x0,y0,ω)ds.
In (2.6), the superscript 2 corresponds to the domain inside Inclusion 1, Medium 2. Gf(2)(x,y,x0,y0,ω) and Hf(2)(x,y,-nn1,x0,y0,ω) are Green’s functions for pressure and the gradient of pressure on the boundary S1 at (x,y), caused by a unit point pressure at the nodal point, (x0,y0) (see the appendix). p(2)(x,y,ω) corresponds to the pressure at (x,y), q(2)(x,y,nn1,ω) 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 loadsanj,n_ext(2) and anj,n_int⁡(2), it is also necessary to impose the continuity of displacements and tractions at interface S2, which is the boundary of Inclusion 2, along NS collocation points (xcol,ycol). This must be done so as to take into account the scattered field generated at Inclusion 1. The following two equations are thus defined:
(2.7)∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,xcol,ycol,ω)ds+uiinc(xcol,ycol,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)Gji(1)(xcol,ycol,xn_ext,yn_ext,ω)=∑n=1NS∑j=13anj,n_int(2)Gji(3)(xcol,ycol,xn_int⁡,yn_int⁡,ω),(2.8)∫S1t1(1)(x,y,nn1,ω)G̅i1(1)(x,y,nn2,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)H̅ij(1)(x,y,nn1,nn2,xcol,ycol,ω)ds+u̅iinc(xcol,ycol,nn2,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)G̅ji(1)(xcol,ycol,nn2,xn_ext,yn_ext,ω)=∑n=1NS∑j=13anj,n_int⁡(2)G̅ji(3)(xcol,ycol,nn2,xn_int⁡,yn_int⁡,ω).
Green’s functionsH̅ij(1)(x,y,nn1,nn2,xcol,ycol,ω) and G̅ij(1)(x,y,nn2,xcol,ycol,ω) are defined by applying the traction operator to Hij(1)(x,y,nn1,xcol,ycol,ω)and Gij(1)(x,y,xcol,ycol,ω), 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, nn2 is the unit outward normal to the boundary S2 at the collocation points (xcol,ycol).(d) Final System of Equations.
The global solution is obtained by solving (2.5)–(2.8). This requires the discretization of the interface S1, 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 u1(1)=-(1/ρ2)(∂p(2)/∂nn1) 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 boundaryS1 and the unknown virtual load amplitudes, anj,n_ext(2) and anj,n_int⁡(2), 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 boundaryS2. Thus, (2.5) and (2.6) are kept as before and (2.8) is simplified to
(2.9)∫S1t1(1)(x,y,nn1,ω)G̅i1(1)(x,y,nn2,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)H̅ij(1)(x,y,nn1,nn2,xcol,ycol,ω)ds+u̅iinc(xcol,ycol,nn2,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)G¯ji(1)(xcol,ycol,nn2,xn_ext,yn_ext,ω)=0.The solution of the boundary integral along the surface (S1) 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.


## 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:(2.10)∫S1t1(1)(x,y,nn1,ω)Gi1(1)(x,y,xcol,ycol,ω)ds-∑j=13∫S1uj(1)(x,y,ω)Hij(1)(x,y,nn1,xcol,ycol,ω)ds+uiinc(xcol,ycol,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)Gji(1)(xcol,ycol,xn_ext,yn_ext,ω)=0.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.


## 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:(2.11)aijui(1)(x0,y0,ω)+ci1t1(1)(x0,y0,nn1,ω)=∫S1t1(1)(x,y,nn1,ω)G¯i1(1)(x,y,nn2,x0,y0,ω)ds-∑j=13∫S1uj(1)(x,y,ω)H¯ij(1)(x,y,nn1,nn2,x0,y0,ω)ds+u¯iinc(x0,y0,nn2,xs,ys,ω)+∑n=1NS∑j=13anj,n_ext(2)G¯ji(1)(x0,y0,nn2,xn_ext,yn_ext,ω),ap(2)(x0,y0,ω)+cq(2)(x0,y0,-nn1,ω)=∫S1q(2)(x,y,-nn1,ω)G¯f(2)(x,y,nn2,x0,y0,ω)ds-∫S1H¯f(2)(x,y,-nn1,nn2,x0,y0,ω)p(1)(x,y,ω)ds.
As noted by Guiggiani [49], the coefficients aij and a are zero for piecewise straight boundary elements. The factors cij 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 (S1) 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.


## 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.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 kz=0.2 rad/m.The responses are computed at receiverR1 and R2, 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 distances0.9×r and 1.1×r from the centre of the inclusion, with r being radii.Figures4, 5, and 6  present the real (solid line) and imaginary (dashed line) parts of the displacements ux, uy, and uz (receiver R1) and pressure response (receiver R2) 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.Case  1: BEM, TBEM, MFS, and coupling formulations’ results: displacements at receiver (R1) and pressure at receiver (R2) when the system is excited by a blast load.
(a)(b)(c)(d)Case 2: BEM, TBEM, MFS, and coupling formulations’ results: displacements at receiver (R1) and pressure at receiver (R2) when the system is excited by a blast load.
(a)(b)(c)(d)Case 3: BEM, TBEM, MFS, and coupling formulations’ results: displacements at receiver (R1) and pressure at receiver (R2) when the system is excited by a blast load: BEM.
(a)(b)(c)(d)


## 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 (Figure7(a)).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.
(a)(b)(c)(d)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 placed0.05m (kz=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 x and y directions at equal intervals and placed from x=-0.10m to x=0.25m and from y=-0.15m to y=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 Figure7(b)).Figures8 and 9 illustrate the real and imaginary part of the reference solutions for both excitation frequencies.Displacements and pressures solutions for the excitation frequency 140 Hz.
(a)(b)(c)(d)(e)(f)Displacements and pressures solutions for the excitation frequency 9000 Hz.
(a)(b)(c)(d)(e)(f)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 Figure7(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 IntelCore Duo CPU E6750.Figure10 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.Global domain error versus CPU time for solving the system composed of a fluid-filled borehole placed in the vicinity of a crack.
(a)(b)(c)The global domain errors shown in Figure10 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 of0.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 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.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.Snapshots illustrating the displacement,ux and uy, and the pressure generated by a line blast load, modelled as a Ricker pulse with a characteristic frequency of 30000 Hz.
(a)(b)(c)(d)(e)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(5.1)u(τ)=A(1-2τ2)e-τ2,
where A represents the amplitude; τ=(t-ts)/t0, t corresponds to the time, ts is the time when the wavelet takes its maximum value, and πt0 is the characteristic (dominant) period of the Ricker wavelet. In the frequency domain, this pulse is written as
(5.2)U(ω)=A[2t0πe-iωts]Ω2e-Ω2,
where Ω=ωt0/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 from140Hz 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. Figure11 illustrates the position of the virtual sources, collocation points, and boundary elements.TheP-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 thex and y directions at equal intervals and placed in the vicinity of the inclusion and crack from x=-0.4m to x=0.4m and from y=-0.4m to y=0.4m.A set of snapshots taken from computer animations is presented in Figure12 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 (ux) and the vertical displacements (uy) 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, att=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.Att=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 t=0.07ms, but it is not easy to distinguish the two types of waves since they are almost coincident at this early stage. At t=0.10ms, the reflected P- and S-waves are very well developed as they spread away from the crack. At t=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.Byt=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 (t=0.10ms and t=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.



---

*Source: 159389-2011-11-22.xml*