Abstract
A volume integral equation method (VIEM) is applied for the effective analysis of elastic wave scattering problems in unbounded solids containing general anisotropic inclusions. It should be noted that this numerical method does not require use of Green's function for anisotropic inclusions to solve this class of problems since only Green's function for the unbounded isotropic matrix is necessary for the analysis. This new method can also be applied to general two-dimensional elastodynamic problems involving arbitrary shapes and numbers of anisotropic inclusions. A detailed analysis of SH wave scattering problems is developed for an unbounded isotropic matrix containing multiple orthotropic elliptical inclusions. Numerical results are presented for the displacement fields at the interfaces of the inclusions in a broad frequency range of practical interest. Through the analysis of plane elastodynamic problems in an unbounded isotropic matrix with multiple orthotropic elliptical inclusions, it is established that this new method is very accurate and effective for solving plane elastic problems in unbounded solids containing general anisotropic inclusions of arbitrary shapes.
1. Introduction
A micrograph of the cross-section of a phosphate glass fiber/polymer composite is shown in Figure 1 [1]. Figure 1 indicates that the fibers are close to ellipses and the major axis of the elliptical fibers is not aligned in any one direction. Furthermore, the micrograph cross section of the glass fiber composite in Figure 2 [2] shows that the fibers are triangular in shape.

Phosphate glass fiber/polymer composite cross section [1]. Used with permission from Journal of Materials Science: Materials in Medicine.

Micrograph cross-section of a unidirectional 60 μm triangular glass fiber composite (V f = 0.5) [2]. Used with permission from Composites Science and Technology.
A number of analytical techniques are available for solving stress analysis of isotropic inclusion problems when the geometry of the inclusions is simple (i.e., cylindrical, spherical, or ellipsoidal) and when they are well separated [3–6]. However, these approaches cannot be applied to more general problems where the inclusions are both anisotropic and arbitrary in shape when their concentration is high. Thus, analysis of elastic wave scattering problems in heterogeneous solids often requires the use of numerical techniques based on the finite element method (FEM) or boundary element method (BIEM). Unfortunately, both methods encounter limitations in dealing with elastic wave scattering problems involving multiple anisotropic inclusions of arbitrary shapes. However, it has been demonstrated that a newly developed numerical method based on a volume integral equation formulation can overcome such difficulties in solving this class of inclusion problems [7, 8]. In contrast to the conventional boundary integral equation method (BIEM), where the infinite medium Green's functions for both the matrix material and the anisotropic inclusion material are needed, the present method does not require the latter. Since elastodynamic Green's functions for anisotropic media are extremely difficult to calculate, the present method offers a definite advantage over the boundary integral equation method. In addition, the VIEM is not sensitive to the geometry or concentration of the inclusions. Moreover, in contrast to the finite element method (FEM), where the full domain needs to be discretized, the VIEM requires discretization of the inclusions only.
Therefore, in order to investigate the influence of orthotropic elliptical inclusions on the interfacial field, a detailed analysis of the displacement field at the interface between the matrix and the central inclusion is first carried out for an unbounded isotropic matrix containing five orthotropic elliptical inclusions. The major axis of each elliptical inclusion is 0° or 90° with the x-axis and the aspect ratio is assumed to be 1.0 (circular) and 0.8. The inclusion separation distance (s) is assumed to be 3.24a and 2.51a (a: the radius of each circular inclusion or the major radius of each elliptical inclusion). It should be noted that the central inclusion is selected since the interaction effects of the multiple inclusions are likely to be highest in value there. The incident wave represents SH waves propagating parallel to the x-axis with particle motion being parallel to the z-axis.
It is demonstrated that this new method is very accurate and effective for solving plane elastodynamic problems in unbounded solids containing general anisotropic inclusions of arbitrary shapes. In forthcoming papers, in order to accurately predict the failure and damage mechanisms in real composites, the formulations developed herein will be used to calculate other quantities of practical interest in realistic models of composites containing multiple strong anisotropic heterogeneities [9].
2. Volume Integral Equation Method (VIEM)
The geometry of the general elastodynamic problem considered here is shown in Figure 3 where an unbounded isotropic elastic solid containing a number of isotropic or anisotropic inclusions of arbitrary shape are subjected to prescribed dynamic loading at infinity. The symbols

