Three-Dimensional Volume Integral Equation Method for Solving Isotropic/Anisotropic Inhomogeneity Problems

: In this paper, the volume integral equation method (VIEM) is introduced for the analysis of an unbounded isotropic solid composed of multiple isotropic/anisotropic inhomogeneities. A comprehensive examination of a three-dimensional elastostatic VIEM is introduced for the analysis of an unbounded isotropic solid composed of isotropic/anisotropic inhomogeneity of arbitrary shape. The authors hope that the volume integral equation method can be used to compute critical values of practical interest in realistic models of composites composed of strong anisotropic and/or heterogeneous inhomogeneities of arbitrary shapes.


Introduction
The fibers and matrix in composites are usually composed of isotropic material. However, to satisfy advanced composites, a portion of the constituents can be anisotropic. For example, in SiC/Ti (silicon carbide/titanium) metal matrix composites, the matrix is almost isotropic, whereas the SiC fiber has strong anisotropy and consists of an interphase and a core. A transverse cross section of a SiC/Ti-15-3 composite is shown in Figure 1. Various analytical methods are accessible for solving inhomogeneity problems when the geometry of the inhomogeneities is straightforward (e.g., ellipsoidal, cylindrical and spherical) or when the inhomogeneities are separated well from one another [1][2][3]. We cannot apply these methods to realistic models for a general problem when the inhomogeneities or voids are of arbitrary shape and the density of the voids or inhomogeneities is high. Therefore, stress analysis of heterogeneous solids frequently requires the utilization of numerical techniques in view of the finite element method (FEM) or the boundary integral equation method (BIEM or BEM). These techniques experience challenges in dealing with problems which include infinite media or multiple anisotropic inhomogeneities. In response to this concern, it has been shown that the volume integral equation method (VIEM) can eliminate both of these limitations in heterogeneous problems which include infinite media [4,5].
By contrast with the BIEM, the VIEM does not require utilization of Green's function for anisotropic inhomogeneities and is not sensitive to the geometry of the inhomogeneities. Moreover, as opposed to the FEM, where the whole domain needs to be discretized, the VIEM needs discretization of the multiple inhomogeneities only.
Problems associated with multiple inhomogeneities have been examined by several authors [6][7][8]. In this paper, the stress field in an unbounded isotropic elastic matrix, composed of multiple isotropic inhomogeneities and whose number and shape are arbitrary under uniform loading at infinity could be evaluated. The efficiency, accuracy and capability of the volume integral equation method will also be examined using these results.
The VIEM is applied to a class of three-dimensional elastostatic inhomogeneity problems. The details of the numerical treatment, especially the evaluation of singular integrals, for resolving three-dimensional problems in view of the VIEM are introduced in this paper. The accuracy of the VIEM is examined by comparing the results obtained from analytical solutions.
The purpose of this paper is to introduce the volume integral equation method (VIEM) as an efficient numerical method for solving multiple anisotropic inhomogeneity problems.

