Research Article  Open Access
Jungki Lee, Hyechang Lee, Hogwan Jeong, "Multiple Scattering Using Parallel Volume Integral Equation Method: Interaction of SH Waves with Multiple Multilayered Anisotropic Elliptical Inclusions", Mathematical Problems in Engineering, vol. 2015, Article ID 809320, 48 pages, 2015. https://doi.org/10.1155/2015/809320
Multiple Scattering Using Parallel Volume Integral Equation Method: Interaction of SH Waves with Multiple Multilayered Anisotropic Elliptical Inclusions
Abstract
The parallel volume integral equation method (PVIEM) is applied for the analysis of elastic wave scattering problems in an unbounded isotropic solid containing multiple multilayered anisotropic elliptical inclusions. This recently developed numerical method does not require the use of Green’s function for the multilayered anisotropic inclusions; only Green’s function for the unbounded isotropic matrix is needed. This method can also be applied to solve general two and threedimensional elastodynamic problems involving inhomogeneous and/or multilayered anisotropic inclusions whose shape and number are arbitrary. A detailed analysis of the SH wave scattering is presented for multiple triplelayered orthotropic elliptical inclusions. Numerical results are presented for the displacement fields at the interfaces for square and hexagonal packing arrays of triplelayered elliptical inclusions in a broad frequency range of practical interest. It is necessary to use standard parallel programming, such as MPI (message passing interface), to speed up computation in the volume integral equation method (VIEM). Parallel volume integral equation method as a pioneer of numerical analysis enables us to investigate the effects of single/multiple scattering, fiber packing type, fiber volume fraction, single/multiple layer(s), multilayer’s shape and geometry, isotropy/anisotropy, and softness/hardness of the multiple multilayered anisotropic elliptical inclusions on displacements at the interfaces of the inclusions.
1. Introduction
A number of analytical techniques are available for the solution of isotropic inclusion problems when geometries of the inclusions are simple (e.g., cylindrical, spherical, or ellipsoidal) and when they are well separated [1, 2]. However, none of these techniques can be applied to more general problems, where the inclusions are multilayered, anisotropic, and arbitrary in shape, particularly when their concentration is high. For example, a micrographic cross section of a phosphate glass fiber/polymer composite is shown in Figure 1 [3]. 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.
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). However, both methods have limitations in dealing with elastic wave scattering problems in unbounded media containing anisotropic and/or heterogeneous inclusions of arbitrary shapes. It has been demonstrated that a numerical method based on the volume integral equation formulation can overcome such difficulties in solving this class of inclusion problems [2, 4–8]. In contrast to the conventional boundary integral equation method (BIEM), where infinite medium Green’s functions for both the matrix and the multilayered anisotropic inclusion are needed, the volume integral equation method (VIEM) does not require the use of Green’s function for the multilayered anisotropic inclusion. Since elastodynamic Green’s functions for multilayered anisotropic media are extremely difficult to calculate, the volume integral equation method offers a clear advantage over the boundary integral equation method. In addition, the VIEM is not sensitive to the geometry and concentration of the inclusions, and in contrast to the finite element method (FEM), where the full domain needs to be discretized, the VIEM requires discretization of the inclusions only.
The volume integral equation method has been proposed as a numerical solution scheme for twodimensional plane wave scattering problems with multiple isotropic inclusions (P, SV, and SH waves) [4], an orthotropic inclusion and a void (P and SV waves) [5], and five orthotropic inclusions (SH waves) [6].
The matrix and the fibers in composites are usually made of isotropic material. However, in order to satisfy advanced composites, some of the constituents can be anisotropic. As an example, in SiC/Ti metal matrix composites, the matrix is nearly isotropic, but the SiC fiber has strong anisotropy and a multilayered structure including an interphase and a core. As another example, Figure 2 shows schematically a cross section of a fiber containing multilayers [9].
In order to more quickly facilitate VIEM computations, parallel volume integral equation method (PVIEM) with message passing interface (MPI) and optimization [10, 11] using SUN computer systems (Tachyon system) from Korea Institute of Science and Technology Information (KISTI) has been developed recently. Thus, the PVIEM allows for a substantial speedup of VIEM calculations.
In this paper, the scattering problem of elastic waves by multiple threelayered orthotropic elliptical inclusions embedded in an unbounded isotropic elastic solid is solved using the parallel volume integral equation method (PVIEM). In order to investigate the influence of multiple threelayered orthotropic elliptical inclusions on the interfacial field in composites, a detailed analysis of the displacement field at the interface between the matrix and the outermost layer and the outermost layer and the middle layer of the central inclusion is first carried out for an unbounded isotropic matrix containing square and hexagonal packing arrays of threelayered orthotropic elliptical inclusions. The major axis of each elliptical inclusion is assumed to be oriented at 0° or 90° with axis and the aspect ratio is assumed to be 0.75. The fiber volume fraction is assumed to be 0.30 and 0.50. Here, the fiber volume fractions of the elliptical inclusions are defined as those of the circular inclusions circumscribing the elliptical inclusions [12, 13]. It should be noted that the central inclusion is selected since the interaction effects of multiple threelayered orthotropic elliptical inclusions are likely to be the highest there. The incident waves are assumed to be timeharmonic SH type waves propagating parallel to axis with particle motion parallel to axis. The results are presented for a broad range of frequencies.
Parallel volume integral equation method as a pioneer of numerical analysis makes it possible to investigate the effects of single/multiple scattering, fiber packing type, fiber volume fraction, single/multiple layer(s), multilayer’s shape and geometry, isotropy/anisotropy, and softness/hardness of the multiple multilayered anisotropic elliptical inclusions on displacements at the interfaces of the inclusions. Therefore, a precise knowledge of the motion at the interfaces of multiple multilayered anisotropic elliptical inclusions under remote dynamic loading can be used as one of the most important influential factors on predicting the failure mechanism in composites.
It should be noted that, to the best of the authors’ knowledge, the solution of the elastodynamic problem considered in this paper is not available in the literature [2, 14–19].
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 is subjected to prescribed dynamic loading at infinity. In Figure 3, and represent the volume and surface of the inclusion, respectively, and is the outward unit normal to , while and represent the infinite volume and surface, respectively. The symbols and denote the density and the elastic stiffness tensor of the inclusion. and denote the density and the elastic stiffness tensor of the unbounded matrix material. The matrix is assumed to be homogeneous and isotropic so that is a constant isotropic tensor, while can be arbitrary; that is, the inclusions may, in general, be inhomogeneous and anisotropic. The interfaces between the inclusions and the matrix are assumed to be perfectly bonded ensuring continuity of the displacement and stress vectors.
Let denote the th component of the displacement vector due to the incident field at in the absence of the inclusions. Additionally, let denote the same in the presence of the inclusions, where is the circular frequency of the waves. In what follows, the common time factor and the explicit dependence on for all field quantities will be suppressed.
It has been shown in Mal and Knopoff [20] that the elastodynamic displacement, , in the composite satisfies the volume integral equationwhere the integral is over the domain occupied by the inclusions, and represent contrasts in the density and elastic tensor of the inclusions and the matrix (i.e., and ), and is elastodynamic Green’s function for the unbounded homogeneous matrix material; that is, represents the th component of the displacement at due to unit concentrated force, , at in the th direction. In (1), the summation convention and comma notation have been used and the differentiations are with respect to the integration variable, . It should also be noted that the integrand is nonzero only within the inclusions, since outside the inclusions.
If lies in the domain occupied by the inclusions, then (1) is an integrodifferential equation for the unknown displacement vector within the inclusions, which can, in principle, be determined through the solution of (1). It is extremely difficult if not impossible to solve this coupled system of equations analytically even in the presence of a single inclusion of arbitrary shape. An algorithm based on the discretization of (1) was developed by Lee and Mal [4] to calculate numerically the unknown displacement vector by discretizing the inclusions using standard finite elements. Once within the inclusions is determined, the displacement field outside the inclusions can be obtained from (1) by evaluating the integral, and the stress field within and outside the inclusions can also be determined without any difficulty. Details of the numerical treatment of (1) for plane elastodynamic problems can be found in [4–6]. Further mathematical formulation of the elastostatic volume integral equation method can also be found in Section of the book Volume Integral Equation Method by Buryachenko [21].
2.1. ThreeDimensional Elastodynamic Problems Using the Parallel Volume Integral Equation Method
For threedimensional elastodynamic problems, the displacement components for SH waves in the volume integral equation (1) for fully anisotropic inclusions can be expressed in the following form:where , , and are the threedimensional displacement components, is the circular frequency of the waves, and , where denotes the density of the inclusions while represents that of the matrix material. In addition, , where denote the elastic stiffness constants of the inclusions while represent those of the matrix material. In (2), () represents the threedimensional elastodynamic Green’s function for the unbounded isotropic matrix material and is given by Pao and Varatharajulu [22] aswhere , , and , where is wave speed while is wave speed.
2.2. Parallel Computing with MPI and Optimization
The SUN computer system (Tachyon system) of KISTI (Korea Institute of Science and Technology Information) was used to conduct parallel programming and optimization. Each node has Intel Xeon X5570 Nehalem 2.93 GHz. In a representative MPI parallelization approach, global stiffness matrices were formed using MPI_ALLGATHERV (MPI_ALLGATHERV gathers data from all tasks and delivers the combined data to all tasks). Then, MPI_ALLREDUCE (MPI_ALLREDUCE combines values from all processes and distributes the result back to all processes) was applied to vectors in the conjugate gradient solver. Figure 4 shows procedures of the representative MPI parallelization approach.
3. Scattering of SH Waves by Multiple SingleLayered Orthotropic Elliptical Inclusions
3.1. Scattering of SH Waves in an Unbounded Isotropic Matrix Containing an Isotropic Inclusion
We first consider an isotropic inclusion in an unbounded matrix. The incident wave is an SH wave with particle motion along axis (Figure 5).
For SH waves, the volume integral equation (2) reduces towhere is the antiplane displacement component (i.e., the displacement component in axis), , , and is the shear modulus.
Green’s function for the SH problem is given by Beer [23] aswhere , is wavenumber, is the shear wave speed in the matrix material, and is the Hankel function of the first kind of the zeroth order.
Finite element discretization of the inclusion in (4) results in a system of linear algebraic equations for the unknown nodal displacements inside the inclusion, which can be solved for the unknowns. Once the displacement field, , within the inclusion is determined, that outside the inclusion can be obtained from (4) by evaluating the integral. The stress fields within and outside the isotropic inclusion can also be readily determined. Further details of this numerical treatment can be found in Lee and Mal [4] and in Lee et al. [5, 6].
3.2. Scattering of SH Waves in an Unbounded Isotropic Matrix Containing an Orthotropic Inclusion
We next consider an orthotropic inclusion. Let the coordinate axes , , and be taken parallel to the symmetry axes of the orthotropic material; and , denote its density and elastic constants.
For SH waves, the volume integral equation (2) reduces towhere is the antiplane displacement component (i.e., the displacement component in axis), , , and .
In (6), is Green’s function for the unbounded isotropic matrix material. Thus, the parallel volume integral equation method does not require use of Green’s function for the orthotropic material of the inclusion.
3.3. Numerical Formulation
The integrands in the volume integral equations (1), (2), (4), and (6) contain singularities of different orders due to the singular nature of Green’s function at (i.e., ), and the evaluation of the integrals requires special attention. The order of the singularity for the isotropic matrix is in Green’s function and in its derivatives. It should be noted that the volume integral equation method involves only for the isotropic matrix and its derivatives, while the boundary element method involves Green’s function for the multilayered anisotropic inclusions and its derivatives in addition to these. However, the closed form solution for the twodimensional timeharmonic elastodynamic Green’s function for the anisotropic material is not available in the literature [2]. Therefore, the numerical implementation of the boundary element method for the solving wave scattering by anisotropic inclusions is an extremely difficult if not impossible task. Furthermore, the singularities in the VIEM are weaker (integrable) than those in the BEM, where they are of the Cauchy type. We have used the direct integration scheme introduced by Cerrolaza and Alarcon [24], Li et al. [25], and Lu and Ye [26] after suitable modifications to address these singularities in the integrands. A description of evaluation of the singular integrals used in the volume integral equation is as follows [4, 27]. Triangular polar coordinates are used to reduce the order of singularity of the singular elements by one degree. This strategy results in good numerical approximations of singular integrals and converts weakly singular integrals into integrals over smooth functions. As an example, the application of the special polar coordinates to quadratic, quadrilateral, and isoparametric singular elements is explained here. First, the quadratic, quadrilateral element , representing the part of the surface of the considered body, is mapped onto a square of sidelength 2 as shown in Figure 6(a). While is embedded in a global twodimensional Cartesian space, is associated with a local twodimensional Cartesian space, characterized by the natural coordinates and . The second step of the concept consists of a subdivision of into two or three triangular subelements depending on whether the singular point, designated locally as point 1, is coinciding with a corner or a midpoint node as in Figure 6(b). The third step involves the mapping of each one of the subelements onto a unit square, employing the triangular polar coordinates. Figure 6(c) illustrates the mapping of the subelement of the lefthand part of Figure 6(b) onto and subsequently the mapping of onto . The last mapping represents a linear transformation for the purpose of carrying out numerical integrations in a standard way. Similarly, all other subdomains are to be mapped onto unit squares which in turn must be mapped onto squares .
(a)
(b)
(c)
3.4. Verification of the VIEM Solution
To the best of the authors’ knowledge, neither the analytical solution to this problem nor the closed form solution for the twodimensional timeharmonic elastodynamic Green’s function for the orthotropic material is available in the literature. Thus, the verification of the VIEM is replaced with the comparison between the VIEM solution and the analytical results for simple packing sequences (hexagonal and square) and isotropic circular inclusions in an unbounded isotropic elastic solid (see Figure 5).
The incident wave is given aswhere is wavenumber in the matrix material.
Figure 7 shows discretized models used in the volume integral equation method. The total number of standard eightnode quadrilateral and sixnode triangular elements for each inclusion used in the VIEM was 256. Standard graphite/epoxy composites with volume concentration of 0.6 and a frequency of 10 MHz are used and their material properties are given in Table 1. Table 2 shows real and imaginary parts of the calculated average strains in the central isotropic fiber by the generalized selfconsistent method (GSCM) [2] 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. It should be noted that the percentage differences in the two sets of results are less than 1% in all cases [6].