Geometry of the general elastodynamic problem.
Let u
m
o
(
It has been shown in Mal and Knopoff [10] that the elastodynamic displacement, u
m
(
where the integral is over the domain occupied by the inclusions, δρ and δc
ijkl
represent contrasts in the density and elastic tensor of the inclusions and the matrix (i.e., δρ = ρ(1) – ρ(2) and δc
ijkl
= c
ijkl
(1) – c
ijkl
(2)), and g
i
m
(
If
3. Scattering of SH Waves
3.1. Scattering of SH Waves in an Unbounded Isotropic Matrix Containing Multiple Orthotropic Circular Inclusions
We first consider multiple orthotropic circular inclusions in an unbounded isotropic matrix. The incident wave is represented as an SH wave with particle motion along the z-axis (see Figure 4). Let the coordinate axes x1(x), x2(y), x3(z) be taken parallel to the symmetry axes of the orthotropic material, let ρ1 and c44, c55 denote the density and elastic constants, respectively, of the orthotropic inclusion, and let ρ2 and μ2 denote the density and shear modulus of the isotropic matrix, respectively.

SH wave interaction with five orthotropic elliptical inclusions in an unbounded isotropic matrix.
For SH waves, the volume integral equation (1) reduces to
where u3(
In (2), g
i
m
(
where r = |
Finite element discretization of the inclusions in (2) results in a system of linear algebraic equations for the unknown nodal displacements inside the orthotropic inclusions. Once the displacement field, u3(
It should be noted that the singularities in VIEM are weaker (integrable) than those in BIEM, where they are of the Cauchy type. In response, we have used the direct integration scheme introduced by Cerrolaza and Alarcon [13], Li et al. [14], and Lu and Ye [15] after suitable modifications to address these singularity concerns. A description of the modified method used in the discretization of the volume integral equation is given by Lee and Mal [7, 16].
In the SH wave case, the incident wave is given as
where k2 is the S wavenumber in the matrix material. The elastic constants for the materials of the isotropic matrix and the orthotropic inclusions are listed in Table 1. Two different elastic constants for the orthotropic inclusions are considered: for “orthotropic inclusion #1,” the density values, c44 and c55 in the orthotropic inclusions, are smaller than those in the matrix, and for “orthotropic inclusion #2,” the density values, c44 and c55 in the orthotropic inclusions, are greater than those in the matrix. The calculations are carried out for three different normalized wavenumbers (i.e., k2a = 3.14, 5.03, and 7.85). The ratios of the corresponding wavelengths to the radius of each circular inclusion are approximated to be between 0.8 and 2.0. This represents a useful range of intermediate frequencies for dynamic loading and ultrasonic testing.
Material properties of the isotropic matrix and the orthotropic inclusions.
Figure 5 shows a typical discretized model used in the volume integral equation method [17]. Standard eight-node quadrilateral and six-node triangular elements were used. The number of elements in each circular inclusion used in the VIEM was 3,600 and was determined based on a convergence test. Figure 6 shows real and imaginary parts of the displacement component in the z-axis, u3, at three separate interfaces: (a) for the single orthotropic circular inclusion (#1), (b) for the central inclusion with five orthotropic circular inclusions (#1) where s = 3.24a, and (c) for the central inclusion with five orthotropic circular inclusions (#1) where s = 2.51a, all of which use the VIEM for three different frequencies (i.e., k2a = 3.14, 5.03, and 7.85). Figure 7 shows real and imaginary parts of the displacement component in the z-axis, u3, at three separate interfaces: (a) for the single orthotropic circular inclusion (#2), (b) for the central inclusion with five orthotropic circular inclusions (#2) where s = 3.24a, and (c) for the central inclusion with five orthotropic circular inclusions (#2) where s = 2.51a, all of which use the VIEM for three different frequencies (i.e., k2a = 3.14, 5.03, and 7.85).

A typical discretized model in the volume integral equation method.

Real and imaginary parts of the displacement component in the z-axis, u3, at the interface for: (a) the single orthotropic circular inclusion (#1), (b) the central inclusion with five orthotropic circular inclusions (#1) where s = 3.24a, and (c) the central inclusion with five orthotropic circular inclusions (#1) where s = 2.51a. The frequencies used are k2a = 3.14, 5.03, and 7.85 where k2 is the S wavenumber in the unbounded isotropic matrix and a is the radius of the circular inclusion. The inclusion cited represents orthotropic inclusion #1 in Table 1.

Real and imaginary parts of the displacement component in the z-axis, u3, at the interface for: (a) the single orthotropic circular inclusion (#2), (b) the central inclusion with five orthotropic circular inclusions (#2) where s = 3.24a, and (c) the central inclusion with five orthotropic circular inclusions (#2) where s = 2.51a. The frequencies used are k2a = 3.14, 5.03, and 7.85 where k2 is the S wavenumber in the unbounded isotropic matrix and a is the radius of the circular inclusion. The inclusion cited represents orthotropic inclusion #2 in Table 1.
3.2. Verification of the Viem Solution
To the best of the authors' knowledge, neither an analytical solution to this problem nor a closed form solution for the two-dimensional time-harmonic elastodynamic Green's function for the orthotropic material is presently available in the literature. Thus, in order to compare the calculated results with available analytical solutions, simple packing sequences (hexagonal and square) and isotropic circular inclusions in an unbounded isotropic elastic solid are considered (see Figure 8). Let ρ1 and μ1 denote the density and elastic constants, respectively, of the isotropic inclusion and let ρ2 and μ2 denote the density and shear modulus of the isotropic matrix, respectively. The material properties used are those for typical graphite/epoxy composites and are given in Table 2.
Material properties of graphite/epoxy composites.

SH wave interaction with multiple isotropic inclusions.
For SH waves, the volume integral equation (1) reduces to
where u3(
In (5), g
i
m
(
The incident wave is given as
where k2 is the S wavenumber in the matrix material.
Figure 9 shows discretized models used in the VIEM. The standard eight-node quadrilateral and six-node triangular elements were used in the VIEM. The total number of elements for each inclusion used in the VIEM was 256 and was determined based on a convergence test. Tables 3 and 4 show calculated average strains in the central isotropic fiber for real and imaginary parts as determined by an analytical method using the generalized self-consistent model (GSCM) [5] and the volume integral equation method for both (i) the square packing of 9 and 25 isotropic circular inclusions and (ii) the hexagonal packing of 7 and 19 isotropic circular inclusions. The material used was standard graphite/epoxy composites with an inclusion volume fraction of 0.6 and a frequency of 10 MHz. The percentage differences in the two sets of results were less than 1% in all cases [7].
Calculated average strains in the central fiber for the real parts using both analytical and volume integral equation methods. The material is graphite/epoxy with a volume concentration of 0.6 and a frequency of 10 MHz.
Calculated average strains in the central fiber for the imaginary parts using both analytical and volume integral equation methods. The material is graphite/epoxy with a volume concentration of 0.6 and a frequency of 10 MHz.

Discretized models for multiple inclusion problems used in the volume integral equation method.
3.3. Scattering of SH Waves in an Unbounded Isotropic Matrix Containing Multiple Orthotropic Elliptical Inclusions
We next consider multiple orthotropic elliptical inclusions in an unbounded isotropic matrix (see Figure 4). For elliptical cylindrical inclusions, the value of the aspect ratio for each elliptical inclusion was chosen to be 0.8 (see Figures (4) and (5)). Two different types of elliptical inclusions are considered where (1) the major axis of each inclusion coincides with the x axis (see Figure 4) and (2) the major axis of each inclusion coincides with the y axis (see Figure 5). The general features of the volume integral equation method for multiple orthotropic elliptical inclusions are similar to those for multiple orthotropic circular inclusions.
In the SH wave case, the incident wave is given as
where k2 is the S wavenumber in the matrix material. The elastic constants for the materials of the isotropic matrix and the orthotropic inclusions are listed in Table 1. The calculations are carried out for three different normalized wavenumbers (i.e., k2a = 3.14, 5.03, and 7.85). The ratios of the corresponding wavelengths to the major radius of each elliptical inclusion are approximated to be between 0.8 and 2.0. This is a useful range of intermediate frequencies for dynamic loading and ultrasonic testing.
Figure 5 shows a typical discretized model used in the volume integral equation method [17]. Standard eight-node quadrilateral and six-node triangular elements were used. The number of elements in each elliptical inclusion used in the VIEM was 3,600 and was determined based on a convergence test.
When the major axis of each elliptical inclusion is parallel to the x-axis, Figure 10 shows real and imaginary parts of the displacement component in the z-axis, u3, at three separate interfaces: (a) for the single orthotropic elliptical inclusion (#1), (b) for the central inclusion with five orthotropic elliptical inclusions (#1) where s = 3.24a, and (c) for the central inclusion with five orthotropic elliptical inclusions (#1) where s = 2.51a, all of which use the volume integral equation method for three different frequencies (i. e., k2a = 3.14, 5.03, and 7.85). Figure 11 shows real and imaginary parts of the displacement component in the z-axis, u3, at three separate interfaces: (a) for the single orthotropic elliptical inclusion (#2), (b) for the central inclusion with five orthotropic elliptical inclusions (#2) where s = 3.24a, and (c) for the central inclusion with five orthotropic elliptical inclusions (#2) where s = 2.51a, all of which use the VIEM for three different frequencies (i.e., k2a = 3.14, 5.03, and 7.85).

Real and imaginary parts of the displacement component in the z-axis, u3, at the interface for: (a) the single orthotropic elliptical inclusion (#1), (b) the central inclusion with five orthotropic elliptical inclusions (#1) where s = 3.24a, and (c) the central inclusion with five orthotropic elliptical inclusions (#1) where s = 2.51a. The frequencies used are k2a = 3.14, 5.03, and 7.85 where k2 is the S wavenumber in the unbounded isotropic matrix and a is the major radius of the elliptical inclusion. The inclusion cited represents orthotropic inclusion #1 in Table 1.

Real and imaginary parts of the displacement component in the z-axis, u3, at the interface for: (a) the single orthotropic elliptical inclusion, (b) the central inclusion with five orthotropic elliptical inclusions (#2) where s = 3.24a, and (c) the central inclusion with five orthotropic elliptical inclusions (#2) where s = 2.51a. The frequencies used are k2a = 3.14, 5.03, and 7.85 where k2 is the S wavenumber in the unbounded isotropic matrix and a is the major radius of the elliptical inclusion. The inclusion cited represents orthotropic inclusion #2 in Table 1.
When the major axis of each elliptical inclusion is parallel to the y-axis, Figure 12 shows real and imaginary parts of the displacement component in the z-axis, u3, at three separate interfaces: (a) for the single orthotropic elliptical inclusion (#1), (b) for the central inclusion with five orthotropic elliptical inclusions (#1) where s = 3.24a, and (c) for the central inclusion with five orthotropic elliptical inclusions (#1) where s = 2.51a, all of which use the VIEM for three different frequencies (i.e., k2a = 3.14, 5.03, and 7.85). Figure 13 shows real and imaginary parts of the displacement component in the z-axis, u3, at three separate interfaces: (a) for the single orthotropic elliptical inclusion (#2), (b) for the central inclusion with five orthotropic elliptical inclusions (#2) where s = 3.24a, and (c) for the central inclusion with five orthotropic elliptical inclusions (#2) where s = 2.51a, all of which use the VIEM for three different frequencies.

Real and imaginary parts of the displacement component in the z-axis, u3, at the interface for: (a) the single orthotropic elliptical inclusion, (b) the central inclusion with five orthotropic elliptical inclusions (#1) where s = 3.24a, and (c) the central inclusion with five orthotropic elliptical inclusions (#1) where s = 2.51a. The frequencies used are k2a = 3.14, 5.03, and 7.85 where k2 is the S wavenumber in the unbounded isotropic matrix and a is the major radius of the elliptical inclusion. The inclusion cited represents orthotropic inclusion #1 in Table 1.

Real and imaginary parts of the displacement component in the z-axis, u3, at the interface for: (a) the single orthotropic elliptical inclusion, (b) the central inclusion with five orthotropic elliptical inclusions (#2) where s = 3.24a, and (c) the central inclusion with five orthotropic elliptical inclusions (#2) where s = 2.51a. The frequencies used are k2a = 3.14, 5.03, and 7.85 where k2 is the S wavenumber in the unbounded isotropic matrix and a is the major radius of the elliptical inclusion. The inclusion cited represents orthotropic inclusion #2 in Table 1.
For different shapes of orthotropic elliptical inclusions, the interaction effect of the inclusions on the displacement component in the z-axis, u3, at the interface of the central inclusion appeared to be significant for orthotropic elliptical inclusion #1. In this case, the density values, c44 and c55, in the inclusions were determined to be smaller than those in the matrix. However, the interaction effect of the inclusions appeared to be minimal for orthotropic elliptical inclusion #2. In that case, the density values, c44 and c55, in the inclusions were found to be greater than those in the matrix.
4. Concluding Remarks
The volume integral equation method was applied to the solution of SH wave scattering problems with multiple orthotropic inclusions involving variations in orientation angle and inclusion separation distances for three different frequencies in an unbounded isotropic matrix. The aspect ratio of each elliptical inclusion was given as 1.0 (circular) and 0.8. It was determined that the displacement component in the z-axis, u3, at the interface of the central inclusion for orthotropic elliptical inclusion #1 (soft inclusion) differed significantly. However, for orthotropic elliptical inclusion #2 (hard inclusion), the interaction effect of orthotropic elliptical inclusions on the displacement component, u3 (in the z-axis), at the interface of the central inclusion appeared to be minimal.
The main advantage of this technique over those based on finite elements is that it requires discretization of the inclusions only as opposed to having to discretize the entire domain. This approach is similar to the boundary integral equation method except for the presence of the volume integral over the inclusions instead of the surface integrals over the two sides of the interface. It should be noted that if the medium contains a small number of isotropic inclusions, this method may not be as advantageous when compared to the BIEM. However, in the presence of multiple nonsmooth inclusions, the BIEM numerical treatment becomes noticeably cumbersome. Specifically, in elastodynamic problems involving multiple anisotropic inclusions, BIEM numerical treatment becomes extremely difficult since closed-form expressions for elastodynamic Green's functions for anisotropic media are currently not available. However, the volume integral equation method does not require the use of Green's functions for anisotropic inclusions. Furthermore, since standard finite elements are used in the VIEM, it is easier and more convenient to handle multiple nonsmooth anisotropic inclusions. Therefore, the formulations developed in this paper in conjunction with other useful methods [18, 19] can be used to calculate the dynamic stress intensity factors and other quantities of practical interest in realistic models of materials containing strong heterogeneities.
Through the analysis of plane elastodynamic problems in an unbounded isotropic matrix with multiple orthotropic elliptical inclusions, it is established that this new method is very accurate and effective for solving plane elastodynamic problems in unbounded solids containing general anisotropic inclusions of arbitrary shapes. In forthcoming papers, the formulations developed in this research will be used to calculate other quantities of practical interest in realistic models of composites containing multiple strong anisotropic heterogeneities [9].
Footnotes
Appendix
For three-dimensional elastodynamic problems, u3(
where u1(
where
Acknowledgments
This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology (Grant no. 2012-0007636) and Business for Cooperative R&D between Industry, Academy, and Research Institute funded by Korea Small and Medium Business Administration in 2013 (Grant no. C0096465). The authors would also like to express their sincere appreciation to Professor Jason Buschman for his excellent English proofreading and to acknowledge the support from Korea Institute of Science and Technology Information (KISTI) supercomputing center through the strategic support Program for the supercomputing application research (no. KSC-2012-C1-06).