Volume Integral Equation Method (VIEM)
The geometry of the overall elastodynamic problem considered herein is shown in Figure 2a. It is assumed that an unbounded isotropic elastic solid, containing a number of isotropic or anisotropic inhomogeneities of arbitrary shape, is subjected to prescribed dynamic loading at infinity. In Figure 2a, V and S show the volume and surface of the inhomogeneity and the symbol n indicates the outward unit normal to S. The symbols ρ (1) and cijkl (1) denote the density and fourth-order elasticity tensors of the inhomogeneity whereas ρ (2) and cijkl (2) denote the density and the fourth-order elasticity tensors of the unbounded matrix material. The unbounded matrix material is presumed to be isotropic and homogeneous, so that cijkl (2) is a constant isotropic tensor, whereas cijkl (1) can be arbitrary. Therefore, the inhomogeneities may, in general, be anisotropic and inhomogeneous. The interfaces between the inhomogeneities as well as the unbounded matrix material are assumed to be perfectly bonded, ensuring continuity of the stress and displacement vectors.
Let um o (x, ω)e -iωt signify the mth component of the displacement vector as a result of the incident field at x in the absence of the inhomogeneities. Let um(x, ω)e -iωt signify the same within the inhomogeneities, where ω is the circular frequency of the waves. In the following example, we suppress the common time factor e -iωt and the explicit dependence of ω for all field quantities.
Mal and Knopoff [9] showed that the elastodynamic displacement um(x), within the composite, fulfills the volume integral equation: where the integral is over entire space R, δρ and δcijkl show differences in the densities and the fourth-order elasticity tensors of the inhomogeneities and the unbounded matrix material: δρ = ρ (1) − ρ (2) and δcijkl = cijkl (1) − cijkl (2) . Since cijkl = cjikl = cijlk = cklij (i, j, k, l = 1, 2, 3), cijkl can be reduced as cαβ (α, β = 1, 6). The value gi m (ξ,x) denotes the elastodynamic Green's function for the unbounded homogeneous matrix material. The value gi m (ξ,x) stands for the ith displacement component at the point ξ as a result of unit concentrated force, eme -iωt , at the point x in the mth direction. It should be noted that the symbol em represents a unit vector in the mth direction. uk,l(ξ) represents the strain field inside the inhomogeneities. A detailed expression of three-dimensional elastostatic gi,j m (ξ,x) will be in a later section.
(a) (b) The geometry of the general elastostatic problem is shown in Figure 2b. Lee and Mal [4] showed that the elastostatic displacement, um(x), within the composite, fulfills the volume integral equation: where the integral is over entire space R and δcijkl = cijkl (1) − cijkl (2) . The value gi m (ξ,x) is the elastostatic Green's function (or Kelvin's solution) for the unbounded homogeneous matrix material. The value gi m (ξ,x) stands for the ith displacement component at the point ξ as a result of unit concentrated force at the point x in the mth direction. The elastodynamic volume integral Equation (1) requires gi m (ξ,x) and gi,j m (ξ,x) while the elastostatic volume integral Equation (2) requires only gi,j m (ξ,x). This is in contrast to the BEM. In Equations (1) and (2), the summation convention and comma notation have been utilized and we differentiate them with respect to the integration variable ξi. It is evident that the integrand is non-zero within the inhomogeneities only, since δcijkl = 0, outside the inhomogeneities.
In Equations (1) and (2), x represents a field point while ξ represents a point inside and on the boundary of the inhomogeneities (see Figure 3). Here, the field point is an exterior point or a point inside and on the boundary of the inhomogeneities. The field point (x) and the point inside and on the boundary of the inhomogeneities (ξ) are independent of each other. If x lies in the domain occupied by the inhomogeneities, then Equations (1) and (2) are integro-differential equations for the unknown displacement vector u(x) within the inhomogeneities, which can, in principle, be resolved through the solution of Equations (1) and (2). It is extremely challenging and sometimes impossible to solve systems of linear equations analytically even within the presence of a single inhomogeneity of arbitrary shape. An algorithm based on the discretization of Equations (1) and (2) was developed by Lee and Mal [4,5] to compute numerically the unknown displacement vector u(x) by discretizing the inhomogeneities utilizing standard finite elements. Once u(x) inside the inhomogeneities is determined, the displacement field outside the inhomogeneities can be acquired from Equations (1) and (2) by assessing the integral. The stress field inside and outside the inhomogeneities can also be resolved without difficulty. Details of the numerical treatment of Equations (1) and (2) can be seen in references [5,10,11] for plane elastodynamic problems and in Lee and Mal [4] for plane elastostatic problems. Further mathematical detailing of the elastostatic VIEM can also be seen in Section 4.3 from the book "Volume Integral Equation Method" by Buryachenko [12]. A general description of the volume integral equation method can be viewed in Chapter 4 entitled "Volume Integral Equation Method (VIEM)" written by the first author of this paper, contained in the book "Advances in Computers and Information in Engineering Research, Vol. 2", edited by Michopoulos et al. (eds.) [13].
In Section 3.1 of the reference [16], Buryachenko pointed out that "the VIEM was quite time-consuming and no optimized commercial software existed for its application." Standard parallel programming, such as MPI (message passing interface), has been utilized to speed up computation in the VIEM. The parallel volume integral equation method (PVIEM) enabled us to investigate more complicated multiple anisotropic inhomogeneity problems elastostatically or elastodynamically. Since there is no commercial software for the VIEM, we are developing a VIEMAP (volume integral equation method application program). This VIEM modeling software includes a pre-processor, a solver and a post-processor, which are adapted to solve multiple isotropic/anisotropic inhomogeneity problems more easily and efficiently.
The authors intend to help university researchers and college students in undergraduate degree programs make VIEM models using the VIEMAP more conveniently than the finite element method and solve multiple isotropic/anisotropic inhomogeneity problems in an unbounded isotropic media more accurately and easily than the boundary element method.