(a) 19 inclusions in a hexagonal packing sequence
(b) 25 inclusions in a square packing sequence
4. Scattering of SH Waves by Multiple TripleLayered Orthotropic Elliptical Inclusions
In this section, a detailed analysis of the displacement field at the interface between the matrix and the outermost layer and the outermost layer and the middle layer of the central inclusion is carried out for an unbounded isotropic matrix containing square and hexagonal packing arrays of threelayered orthotropic elliptical inclusions. The major axis of each elliptical inclusion is assumed to be oriented at 0° or 90° with axis and the aspect ratio is assumed to be 0.75. As a pioneer of numerical analysis, the influence of multiple threelayered orthotropic elliptical inclusions on the interfacial field in composites will be quantitatively analyzed.
4.1. Scattering of SH Waves in an Unbounded Isotropic Matrix Containing a Square Packing of 9 TripleLayered Orthotropic Elliptical (0°) Inclusions
Using the parallel volume integral equation method (PVIEM) developed in Section 2.2, we first consider a square packing of 9 triplelayered orthotropic elliptical (0°) inclusions in an unbounded isotropic matrix (see Figures 8(a) and 8(b)). The value of the aspect ratio for each elliptical inclusion was chosen to be 0.75 when the major axis coincides with axis (see Figures 8(a) and 8(b)). The incident wave is represented as an SH wave with particle motion along axis (see Figure 8). Let the coordinate axes , , and be taken parallel to the symmetry axes of the orthotropic material. Four different kinds of 9 triplelayered orthotropic elliptical inclusions made of two different orthotropic materials and one isotropic material are considered. The elastic constants for the materials of the matrix, the isotropic inclusion, the orthotropic inclusion #1, and the orthotropic inclusion #2 are listed in Table 3. For “orthotropic inclusion #1 (red color),” the density and the elastic constants and in the orthotropic layer are greater than those in the matrix. For “orthotropic inclusion #2 (blue color),” the density and and in the orthotropic layer are smaller than those in the matrix. Thus, for “triplelayered inclusion #1 (Inc. #1),” the outermost and innermost layers represent orthotropic inclusion #1, while the middle layer consists of orthotropic inclusion #2. For “triplelayered inclusion #2 (Inc. #2),” the outermost and innermost layers represent orthotropic inclusion #1 and the middle layer consists of “isotropic inclusion (yellow color).” For “triplelayered inclusion #3 (Inc. #3),” the outermost and innermost layers represent orthotropic inclusion #2 and the middle layer consists of orthotropic inclusion #1. Finally, for “triplelayered inclusion #4 (Inc. #4),” the outermost and innermost layers represent orthotropic inclusion #2 and the middle layer consists of the “isotropic inclusion.” The thickness ratio of the layers is assumed to be 3 (outermost layer), 4 (middle layer), and 3 (innermost layer).

