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

Abstract

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

1. Introduction

Various numerical methods have been proposed to simulate the propagation of waves in elastic and acoustic media, since analytical solutions are only known for simple and regular geometries (e.g., [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 𝜌₁ (Figure 1) and allowing longitudinal (P-wave) and shear waves (S-wave) to travel at velocities 𝛼₁ and 𝛽₁, respectively. Medium 2, inside Inclusion 1, is fluid, has density 𝜌₂, and permits pressure waves (P-wave) to travel at velocity 𝛼₂. Inclusion 2 is elastic (Medium 3), has density 𝜌₃, and allows longitudinal and shear waves to travel at velocities 𝛼₃ and 𝛽₃, respectively.

Figure 1 

The geometry of the problem.

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

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

Then, the displacement field can be expressed as𝑢^inc_𝑥𝑥,𝑦,𝑥_𝑠,𝑦_𝑠,𝜔=i𝐴2𝑘^𝛼1𝐻₁𝑘_𝛼1𝑟𝜕𝑟𝑢𝜕𝑥_inc𝑦,𝑢𝑥,𝑦,𝑥_𝑠,𝑦_𝑠,𝜔=i2𝐴𝑘^𝛼1𝐻₁𝑘_𝛼1𝑟𝜕𝜕𝑦_inc𝑧𝑟,𝑥,𝑦,𝑥_𝑠,𝑦_𝑠,𝜔=𝐴2𝑘^𝑧𝐻₀𝑘_𝛼1𝑟.(2.2)

2.1. BEM/MFS Coupling Formulation

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

Figure 2 

Discretization of the system: position of virtual loads, collocation points and boundary elements.

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

2.1.1. Fluid Inclusion 1 and Elastic Inclusion 2

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

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

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

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

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

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

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

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

In this case, null tractions are prescribed along the boundary 𝑆₂. Thus, (2.5) and (2.6) are kept as before and (2.8) is simplified to _𝑆1𝑡−¹⁽¹⁾𝑥,𝑦,𝐧_𝑛1,𝜔𝐺^(1)𝑖1𝑥,𝑦,𝐧_𝑛2,𝑥_col,𝑦_col,𝜔d𝑠₃_𝑗=1_𝑆1+𝑢^𝑗(1)(𝑥,𝑦,𝜔)𝐻^(1)𝑖𝑗𝑥,𝑦,𝐧_𝑛1,𝐧_𝑛2,𝑥_col,𝑦_col,𝜔d𝑠𝑢^inc_𝑖𝑥_col,𝑦_col,𝐧_𝑛2,𝑥_𝑠,𝑦_𝑠,𝜔+^𝑁𝑆_3𝑛=1_𝑗=1𝑎^{(2)𝑛𝑗,𝑛}_{_ext}𝐺^(1)𝑗𝑖𝑥_col,𝑦_col,𝐧_𝑛2,𝑥_{𝑛_ext},𝑦_{𝑛_ext},𝜔=0.(2.9)

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

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

Null displacements on the surface of Inclusion 2 are now prescribed, which leads to (2.5) and (2.6) and to the following equation:_𝑆1𝑡−¹⁽¹⁾𝑥,𝑦,𝐧_𝑛1,𝜔𝐺^(1)𝑖1𝑥,𝑦,𝑥_col,𝑦_col,𝜔d𝑠₃_𝑗=1_𝑆1+𝑢_inc𝑖𝑢^𝑗(1)(𝑥,𝑦,𝜔)𝐻^(1)𝑖𝑗𝑥,𝑦,𝐧_𝑛1,𝑥_col,𝑦_col,𝜔d𝑠𝑥_col,𝑦_col,𝑥_𝑠,𝑦_𝑠,𝜔+^𝑁𝑆_3𝑛=1_𝑗=1𝑎^{(2)𝑛𝑗,𝑛}_{_ext}𝐺^(1)𝑗𝑖𝑥_col,𝑦_col,𝑥_{𝑛_ext},𝑦_{𝑛_ext},𝜔=0.(2.10)

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

2.2. TBEM/MFS Coupling Formulation

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

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

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

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

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

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

3. Verification of the Coupling Algorithms

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

Figure 3 

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

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

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

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

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

(a)

(b)

(c)

(d)

Figure 4 

Case  1: BEM, TBEM, MFS, and coupling formulations’ results: displacements at receiver (𝑅1) and pressure at receiver (𝑅2) when the system is excited by a blast load.

(a)

(b)

(c)

(d)

Figure 5 

Case 2: BEM, TBEM, MFS, and coupling formulations’ results: displacements at receiver (𝑅1) and pressure at receiver (𝑅2) when the system is excited by a blast load.

(a)

(b)

(c)

(d)

Figure 6 

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

4. Computational Efficiency of the Coupling Algorithms

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

(a)

(b)

(c)

(d)

Figure 7 

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

The host medium, with a density of 2250 kg/m³, 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/m³ and permits a P-wave speed of 1500 m/s. A null-thickness arc-shaped crack centred at (0.0m,0.0m) has a radius of 0.10m and a length of (1.2𝜋/32)m.

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

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

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

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

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8 

Displacements and pressures solutions for the excitation frequency 140 Hz.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 9 

Displacements and pressures solutions for the excitation frequency 9000 Hz.

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

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

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

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

(a)

(b)

(c)

Figure 10 

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

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

5. Numerical Application

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

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.

(a)

(b)

(c)

(d)

(e)

Figure 12 

Snapshots illustrating the displacement, 𝑢𝑥 and 𝑢𝑦, and the pressure generated by a line blast load, modelled as a Ricker pulse with a characteristic frequency of 30000 Hz.

The responses in the time domain are computed by applying an inverse (fast) Fourier transform to the responses in the frequency domain 𝜔. The Ricker pulse modelled is expressed in the time domain by𝑢(𝜏)=𝐴1−2𝜏²e^−𝜏2,(5.1) where 𝐴 represents the amplitude; 𝜏=(𝑡−𝑡_𝑠)/𝑡₀, 𝑡 corresponds to the time, 𝑡_𝑠 is the time when the wavelet takes its maximum value, and 𝜋𝑡₀ is the characteristic (dominant) period of the Ricker wavelet. In the frequency domain, this pulse is written as 𝑈(𝜔)=𝐴2𝑡₀√𝜋e^{−i𝜔𝑡𝑠}Ω²e^−Ω2,(5.2) where Ω=𝜔𝑡₀/2.

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

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

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

The P-wave and S-wave velocities allowed in the host medium and its density remain constant at 4208m/s, 2656m/s, and 2140kg/m³, 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/m³, a P-wave velocity of 1500m/s. A null-thickness crack is embedded in the vicinity of this elastic inclusion.

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

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

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

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

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

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

6. Conclusions

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

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

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

Appendix

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

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

Solid Media Green’s Functions
One has𝐺_𝑥𝑥=14i𝜌𝜔²𝜔𝛽²²𝐻_0𝛽−1𝑟𝐵¹+𝑥−𝑥⁰𝑟²𝐵₂𝐺_𝑦𝑦,=14i𝜌𝜔²𝜔𝛽²²𝐻_0𝛽−1𝑟𝐵¹+𝑦−𝑦⁰𝑟²𝐵₂𝐺_𝑧𝑧,=14i𝜌𝜔²𝜔𝛽²²𝐻_0𝛽−𝑘^2𝑧𝐵₀𝐺_𝑥𝑦,=𝐺_𝑦𝑥=14i𝜌𝜔²𝑥−𝑟𝑥⁰𝑦−𝑟𝑦⁰𝐵²𝐺_𝑥𝑧,=𝐺_𝑧𝑥=i𝑘_𝑧14i𝜌𝜔²𝑥−𝑟𝑥⁰𝐵¹𝐺_𝑦𝑧,=𝐺_𝑧𝑦=i𝑘_𝑧14i𝜌𝜔²𝑦−𝑟𝑦⁰𝐵¹.(A.1)
The derivatives of the above Green’s functions give the following tractions along the 𝑥, 𝑦, and 𝑧 directions, in the solid medium, 𝐻_𝑟𝑥=2𝜇𝛼²2𝛽²𝜕𝐺_𝑟𝑥𝜕𝑥+2𝛼𝛽²²−1𝜕𝐺^𝑟𝑦𝜕𝑦+𝜕𝜕𝐺𝑧^𝑟𝑧cos𝜃^𝑛1+𝜇𝜕𝐺^𝑟𝑦𝜕𝑥+𝜕𝜕𝐺𝑦^𝑟𝑥sin𝜃^𝑛1𝐻_𝑟𝑦,=2𝜇𝛼²2𝛽²−1𝜕𝐺^𝑟𝑥𝜕𝑥+𝜕𝜕𝐺𝑧^𝑟𝑧+2𝛼𝛽²²𝜕𝐺_𝑟𝑦𝜕𝑦sin𝜃^𝑛1+𝜇𝜕𝐺^𝑟𝑦𝜕𝑥+𝜕𝜕𝐺𝑦^𝑟𝑥cos𝜃^𝑛1𝐻_𝑟𝑧,=𝜇𝜕𝐺^𝑟𝑥𝜕𝑧+𝜕𝜕𝐺𝑥^𝑟𝑧cos𝜃^𝑛1+𝜇𝜕𝐺^𝑟𝑦𝜕𝑧+𝜕𝜕𝐺𝑦^𝑟𝑧sin𝜃^𝑛1,(A.2) with 𝐧_𝑛1=(cos𝜃_𝑛1,sin𝜃_𝑛1), 𝐻_𝑟𝑡=𝐻_𝑟𝑡(𝑥,𝑦,𝐧_𝑛1,𝑥₀,𝑦₀,𝜔), 𝐺_𝑟𝑡=𝐺_𝑟𝑡(𝑥,𝑦,𝑥₀,𝑦₀,𝜔), and 𝑟, 𝑡=𝑥,𝑦,𝑧. These expressions can be combined to obtain 𝐻_𝑖𝑗(𝑥,𝑦,𝐧_𝑛1,𝑥₀,𝑦₀,𝜔) in the normal and tangential directions. In these equations, 𝜇=𝜌𝛽².

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

Fluid Media Green’s Functions
One has𝐺_𝑓=−i4𝐻⁰𝑘_𝛼𝐻_𝑓𝑟,=i4𝑘^𝛼𝐻₀𝑘_𝛼𝑟𝜕𝑟𝜕𝐧_𝑛1,(A.4) where 𝐺_𝑓=𝐺_𝑓(𝑥,𝑦,𝑥₀,𝑦₀,𝜔) and 𝐻_𝑓=𝐻_𝑓(𝑥,𝑦,𝐧_𝑛1,𝑥₀,𝑦₀,𝜔).

Fluid Media Traction Green’s Functions
One has𝐺_𝑓=i4𝑘^𝛼𝐻₁𝑘_𝛼𝑟𝜕𝑟𝜕𝐧_𝑛2,𝐻_𝑓=i4𝑘^𝛼−𝑘^𝛼𝐻₂𝑘_𝛼𝑟𝜕𝑟𝜕𝑥²𝜕𝑥𝜕𝐧_𝑛1+𝜕𝑟𝜕𝑥𝜕𝜕𝑦𝑟𝜕𝜕𝐧𝑦_𝑛1+𝐻¹𝑘_𝛼𝑟𝑟𝜕𝜕𝐧𝑥_𝑛1𝜕𝜕𝐧𝑥_𝑛2+i4𝑘^𝛼−𝑘^𝛼𝐻₂𝑘_𝛼𝑟𝜕𝑟𝜕𝑥𝜕𝜕𝑦𝑟𝜕𝜕𝐧𝑥_𝑛1+𝜕𝑟𝜕𝑦²𝜕𝑦𝜕𝐧_𝑛1+𝐻¹𝑘_𝛼𝑟𝑟𝜕𝜕𝐧𝑦_𝑛1𝜕𝜕𝐧𝑦_𝑛2,(A.5) where 𝐺_𝑓=𝐺_𝑓(𝑥,𝑦,𝐧_𝑛2,𝑥₀,𝑦₀,𝜔) and 𝐻_𝑓=𝐻_𝑓(𝑥,𝑦,𝐧_𝑛1,𝐧_𝑛2,𝑥₀,𝑦₀,𝜔).

Acknowledgment

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