The scattering of harmonic waves has been studied extensively for problems in quantum mechanics, acoustics, electromagnetics, and elasticity. Solutions to elastodynamic problems are the basis for ultrasonic non-destructive evaluation measurement models. Therefore, in this study, we investigate the use of the boundary element method (BEM) in the frequency domain using an off-boundary technique in which the observation points are taken inside the scattering object. This methodology removes both non-integrable singularities from the domain of integration along with avoiding ill-conditioning effects that occur at fictitious eigenfrequencies of which both require using special computationally demanding procedures to obtain solutions. Additionally, we employ both free and half-space fundamental solutions (Green’s displacement tensors) to investigate the elastodynamic scattering of an incident, plane, time-harmonic longitudinal wave in the frequency domain of a homogeneous, isotropic, and linear elastic solid with one or two spherical cavities. The half-space fundamental solution reduces the number of required boundary elements in half, which significantly reduces computational resource requirements. We only consider spherical cavities in this paper to illustrate the full and half-space off-boundary BEM and analyze the interaction effects associated with the elastodynamic scattering by two symmetric spherical cavities. To verify the validity of the free and half-space off-boundary BEM formulations, surface displacement results are compared with existing surface displacement results and show good agreement. Finally, the half-space off-boundary BEM is used to illustrate the interaction effects of the back and forward scattered displacement fields as a function of the distance from the spherical cavities.
INTRODUCTION
The use of ultrasonic testing has become an important non-destructive evaluation (NDE) method for determining the presence and character of imperfections that exist in materials and fabricated parts. As pointed out by Achenbach,1 the solutions to elastodynamic problems form the basis for ultrasonic measurement models used for quantitative non-destructive evaluation (QNDE). These methods are used to assess the deterioration of structures and materials by detecting and characterizing inherent and developing discrete flaws. The results are then utilized to evaluate the requirements for the reliability of structures and their components along with their constituent materials. Therefore, these methods currently provide an increasing role in the quality control of material processing, pre-assembly validation, in service assessments, and the maintenance of existing structures and machinery by analyzing the scattering profiles of elastic waves by the inherent imperfections.
Numerous methods have been employed over the years for solving both acoustic and elastodynamic wave propagation problems. Wave scattering problems in mathematical physics concentrate on finding the solution to the Helmholtz equation2 in an exterior domain. In the process of this solution, the integral equation method is very convenient. Helmholtz equation arises when steady-state solutions of the scalar wave equation for single frequency sources are required. For acoustic scattering problems, an incident wave with harmonic time dependence is intercepted by a bounded scatterer producing reflected and diffracted waves that travel outward from the scattering region. These acoustic waves are described using a linear wave equation in terms of a reduced scalar velocity potential wave function φ (here, the time factor e−iωt has been omitted, where , and t is time), sound velocity c, angular velocity ω = kc, wave number k = 2π/λ, and the source wavelength λ that define Helmholtz equation as ∇2φ + k2φ = 0. Therefore, it is important to know how the incident wave field interacts with the scattering objects. For acoustic scattering problems, it is generally convenient to regard the total scalar wave function φt = φinc + φsc at any point as consisting of a known incident wave function φinc that would exist in the absence of the scattering objects, and a scattered wave function φsc that is to be determined, and diverges from the scattering regions, and satisfies a radiation condition at infinity to ensure it represents a diverging wave function. At the boundary of the scatterer, conditions are imposed on the total wave function, or its normal derivative, or a linear combination of both in which the mathematical problem is to obtain an exterior wave function that satisfies the boundary conditions, and the scattered field satisfies the radiation condition at infinity. In general, this problem requires a numerical solution; therefore, it is convenient to reformulate the problem in terms of a boundary integral equation. The fundamental feature of the boundary integral formulation is that it shifts the problem of solving a domain equation to a boundary equation that is generally represented by a boundary integral equation. For time-harmonic acoustic scattering problems the resulting boundary integral equations cannot usually be solved exactly in closed form; however, the boundary element method (BEM) has been shown to provide fairly accurate approximate solutions for acoustic scattering problems defined in terms of boundary integral equations. A comprehensive review of acoustic scattering problems is given by Zaman.3
More recently, a semi-analytical approach has been proposed for solving three-dimensional acoustic pressure wave scattering problems with one, two, and three sound-hard spheres of radius a in an unbounded acoustic medium.4 The sound-hard spheres are subjected to an incident time-harmonic plane sound pressure wave with unit amplitude, and a truncated collocation multipole expansion is obtained in terms of Legendre and spherical Hankel functions that give a finite linear algebraic system of equations that are solved giving the acoustic pressure wave scatter field that satisfies the radiation condition at infinity. Results are obtained for the near-field pressure intensities and far field scattering patterns. Solutions for the single sound-hard sphere were compared with both existing analytical and BEM solutions, and the multiple sphere solutions were compared only with existing BEM results. First, a convergence analysis was conducted to determine the number of expansion terms (truncation number) required for accurate pressure intensity solutions relative to a root mean square error. Next, the effects of sphere separation distance and the incident dimensionless wave numbers were analyzed. It is shown that the collocation multipole expansion truncation number for a single sound-hard sphere increases with increasing dimensionless wave numbers. However, for multiple sound-hard spheres, the collocation multipole expansion truncation number depends on the specified relative root mean square errors, separation distances of the scatters, and the dimensionless incident wave numbers. It is concluded from numerical simulations that fictitious eigenfrequencies do not occur with the proposed truncated collocation multipole expansion formulation. Furthermore, as the separation distance between sound-hard spheres is decreased, the near-field acoustic pressure intensity increases; conversely, the far-field acoustic pressure field becomes more fluctuated along the azimuthal angle as the separation distance is increased. Therefore, the effect of distance between scatters on the near-field pressure intensity is opposite from that at the far-field.
However, scalar wave models have been found to give poor approximations for elastic wave descriptions;5 therefore, by utilizing a complete elastodynamic scattering theory5 provides significantly more information about the different types of waves that propagate in elastic solids due to a solids resistance to shear in which more than one type of wave generally propagates.2 For the elastodynamic problems considered in this study, both longitudinal waves that are similar to the scalar acoustic waves also coexist with transverse waves that are dependent on the solid materials resistance to shear deformations and are not considered with acoustic wave scattering problems. Therefore, elastodynamic wave scattering in solid materials is significantly more complicated than the scattering of acoustic waves. Furthermore, elastodynamic boundary integral representations are formulated in terms of tensorial functions5 as opposed to scalar functions.3
Solutions to elastodynamic scattering problems provide results for different scattering source distributions in diverse media that are of fundamental importance. The quantification of elastic wave transmission and reflection at boundaries and interfaces along with the scattering associated by assorted imperfections is a mandatory requirement for QNDE. However, except for very simple geometries, along with homogeneous, isotropic, and linear elastic material properties,6 closed form analytical solutions are generally not possible for most geometric configurations under consideration. In some cases, the displacement and traction fields produced by scattering objects can be approximated; however, generally these approximations are only accurate for either very low or very high frequencies. For low frequency problems, it has been shown5,7 that by using volume or surface integral equation formulations based on elastostatic solutions and resulting static elastic fields provide a starting approximation for use with a systematic perturbation analysis to produce a Born series approximation or a regular non-singular power series that gives approximate elastodynamic solutions that are applicable for low frequencies. Additionally, a class of approximate solutions also exists for many high frequency elastodynamic scattering problems through the use of ray theory.2 Furthermore, the use of finer surface mesh discretizations with smaller elements to adequately represent rapidly changing spatial boundary variables associated with surface elements exhibiting high frequency responses along with special arrangements of the elements and also using more element Gaussian quadrature integration points all increase the accuracy of BEM solutions with increasing frequencies.8 However, there is a cost due to the increase in computational requirements that are required for evaluating finer surface mesh discretizations along with an increased number of Gaussian quadrature integration points.
For the frequencies of the elastodynamic scattering problems being considered in this study, a rigorous approach is employed for the formulation of solutions through the use of boundary integral representations. With this approach, the elastodynamic problem is reduced to an integral representation over the surfaces of the scattering objects in terms of fundamental solutions1,5,7,9 (Green’s displacement tensor), double layer kernel tensor (obtained from substituting the fundamental solution into Hooke’s law), and displacement and/or traction fields on the scattering objects. After the resulting integral representations are solved for the surface displacements and/or tractions, the displacement and stress fields everywhere in the exterior domain V can be obtained using the corresponding boundary integral representations. Unfortunately, these integral equations generally cannot be solved analytically, and require using numerical techniques, such as the BEM, which can be utilized in either the time or frequency domains. Additionally,10 domain formulations, such as the finite element method (FEM) and finite difference method (FDM), can be used for some problems along with hybrid schemes using a combination of FEM and BEM for an optimal solution. However, the problem with domain based methods is that they require discretization of not only surfaces of the scattering domains as with the BEM but also all other domains except voids making computer resource requirements unacceptable for many problems. Indirect BEM and other boundary methods5,10–12 have also been extensively used for many elastodynamic scattering problems, and especially for those involving multiple inhomogeneities. Some of the most utilized indirect boundary methods include the T-matrix, Born series approximations, and regular non-singular power series approximations.
Another method based on the theory of elastodynamics for multiple scatters is the addition theorem for spherical Bessel functions.13,14 The multiple scattering of elastic waves and corresponding dynamic stresses are analyzed for two-particle reinforced composite systems,13 and analytical expressions for the wave fields in each medium are obtained in terms of spherical Bessel functions of the first kind, Lagendre polynomials, and a translational coefficient of coordinate systems. The primary objective in this study was to analyze the multiple elastodynamic scattering of elastic waves and the dynamic stresses in two-particle reinforced metal matrix composites resulting from an incident plane longitudinal wave with unit amplitude propagating in the positive z direction. Due to the continuous boundary conditions around the spherical particles, the expansion mode coefficients are determined and an analytical expression for the dynamic stress concentration factor around the interfaces between particles and matrix using different parameters is analyzed. Results show that, when the frequencies are different, effects of the center to center distance between particles varies the dynamic stress values significantly, and only when this distance between particles is greater than a specific value can the dynamic stress values be ignored. Additionally, for different distances between particles, the effect of material properties of both particles and matrix materials on dynamic stress values is examined. Results illustrate that the effect of the distance between particles on dynamic stress concentration factors around the particle matrix interface is closely related to frequency of the incident wave. It was concluded that the dynamic stress concentration factors around the particles is strongly dependent on the incident wave frequency, material properties of both particles and matrix, and the center to center distance between particles that has the largest effect on the dynamic stress concentrations. In addition, the effect of material properties on dynamic stress concentrations also varies significantly with the center distance between particles.
The scattering of waves by two elastic spheres of the same material separated by a distance d between their origins in an elastic sold of a different material has also been analyzed14 using the addition theorem for spherical Bessel functions. For this problem, the incident wave, the scattered waves in the external domain of elastic spheres, and the transmitted waves within the spheres are all expanded in series forms of spherical Bessel functions of the first and third kinds along with Legendre polynomials defining local coordinates with origins at the center of the spheres of radius a using an addition theorem method similar to that proposed with the two-particle reinforced composite system.13 For any point in the external domain of the spheres, the total wave function is represented by the incident wave function and the scattered wave functions from both of the individual spheres. Coordinate transformations for each sphere are employed through the application of the addition theorem of the wave functions of spherical Bessel functions. The interface conditions at the surface of the spheres are obtained from the continuity of the radial and circumferential displacements and traction components in two dimensions that are defined by two potential functions. Therefore, each sphere has four boundary conditions and four unknown expansion coefficients giving a total of eight boundary conditions and eight unknown expansion coefficients for this problem. The expansion coefficients for the scattering waves are subsequently obtained by solving for the scattering cross section using the boundary conditions at each spherical interface that represents the ratio of energy associated with the incident and scattering waves per unit area. The authors then provide numerical examples for two aluminum spheres embedded in a germanium matrix. Solutions are obtained by first obtaining the scatter fields from each sphere, and then the total scatter field from both spheres is obtained using the addition theorem from which the scattering cross section is obtained. Numerical results illustrate that the interaction effects associated with the scattering cross section produce oscillations between the spheres that decreases with both increasing dimensionless distance D = d/a and increasing dimensionless wave number kLa. However, for both one and two spheres, the magnitude of the total scattering cross section increases with the dimensionless longitudinal wave number kLa.
More recently,15 the classical BEM method for elastodynamic scattering has been reformulated giving what the authors refer to as the energetic space–time weak BEM formulation (it is considered to be energetic due to the fact it takes advantage of some properties of the system’s energy) of the boundary integral equations. This new method has been employed to solve general two-dimensional elastodynamic scattering problems by open arcs, with vanishing initial and Dirichlet boundary conditions, and three-dimensional scalar wave propagation problems. However, elastodynamic solutions using this new method are currently significantly restricted to a very limited number of problems, and are not applicable for the problems considered in this paper. Current research associated with this new method is focusing on the development of an elastodynamic space–time double layer potential along with a hypersingular (strongly singular) boundary integral operator that could make this method applicable to the problems considered in this paper at some point.
As opposed to the acoustic and elastodynamic scattering problems previously discussed, the objective of this present study is to quantitatively investigate the interaction effects from both surface and field point displacements obtained through the rigorous formulation of the boundary integral representation for the surface and field point displacements associated with the elastodynamic scattering by two symmetric spherical cavities subjected to a plane, time-harmonic longitudinal displacement wave with a constant dimensionless wave number kLa over a range of dimensionless separation distances d/a and field point locations (as opposed to stress concentrations13 or total cross sections14 as a function dimensionless wave numbers, of which both of these results can be obtained from the surface and field point displacements obtained from the off-boundary BEM method developed in this study) through the use of a half-space off-boundary BEM in the frequency domain.16 If time domain results are required, a fast Fourier transform8 can subsequently be applied with the frequency domain results to construct the solution in the time domain. However, in this study, only frequency domain surface field point displacements are of interest and are used to obtain theoretical scattering signatures by frequency-averaging the scattered power obtained through the use of field point displacement values. These theoretical scattering signatures can then be employed for comparison with frequency domain scattering signatures obtained from NDE ultrasonic transducers for the identification of specific inhomogenities.1,5 The traction free cavity surfaces being integrated in this study are subjected to an incident, plane, time-harmonic longitudinal displacement wave with unit amplitude in a homogeneous, isotropic, and linear elastic solid material.
The first modification in this paper to the standard BEM is the implementation of the off-boundary technique.16 With this modification, the elements and source points are located on the cavity surface; however, the observation points (field points) are taken inside the cavity. Furthermore, the constant element approximation is employed; therefore, as with the standard BEM7,16,17 using the constant element approximation, the number of observation points is equal to the number of source points and elements. Thus, at each observation point, integral equations associated with each element source point are evaluated generating a set of linear complex algebraic equations that are solved using complex Gaussian elimination giving the constant complex displacements associated with each surface element. Once the surface displacements and/or tractions are known for each element, they can be used with the corresponding boundary integral representation to give the total and scattered displacement fields anywhere in the exterior domain V of the solid material. Using Hooke’s law with the displacement fields gives the stress and strain fields anywhere in the exterior domain V of the solid material. For illustrative purposes, and comparisons with existing results, we will only consider elastodynamic scattering by spherical cavities in this paper; however, this method can be extended for the scattering by both cavities and/or fixed rigid bodies16 of general shape (unfortunately, this method is not applicable to cracks due to the fact there is no interior domain for the observation points). It is shown for a single cavity subjected to a plane, time-harmonic longitudinal displacement wave16 that using the off-boundary technique with the BEM removes complicated nonintegrable singularities from the domain of integration along with eliminating ill conditioned problems encountered at fictitious eigenfrequencies8,16 of which both require the use of special computationally demanding procedures to obtain solutions. Therefore, the off-boundary BEM reduces computational resource requirements making the method more efficient and applicable for a larger number of problems.
The second modification to the standard BEM in this study is based on considering symmetry of the spherical cavities about the plane contained within the inertial Cartesian coordinate system shown in Fig. 1. With this modification, a half-space Green’s displacement tensor is constructed using the method of images18 by the superposition of two free-space Greens displacement tensors. Then, substituting the half-space Green’s displacement tensor into Hooke’s law gives the half-space double layer kernel tensor, in which both tensors satisfy the required boundary conditions on the plane of symmetry. The half-space double layer kernel tensor is then used with the boundary integral representation, where only one cavity requires integration to obtain the elastodynamic scattering solution for the cavity being integrated along with its exact geometrical image. This then reduces the number of required elements in half, which significantly reduces computational resource requirements.
Incident field, scattered field, and Cartesian coordinate system for two. spherical cavities symmetric about the plane.
Incident field, scattered field, and Cartesian coordinate system for two. spherical cavities symmetric about the plane.
Results using the off-boundary BEM are initially verified using only a single elastodynamic free-space double layer kernel tensor, and comparing the magnitudes of the complex radial and polar constant displacements on the surface of a single cavity with corresponding existing closed form analytical and constant element approximation BEM results,6,7 and are shown to be in good agreement. Next, solutions for both the free-space and half-space Green’s displacement tensors using their corresponding double layer kernel tensors are verified by comparing magnitudes of polar displacements on the inner and outer surfaces of two symmetric spherical cavities as shown in Fig. 1. These results are compared with existing results for the same problem using a generalized Born series approximation19 and are also shown to be in reasonable agreement. However, it must be noted that generalized Born series approximations19 are expected to lose accuracy with strong scatters, such as cavities,5,11 as considered in this study.
Finally, the main objective of this study is to quantitatively analyze the interaction effects associated with the scattering of elastic waves by two symmetric spherical cavities using the half space off-boundary BEM. Results are obtained for the magnitude of the complex back and forward scattered displacement fields in the X3 direction along the X3 axis as illustrated in Fig. 1 as a function of the dimensionless separation distance d/a, between the inner surfaces of the two symmetric spherical cavities with radius a.
BOUNDARY INTEGRAL REPRESENTATION FOR A DOMAIN OF INFINITE EXTENT
The symmetric spherical cavity scatterers, complex incident displacement field, and complex scattered displacement field are shown in Fig. 1. The surfaces of the two symmetric spherical cavity scatterers are denoted by S1 and S2. The void regions inside the scatterers are denoted by and , and the region outside the scatterers is a homogeneous, isotropic, and linear elastic solid material denoted by V. In the Fourier transformed frequency domain, the displacement field is time-harmonic; therefore, for simplicity, the time factor e−iωt will be omitted leaving only the spatial dependent terms and removing all time dependence. This off-boundary BEM is applicable for both cavities and fixed rigid inclusions;16 however, in this study, we only consider spherical cavities.
The total exterior complex displacement field at a field point x in V is given by
where is the complex incident displacement field in a homogeneous, isotropic, and linear elastic solid in the absence of cavities, and is the complex scattered displacement field, where the longitudinal and transverse components satisfy the elastodynamic equivalent of the Sommerfeld radiation condition.20 This radiation condition requires that at large distances from the sources of radiation (scatters), vanishes; however, does not. The integral representation for the total exterior displacement field is contingent on the use of a free-space elastodynamic fundamental solution (free-space Green’s displacement tensor)21
where , y is the source point, x is the field point, kL = ω/cL is the longitudinal wave number with wave velocity , kT = ω/cT is the transverse wave number with wave velocity , λ and μ are the Lamé elastic constants of the solid material in V, ρ is the density of the solid material in V, and δij is the Kronecker delta. The expression in Eq. (2) physically represents21 the elastodynamic displacement in the j direction at point y of a solid material due to a concentrated unit force acting in the i direction at a point x, and satisfies the displacement equation of motion17
where is the Dirac delta function.
Now, we consider the elastodynamic scattering by two symmetric spherical cavities (this formulation can also be used for a single cavity) whose surfaces are free of tractions. Using Eq. (2) with the dynamic Betti–Rayleigh reciprocal theorem20 along with elastodynamic equivalent of the Sommerfeld radiation condition20 allows the derivation of the boundary integral representation for the total displacement field, and is given by
where in this formulation are the unit outward normals on S1 and S2 into V, and the double layer kernel tensor in Eq. (4) is obtained by substituting the fundamental solution from Eq. (2) into Hooke’s law, giving21
BOUNDARY INTEGRAL REPRESENTATION FOR A HALF-SPACE
A substantial decrease in computer time and memory requirements can be attained for many elastodynamic problems that exhibit symmetry. This is achieved by constructing a half-space Green’s displacement tensor using the method of images18 by the superposition of two free-space Green’s displacement tensors. The half-space Green’s displacement and double layer kernel tensors are defined by and , respectively, for two symmetric spherical cavities, where only one cavity requires integration to obtain the elastodynamic scattering solution for the cavity being integrated along with its exact geometrical image. With this formulation, the half-space Green’s displacement and double layer kernel tensors must satisfy the displacement and shear stress boundary conditions on the plane of symmetry when y2 = 0 as defined by
and
From the superposition of two free-space Green’s displacement tensors, the half-space Green’s displacement tensor is given by
where the second free-space Green’s displacement tensor on the right is subtracted only when i = 2, and the vector is just a function of y. The half-space double layer kernel tensor is then obtained by substituting Eq. (7) into Eq. (5) such that only one cavity requires integration and reduces the integration surface area in half.
The boundary integral representation for the total displacement field for two symmetric spherical cavities with the spherical cavity located on the negative X2 axis being the exact geometrical image of the spherical cavity on the positive X2 axis, giving the boundary integral representation
where in this formulation, integration is done only over surface S1, and can only have non-zero components in the i = 1 and 3 directions.
NUMERICAL PROCEDURE
The boundary integral representations in Eqs. (4) and (8) are solved using the off-boundary BEM with the free-space , and half space double layer kernel tensors, respectively. Additionally, we define to represent either or depending on the specific problem. Numerical integration is achieved by discretizing the surface of S1 or surfaces S1 + S2 into M elements that are denoted by Sq (q = 1,2, … M). In this study, we use the constant element approximation7,17 for the unknown surface element displacements, in which four node two-dimensional elements are employed. However, higher order shape functions and elements can also be easily implemented.7 The constant surface elements used in this study consist of four corner nodes defined by yβ (β = 1–4) with shape functions
where −1 ≤ ξ ≤ 1 and −1 ≤ η ≤ 1, and ξ and η are the serendipity coordinates. An arbitrary source point position vector yq on Sq is defined by the shape functions from Eq. (9) as
where yβ is the position vector of the βth node of the qth element. Using the constant element approximation with the integrals in Eqs. (4) or (8) over the qth element gives
where J denotes the absolute value of the determinate of the Jacobian matrix
and is the unit outward normal vector of the qth element in the direction of the exterior V of the cavity surface and is given by
We utilize the off-boundary BEM technique16 to obtain values for the three constant complex displacement terms of each element. This is done by taking the off boundary field points x in Eqs. (4) or (8) as xα (α = 1, 2, … M), where xα lies inside the cavity under the center of each element, and integrating over all the elements (q = 1, 2, … M). This process generates the discretized version of the boundary integral representation for x ∈ Vc as
where are the complex surface displacements of the qth element, is the complex incident displacement field at xα, and
where x = xα at the αth off-boundary field point with the element source point taken on the qth element. The integrations in Eq. (15) are carried out over the serendipity coordinates ξ and η using 16 and 64 Gaussian quadrature integration points in which the latter is used for elements in which the field points xα correspond to the α = q source points [these are dominant diagonal terms, and with α = q, xα is close to ; therefore, more Gaussian quadrature integration points are generally required for the integration value to be accurate22]. The numerical values of the constant complex element displacements are then obtained by applying complex Gaussian elimination to the set of 3M complex linear algebraic equations given by Eq. (14).
After the surface displacements of all elements are obtained, they are substituted into discretized versions of Eqs. (4) or (8) to obtain values of the total exterior displacement fields at any desired field point x ∈ V. The discretized version of the integral representation for total exterior displacement field with is
where is the number of field points in V being evaluated, are the total exterior complex displacements at field points , are the incident complex displacements at the field points , are the cavity surface displacements at the qth element, and is the same as Eq. (15) except it is now evaluated with the field points corresponding to .
NUMERICAL RESULTS
Numerical results are presented that first validate the off-boundary BEM formulation used in this study using both the free-space and half-space double layer kernel tensors. First, a single spherical cavity is analyzed, and then two spherically symmetric spherical cavities are analyzed using both a free-space double layer kernel tensor, and a half-space double layer kernel tensor that represents a single spherical cavity and its exact geometrical image. Next, results are obtained that illustrate the interaction effects as function of the non-dimensional separation distance d/a between the two symmetric spherical cavities along the X3 axis as shown in Fig. 1. For all of the simulations, an incident, plane, time-harmonic longitudinal displacement wave traveling in the X3 direction with unit amplitude is employed and given by7,17
where e3 is a unit vector in the direction of the X3 coordinate axis as shown in Fig. 1. Additionally, in all simulations, the dimensionless longitudinal wave numbers used were either kLa = 1.0 or 1.3 with a Poisson’s ratio ν = 0.25, cavity radius a = 1, elastic modulus E = 1, and density ρ = 1.2. Therefore, the dimensionless longitudinal wave number is always equal to the angular frequency such that kLa = ω. All of the spherical cavities considered in this study were discretized using 288 four node rectangular surface elements with constant complex displacements and defined using the shape functions from Eq. (9) with 12 and 24 elements in the polar θ and azimuthal ϕ directions, respectively. The off-boundary field points xαwere taken at a radius of 0.9a under the center of each α = q elements. It was determined16 that the precise location of the field points under each of the elements inside the cavity is not extremely critical; however, ill conditioning may result if these points are taken too far away from the cavity surface. At the top and bottom of the spherical cavities, triangular elements are required and were constructed by using the four node elements with the shape functions in Eq. (9) and collapsing one side of the four node element.23 Therefore, considering that each pole has 1 node common to all triangular elements, the total number of nodes for the 288 element cavity mesh is 266.
To initially validate our implementation of the off-boundary BEM developed in this study, we first considered the elastodynamic scattering by a single spherical cavity located at the origin of the coordinate system shown in Fig. 1 and uses the free-space double layer kernel tensor in Eq. (5). For this problem, we subjected the cavity to an incident, plane, time-harmonic longitudinal wave given by Eq. (17) with a dimensionless wave number of kLa = 1.3. The magnitudes of the complex radial and polar surface displacements as a function of the polar angle θ in the frequency domain were compared with results for the same problem using the closed form analytical solution by Pao and Mow6 along with the standard constant element approximation BEM by Kitahara and Nakagawa,7 as shown in Figs. 2(a) and 2(b), respectively. Furthermore, results using the off-boundary BEM16 implemented in this study using Eq. (14) with the double layer kernel tensor from Eq. (5) for a single spherical cavity in a domain of infinite extent can be considered to provide reasonably good results for the magnitudes of both radial and polar complex displacements at the cavity surface.
Magnitude of the complex (a) radial and (b) polar displacements on the surface. of a single spherical cavity with KLa = 1.3.
Magnitude of the complex (a) radial and (b) polar displacements on the surface. of a single spherical cavity with KLa = 1.3.
The accuracy of the off-boundary16 and standard7 constant element approximation BEM techniques are evaluated using root mean-square deviation percentage error4 with respect to the closed form analytical complex displacement magnitude results6 given by
where for a single spherical cavity is the root mean-square deviation percentage error associated with an individual complex displacement magnitude component, n is the number of data points, is the corresponding individual complex displacement component magnitude components associated with the known closed form analytical solution,6 and is the individual complex displacement magnitude components associated with the BEM used. Due to the axial symmetry associated with a single cavity subjected to a plane, time-harmonic longitudinal displacement wave traveling in the X3 direction with unit amplitude, the magnitude of the azimuthal components of the surface displacements are , and there are only displacements in the radial and polar directions. For the off-boundary BEM,16 the root mean-square deviation percentage error of the magnitudes of the complex radial and polar displacements are = 3.52%, and = 2.68%, respectively. For the standard constant element BEM,7 the number of boundary elements used was only 110 giving root mean-square deviation percentage errors of the radial and polar complex displacement magnitudes of = 4.82%, and = 2.86%, respectively. Here, the root mean-square deviation percentage errors with the standard BEM using the constant element approximation is slightly greater than that using the off-boundary BEM with the constant element approximation. It is conjectured that the difference in errors with the standard BEM7 with the constant element approximation is due to several factors. First, only 110 constant elements were used with the standard BEM as opposed to 288 constant elements with the off-boundary BEM, second, with the off-boundary BEM, the free term7 (for a smooth surface) is not required as when x = y with the standard BEM, and third, a different Gaussian quadrature integration scheme was employed with the off-boundary BEM in this study.
It is shown7,17 that the most accurate values for radial and polar complex surface displacement magnitudes were obtained using only 30 quadratic isoparametric elements with the standard BEM (quadratic elements can also be used with the off-boundary technique provided observation points are placed under each individual node in the mesh) with the root mean-square deviation percentage errors of = 1.30%, and = 0.64%. However, the errors associated with the off-boundary constant element approximation can still be considered small enough to provide an acceptable approximation for the magnitudes of the complex displacements in the frequency domain. Additionally, it appears that these errors can be reduced with further mesh refinement in view of the error differences between the standard and off-boundary constant element BEM approximations.
Next, we validate the half-space double layer kernel tensor used with the half-space off-boundary BEM developed using the method of images.18 For this problem, we considered the elastodynamic scattering by two symmetric spherical cavities with a dimensionless inner surface separation distance of d/a = 1/5. In this case, we subjected both two symmetric spherical cavities and a single spherical cavity with its exact geometrical image to the incident, plane, time-harmonic longitudinal wave given by Eq. (17) with a dimensionless wave number of kLa = 1.0. For the case of two symmetric spherical cavities using the free-space double layer kernel tensor , we integrated over both cavity surfaces S1 and S2 with a total of M = 576 elements (giving a total of 3M = 1728 complex linear algebraic equations that are solved using complex Gaussian elimination), and for the case of a single spherical cavity and its exact geometrical image using the half-space double layer kernel tensor we integrated over S1 with a total of M = 288 elements. Results for the inner and outer cavity complex polar surface displacement magnitudes as a function of the polar angle θ in the frequency domain using these two off-boundary BEM methods are compared with results for the same problem obtained by Kitahara and Nakagawa19 using a generalized Born series approximation24,25 to model the elastodynamic scattering of the complex polar displacement magnitudes with interaction effects from two symmetric spherical cavities, as shown in Figs. 3(a) and 3(b), at azimuthal angles of ϕ = 270° and 90°, respectively. From these numerical results, it is observed that both the free-space (integrating over both symmetric spherical cavities) and half-space (integrating over a single cavity) double layer kernels used with the off-boundary BEM give essentially the exact same results for the magnitudes of the polar surface displacements as a function of the polar angle θ for two symmetric spherical cavities. Additionally, it is observed that both the free and half-space off-boundary BEM results are in reasonably good agreement with the magnitudes of the complex polar surface displacement magnitudes obtained from the generalized Born series approximation19 for both the inner and outer surfaces of the symmetric spherical cavities (the differences observed between the displacement magnitudes between both of the off-boundary BEM results, and the generalized Born series approximation results are not unexpected due to the fact the generalized Born series approximation is known to lose accuracy, especially with multiple strong scatters such as voids25). However, we can determine the accuracy of using the half-space double layer kernel tensor using Eq. (18) for the root mean-square deviation percentage error with respect to the free-space double layer kernel tensor, which was shown to give fairly accurate complex displacement magnitudes for the single cavity problem using the off-boundary technique with constant element displacements. Therefore, the root mean-square deviation percentage error obtained for the magnitude of the complex polar surface displacements as a function of the polar angle θ is = 0.29% and = 0.00% for the inner and outer cavity surfaces, respectively. Therefore, it is concluded that the half-space double layer kernel tensor gives extremely accurate results for the magnitude of the complex polar surface displacements as a function of the polar angle θ with respect to the free-space double layer kernel tensor using the off-boundary technique with the constant element approximation.
Magnitude of the complex polar surface displacements at the (a) inner cavity. surface (ϕ = 270° on S1), and (b) outer cavity surface (ϕ = 90° on S1) for two. symmetric spherical cavities with KLa = 1.0, and d = a/5.
Magnitude of the complex polar surface displacements at the (a) inner cavity. surface (ϕ = 270° on S1), and (b) outer cavity surface (ϕ = 90° on S1) for two. symmetric spherical cavities with KLa = 1.0, and d = a/5.
Based on the results from the previous spherical cavity elastodynamic scattering simulations, it was shown that both the free-space and half-space off-boundary BEM developed in this study provide acceptable approximate results for the magnitudes of the constant complex displacements of the elements at the surfaces of the cavities. Therefore, we now use the three components of the complex surface element displacement vectors obtained using the half-space off-boundary BEM to analyze the interaction effects from two symmetric spherical cavities in the frequency domain. The interaction effects are determined as a function of the dimensionless separation distance d/a between the inner surfaces of the two symmetric spherical cavities with a dimensionless wave number of kLa = 1.3, and subjected to an incident, plane, time-harmonic longitudinal displacement wave given by Eq. (17). The interaction effects are analyzed using the half-space off-boundary BEM formulation in Eq. (8) using the constant element displacements from Eqs. (8) and (14) with the discrete boundary integral representation given by Eq. (16). Here, the magnitudes of the complex scattered field displacements in the X3 direction along the X3 axis are analyzed.
In Fig. 4, we evaluate the back scattered complex displacement magnitudes with and without interaction effects [for no interaction effects we integrate over the single spherical cavity S1 using Eqs. (4) and (14) with the discrete boundary integral representation given by Eq. (16) and multiply the resulting values of by two] at locations, (a) −25a, (b) −50a, and (c) −75a, on the X3 axis. It must be noted that a minimum separation distance of d/a = 0.05 was used for the results in both Figs. 4 and 5; therefore, the cavity being integrated is never actually in contact with the plane of symmetry (the cavities never actually come in contact with each other), and problems were not encountered with the spherical cavities overlapping. The magnitude of the complex back scatter displacement field results in Fig. 4 indicate the displacement interaction effects produce constructive and destructive interference between the two symmetric spherical cavities as a function of d/a (analyzing the interaction effects associated with the magnitudes of the complex scattered field displacements at several external field points appears to be a novel result). This interference causes the magnitude of the displacement field to go in and out of phase, and produces the observed small amplitude oscillatory behavior about the magnitude of the complex displacement scatter field solution for that neglects interaction effects. Additionally, at all field point locations considered, the small amplitude oscillatory interference behavior decays significantly for d/a > 30. Furthermore, the maximum overall magnitude of the complex scattered displacement field (with and without interaction effects) along the negative X3 axis increases with decreasing distances away from the origin and also exhibits the least overall fluctuations in amplitude as a function of d/a.
Magnitude of the complex back-scattered displacement field in the X3. Direction as a function of the dimensionless cavity separation distance with and without interaction: (a) X3 = −25a, (b) X3 = −50a, and (c) X3 = −75a.
Magnitude of the complex back-scattered displacement field in the X3. Direction as a function of the dimensionless cavity separation distance with and without interaction: (a) X3 = −25a, (b) X3 = −50a, and (c) X3 = −75a.
Magnitude of the forward scattered displacement field in the X3 direction as a function of the dimensionless cavity separation distance with and without interaction: (a) X3 = 25a, (b) X3 = 50a, and (c)X3 = 75a.
Magnitude of the forward scattered displacement field in the X3 direction as a function of the dimensionless cavity separation distance with and without interaction: (a) X3 = 25a, (b) X3 = 50a, and (c)X3 = 75a.
In Fig. 5, we analyze the magnitude of the forward complex scattered displacement field magnitudes with and without interaction effects at locations, (a) 25a, (b) 50a, and (c) 75a, on the X3 axis. As with the back scattered results, the forward scattered results also exhibit interaction effects producing constructive and destructive interference between the two symmetric spherical cavities. This interference also causes the magnitude of the complex scattered displacement field to go in and out of phase and produces small amplitude oscillatory behavior about the magnitude of the complex scattered displacement field that neglects interaction effects. However, the forward scattered interference has slightly less amplitude than the back scattered interference, but also decays significantly for d/a > 30 Additionally, the maximum overall magnitude of the forward complex scattered displacement field along the positive X3 axis also increases with decreasing distances away from the origin along with exhibiting the largest overall changes in amplitude as a function of d/a.
CONCLUSIONS
The elastodynamic scattering of an incident, plane, time-harmonic longitudinal displacement wave given by Eq. (17) striking spherical cavities and generating scattered complex displacement fields are analyzed in this study. Figure 1 illustrates the solid material of infinite extent under consideration with two symmetric spherical cavities along with the complex incident and scattered displacement fields. For all of the problems considered, solutions for the complex scattered and total displacement fields are solved in the frequency domain that eliminates the time variable t. If time domain results are required, a fast Fourier transform can subsequently be applied8 with the frequency domain results to construct the solution in the time domain. In this study, we employ a rigorous approach for the formulation of elastodynamic scattering solutions through the use of the boundary integral representations given by Eqs. (4) and (8) using free and half-space double layer kernel tensors, respectively (this is a strong formulation of the problem as opposed to a weak formulation15). These double layer kernel tensors are obtained by substituting the Green’s displacement tensors into Hooke’s law given by Eq. (5).
The half-space double layer kernel tensor is obtained from the half-space Green’s displacement tensor developed using method of images18 and substituting it into Hooke’s law. It is then used with the boundary integral representation in Eq. (8) and only requires integration over the surface of the single cavity S1 to obtain the complex elastodynamic scattering surface displacements for the single cavity being integrated along with its exact geometrical image that reduces the total integration surface area in half. Furthermore, reducing the surface area of integration in half allows for obtaining complex surface and exterior field point displacements for elastodynamic scattering problems consisting of voids and/or fixed rigid bodies of arbitrary geometries that can be reasonably approximated as being symmetric with increasingly higher frequencies. Thus, by reducing the total surface area of integration provides a significant reduction of computational resource requirements allowing for the use of finer surface mesh discretizations with smaller boundary elements (mesh refinement). This, along with the construction of special optimized boundary element arrangements for the scattering objects provides the capability of representing rapidly changing spatial boundary variables along with the use of more element Gaussian quadrature integration points all of which increase the accuracy of the BEM solutions with increasing frequencies.8
The boundary integral representations for the surface and exterior field point displacements given by Eqs. (4) and (8) generally cannot be solved analytically, and require the use of numerical methods to obtain solutions. In this paper, the off-boundary BEM16 is employed to solve the discrete boundary integral representations in Eqs. (14) and (16) for the complex surface and exterior field point displacements, respectively. With this technique, the boundary elements and source points are located on the cavity surface; however, the field points are taken inside the cavities for obtaining the complex displacements of the surface elements using Eq. (14). By taking the field points inside the cavity eliminates both the problem with integrating the strongly singular and non-integrable double layer kernel tensors when x = y on the cavity surface as with the standard BIE,7,17 along with ill-conditioning problems that occur at fictitious eigenfrequencies8,16 of which both require the use of special computational demanding procedures to obtain solutions. Additionally, the constant element approximation is employed; therefore, as with the standard BEM,7 the number of field points is equal to the number of source points and elements. Thus, at each off-boundary field point, integral equations for a cavity’s interior associated with each element source point are evaluated generating a set of 3M complex linear algebraic equations that are solved using complex Gaussian elimination providing the complex constant surface displacements associated with each element. Once the complex surface displacements are known for each constant element, they are used with the appropriate boundary integral representation for the exterior domain V, giving the total and scattered complex displacement fields anywhere in the exterior domain V.
In this study, we initially verify that both our free-space and half-space off-boundary BEM formulations are able to provide accurate results. Therefore, we initially consider the elastodynamic scattering of an incident, plane, time-harmonic longitudinal given by Eq. (17), striking a single spherical cavity.7,17 In this case, we take the dimensionless wave number as kLa = 1.3 and use the free-space double layer kernel tensor boundary integral representation given by Eq. (4). Next, we discretize the cavity into 288 constant displacement surface elements, with the corresponding equal number of field points located at a radius of 0.9a under the center of each element. The precise location of the field points under each element does not appear to be critical; however, ill-conditioning may result if taken to far from the surface.16 Solving the system of resulting complex linear algebraic equations in Eq. (14) using complex Gaussian elimination gives the three resulting complex displacement components for each of the surface elements. Results for the magnitude of the complex radial and polar surface displacements as a function of the polar angle θ in the frequency domain from our off-boundary BEM are compared with existing results using two different methods6,7 and are shown to be in good agreement, as shown in Figs. 2(a) and 2(b). Furthermore, we compare both the off-boundary BEM surface displacement magnitude results with the corresponding closed form analytical solutions obtained by Pow and Mao6 using the root mean-square deviation percentage error of the magnitude of the complex radial and polar surface displacements that were calculated to be = 3.52%, and = 2.68%, respectively, and are considered to provide acceptable approximations for the elastodynamic cavity surface displacement magnitudes under consideration.
Next, we use both our free-space and half-pace double layer kernel tensors for the elastodynamic scattering by two symmetric spherical cavities giving the magnitudes of the complex polar surface displacements . We first use the free-space double layer kernel tensor and integrate over the two spherical cavity surfaces S1 + S2. We then use the half-space double layer kernel tensor that represents the single spherical cavity S1 and its exact geometrical image, and only requires integrating over the single spherical cavity S1. In both cases, the cavities are subjected to an incident, plane, time-harmonic longitudinal displacement wave given by Eq. (17). For this problem, the dimensionless wave number was taken as kLa = 1.0 with a separation distance between the inner cavity surfaces of d/a = 1/5. The same cavity discretization was used for each cavity with the field points located at 0.9a under the center of each constant element with both methods. The magnitude of the complex polar surface displacements as a function of the polar angle θ using both BEM solutions are compared with the solution to the same problem using a generalized Born series approximation19 for the inner and outer cavity surfaces, as shown in Figs. 3(a) and 3(b), and are also in agreement. However, the free-space and half-space off-boundary BEM tend to give slightly larger values for in most cases, but are still reasonably close in view of the fact Born series approximation results are known to lose accuracy with strong scatters, such as voids.25 Additionally, the root mean-square deviation percentage error obtained for the magnitude of the complex polar surface displacements as a function of the polar angle θ is = 0.29% and = 0.00% for the inner and outer cavity surfaces, respectively. Therefore, it is concluded that the half-space double layer kernel tensor gives extremely accurate results for the magnitude of the complex polar surface displacements as a function of the polar angle θ with respect to the free-pace double layer kernel tensor using the off-boundary technique with the constant element approximation.
However, the main objective of this study was to specifically analyze the interaction effects associated with the magnitudes of the exterior field point complex displacements along the positive and negative X3 axis produced by an incident time-harmonic plane longitudinal displacement wave given by Eq. (17). Using the half-space double layer kernel tensor represents the single spherical cavity S1 and its exact geometrical image S2 in the boundary integral representation for the complex surface displacements for both symmetric spherical cavities and are a function of the dimensionless separation distance d/a and the field point locations on the X3 axis in between the inner cavity surfaces at several locations. For these simulations, the dimensionless wave number was taken as kLa = 1.3, and the exterior field points were taken along the positive and negative X3 axis at distances of 25a, 50a, and 75a. Results for the magnitudes of the complex back and forward displacement scatter fields are shown in Figs. 4 and 5, respectively. In both cases, the interaction effects generate constructive and destructive interference between the two symmetric spherical cavities causing the displacement field to go in and out of phase, producing a small amplitude oscillatory behavior about the magnitude of the corresponding displacement scatter field solutions that neglect interaction effects. From Figs. 4 and 5, it is observed that the back scattered interference has a slightly larger amplitude than the forward scattered interference. Furthermore, the oscillatory interference effects are observed to decay significantly for d/a > 30, and the maximum overall magnitude of the far field displacements along both the positive and negative X3 axis increase with decreasing distances from the origin along with exhibiting the largest overall fluctuations in amplitude as a function of d/a.
Of course, the assumption that both spherical cavities are exactly symmetric is an over simplification of any naturally occurring configurations of two cavities in a solid material. However, even in view of this idealization, the magnitudes of the complex displacement scatter fields on the X3 axis obtained from the boundary integral representation given by Eq. (8) provide significant insight into what can be expected from experimental measurements obtained from ultrasonic sensors used with NDE. As observed from the magnitudes of the complex back and forward scattered displacements , the elastodynamic interaction effects from the two symmetric spherical cavities are essentially an order of magnitude less than the effects from the magnitudes of the scatted displacements without interaction effects. Additionally, the interaction effects are shown to attenuate with increasing separation distances between the inner spherical cavity surfaces, especially for d/a > 30 at all field point locations along the X3 axis. This result is also consistent with multiple scattering results obtained using indirect BEM methods that include T-matrix,10–12 Born series,5,10,11 generalized Born series,18,24,25 and far-field25 approximations. Therefore, due to the small amplitude interaction effects observed in this study for spherical cavities that are considered to be strong scatters (i.e., when there is a very high impedance gradient between the interior and exterior of the scatters) indicates that a reasonable approximation may be obtained for many other multiple scatterer problems (except in some cases in which the scattering objects are extremely close) by just summing the scatter field results from each individual scatter and completely neglecting any interaction effects and applying appropriate position dependent phase corrections to each individual scatterer.25 Some of these other problems include a subsurface cavity with a free planar boundary,11 along with numerous other non-cavity weak scattering problems. Furthermore, this author is unaware of any previous use of a half-space double layer kernel tensor (obtained from a half-space Green’s displacement tensor) employed with the off-boundary BEM derived from the rigorous derivation of a boundary integral representation for the complex displacements on the cavity surface and at exterior field points. Additionally, these results also appear to be relevant to numerous other problems associated with both quantum and Newtonian mechanics for time-harmonic wave scattering problems.
ACKNOWLEDGMENTS
The author greatly appreciates his opportunity to work on this subject under the direction of Walter P Murphy Professor Jan D. Achenbach at the Center for Quality Engineering and Failure Prevention, the Technological Institute, Northwestern University, Evanston, Illinois.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Thomas L. Warren: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.