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 three-dimensional 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 triple-layered orthotropic elliptical inclusions. Numerical results are presented for the displacement fields at the interfaces for square and hexagonal packing arrays of triple-layered 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 two-dimensional 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 three-layered 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 three-layered 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 three-layered 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 three-layered orthotropic elliptical inclusions are likely to be the highest there. The incident waves are assumed to be time-harmonic 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. Three-Dimensional Elastodynamic Problems Using the Parallel Volume Integral Equation Method
For three-dimensional 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 three-dimensional 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 three-dimensional 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 Single-Layered 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 two-dimensional time-harmonic 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 side-length 2 as shown in Figure 6(a). While is embedded in a global two-dimensional Cartesian space, is associated with a local two-dimensional 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 left-hand 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 two-dimensional time-harmonic 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 eight-node quadrilateral and six-node 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 self-consistent 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 Triple-Layered 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 three-layered 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 three-layered 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 Triple-Layered 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 triple-layered 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 triple-layered 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 “triple-layered inclusion #1 (Inc. #1),” the outermost and innermost layers represent orthotropic inclusion #1, while the middle layer consists of orthotropic inclusion #2. For “triple-layered inclusion #2 (Inc. #2),” the outermost and innermost layers represent orthotropic inclusion #1 and the middle layer consists of “isotropic inclusion (yellow color).” For “triple-layered inclusion #3 (Inc. #3),” the outermost and innermost layers represent orthotropic inclusion #2 and the middle layer consists of orthotropic inclusion #1. Finally, for “triple-layered 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 triple-layered elliptical (0°) inclusions #1 (Inc. #1)
(b) Nine triple-layered elliptical (0°) inclusions #2 (Inc. #2)
(c) Seven triple-layered elliptical (90°) inclusions #3 (Inc. #3)
(d) Seven triple-layered elliptical (90°) inclusions #4 (Inc. #4)
In each triple-layered 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 triple-layered 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 triple-layered 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 one-quarter symmetry model used in the volume integral equation method [28]. Standard eight-node quadrilateral and six-node 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 triple-layered 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 triple-layered orthotropic elliptical (0°) inclusion #1 (Inc. #1), (b) the central inclusion in a square packing of 9 triple-layered orthotropic elliptical (0°) inclusions #1 (Inc. #1), where , and (c) the central inclusion in a square packing of 9 triple-layered 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 triple-layered 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 triple-layered 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 triple-layered 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 Triple-Layered and 9 Single-Layered 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 single-layered orthotropic elliptical (0°) inclusion made of the orthotropic inclusion #1 in Table 3, (b) the central inclusion in a square packing of 9 single-layered 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 single-layered 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 single-layered 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 triple-layered 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 single-layered 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 triple-layered 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 single-layered 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 triple-layered orthotropic elliptical (0°) inclusion #1 (Inc. #1), (b) the central inclusion in a square packing of 9 triple-layered orthotropic elliptical (0°) inclusions #1 (Inc. #1), where , and (c) the central inclusion in a square packing of 9 triple-layered 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 triple-layered 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 triple-layered 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 triple-layered 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 triple-layered 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 Triple-Layered Orthotropic Elliptical (90°) Inclusions
We next consider a hexagonal packing of 7 triple-layered 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 triple-layered orthotropic elliptical (90°) inclusions are similar to those for the square packing of 9 triple-layered 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 one-quarter symmetry model used in the volume integral equation method [28]. Standard eight-node quadrilateral and six-node 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 triple-layered 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 triple-layered orthotropic elliptical (90°) inclusion #1 (Inc. #1), (b) the central inclusion in a hexagonal packing of 7 triple-layered orthotropic elliptical (90°) inclusions #1 (Inc. #1), where , and (c) the central inclusion in a hexagonal packing of 7 triple-layered 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 triple-layered 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 triple-layered 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 triple-layered 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 Triple-Layered and 7 Single-Layered 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 single-layered orthotropic elliptical (90°) inclusion made of the orthotropic inclusion #1 in Table 3, (b) the central inclusion in a hexagonal packing of 7 single-layered 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 single-layered 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 single-layered 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 triple-layered 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 single-layered 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 triple-layered 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 single-layered 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 triple-layered orthotropic elliptical (90°) inclusion #1 (Inc. #1), (b) the central inclusion in a hexagonal packing of 7 triple-layered orthotropic elliptical (90°) inclusions #1 (Inc. #1), where , and (c) the central inclusion in a hexagonal packing of 7 triple-layered 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 triple-layered 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 triple-layered 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 triple-layered 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 triple-layered 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 triple-layered elliptical (0°) inclusions are similar to those for the hexagonal packing of 7 triple-layered 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 triple-layered 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 triple-layered 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 triple-layered orthotropic elliptical inclusions #3 and #4 (basically soft inclusions), where . However, for the triple-layered 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. 2014-044104) 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. KSC-2014-C1-027) 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.