A Single Cubic Inhomogeneity Problem
In this section, in order to introduce the volume integral equation method for three-dimensional elastostatic problems, we consider a single isotropic or orthotropic cubic inhomogeneity in an unbounded isotropic matrix subject to uniform remote tensile loading, σ o xx, as shown in Figure 4 (see also Figure 2b).
Next we consider a single orthotropic cubic inhomogeneity (−4 mm ≤ x, y, z ≤ 4 mm) in the unbounded isotropic matrix. The same four models for the isotropic inhomogeneity are used for the orthotropic inhomogeneity.  The elastic constants for the isotropic matrix, the isotropic inhomogeneity and three different kinds of orthotropic inhomogeneity are listed in Tables 1 and 2. In order to distinguish different material properties easily, we assign a distinct material name (mat_01, mat_02, mat_03, ---) to each material. The coordinate axes are x1, x2, x3. The stress-strain relationships (σα = cαβεβ (α, β = 1, 6)) for the isotropic inhomogeneities can be expressed in the form: It should be noted that the stress-strain relationships are also referred to as constitutive equations.
In Equations (4)- (6), gi m (ξ,x) shows the Green's function for the unbounded isotropic matrix material and is presented by Banerjee [18] and Pao and Varatharajulu [19] as: , ν is Poisson's ratio, and μ is the shear modulus. Also, gi,j m (ξ,x) is presented as: It should be noted that for the applied uniformly remote stress, the displacement vector u° is of the form: where the constants C1-C3 are related to the tensile and shear components of the applied uniformly remote stress.

Volume Integral Equation for an Unbounded Isotropic Matrix Containing a Single Orthotropic Cubic Inhomogeneity
Let the coordinate axes x1, x2, x3 be taken parallel to the symmetry axes of the orthotropic material. The constitutive equations for the orthotropic inhomogeneities can be expressed in the form: .

Numerical Formulation
The integrands in the volume integral Equations (4)-(6) and Equations (13)-(15) contain singularities with different orders due to the singular characteristics of Green's function at x = ξ (i.e., r = 0), and the evaluation of the singular integrals requires special attention. In general, gi m (ξ,x) behaves as 1/r while its derivatives behave as 1/r 2 as r → 0. It ought to be noted that the VIEM involves only gi m (ξ,x) for the isotropic matrix and its derivatives. By contrast, the BEM involves Green's function for these and for the anisotropic inhomogeneities and their derivatives. The singularities in the VIEM are weaker (integrable) when compared with those in the BEM. Therefore, we have utilized the direct integration scheme introduced by Li et al. [20] after reasonable adjustments to address these singularities in the integrands.
The order of singularity of the singular elements decreases by one degree with tetrahedron polar co-ordinates [20] in the VIEM. This strategy also converts weakly singular integrals into integrals over smooth functions. The tetrahedron polar co-ordinates are applied to quadratic, hexahedral, isoparametric singular elements as follows: (1) In Figure 5a, Ω f stands for the quadratic hexahedral element indicating the interior of the considered body. Ω f is associated with a global three-dimensional Cartesian coordinate system. Ω f is mapped onto a cube Ω f′ of side-length 2, characterized by the local three-dimensional Cartesian coordinates, ξ1, ξ2 and ξ3. (2) If the singular point is located at a corner node, Ω f′ is divided into two triangular prisms, Ω1 f′ and Ω2 f′ , as shown in  . If the singular point is located at a midside node, then Ω f′ is divided into three triangular prisms, Ω1 f′ , Ω2 f′ and Ω3 f′ , as shown in Figure 5b. (3) A subdivision of each triangular prism is divided into three tetrahedral subcells as shown in Figure 5c. (4) Each subcell maps onto a unit cube, applying the tetrahedron polar coordinates [20]; the subcell Ω11 f ′ maps onto Ω11 f ″ .