(a) Nine triplelayered elliptical (0°) inclusions #1 (Inc. #1)
(b) Nine triplelayered elliptical (0°) inclusions #2 (Inc. #2)
(c) Seven triplelayered elliptical (90°) inclusions #3 (Inc. #3)
(d) Seven triplelayered elliptical (90°) inclusions #4 (Inc. #4)
In each triplelayered inclusion, for the orthotropic layer, , , and are used in the volume integral equation (6), whereas, for the isotropic layer, and are used in the volume integral equation (4). In (4) and (6), is Green’s function for the unbounded isotropic matrix material. Thus, the parallel volume integral equation method does not require use of Green’s function for the orthotropic material of the triplelayered orthotropic elliptical inclusions. This is in contrast to the boundary integral equation method, where the infinite medium Green’s functions for both the isotropic matrix and the orthotropic material of the triplelayered orthotropic elliptical inclusions are involved in the formulation of the equations.
Calculations are carried out for three different normalized wavenumbers (i.e., = 3.14, 5.03, and 7.85). The ratios of the corresponding wavelengths to the major radius of each elliptical inclusion are assumed to be between 0.8 and 2.0. This represents the usual range of frequencies for dynamic loading and ultrasonic testing.
The incident wave is given aswhere is wavenumber in the matrix material.
Figure 9(a) shows a typical discretized onequarter symmetry model used in the volume integral equation method [28]. Standard eightnode quadrilateral and sixnode triangular elements were used. The number of elements in each elliptical (0°) inclusion used in the VIEM was 3,600 and was determined based on a convergence test. The number of elements in the layers was 1,080 (outermost layer), 1,440 (middle layer), and 1,080 (innermost layer). Figures 9(b), 9(c), 9(d), and 9(e) show discretized VIEM models for the triplelayered orthotropic elliptical inclusion #1 (Inc. #1), inclusion #2 (Inc. #2), inclusion #3 (Inc. #3), and inclusion #4 (Inc. #4), respectively [28].
(a)
(b)
(c)
(d)
(e)
4.1.1. Displacement Field at the Interface between the Matrix and the Outermost Layer of the Central Inclusion
Figure 10 shows real and imaginary parts of the displacement component in axis, , at the interface between the isotropic matrix and the outermost layer for (a) the single triplelayered orthotropic elliptical (0°) inclusion #1 (Inc. #1), (b) the central inclusion in a square packing of 9 triplelayered orthotropic elliptical (0°) inclusions #1 (Inc. #1), where , and (c) the central inclusion in a square packing of 9 triplelayered orthotropic elliptical (0°) inclusions #1 (Inc. #1), where . Table 4 shows fiber separation distances for different fiber volume fractions. It should be noted that the fiber volume fractions of the elliptical inclusions are defined as those of the circular inclusions circumscribing the elliptical inclusions [12, 13]. The frequencies used are = 3.14, 5.03, and 7.85, where is wavenumber in the unbounded isotropic matrix and is the major radius of each elliptical inclusion.