Numerical Results
In the VIEM, the displacements and stresses inside the cubic inhomogeneity are first calculated using the models in Figure 4. We calculate the displacements and stresses outside the cubic inhomogeneity using Equation (2), by adding standard finite elements in appropriate locations. Figure 8 shows typical VIEM models for the single cubic inhomogeneity (−4 mm ≤ x, XX ≤ 4 mm) (between two green lines) including the outside of the cubic inhomogeneity (−8 mm ≤ x ≤ −4 mm and 4 mm ≤ x ≤ 8 mm). Since x ≠ ξ outside the inhomogeneity, singularities do not exist while calculating the displacements and stresses outside the inhomogeneity in the VIEM. Therefore, the displacement and stress fields inside and outside the inhomogeneities can be solved without difficulty. Figure 9 shows normalized tensile stress component (σxx/σ o xx) along the x-axis and the XX-axis (−8 mm ≤ x, XX ≤ 8 mm) for Model_8 × 8 × 8 (512 elements), Model_10 × 10 × 10 (1000 elements), Model_14 × 14 × 14 (2744 elements) and Model_16 × 16 × 16 (4096 elements) under uniform remote tensile loading for the isotropic cubic inhomogeneity. The remote applied load was assumed to be σ o xx = 143.1 GPa.
Standard 20-node quadratic hexahedral elements in Figure 6 were used in the VIEM models. In the Model_nxnxn, the symbol n represents the number of line segments split in each side of a cube. A cubic inhomogeneity (−4 mm ≤ x, y, z ≤ 4 mm) is divided into n 3 hexahedral elements in the Model_nxnxn. For example, in the Model_14 × 14 × 14, a cubic inhomogeneity (−4 mm ≤ x, y, z ≤ 4 mm) is divided into 2744 (=14 3 ) elements. Figure 9a shows the normalized tensile stress component (σxx/σ o xx) along the x-axis inside the isotropic cubic inhomogeneity (−4 mm ≤ x ≤ 4 mm) and outside the isotropic cubic inhomogeneity (−8 mm ≤ x ≤ −4 mm and 4 mm ≤ x ≤ 8 mm), respectively. Figure 9b shows the normalized tensile stress component (σxx/σ o xx) along the XX-axis inside the isotropic cubic inhomogeneity (−4 mm ≤ XX ≤ 4 mm) and outside the isotropic cubic inhomogeneity (−8 mm ≤ XX ≤ −4 mm and 4 mm ≤ XX ≤ 8 mm), respectively. The results using the VIEM for the different models converged very well within this range of hexahedral elements.  Table 2) (−4 mm ≤ x ≤ 4 mm) and outside the cubic inhomogeneity (−8 mm ≤ x ≤ −4 mm and 4 mm ≤ x ≤ 8 mm), respectively. Figure 10b shows the normalized tensile stress component (σxx/σ o xx) along the XX-axis inside the orthotropic cubic inhomogeneity #1 (−4 mm ≤ XX ≤ 4 mm) and outside the cubic inhomogeneity (−8 mm ≤ XX ≤ −4 mm and 4 mm ≤ x ≤ 8 mm), respectively. The VIEM solutions using the different models converged very well within this range of hexahedral elements. Figures 11a and 12a show the same as Figure 10a for the orthotropic cubic inhomogeneity #2 and #3 (mat_03 and mat_04 in Table 2) (−8 mm ≤ x ≤ 8 mm), respectively. Figures 11b and 12b show the same as Figure 10b for the orthotropic cubic inhomogeneity #2 and #3 (−8 mm ≤ XX ≤ 8 mm), respectively. The results using the VIEM for the different models converged very well within this range of hexahedral elements.
The normalized tensile stress component (σxx/σ o xx) for the isotropic cubic inhomogeneity in Figure 9 appears to be considerably different from those of the orthotropic cubic inhomogeneities in