(a)
(b)
(c)
Figures 11, 12, and 13 show the same quantities as those in Figure 10 for the triplelayered orthotropic elliptical (0°) inclusion #2 (Inc. #2), inclusion #3 (Inc. #3), and inclusion #4 (Inc. #4), respectively.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
The interaction effect of the square packing of 9 inclusions on the displacement component at the interface between the matrix and the outermost layer of the central inclusion appeared to be significant for the triplelayered orthotropic elliptical (0°) inclusions #3 and #4, where . In this case, the density value and and of orthotropic inclusion #2 in the outermost and innermost layers were assumed to be smaller than those in the matrix. In addition, the inclusions are closer to each other. However, the interaction effect of the square packing of 9 inclusions on the displacement component appeared to be minimal for the triplelayered orthotropic elliptical (0°) inclusions #1 and #2, where . In that case, the density value and and of orthotropic inclusion #1 in the outermost and innermost layers were assumed to be greater than those in the matrix. In addition, the inclusions are well separated from each other.
4.1.2. Comparison between Results for Square Packing of 9 TripleLayered and 9 SingleLayered Orthotropic Elliptical (0°) Inclusions
Figure 14 shows the real and imaginary parts of the displacement component in axis, , at the interface between the isotropic matrix and the outermost layer for (a) the single singlelayered orthotropic elliptical (0°) inclusion made of the orthotropic inclusion #1 in Table 3, (b) the central inclusion in a square packing of 9 singlelayered orthotropic elliptical (0°) inclusions made of the orthotropic inclusion #1 in Table 3, where , and (c) the central inclusion in a square packing of 9 singlelayered orthotropic elliptical (0°) inclusions made of the orthotropic inclusion #1 in Table 3, where . It should be noted that the fiber volume fractions of the elliptical inclusions are defined as those of the circular inclusions circumscribing the elliptical inclusions [12, 13]. The frequencies used are , 5.03, and 7.85, where is wavenumber in the unbounded isotropic matrix and is the major radius of each elliptical inclusion. Figure 15 shows the same quantities as those in Figure 14 for the square packing of 9 singlelayered orthotropic elliptical (0°) inclusions made of the orthotropic inclusion #2 in Table 3. It was observed that the motion at the interface between the isotropic matrix and the outermost layer for the central inclusion in a square packing of 9 triplelayered orthotropic elliptical (0°) inclusions (Inc. #1 or Inc. #2) is significantly different from the motion at the interface between the isotropic matrix and the outermost layer for the central inclusion in a square packing of 9 singlelayered orthotropic elliptical (0°) inclusions made of the orthotropic inclusion #1 in Table 3. It was also observed that the motion at the interface between the isotropic matrix and the outermost layer for the central inclusion in a square packing of 9 triplelayered orthotropic elliptical (0°) inclusions (Inc. #3 or Inc. #4) is drastically different from the motion at the interface between the isotropic matrix and the outermost layer for the central inclusion in a square packing of 9 singlelayered orthotropic elliptical (0°) inclusions made of the orthotropic inclusion #2 in Table 3.
(a)
(b)
(c)
(a)
(b)
(c)
4.1.3. Displacement Field at the Interface between the Outermost Layer and the Middle Layer of the Central Inclusion
Figure 16 shows real and imaginary parts of the displacement component in axis, , at the interface between the outermost layer and the middle layer for (a) the single triplelayered orthotropic elliptical (0°) inclusion #1 (Inc. #1), (b) the central inclusion in a square packing of 9 triplelayered orthotropic elliptical (0°) inclusions #1 (Inc. #1), where , and (c) the central inclusion in a square packing of 9 triplelayered orthotropic elliptical (0°) inclusions #1 (Inc. #1), where . The frequencies used are = 3.14, 5.03, and 7.85, where is wavenumber in the unbounded isotropic matrix and is the major radius of each elliptical inclusion.
(a)
(b)
(c)
Figures 17, 18, and 19 show the same quantities as those in Figure 16 for the triplelayered orthotropic elliptical (0°) inclusion #2 (Inc. #2), inclusion #3 (Inc. #3), and inclusion #4 (Inc. #4), respectively.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
The interaction effect of the square packing of 9 inclusions on the displacement component at the interface between the outermost layer and the middle layer of the central inclusion appeared to be significant for the triplelayered orthotropic elliptical (0°) inclusions #3 and #4, where . However, the interaction effect of the square packing of 9 inclusions on the displacement component appeared to be minimal for the triplelayered orthotropic elliptical (0°) inclusions #1 and #2, where .
It should be noted that the interaction effect on the displacement component at the interfaces of the central inclusion for the hexagonal packing of 7 triplelayered elliptical (0°) inclusions was omitted here due to lack of space.
4.2. Scattering of SH Waves in an Unbounded Isotropic Matrix Containing a Hexagonal Packing of 7 TripleLayered Orthotropic Elliptical (90°) Inclusions
We next consider a hexagonal packing of 7 triplelayered orthotropic elliptical (90°) inclusions in an unbounded isotropic matrix (see Figures 8(c) and 8(d)). The value of the aspect ratio for each elliptical (90°) inclusion was chosen to be 0.75 when the major axis coincides with axis (see Figures 8(c) and 8(d)). The incident wave is represented as an SH wave with particle motion along axis (see Figure 8). The general features of the parallel volume integral equation method for a hexagonal packing of 7 triplelayered orthotropic elliptical (90°) inclusions are similar to those for the square packing of 9 triplelayered orthotropic elliptical (0°) inclusions.
In the SH wave case, the incident wave is given aswhere is wavenumber in the matrix material. The calculations are carried out for three different normalized wavenumbers (i.e., = 3.14, 5.03, and 7.85). The ratios of the corresponding wavelengths to the major radius of each elliptical (90°) inclusion are assumed to be between 0.8 and 2.0. This represents the usual range of frequencies for dynamic loading and ultrasonic testing.
Figure 20(a) shows a typical discretized onequarter symmetry model used in the volume integral equation method [28]. Standard eightnode quadrilateral and sixnode triangular elements were used. The number of elements in each elliptical (90°) inclusion used in the VIEM was 3,600 and was determined based on a convergence test. Figures 20(b), 20(c), 20(d), and 20(e) show discretized VIEM models for the triplelayered orthotropic elliptical (90°) inclusion #1 (Inc. #1), inclusion #2 (Inc. #2), inclusion #3 (Inc. #3), and inclusion #4 (Inc. #4), respectively [28].
(a)
(b)
(c)
(d)
(e)
4.2.1. Displacement Field at the Interface between the Matrix and the Outermost Layer of the Central Inclusion
Figure 21 shows real and imaginary parts of the displacement component in axis, , at the interface between the isotropic matrix and the outermost layer for (a) the single triplelayered orthotropic elliptical (90°) inclusion #1 (Inc. #1), (b) the central inclusion in a hexagonal packing of 7 triplelayered orthotropic elliptical (90°) inclusions #1 (Inc. #1), where , and (c) the central inclusion in a hexagonal packing of 7 triplelayered orthotropic elliptical (90°) inclusions #1 (Inc. #1), where . Table 4 shows fiber separation distances for different fiber volume fractions. It should be noted that the fiber volume fractions of the elliptical (90°) inclusions are defined as those of the circular inclusions circumscribing the elliptical inclusions [12, 13]. The frequencies used are = 3.14, 5.03, and 7.85, where is wavenumber in the unbounded isotropic matrix and is the major radius of each elliptical (90°) inclusion.
(a)
(b)
(c)
Figures 22, 23, and 24 show the same quantities as those in Figure 21 for the triplelayered orthotropic elliptical (90°) inclusion #2 (Inc. #2), inclusion #3 (Inc. #3), and inclusion #4 (Inc. #4), respectively.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
The interaction effect of the hexagonal packing of 7 inclusions on the displacement component at the interface between the matrix and the outermost layer of the central inclusion appeared to be significant for the triplelayered orthotropic elliptical (90°) inclusions #3 and #4, where . In this case, the density value and and of orthotropic inclusion #2 in the outermost and innermost layers were assumed to be smaller than those in the matrix. In addition, the inclusions are closer to each other. However, the interaction effect of the hexagonal packing of 7 inclusions on the displacement component appeared to be minimal for the triplelayered orthotropic elliptical (90°) inclusions #1 and #2, where . In that case, the density value and and of orthotropic inclusion #1 in the outermost and innermost layers were assumed be greater than those in the matrix. In addition, the inclusions are well separated from each other.
4.2.2. Comparison between Results for Hexagonal Packing of 7 TripleLayered and 7 SingleLayered Orthotropic Elliptical (90°) Inclusions
Figure 25 shows the real and imaginary parts of the displacement component in axis, , at the interface between the isotropic matrix and the outermost layer for (a) the single singlelayered orthotropic elliptical (90°) inclusion made of the orthotropic inclusion #1 in Table 3, (b) the central inclusion in a hexagonal packing of 7 singlelayered orthotropic elliptical (90°) inclusions made of the orthotropic inclusion #1 in Table 3, where , and (c) the central inclusion in a hexagonal packing of 7 singlelayered orthotropic elliptical (90°) inclusions made of the orthotropic inclusion #1 in Table 3, where . It should be noted that the fiber volume fractions of the elliptical inclusions are defined as those of the circular inclusions circumscribing the elliptical inclusions [12, 13]. The frequencies used are = 3.14, 5.03, and 7.85, where is wavenumber in the unbounded isotropic matrix and is the major radius of each elliptical inclusion. Figure 26 shows the same quantities as those in Figure 25 for the hexagonal packing of 7 singlelayered orthotropic elliptical (90°) inclusions made of the orthotropic inclusion #2 in Table 3. It was observed that the motion at the interface between the isotropic matrix and the outermost layer for the central inclusion in a hexagonal packing of 7 triplelayered orthotropic elliptical (90°) inclusions (Inc. #1 or Inc. #2) is significantly different from the motion at the interface between the isotropic matrix and the outermost layer for the central inclusion in a hexagonal packing of 7 singlelayered orthotropic elliptical (90°) inclusions made of the orthotropic inclusion #1 in Table 3. It was also observed that the motion at the interface between the isotropic matrix and the outermost layer for the central inclusion in a hexagonal packing of 7 triplelayered orthotropic elliptical (90°) inclusions (Inc. #3 or Inc. #4) is drastically different from the motion at the interface between the isotropic matrix and the outermost layer for the central inclusion in a hexagonal packing of 7 singlelayered orthotropic elliptical (90°) inclusions made of the orthotropic inclusion #2 in Table 3.
(a)
(b)
(c)
(a)
(b)
(c)
4.2.3. Displacement Field at the Interface between the Outermost Layer and the Middle Layer of the Central Inclusion
Figure 27 shows real and imaginary parts of the displacement component in axis, , at the interface between the outermost layer and the middle layer for (a) the single triplelayered orthotropic elliptical (90°) inclusion #1 (Inc. #1), (b) the central inclusion in a hexagonal packing of 7 triplelayered orthotropic elliptical (90°) inclusions #1 (Inc. #1), where , and (c) the central inclusion in a hexagonal packing of 7 triplelayered orthotropic elliptical (90°) inclusions #1 (Inc. #1), where . It should be noted that the fiber volume fractions of the elliptical (90°) inclusions are defined as those of the circular inclusions circumscribing the elliptical inclusions [12, 13]. The frequencies used are = 3.14, 5.03, and 7.85, where is wavenumber in the unbounded isotropic matrix and is the major radius of each elliptical (90°) inclusion.
(a)
(b)
(c)
Figures 28, 29, and 30 show the same quantities as those in Figure 27 for the triplelayered orthotropic elliptical (90°) inclusion #2 (Inc. #2), inclusion #3 (Inc. #3), and inclusion #4 (Inc. #4), respectively.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
The interaction effect of the hexagonal packing of 7 inclusions on the displacement component at the interface between the outermost layer and the middle layer of the central inclusion appeared to be significant for the triplelayered orthotropic elliptical (90°) inclusions #3 and #4, where . However, the interaction effect of the hexagonal packing of 7 inclusions on the displacement component appeared to be minimal for the triplelayered orthotropic elliptical (90°) inclusions #1 and #2, where .
It should be noted that the interaction effect on the displacement component at the interfaces of the central inclusion for the square packing of 9 triplelayered elliptical (90°) inclusions was omitted in this section due to lack of space.
The general features of the interaction effect on the displacement component at the interfaces of the central inclusion for the square packing of 9 triplelayered elliptical (0°) inclusions are similar to those for the hexagonal packing of 7 triplelayered elliptical (90°) inclusions.
5. Concluding Remarks
The parallel volume integral equation method (PVIEM) was applied to the solution of SH wave scattering problems with square and hexagonal packing arrays of triplelayered orthotropic elliptical inclusions made of two different orthotropic materials and one isotropic material for three different frequencies in an unbounded isotropic matrix. The value of the aspect ratio for each elliptical inclusion was chosen to be 0.75 when the major axis coincides with axis or axis. It was determined that, for different packing arrays of triplelayered orthotropic elliptical inclusions, the interaction effect of the inclusions on the displacement component in axis, , at the interface between the matrix and the outermost layer or between the outermost layer and the middle layer of the central inclusion was significant for the triplelayered orthotropic elliptical inclusions #3 and #4 (basically soft inclusions), where . However, for the triplelayered orthotropic elliptical inclusions #1 and #2 (basically hard inclusions), where , the interaction effect of the inclusions on the displacement component at the interface between the matrix and the outermost layer or between the outermost layer and the middle layer of the central inclusion appeared to be minimal.
The main advantage of parallel volume integral equation method (PVIEM) over the finite element method (FEM) is the fact that it requires discretization of the inclusions only in contrast to the need to discretize the entire domain. In the presence of multiple multilayered anisotropic inclusions, the boundary integral equation method (BIEM) numerical treatment becomes cumbersome. Specifically, in elastodynamic problems involving multiple multilayered anisotropic inclusions, BIEM numerical treatment becomes extremely difficult since closed form expressions for the elastodynamic Green’s functions for anisotropic media are currently not available. However, the parallel volume integral equation method does not require use of Green’s functions for anisotropic inclusions. Furthermore, since standard finite elements are used in the PVIEM, it is easier and more convenient to handle multiple nonsmooth multilayered anisotropic inclusions.
The FORTRAN source code for the VIEM was parallelized and optimized with support from the Korea Institute of Science and Technology Information (KISTI). As a result, the program execution time has been greatly reduced. Hence, the parallel volume integral equation method (PVIEM) is now generally more efficient to apply and execute than the finite element method or boundary element method. The detailed procedures of the PVIEM were omitted in this paper to reduce space.
In this paper, the effects of single/multiple scattering, fiber packing type, fiber volume fraction, single/multiple layer(s), the multilayer’s shape and geometry, isotropy/anisotropy, and softness/hardness of the multiple multilayered anisotropic elliptical inclusions quantitatively could be investigated by calculation of displacements at interfaces of multiple multilayered anisotropic elliptical inclusions using parallel volume integral equation method as a pioneer of numerical analysis. Then, a precise knowledge of the displacements on interfaces of multiple multilayered anisotropic elliptical inclusions under remote dynamic loading can be used as one of the most important influential factors on predicting the failure and damage mechanisms in composites [29]. Therefore, the formulations developed in this paper in conjunction with other useful methods [8, 30] can be used to calculate other quantities of practical interest in realistic models of materials containing strong anisotropic and/or heterogeneous inclusions of arbitrary shapes.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
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. 2014044104) and 2014 Hongik University Research Fund. The authors would also like 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. KSC2014C1027) and express their sincere appreciation to Dr. Kyunghun Lim for parallelizing and optimizing the VIEM code and Professor Daisuke Takahashi at Center for Computational Sciences, University of Tsukuba, for his constructive suggestions.
References
 R.B. Yang and A. K. Mal, “The effective transverse moduli of a composite with degraded fibermatrix interfaces,” International Journal of Engineering Science, vol. 33, no. 11, pp. 1623–1632, 1995. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 P. A. Martin, Multiple Scattering: Interaction of TimeHarmonic Waves with N Obstacles, Cambridge University Press, London, UK, 2006.
 D. S. Brauer, C. Rüssel, S. Vogt, J. Weisser, and M. Schnabelrauch, “Degradable phosphate glass fiber reinforced polymer matrices: mechanical properties and cell response,” Journal of Materials Science: Materials in Medicine, vol. 19, no. 1, pp. 121–127, 2008. View at: Publisher Site  Google Scholar
 J. K. Lee and A. K. Mal, “A volume integral equation technique for multiple scattering problems in elastodynamics,” Applied Mathematics and Computation, vol. 67, no. 1–3, pp. 135–159, 1995. View at: Publisher Site  Google Scholar
 J. K. Lee, H. M. Lee, and A. Mal, “A mixed volume and boundary integral equation technique for elastic wave field calculations in heterogeneous materials,” Wave Motion, vol. 39, no. 1, pp. 1–19, 2004. View at: Publisher Site  Google Scholar
 J.K. Lee, Y.B. Han, and Y.J. Ahn, “SH wave scattering problems for multiple orthotropic elliptical inclusions,” Advances in Mechanical Engineering, vol. 2013, Article ID 370893, 14 pages, 2013. View at: Publisher Site  Google Scholar
 J.T. Chen, J.W. Lee, C.F. Wu, and I.L. Chen, “SHwave diffraction by a semicircular hill revisited: a nullfield boundary integral equation method using degenerate kernels,” Soil Dynamics and Earthquake Engineering, vol. 31, no. 56, pp. 729–736, 2011. View at: Publisher Site  Google Scholar
 J.T. Chen, J.W. Lee, and W.S. Shyu, “SHwave scattering by a semielliptical hill using a nullfield boundary integral equation method and a hybrid method,” Geophysical Journal International, vol. 188, no. 1, pp. 177–194, 2012. View at: Publisher Site  Google Scholar
 S. Nair, S. Stankovich, Y. Wang, and V. Raghavendran, “Multilayered fiber,” US Patent No. 7960024, 2011. View at: Google Scholar
 P. Pacheco, Parallel Programming with MPI, Morgan Kaufmann Publishers, San Francisco, Calif, USA, 1996.
 Message Passing Interface Forum. MPI: A MessagePassing Interface Standard, Version 3.0, 2012, http://www.mpiforum.org/docs/mpi3.0/mpi30report.pdf.
 J. Lee and H.R. Kim, “Volume integral equation method for multiple circular and elliptical inclusion problems in antiplane elastostatics,” Composites Part B: Engineering, vol. 43, no. 3, pp. 1224–1243, 2012. View at: Publisher Site  Google Scholar
 J. K. Lee, S. M. Oh, and A. Mal, “Calculation of interfacial stresses in composites containing elliptical inclusions of various types,” European Journal of Mechanics ASolids, vol. 44, pp. 17–40, 2014. View at: Publisher Site  Google Scholar
 J. A. DeSanto, “Theory of scattering from multilayered bodies of arbitrary shape,” Wave Motion, vol. 2, no. 1, pp. 63–73, 1980. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 S. K. Kanaun and V. M. Levin, “Effective medium method in the problem of axial elastic shear wave propagation through fiber composites,” International Journal of Solids and Structures, vol. 40, no. 18, pp. 4859–4878, 2003. View at: Publisher Site  Google Scholar
 S. Biwa, S. Yamamoto, F. Kobayashi, and N. Ohno, “Computational multiple scattering analysis for shear wave propagation in unidirectional composites,” International Journal of Solids and Structures, vol. 41, no. 2, pp. 435–457, 2004. View at: Publisher Site  Google Scholar
 M. Dravinski and M. C. Yu, “Scattering of plane harmonic SH waves by multiple inclusions,” Geophysical Journal International, vol. 186, no. 3, pp. 1331–1346, 2011. View at: Publisher Site  Google Scholar
 T. Sumiya, S. Biwa, and G. Haïat, “Computational multiple scattering analysis of elastic waves in unidirectional composites,” Wave Motion, vol. 50, no. 2, pp. 253–270, 2013. View at: Publisher Site  Google Scholar
 R. Sheikhhassani and M. Dravinski, “Scattering of a plane harmonic SH wave by multiple layered inclusions,” Wave Motion, vol. 51, no. 3, pp. 517–532, 2014. View at: Publisher Site  Google Scholar
 A. K. Mal and L. Knopoff, “Elastic wave velocities in twocomponent systems,” IMA Journal of Applied Mathematics, vol. 3, no. 4, pp. 376–387, 1967. View at: Publisher Site  Google Scholar
 V. A. Buryachenko, Micromechanics of Heterogeneous Materials, Springer, New York, NY, USA, 2007.
 Y.H. Pao and V. Varatharajulu, “Huygens’ principle, radiation conditions, and integral formulas for the scattering of elastic waves,” Journal of the Acoustical Society of America, vol. 59, no. 6, pp. 1361–1369, 1976. View at: Publisher Site  Google Scholar
 G. Beer, Programming the Boundary Element Method: An Introduction for Engineers, Wiley, Chichester, UK, 2001.
 M. Cerrolaza and E. Alarcon, “A bicubic transformation for the numerical evaluation of the Cauchy principal value integrals in boundary methods,” International Journal for Numerical Methods in Engineering, vol. 28, no. 5, pp. 987–999, 1989. View at: Publisher Site  Google Scholar
 H. B. Li, G. M. Han, and H. A. Mang, “A new method for evaluating singular integrals in stress analysis of solids by the direct boundary element method,” International Journal for Numerical Methods in Engineering, vol. 21, no. 11, pp. 2071–2098, 1985. View at: Publisher Site  Google Scholar  MathSciNet
 S. Lu and T. Q. Ye, “Direct evaluation of singular integrals in elastoplastic analysis by the boundary element method,” International Journal for Numerical Methods in Engineering, vol. 32, no. 2, pp. 295–311, 1991. View at: Publisher Site  Google Scholar
 J. Lee and A. Mal, “A volume integral equation technique for multiple inclusion and crack interaction problems,” Journal of Applied Mechanics, vol. 64, no. 1, pp. 23–30, 1997. View at: Publisher Site  Google Scholar
 PATRAN User's Manual, Version 7.0, MSC/PATRAN, 1998.
 J. K. Lee, S. J. Choi, and A. Mal, “Stress analysis of an unbounded elastic solid with orthotropic inclusions and voids using a new integral equation technique,” International Journal of Solids and Structures, vol. 38, no. 16, pp. 2789–2802, 2001. View at: Publisher Site  Google Scholar
 M. Dravinski and R. Sheikhhassani, “Scattering of a plane harmonic SH wave by a rough multilayered inclusion of arbitrary shape,” Wave Motion, vol. 50, no. 4, pp. 836–851, 2013. View at: Publisher Site  Google Scholar  MathSciNet
Copyright
Copyright © 2015 Jungki Lee et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.