A Single Isotropic Spherical Inhomogeneity
In order to examine the accuracy of the numerical results using the VIEM, the numerical results using the VIEM for a single isotropic spherical inhomogeneity were first compared to the analytical solutions [21]. We considered a single isotropic spherical inhomogeneity with a radius of 7 mm in an unbounded isotropic matrix subject to uniform remote tensile loading, σxx o , as shown in Figure 13 (also see Figure 2b). In Figure 13, standard 20-node quadratic hexahedral elements in Figure 6 were used in the discretization [17]. The number of hexahedral elements was 4320. In order to investigate the accuracy of the hexahedral elements, only hexahedral elements without tetrahedral elements were used in Figure 13. The elastic properties of the isotropic matrix and the isotropic inhomogeneity for the calculations are given in Table 5.   Figure 14 shows a comparison between the numerical results using the volume integral equation method (VIEM) and the analytical solutions [21] for the normalized tensile stress component (σxx/σ o xx) along the x-axis inside the isotropic spherical inhomogeneities with a radius of 7 mm under uniform remote tensile loading. The remote applied load was assumed to be σ o xx = 150.0 GPa. It should be noted that the normalized tensile stress component (σxx/σ o xx) inside the isotropic spherical inhomogeneities was found to be constant [21,22].
Eshelby [22] solved the inhomogeneity problem of an ellipsoidal inclusion using the equivalent inclusion method. The stresses inside the inclusion are determined by: ijkl kl kl ijkl klpq pq kl (16) where cijkl (2) represents the fourth-order elasticity tensors of the matrix and εkl (c) stands for the constraint strain (εkl (c) = Sklpqεpq*). εkl* represents the eigenstrain of the equivalent inclusion and Sklpq represents Eshelby's fourth-order tensor. It should be noted that the Eshelby's tensor satisfies minor symmetries, Sijkl = Sjikl = Sijlk. In order to obtain the analytical solution, we used Equations (5) and (6) in Shibata and Ono [21] and the expressions in Equation (17) In Equation (17), ν is Poisson's ratio for the matrix material. For the isotropic spherical inhomogeneity (mat_05), the analytical solution was 1.309091. For the isotropic spherical inhomogeneity (mat_06), the analytical solution was 1.617336. An excellent agreement was found between the analytical and numerical solutions using the VIEM for both cases considered. Figure 14 shows that the percentage differences for the two sets of results are less than 1% in both cases.

A Single Orthotropic Spherical Inhomogeneity
We consider a single orthotropic spherical inhomogeneity in an unbounded isotropic matrix under uniform remote tensile loading, σxx o , as shown in Figure 13 (also see Figure 2b). The elastic properties of the isotropic matrix and orthotropic inhomogeneities used in these calculations are given in Table 2. Figure 15 shows the numerical results using the volume integral equation method (VIEM) for the normalized tensile stress component (σxx/σ o xx) along the x-axis inside the orthotropic spherical inhomogeneities (the orthotropic inhomogeneity #1 and #2; mat_02 and mat_03 in Table 2) with a radius of 7 mm under uniform remote tensile loading. The remote applied load was assumed to be σ o xx = 143.1 GPa. It should be noted that the normalized tensile stress component (σxx/σ o xx) inside the orthotropic spherical inhomogeneities was found to be constant [23].

Conclusions
The volume integral equation method was applied to a class of three-dimensional elastostatic inhomogeneity problems. First, we introduced the details of the numerical treatment, especially the evaluation of singular integrals, for solving three-dimensional problems using the VIEM. We considered a single isotropic or orthotropic cubic inhomogeneity (−4 mm ≤ x, y, z ≤ 4 mm) in the unbounded isotropic matrix under uniform remote tensile loading.
It was determined that the VIEM solutions, using the different models, converged very well for all cases. Next, we considered a single isotropic/orthotropic spherical inhomogeneity in an unbounded isotropic matrix under uniform remote tensile loading. An excellent agreement was found between the analytical and numerical solutions using the VIEM for single isotropic spherical inhomogeneity problems. The normalized tensile stress component (σxx/σ o xx) inside the spherical inhomogeneities was also found to be constant for the single orthotropic spherical inhomogeneity problems. Figure 15. Numerical results using the volume integral equation method (VIEM) for the normalized tensile stress component (σxx/σoxx) along the x-axis inside the orthotropic spherical inclusions with a radius of 7 mm under uniform remote tensile loading.
The major merit of the volume integral equation method (VIEM), compared to the finite element method (FEM), is that it needs discretization of the inhomogeneities only as opposed to discretization of the whole domain. In elastodynamic problems, involving multiple anisotropic inhomogeneities, numerical treatment of the BIEM is very difficult to achieve, since closed-form analytical expressions of the elastodynamic Green's functions in anisotropic media are currently not available. By contrast, the VIEM does not require the use of Green's functions for anisotropic inhomogeneities. Since the VIEM is designed to use standard finite elements, it is simpler and more advantageous for dealing with multiple non-smooth anisotropic inhomogeneities. Thus, the VIEM is now generally more applicable and executable than the boundary element or finite element methods. Moreover, the volume integral equation method (VIEM) can be used to compute critical values of practical interest in realistic models of composites with strong anisotropic inhomogeneities of arbitrary shapes.