Volume Integral Equation Method Solution for Spheroidal Inclusion Problem

In this paper, the volume integral equation method (VIEM) is introduced for the numerical analysis of an infinite isotropic solid containing a variety of single isotropic/anisotropic spheroidal inclusions. In order to introduce the VIEM as a versatile numerical method for the three-dimensional elastostatic inclusion problem, VIEM results are first presented for a range of single isotropic/orthotropic spherical, prolate and oblate spheroidal inclusions in an infinite isotropic matrix under uniform remote tensile loading. We next considered single isotropic/orthotropic spherical, prolate and oblate spheroidal inclusions in an infinite isotropic matrix under remote shear loading. The authors hope that the results using the VIEM cited in this paper will be established as reference values for verifying the results of similar research using other analytical and numerical methods.


Introduction
The matrix and fibers in composites are usually made of isotropic material. However, in order to have higher strength and stiffness for commercial use, especially in the aerospace and automobile sectors, some constituents of metal matrix composites can be anisotropic. Since anisotropic materials are able to enhance mechanical properties toward orientation, certain mechanical properties (e.g., tensile strength) of anisotropic materials thus depend on orientation. As an example, in titanium-silicon carbide (Ti-SiC) composites, the matrix is nearly isotropic, but the SiC fiber has strong anisotropy and a multilayered structure including an interphase and a core.
A number of analytical techniques for solving inclusion problems are available when the inclusions are simple two-dimensional shapes (cylindrical and elliptical) or simple threedimensional shapes (spherical and ellipsoidal) and when they are well-separated [1][2][3][4][5]. In particular, Eshelby developed a simple and elegant method for solving the inclusion problem in isotropic solids in 1957 [1]. Eshelby first pointed out that the resulting elastic field can be found with the help of a sequence of imaginary cutting, straining and welding operations [1]. Eshelby also found that the strain and stress field inside the ellipsoidal inclusion is uniform and has a closed-form solution, regardless of the material properties and initial eigenstrain [1]. Eshelby's findings significantly influenced the mechanics of composites.
In the micromechanical analysis of composite materials, it is often assumed that the inclusions are periodically distributed in the matrix. Then, the unit-cell model with periodic boundary conditions is used to evaluate the overall, microstructure-insensitive, material properties of the composite. However, in real composites, the distribution of the inclusions is not periodic. Thus, the unit-cell model may not provide accurate estimates of the failure and damage mechanisms in composites [6][7][8].
Therefore, stress analysis of heterogeneous solids often requires the use of numerical approaches based on the standard finite element or boundary element formulations. How-ever, both methods present difficulties in dealing with problems involving infinite media or multiple anisotropic inclusions. In response to this concern, it has been demonstrated that the volume integral formulation can overcome both of these limitations in heterogeneous problems involving infinite media [9][10][11].
In comparison to the boundary element method (BEM), the volume integral equation method (VIEM) does not require the use of the Green's function for anisotropic inclusions and is not sensitive to the geometry of the inclusions. Moreover, as opposed to the standard finite element method (FEM), where it is necessary to discretize the full domain, the multiple inclusions only need to be discretized in the VIEM.
In this paper, three-dimensional elastostatic inclusion problems using the volume integral equation method (VIEM) will be investigated.
In order to introduce the VIEM as a versatile numerical method for the three-dimensional elastostatic inclusion problem, we first examine single isotropic/orthotropic spherical, prolate and oblate spheroidal inclusions in an infinite isotropic matrix subject to uniform remote tensile loading. Two different prolate and oblate spheroidal inclusions with an aspect ratio of 0.5 and 0.75 are considered, respectively. The matrix is assumed to be isotropic. Eight isotropic and five orthotropic inclusions with different characteristics are considered in the numerical calculation. The normalized tensile stress inside the inclusions is investigated in two different directions. Next, we examine single isotropic/orthotropic spherical, prolate and oblate spheroidal inclusions in an infinite isotropic matrix subject to remote shear loading. Two different prolate and oblate spheroidal inclusions with an aspect ratio of 0.5 and 0.75 are considered, respectively. The matrix is assumed to be isotropic. Three isotropic and two orthotropic inclusions with different characteristics are considered in the numerical calculation. The normalized shear stress inside the inclusions is investigated in two different directions.
The authors hope that the present solutions using the parallel volume integral equation method for the single isotropic/orthotropic spherical, prolate and oblate spheroidal inclusions with different material properties under uniform remote tensile loading or remote shear loading will be established as reference values for verifying the results of other analytical and numerical methods.
Since the VIEM is a combination of two powerful general-purpose numerical methods, the standard finite element method (FEM) and the boundary element method (BEM), it is also a highly beneficial tool in the field of numerical analysis and can play a very important role in solving inclusion problems. Subsequently, the purpose of this paper is to introduce the parallel volume integral equation method (PVIEM) as an accessible, versatile and powerful numerical method for solving inclusion problems in the areas of computational mechanics and mechanics of composite materials.

Governing Equations of Volume Integral Equation Formulation
The geometry of the general elastodynamic problem is shown in Figure 1a, where an infinite homogeneous, isotropic and linearly elastic solid containing a number of isotropic or anisotropic inclusions of arbitrary number and shape are subjected to prescribed dynamic loading at infinity.
In Figure 1a, V and S represent the volume and surface of the inclusion respectively, and n is the outward unit normal to S while V o and S o represent the infinite volume and surface, respectively.
The symbols ρ (1) and c ijkl (1) denote the density and the elastic stiffness tensor of the inclusion, while ρ (2) and c ijkl (2) denote the density and the elastic stiffness tensor of the infinite homogeneous, isotropic and linearly elastic matrix material, respectively. Therefore, c ijkl (2) is a constant isotropic tensor, while c ijkl (1) can be arbitrary, i.e., the inclusions may, in general, be inhomogeneous and anisotropic. The isotropic or anisotropic inclusions are assumed to be perfectly bonded to the matrix. In Figure 1a, V and S represent the volume and surface of the inclusion respectively, and n is the outward unit normal to S while Vo and So represent the infinite volume and surface, respectively.
The symbols ρ (1) and cijkl (1) denote the density and the elastic stiffness tensor of the inclusion, while ρ (2) and cijkl (2) denote the density and the elastic stiffness tensor of the infinite homogeneous, isotropic and linearly elastic matrix material, respectively. Therefore, cijkl (2) is a constant isotropic tensor, while cijkl (1) can be arbitrary, i.e., the inclusions may, in general, be inhomogeneous and anisotropic. The isotropic or anisotropic inclusions are assumed to be perfectly bonded to the matrix.
In Equation (1), u m o (x, ω)e −iωt represents the mth component of the displacement vector due to the incident field at x in the absence of the inclusions, while u m (x, ω)e −iωt denotes the same quantity in the presence of the isotropic or anisotropic inclusions, where ω is the circular frequency of the waves. In what follows, the explicit dependence on the circular frequency, and the common time factor, e −iωt , for all field quantities will be suppressed. The geometry of the general elastostatic problem is shown in Figure 1b-e. It has been shown by Lee and Mal [9] that the corresponding elastostatic displacement, u m (x), within the composite, fulfills the volume integral equation as: δc ijkl g m i,j (ξ, x)u k,l (ξ)dξ (2) where the integral is over the space V occupied by the isotropic or anisotropic inclusions and δc ijkl = c ijkl (1) − c ijkl (2) . The value g i m (ξ,x) represents the elastostatic Kelvin's solution (or Green's function) for the infinite homogeneous, isotropic and linearly elastic matrix material.
In Equations (1) and (2), the differentiations are with respect to the integration variable, ξ i , and the summation convention and comma notation have been utilized. The integrand is non-zero within the isotropic or anisotropic inclusions only, since δc ijkl = 0 outside the inclusions.
If x lies inside the inclusions, then Equations (1) and (2) are integro-differential equations for the unknown displacement vector u(x) within the inclusions. It should be noted that an algorithm was developed by Lee and Mal [9,10] to numerically calculate the unknown displacement vector u(x) by discretizing the inclusions only using standard finite elements. Once u(x) within the inclusions is determined, the displacement field outside the inclusions can be obtained from Equations (1) and (2) by evaluating the corresponding integrals respectively, and the stress field within and outside the inclusions can also be readily determined.
Furthermore, Section 4.3 entitled 'Volume Integral Equation Method' of the book "Micromechanics of Heterogeneous Materials" by Buryachenko [18] also explains further mathematical formulation of the elastostatic volume integral equation method. In particular, a general description of the volume integral equation method is presented in Chapter 4 entitled 'Volume Integral Equation Method (VIEM)' of the book "Advances in Computers and Information in Engineering Research, Vol. 2" by Michopoulos et al. (eds.) [22]. In addition, complete descriptions of the fundamental numerical technique of Equation (2) can be found in [17] for three-dimensional elastostatic problems.
Although each numerical method has certain advantages, specific disadvantages have led to further discussion and research. For example, in Section 3.1 of Reference [20], Buryachenko points out that the VIEM is quite time-consuming. Moreover, no optimized commercial software exists for its application.
Firstly, in order to resolve this 'time-consuming' problem, we propose the parallel volume integral equation method and implement MPI-based code. Such method allows us not only to solve the large domain but also to speed up computation in the volume integral equation method. The FORTRAN 90 (Version 1.1, IBM, Armonk, NY, USA) source code containing about 9000 lines for the three-dimensional VIEM of the previous paper [17] was parallelized and optimized for this paper, with support from the Korea Institute of Science and Technology Information (KISTI, Daejeon, Korea). Figure 2 shows the procedures of a representative MPI parallelization approach ("pvi3ds01_sm7560xx.f90") for the sequential three-dimensional VIEM code ("svi3ds01_sm4320xx.f"). As a result, the program execution time has been greatly reduced. Furthermore, we could use more finite elements (31,857 nodes and 7560 elements) in the VIEM model of this paper than those (18,109 nodes and 4320 elements) in the VIEM model of the previous paper [17]. The parallel FORTRAN source code for the three-dimensional VIEM is presently being processed in the KISTI-5. It is referred to as "Nurion", which is a system consisting of compute nodes, CPU-only nodes, Omni-Path interconnect networks, Burst Buffer high-speed storage, a Luster-based parallel file system and a water-cooling device based on a Rear Door Heat Exchanger (RDHx). The CPU-only nodes consist of 132 Intel Xeon 6148 2.4 GHz processors (named "Skylake"). The total theoretical performance is 25.7 petaflops, which ranked 11th in the world in June 2018 (http://www.top500.org, accessed on 3 May 2021). It should be noted that, in order to investigate three-dimensional stress problems with multiple inclusions, in addition to parallelization and optimization of the sequential three-dimensional VIEM code, a domain decomposition method (DDM) was applied to the parallel three-dimensional VIEM code, with support from the Korea Institute of Science and Technology Information (KISTI). The domain decomposition method allows decomposition of large-sized problem solutions to solutions of several smaller-sized problems [23]. Therefore, the parallel volume integral equation method (PVIEM) using the domain decomposition method enables us to investigate more complicated multiple inclusion problems elastostatically or elastodynamically. tial three-dimensional VIEM code ("svi3ds01_sm4320xx.f"). As a result, the program execution time has been greatly reduced. Furthermore, we could use more finite elements (31,857 nodes and 7560 elements) in the VIEM model of this paper than those (18,109 nodes and 4320 elements) in the VIEM model of the previous paper [17]. The parallel FORTRAN source code for the three-dimensional VIEM is presently being processed in the KISTI-5. It is referred to as "Nurion", which is a system consisting of compute nodes, CPU-only nodes, Omni-Path interconnect networks, Burst Buffer high-speed storage, a Luster-based parallel file system and a water-cooling device based on a Rear Door Heat Exchanger (RDHx). The CPU-only nodes consist of 132 Intel Xeon 6148 2.4 GHz processors (named "Skylake"). The total theoretical performance is 25.7 petaflops, which ranked 11th in the world in June 2018 (http://www.top500.org, accessed on 3 May, 2021). It should be noted that, in order to investigate three-dimensional stress problems with multiple inclusions, in addition to parallelization and optimization of the sequential three-dimensional VIEM code, a domain decomposition method (DDM) was applied to the parallel threedimensional VIEM code, with support from the Korea Institute of Science and Technology Information (KISTI). The domain decomposition method allows decomposition of largesized problem solutions to solutions of several smaller-sized problems [23]. Therefore, the parallel volume integral equation method (PVIEM) using the domain decomposition method enables us to investigate more complicated multiple inclusion problems elastostatically or elastodynamically.  Secondly, in order to resolve the 'no optimized commercial software' problem, we plan to develop a semi-commercial VIEM software called the "Volume Integral Equation Method Application Program" (VIEMAP). Table 1 shows the analysis capabilities of VIEMAP including a pre-processor (ViemMesh), a solver (VIEM) and a post-processor (ViemPlot) adapted to solve multiple isotropic/anisotropic inclusion problems in a computationally tractable manner. Figure 3 shows the registered trademark for the VIEMAP. The authors aim to help both university students and researchers create VIEM models using the VIEMAP more easily than using the standard finite element method (FEM), as well as solve multiple isotropic/anisotropic inclusion problems in an unbounded isotropic medium more accurately and conveniently than the boundary element method (BEM). MAP including a pre-processor (ViemMesh), a solver (VIEM) and a post-processor (ViemPlot) adapted to solve multiple isotropic/anisotropic inclusion problems in a computationally tractable manner. Figure 3 shows the registered trademark for the VIEMAP. The authors aim to help both university students and researchers create VIEM models using the VIEMAP more easily than using the standard finite element method (FEM), as well as solve multiple isotropic/anisotropic inclusion problems in an unbounded isotropic medium more accurately and conveniently than the boundary element method (BEM).

Three-Dimensional Elastostatic Problems Using the VIEM
In this section, we first examine single isotropic/orthotropic spherical, prolate and oblate spheroidal inclusions in an infinite isotropic matrix subject to uniform remote tensile loading, σ o xx, as shown in Figure 4 (also see Figures 1b and 5). The remote applied load can be arbitrarily chosen and was assumed to be σ o xx = 143.10 GPa for convenience purposes only. Two different prolate spheroidal inclusions were considered: (a) a/b = c/b = 0.5 and (b) a/b = c/b = 0.75 (see Figure 5). Additionally, two different oblate spheroidal inclusions were considered: (a) b/a = c/a = 0.5 and (b) b/a = c/a = 0.75 (see Figure 5).
(a) A spherical inhomogeneity ( b) A prolate spheroidal inhomogeneity

Three-Dimensional Elastostatic Problems Using the VIEM
In this section, we first examine single isotropic/orthotropic spherical, prolate and oblate spheroidal inclusions in an infinite isotropic matrix subject to uniform remote tensile loading, σ o xx , as shown in Figure 4 (also see Figures 1b and 5). The remote applied load can be arbitrarily chosen and was assumed to be σ o xx = 143.10 GPa for convenience purposes only. Two different prolate spheroidal inclusions were considered: (a) a/b = c/b = 0.5 and (b) a/b = c/b = 0.75 (see Figure 5). Additionally, two different oblate spheroidal inclusions were considered: (a) b/a = c/a = 0.5 and (b) b/a = c/a = 0.75 (see Figure 5).
The elastic constants for the isotropic matrix and the isotropic inclusions are listed in Table 2. The elastic constants for the isotropic matrix and the orthotropic inclusions are listed in Table 3.
We next examine single isotropic/orthotropic spherical, prolate and oblate spheroidal inclusions in an infinite isotropic matrix subject to remote shear loading, σ o xy , σ o xz or σ o yz , as shown in Figure 6 (also see Figures 1c-e and 5) [24]. The remote applied load can be arbitrarily chosen and was assumed to be σ o xy = σ o xz = σ o yz = 75.76 GPa for convenience purposes only. We considered the same geometry of the single spherical, prolate (with an aspect ratio of 0.5 and 0.75) and oblate (with an aspect ratio of 0.5 and 0.75) spheroidal inclusions in an infinite isotropic matrix under remote shear loading (σ o xy , σ o xz and σ o yz ).
Materials 2021, 14, x FOR PEER REVIEW 6 Secondly, in order to resolve the 'no optimized commercial software' problem plan to develop a semi-commercial VIEM software called the "Volume Integral Equa Method Application Program" (VIEMAP). Table 1 shows the analysis capabilities of MAP including a pre-processor (ViemMesh), a solver (VIEM) and a post-proce (ViemPlot) adapted to solve multiple isotropic/anisotropic inclusion problems in a putationally tractable manner. Figure 3 shows the registered trademark for the VIEM The authors aim to help both university students and researchers create VIEM mo using the VIEMAP more easily than using the standard finite element method (FEM well as solve multiple isotropic/anisotropic inclusion problems in an unbounded isotr medium more accurately and conveniently than the boundary element method (BEM

Three-Dimensional Elastostatic Problems Using the VIEM
In this section, we first examine single isotropic/orthotropic spherical, prolate oblate spheroidal inclusions in an infinite isotropic matrix subject to uniform remote sile loading, σ o xx, as shown in Figure 4 (also see Figures 1b and 5). The remote applied can be arbitrarily chosen and was assumed to be σ o xx = 143.10 GPa for convenience poses only. Two different prolate spheroidal inclusions were considered: (a) a/b = c/b and (b) a/b = c/b = 0.75 (see Figure 5). Additionally, two different oblate spheroidal in sions were considered: (a) b/a = c/a = 0.5 and (b) b/a = c/a = 0.75 (see Figure 5).   The elastic constants for the isotropic matrix and the isotropic inclusions are listed in Table 2. The elastic constants for the isotropic matrix and the orthotropic inclusions are listed in Table 3.    The elastic constants for the isotropic matrix and the isotropic inclusions are listed in Table 2. The elastic constants for the isotropic matrix and the orthotropic inclusions are listed in Table 3.  Three different material properties (Iso_01, Iso_05 and Iso_06) in Table 2 were used in the numerical calculation. The elastic constants for the isotropic matrix and the orthotropic inclusions are listed in Table 4. Table 5 shows various characteristics of the material properties used in the numerical calculation. In order to demonstrate the capability of the volume integral equation method for the three-dimensional anisotropic inclusion problem, three independent elastic constants, c 44 (shear modulus in the yz plane), c 55 (shear modulus in the xz plane) and c 66 (shear modulus in the xy plane), were assumed to be different from each other [25]. We next examine single isotropic/orthotropic spherical, prolate and oblate spheroidal inclusions in an infinite isotropic matrix subject to remote shear loading, σ o xy, σ o xz or σ o yz, as shown in Figure 6 (also see Figures 1c-e and 5) [24]. The remote applied load can be arbitrarily chosen and was assumed to be σ o xy = σ o xz = σ o yz = 75.76 GPa for convenience purposes only. We considered the same geometry of the single spherical, prolate (with an aspect ratio of 0.5 and 0.75) and oblate (with an aspect ratio of 0.5 and 0.75) spheroidal inclusions in an infinite isotropic matrix under remote shear loading (σ o xy, σ o xz and σ o yz).
In Equations (7)-(9), g i m (ξ,x) is the Green's function for the infinite isotropic matrix material and is stated in Equation (6). Thus, the VIEM does not require the use of the Green's function for the orthotropic material of the inclusion. In general, Green's function for an anisotropic material is much more complex than that of an isotropic material [28]. Furthermore, a closed form solution of the generalized Green's function for an anisotropic material is not available in the literature.
In contrast, in the BEM, Green's functions for both the isotropic matrix and the anisotropic inclusions must be specified in the formulation. In particular, special emphasis is placed on the fact that Green's function for the anisotropic material of the inclusions is not required in the VIEM.

Numerical Formulations in the VIEM
The integrands in Equations (3)-(8) contain singularities with different orders due to the singular characteristics of the Green's function at x = ξ (i.e., r = 0). Thus, evaluation of the singular integrals requires special attention. In general, g i m (ξ,x) behaves as 1/r, while its derivatives behave as 1/r 2 as r → 0. It should be noted that only g i m (ξ,x) for the isotropic matrix and its derivatives are required in the VIEM. Furthermore, in the BEM, the Green's function for anisotropic inclusions and their derivatives must also be specified. As a result, this may be a critical drawback to the BEM when solving multiple anisotropic inclusion problems.
In contrast to the BEM, the singularities in the VIEM are integrable (weak). Thus, we have decided to utilize the direct integration scheme stated by Li et al. [29]. Finally, after suitable adjustments, we have succeeded in addressing these weak singular integrands in the volume integral equation formulations.
A comprehensive elaboration for the accurate evaluation of singular integrals using the tetrahedron polar co-ordinates shown in [29] was presented in [17].

A Single Isotropic Spherical Inclusion
In order to examine the accuracy of the numerical results using the VIEM, the numerical results using the VIEM for a single isotropic spherical inclusion were first compared to the analytical solutions [21,30]. We considered a single isotropic spherical inclusion with a radius of 6 mm in an infinite isotropic matrix subject to uniform remote tensile loading, σ xx o , as shown in Figure 4a. It should be noted that the length of the radius can be arbitrarily chosen. In Figure 7, standard 20-node quadratic hexahedral elements were used in the discretization [31]. The number of hexahedral elements, 7560, was determined based on a convergence test. For the seven different material properties (Iso_2, Iso_03, Iso_04, Iso_05, Iso_06, Iso_07 and Iso_08) in Table 2, a comparison was made between the numerical results using the volume integral equation method (VIEM) and the analytical solutions. As shown in Table 5, there was no restriction to Poisson's ratio in the inclusions and matrices of Iso_02 and Iso_06. However, Poison's ratio was 1/3 in both the inclusion and matrix of Iso_03, Iso_04, Iso_05, Iso_07 and Iso_08. Furthermore, for Iso_02, Iso_03, Iso_04 and Iso_05, Young's modulus (E) in the isotropic inclusion was greater than that in the isotropic matrix. For Iso_06, Iso_07 and Iso_08, Young's modulus (E) in the isotropic matrix was greater than that in the isotropic inclusion. Thus, seven material properties representing a diversity of materials were chosen. Excellent agreement was found between the analytical and numerical solutions using the VIEM for the seven different materials considered. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figure 7. It should also be noted that the normalized tensile stress (σ xx /σ o xx ) inside the isotropic spherical inclusions was found to be constant [1,30]. Tables 6-8 show that the percentage differences for the two sets of results are less than 0.1% in seven cases. Figure 8 shows numerical solution by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0 • ≤ θ (see Figure 7) ≤ 360 • ) of the isotropic spherical inclusions with a radius of 6 mm under uniform remote tensile loading.
two sets of results are less than 0.1% in seven cases. Figure 8 shows numerical solution by the volume integral equation method for the normalized tensile stress (σxx/σ o xx) along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0° ≤ θ (see Figure 7) ≤ 360°) of the isotropic spherical inclusions with a radius of 6 mm under uniform remote tensile loading.     In most references, the numerical results for this problem were obtained in one direction. Thus, in order to show the VIEM results more thoroughly, the normalized tensile stress (σxx/σ o xx) using the VIEM was presented along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0° ≤ θ (see Figure 7) ≤ 360°) of the isotropic spherical inclusions. It was determined in Figure 8 that the normalized tensile stress (σxx/σ o xx) inside the isotropic spherical inclusions is constant in all directions considered.

A Single Orthotropic Spherical Inclusion
In order to show the advantages of the volume integral equation method (VIEM), we In most references, the numerical results for this problem were obtained in one direction. Thus, in order to show the VIEM results more thoroughly, the normalized tensile stress (σ xx /σ o xx ) using the VIEM was presented along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0 • ≤ θ (see Figure 7) ≤ 360 • ) of the isotropic spherical inclusions. It was determined in Figure 8 that the normalized tensile stress (σ xx /σ o xx ) inside the isotropic spherical inclusions is constant in all directions considered.

A Single Orthotropic Spherical Inclusion
In order to show the advantages of the volume integral equation method (VIEM), we consider a single orthotropic spherical inclusion with a radius of 6 mm in an infinite isotropic matrix subject to uniform remote tensile loading, σ o xx , as shown in Figure 4a. It should be noted that the length of the radius can be arbitrarily chosen. In Figure 7, standard 20-node quadratic hexahedral elements were used in the discretization [31]. The number of hexahedral elements was 7560, determined based on a convergence test. For this problem, in comparison to the boundary element method (BEM), since the VIEM is not sensitive to the anisotropy of the inclusions, it does not require use of the Green's function for the anisotropic inclusions. Moreover, as opposed to the standard FEM, where it is necessary to discretize the full domain, the orthotropic inclusion only needs to be discretized in the VIEM.
Five different material properties (Ort_1, Ort_02, Ort_03, Ort_04 and Ort_05) in Table 5 were used in the numerical calculation. As shown in Table 5, it was assumed that c 11 > c 22 = c 33 for five orthotropic inclusions. Additionally, c 11 of the inclusion in Ort_03 > c 11 of the inclusion in Ort_02 > c 11 of the inclusion in Ort_01. Furthermore, c 11 of the inclusion in Ort_04 < c 11 of the inclusion in Ort_05 < c 11 of the inclusion in Ort_01. Thus, five material properties representing a diversity of materials were chosen. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figure 7. Moreover, the normalized tensile stress (σ xx /σ o xx ) inside the orthotropic spherical inclusions was found to be constant [1,30]. Table 9 shows the numerical solution by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) inside the orthotropic spherical inclusions. For the inclusions in Ort_01, Ort_02 and Ort_03, the normalized tensile stress (σ xx /σ o xx ) inside the inclusion was greater than 1.0. However, for the inclusions in Ort_04 and Ort_05, the normalized tensile stress (σ xx /σ o xx ) inside the inclusion was less than 1.0. Figure 9 shows the numerical solution by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) along (left) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (right) the circumferential direction (0 • ≤ θ (see Figure 7) ≤ 360 • ) of the orthotropic spherical inclusions with a radius of 6 mm under uniform remote tensile loading. It was determined in Figure 9 that the normalized tensile stress (σ xx /σ o xx ) inside the orthotropic spherical inclusions is constant in all directions considered. (σxx/σ o xx) inside the inclusion was less than 1.0. Figure 9 shows the numerical solution by the volume integral equation method for the normalized tensile stress (σxx/σ o xx) along (left) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (right) the circumferential direction (0° ≤ Ө (see Figure 7) ≤ 360°) of the orthotropic spherical inclusions with a radius of 6 mm under uniform remote tensile loading. It was determined in Figure 9 that the normalized tensile stress (σxx/σ o xx) inside the orthotropic spherical inclusions is constant in all directions considered.

A Single Spheroidal Inclusion Problem under Uniform Remote Tensile Loading
In order to introduce the VIEM as a versatile numerical method, we considered a single isotropic/orthotropic spheroidal inclusion in an infinite isotropic matrix subject to uniform remote tensile loading, σ o xx , as shown in Figure 4b,c. Figure 5 shows an orientation of the spheroidal inclusion.

A Single Spheroidal Inclusion Problem under Uniform Remote Tensile Loading
In order to introduce the VIEM as a versatile numerical method, we considered a single isotropic/orthotropic spheroidal inclusion in an infinite isotropic matrix subject to uniform remote tensile loading, σ o xx, as shown in Figure 4b,c. Figure 5 shows an orientation of the spheroidal inclusion.

A Single Isotropic Prolate Spheroidal Inclusion
Two different prolate spheroidal inclusions are considered: (a) a/b = c/b = 0.5, where b = 6 mm, and (b) a/b = c/b = 0.75, where b = 6 mm (see Figure 5). It should be noted that the length of b (= 6 mm) can be arbitrarily chosen. Figures 10 and 11 show a typical discretized model for the single (a) prolate spheroidal inclusion (a/b = c/b = 0.5 where b = 6 mm) and (b) prolate spheroidal inclusion (a/b = c/b = 0.75 where b = 6 mm) used in the VIEM [31], respectively. A total of 7560 standard 20-node quadratic hexahedral elements were used for the single prolate spheroidal inclusion in Figures 10 and 11. The number of elements, 7560, was determined based on a convergence test.
Eight different isotropic inclusions (from Iso_01 to Iso_08) in Table 2 were used in the numerical calculation. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figures 10 and 11. It should also be noted that the normalized tensile stress (σxx/σ o xx) inside the isotropic prolate spheroidal inclusions was found to be constant [1,30].  For the inclusions in Iso_01, Iso_02, Iso_03, Iso_04 and Iso_05, the normalized tensile stress (σxx/σ o xx) inside the inclusion was greater than 1.0. However, for the inclusions in Iso_06, Iso_07 and Iso_08, the normalized tensile stress (σxx/σ o xx) inside the inclusion was less than 1.0. Figure 12 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σxx/σ o xx) along (i) the x-axis inside (−3 mm ≤ x ≤ 3 mm) and (ii) the circumferential direction (0° ≤ Ө (see Figure 10) ≤ 360°) of the isotropic prolate spheroidal inclusions (a/b = c/b = 0.5 where b = 6 mm) under uniform remote tensile loading. Figure 13 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σxx/σ o xx) along (i) the x-axis inside (−4.5 mm ≤ x ≤ 4.5 mm) and (ii) the circumferential direction (0° ≤ Ө (see Figure 11) ≤ 360°) of the isotropic prolate spheroidal inclusions (a/b = c/b = 0.75 where b = 6 mm) under uniform remote tensile loading. It was determined in Figures 12 and 13 that the normalized tensile stress (σxx/σ o xx) inside the isotropic prolate spheroidal inclusions is constant in all directions considered.
(i) x-axis (ii) circumferential direction (a) Eight different isotropic inclusions (from Iso_01 to Iso_08) in Table 2 were used in the numerical calculation. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figures 10 and 11. It should also be noted that the normalized tensile stress (σ xx /σ o xx ) inside the isotropic prolate spheroidal inclusions was found to be constant [1,30].
Tables 10-12 show numerical solutions by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) inside the isotropic prolate spheroidal inclusions. For the inclusions in Iso_01, Iso_02, Iso_03, Iso_04 and Iso_05, the normalized tensile stress (σ xx /σ o xx ) inside the inclusion was greater than 1.0. However, for the inclusions in Iso_06, Iso_07 and Iso_08, the normalized tensile stress (σ xx /σ o xx ) inside the inclusion was less than 1.0. Figure 12 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) along (i) the x-axis inside (−3 mm ≤ x ≤ 3 mm) and (ii) the circumferential direction (0 • ≤ θ (see Figure 10) ≤ 360 • ) of the isotropic prolate spheroidal inclusions (a/b = c/b = 0.5 where b = 6 mm) under uniform remote tensile loading. Figure 13 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) along (i) the x-axis inside (−4.5 mm ≤ x ≤ 4.5 mm) and (ii) the circumferential direction (0 • ≤ θ (see Figure 11) ≤ 360 • ) of the isotropic prolate spheroidal inclusions (a/b = c/b = 0.75 where b = 6 mm) under uniform remote tensile loading. It was determined in Figures 12 and 13 that the normalized tensile stress (σ xx /σ o xx ) inside the isotropic prolate spheroidal inclusions is constant in all directions considered.     Figure 5). It should be noted that the length of b (=6 mm) can be arbitrarily chosen. Figures 10 and 11 show a typical discretized model for the single (a) prolate spheroidal inclusion (a/b = c/b = 0.5 where b = 6 mm) and (b) prolate spheroidal inclusion (a/b = c/b = 0.75 where b = 6 mm) used in the VIEM [31], respectively. A total of 7560 standard 20-node quadratic hexahedral elements were used for the single prolate spheroidal inclusion in Figures 10 and 11. The number of elements, 7560, was determined based on a convergence test.
Five different orthotropic inclusions (from Ort_01 to Ort_05) in Table 3 were used in the numerical calculation. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figures 10 and 11. It should also be noted that the normalized tensile stress (σ xx /σ o xx ) inside the orthotropic prolate spheroidal inclusions was found to be constant [1,30]. Table 13 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) inside the orthotropic prolate spheroidal inclusions. For the inclusions in Ort_01, Ort_02 and Ort_03, the normalized tensile stress (σ xx /σ o xx ) inside the inclusion was greater than 1.0. However, for the inclusions in Ort_04 and Iso_05, the normalized tensile stress (σ xx /σ o xx ) inside the inclusion was less than 1.0. Figure 14 shows numerical solution by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) along (left) the x-axis inside (−3 mm ≤ x ≤ 3 mm) and (right) the circumferential direction (0 • ≤ θ (see Figure 10) ≤ 360 • ) of the orthotropic prolate spheroidal inclusions (a/b = c/b = 0.5 where b = 6 mm) under uniform remote tensile loading.   roidal inclusions (a/b = c/b = 0.5 where b = 6 mm) under uniform remote tensile loading. Figure 13 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σxx/σ o xx) along (i) the x-axis inside (−4.5 mm ≤ x ≤ 4.5 mm) and (ii) the circumferential direction (0° ≤ Ө (see Figure 11) ≤ 360°) of the isotropic prolate spheroidal inclusions (a/b = c/b = 0.75 where b = 6 mm) under uniform remote tensile loading. It was determined in Figures 12 and 13   (i) x-axis (ii) circumferential direction (a)  (i) x-axis (ii) circumferential direction (c)

A Single Isotropic Oblate Spheroidal Inclusion
In this section, two different oblate spheroidal inclusions are considered: (a) b/a = c/a = 0.5, where a = 6 mm, and (b) b/a = c/a = 0.75, where a = 6 mm (see Figure 5). It should be noted that the length of a (= 6 mm) can be arbitrarily chosen.  Eight different isotropic inclusions (from Iso_01 to Iso_08) in Table 2 were used in the numerical calculation. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figures 16 and 17. It should also be noted that the normalized tensile stress (σ xx /σ o xx ) inside the isotropic oblate spheroidal inclusions was found to be constant [1,30]. Tables 14-16 show numerical solutions by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) inside the isotropic oblate spheroidal inclusions. For the inclusions in Iso_01, Iso_02, Iso_03, Iso_04 and Iso_05, the normalized tensile stress (σ xx /σ o xx ) inside the inclusion was greater than 1.0. However, for the inclusions in Iso_06, Iso_07 and Iso_08, the normalized tensile stress (σ xx /σ o xx ) inside the inclusion was less than 1.0. Figure 18 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0 • ≤ θ (see Figure 16) ≤ 360 • ) of the isotropic oblate spheroidal inclusions (b/a = c/a = 0.5 where a = 6 mm) under uniform remote tensile loading. Figure 19 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0 • ≤ θ (see Figure 17) ≤ 360 • ) of the isotropic oblate spheroidal inclusions (b/a = c/a = 0.75 where a = 6 mm) under uniform remote tensile loading. It was determined in Figures 18 and 19 that the normalized tensile stress (σ xx /σ o xx ) inside the isotropic oblate spheroidal inclusions is constant in all directions considered.
In this section, two different oblate spheroidal inclusions are considered: (a) b/a = c/a = 0.5, where a = 6 mm, and (b) b/a = c/a = 0.75, where a = 6 mm (see Figure 5). It should be noted that the length of a (= 6 mm) can be arbitrarily chosen. Figures 16 and 17 show a typical discretized model for the single (a) oblate spheroidal inclusion (b/a = c/a = 0.5 where a = 6 mm) and (b) oblate spheroidal inclusion (b/a = c/a = 0.75 where a = 6 mm) used in the VIEM [31], respectively. A total of 7560 standard 20-node quadratic hexahedral elements were used for the single oblate spheroidal inclusion in Figures 16 and 17. The number of elements, 7560, was determined based on a convergence test.

A Single Orthotropic Oblate Spheroidal Inclusion
In this section, two different oblate spheroidal inclusions are considered: (a) b/a = c/a = 0.5, where a = 6 mm, and (b) b/a = c/a = 0.75, where a = 6 mm (see Figure 5). It should be noted that the length of a (=6 mm) can be arbitrarily chosen.  should also be noted that the normalized tensile stress (σxx/σ xx) inside the isotropic oblate spheroidal inclusions was found to be constant [1,30]. Tables 14-16 show numerical solutions by the volume integral equation method for the normalized tensile stress (σxx/σ o xx) inside the isotropic oblate spheroidal inclusions. For the inclusions in Iso_01, Iso_02, Iso_03, Iso_04 and Iso_05, the normalized tensile stress (σxx/σ o xx) inside the inclusion was greater than 1.0. However, for the inclusions in Iso_06, Iso_07 and Iso_08, the normalized tensile stress (σxx/σ o xx) inside the inclusion was less than 1.0. Figure 18 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σxx/σ o xx) along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0° ≤ Ө (see Figure 16) ≤ 360°) of the isotropic oblate spheroidal inclusions (b/a = c/a = 0.5 where a = 6 mm) under uniform remote tensile loading. Figure 19 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σxx/σ o xx) along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0° ≤ Ө (see Figure 17) ≤ 360°) of the isotropic oblate spheroidal inclusions (b/a = c/a = 0.75 where a = 6 mm) under uniform remote tensile loading. It was determined in Figures  18 and 19 that the normalized tensile stress (σxx/σ o xx) inside the isotropic oblate spheroidal inclusions is constant in all directions considered.
(i) x-axis (ii) circumferential direction (a)   (i) x-axis (ii) circumferential direction (c)  Five different orthotropic inclusions (from Ort_01 to Ort_05) in Table 3 were used in the numerical calculation. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figures 16 and 17. It should also be noted that the normalized tensile stress (σ xx /σ o xx ) inside the orthotropic oblate spheroidal inclusions was found to be constant [1,30]. Table 17 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) inside the orthotropic oblate pheroidal inclusions. For the inclusions in Ort_01, Ort_02 and Ort_03, the normalized tensile stress (σ xx /σ o xx ) inside the inclusion was greater than 1.0. However, for the inclusions in Ort_04 and Iso_05, the normalized tensile stress (σ xx /σ o xx ) inside the inclusion was less than 1.0. Figure 20 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) along (left) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (right) the circumferential direction (0 • ≤ θ (see Figure 16) ≤ 360 • ) of the orthotropic oblate spheroidal inclusions (b/a = c/a = 0.5 where a = 6 mm) under uniform remote tensile loading. Figure 21 shows numerical solutions by the volume integral equation method for the normalized tensile stress (σ xx /σ o xx ) along (left) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (right) the circumferential direction (0 • ≤ θ (see Figure 17) ≤ 360 • ) of the orthotropic oblate spheroidal inclusions (b/a = c/a = 0.75 where a = 6 mm) under uniform remote tensile loading. It was determined in Figures 20 and 21 that the normalized tensile stress (σ xx /σ o xx ) inside the orthotropic oblate spheroidal inclusions is constant in all directions considered.   integral equation method for the normalized tensile stress (σxx/σ o xx) along (left) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (right) the circumferential direction (0° ≤ θ (see Figure 17) ≤ 360°) of the orthotropic oblate spheroidal inclusions (b/a = c/a = 0.75 where a = 6 mm) under uniform remote tensile loading. It was determined in Figures 20 and 21 that the normalized tensile stress (σxx/σ o xx) inside the orthotropic oblate spheroidal inclusions is constant in all directions considered.  Figure 5) b/a = c/a = 0.75 (See Figure 5       Both the standard finite element method (FEM) and the boundary element method (BEM) are powerful general-purpose tools in the field of numerical analysis. Since the VIEM is a combination of these two methods, it is also highly beneficial to the field of numerical analysis and can play a very important role in solving "inclusion problems". The authors hope that the results using the VIEM cited in this paper will be used as benchmarked data for verifying the results of similar research using other analytical and numerical methods.

VIEM Formulation Applied to Isotropic/Orthotropic Inclusion Problems
The displacements for isotropic spherical, prolate and oblate spheroidal inclusions can be determined from volume integral Equations (3)-(5), while the displacements for orthotropic spherical, prolate and oblate spheroidal inclusions can be determined from volume integral Equations (6)-(8).

A Single Isotropic Spherical Inclusion
We considered a single isotropic spherical inclusion with a radius of 6 mm in an infinite isotropic matrix subject to remote shear loading, σ o xy , σ o xz and σ o yz , as shown in Figure 6a [24]. It should be noted that the length of the radius can be arbitrarily chosen. In Figure 7, standard 20-node quadratic hexahedral elements were used in the discretization [31]. The number of hexahedral elements, 7560, was determined based on a convergence test. Three different material properties (Iso_01, Iso_05 and Iso_06) in Table 2 were used in the numerical calculation. It should be noted that the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the isotropic spherical inclusions were found to be constant, respectively [1]. It should also be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figure 7. Table 18 shows numerical solutions by the volume integral equation method for the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the isotropic spherical inclusions. For the inclusions in Iso_01 and Iso_05, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were greater than 1.0, respectively. However, for the inclusion in Iso_06, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were less than 1.0, respectively. Figure 22 shows numerical solutions by the volume integral equation method for the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0 • ≤ θ (see Figure 7) ≤ 360 • ) of the isotropic spherical inclusions with a radius of 6 mm under remote shear loading. In most references, spherical inclusion problems under uniform remote tensile loading were considered. Thus, in order to show the VIEM results more thoroughly, the normalized shear stresses, (a) σ xy /σ o xy , (b) σ xz /σ o xz and (c) σ yz /σ o yz , using the VIEM were presented along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0 • ≤ θ (see Figure 7) ≤ 360 • ) of the isotropic spherical inclusions.
It was determined in Figure 22 that the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the single isotropic spherical inclusions are constant in all directions considered and are identical to each other. Since isotropic materials have an infinite number of planes of symmetry, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the single isotropic spherical inclusions turned out to be identical to each other.

A Single Orthotropic Spherical Inclusion
In order to show the advantages of the volume integral equation method (VIEM), we considered a single orthotropic spherical inclusion with a radius of 6 mm in an infinite isotropic matrix subject to remote shear loading, σ o xy , σ o xz and σ o yz , as shown in Figure 6a. It should be noted that the length of the radius can be arbitrarily chosen. In Figure 7, standard 20-node quadratic hexahedral elements were used in the discretization [31]. The number of hexahedral elements was 7560, determined based on a convergence test. For this problem, in comparison to the boundary element method (BEM), since the VIEM is not sensitive to the anisotropy of the inclusions, it does not require the use of the Green's function for the anisotropic inclusions. Moreover, as opposed to the standard FEM, where it is necessary to discretize the full domain, the orthotropic inclusion only needs to be discretized in the VIEM.
Two different material properties (Ort_06 and Ort_07) in Table 4 were used in the numerical calculation [25]. As shown in Table 5, it was assumed that c 55 > c 66 > c 44 for two orthotropic inclusions. Additionally, c 44 , c 55 and c 66 of the inclusion were assumed be greater than µ of the matrix in the Ort_06 material, while µ of the matrix was assumed to be greater than c 44 , c 55 and c 66 of the inclusion in the Ort_07 material. Thus, two material properties representing different characteristics were chosen. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figure 7. Moreover, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the orthotropic spherical inclusions were found to be constant, respectively [1]. Table 19 shows numerical solutions by the volume integral equation method for the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the orthotropic spherical inclusions. For the inclusion in Ort_06, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were greater than 1.0, respectively. However, for the inclusion in Ort_07, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were less than 1.0, respectively. Figure 23 shows numerical solutions by the volume integral equation method for the normalized shear stresses (a) σ xy /σ o xy , (b) σ xz /σ o xz and (c) σ yz /σ o yz along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0 • ≤ θ (see Figure 7) ≤ 360 • ) of the orthotropic spherical inclusions with a radius of 6 mm under remote shear loading. It was determined in Figure 23 that the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the orthotropic spherical inclusions are constant in all directions considered and are different from each other. Since orthotropic materials have three planes/axes of symmetry and the independent shear moduli in three planes of symmetry are different from each other (c 55 > c 66 > c 44 ), the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the orthotropic spherical inclusions turned out to be different from each other. Furthermore, since c 55 (shear modulus in the xz plane) is greater than c 66 (shear modulus in the xy plane) and c 66 is greater than c 44 (shear modulus in the yz plane) in the orthotropic inclusions of the Ort_06 and Ort_07 materials, it was determined that the normalized shear stress, σ xz /σ o xz , was greater than the normalized shear stress, σ xy /σ o xy . Furthermore, σ xy /σ o xy was found to be greater than the normalized shear stress, σ yz /σ o yz , inside the orthotropic spherical inclusions.

A Single Spheroidal Inclusion Problem under Remote Shear Loading
In order to introduce the VIEM as a versatile numerical method, we considered a single isotropic/orthotropic spheroidal inclusion in an infinite isotropic matrix subject to remote shear loading, σ o xy , σ o xz and σ o yz , as shown in Figure 6b,c. Figure 5 shows the orientation of the spheroidal inclusion.

A Single Isotropic Prolate Spheroidal Inclusion
Two different prolate spheroidal inclusions are considered: (a) a/b = c/b = 0.5, where b = 6 mm, and (b) a/b = c/b = 0.75, where b = 6 mm (see Figure 5). It should be noted that the length of b (=6 mm) can be arbitrarily chosen. Figures 10 and 11 show a typical discretized model for the single (a) prolate spheroidal inclusion (a/b = c/b = 0.5 where b = 6 mm) and (b) prolate spheroidal inclusion (a/b = c/b = 0.75 where b = 6 mm) used in the VIEM [31], respectively. A total of 7560 standard 20-node quadratic hexahedral elements were used for the single prolate spheroidal inclusion in Figures 10 and 11. The number of elements, 7560, was determined based on a convergence test.
Three different isotropic inclusions (Iso_01, Iso_05 and Iso_06) in Table 2 were used in the numerical calculation. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figures 10 and 11. It should also be noted that the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the isotropic prolate spheroidal inclusions were found to be constant, respectively [1]. Table 20 shows numerical solutions by the volume integral equation method for the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the isotropic prolate spheroidal inclusions. For the inclusions in Iso_01 and Iso_05, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were greater than 1.0, respectively. However, for the inclusion in Iso_06, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were less than 1.0, respectively. Figure Figure 26, the cross-section in the xy plane is identical to the cross-section in the yz plane in the prolate spheroidal inclusion, the normalized shear stress, σ xy /σ o xy , was identical to the normalized shear stress, σ yz /σ o yz , inside the isotropic prolate spheroidal inclusion under remote shear loading.

A Single Orthotropic Prolate Spheroidal Inclusion
Two different prolate spheroidal inclusions are considered: (a) a/b = c/b = 0.5, where b = 6 mm, and (b) a/b = c/b = 0.75, where b = 6 mm (see Figure 5). It should be noted that the length of b (=6 mm) can be arbitrarily chosen. Figures 10 and 11 show a typical discretized model for the single (a) prolate spheroidal inclusion (a/b = c/b = 0.5 where b = 6 mm) and (b) prolate spheroidal inclusion (a/b = c/b = 0.75 where b = 6 mm) used in the VIEM [31], respectively. A total of 7560 standard 20-node quadratic hexahedral elements were used for the single prolate spheroidal inclusion in Figures 10 and 11. The number of elements, 7560, was determined based on a convergence test.   Two different orthotropic inclusions (Ort_06 and Ort_07) in Table 4 were used in the numerical calculation. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figures 10 and 11. It should also be noted that the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the orthotropic prolate spheroidal inclusions were found to be constant, respectively [1]. Table 21 shows numerical solutions by the volume integral equation method for the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the orthotropic prolate spheroidal inclusions. For the inclusion in Ort_06, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were greater than 1.0, respectively. However, for the inclusion in Ort_07, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were less than 1.0, respectively. Figure 27 shows numerical solutions by the volume integral equation method for the normalized shear stresses Furthermore, even though, as shown in Figure 26, the cross-section in the xy plane is identical to the cross-section in the yz plane in the prolate spheroidal inclusion, since c 55 (shear modulus in the xz plane) is greater than c 66 (shear modulus in the xy plane) and c 66 is greater than c 44 (shear modulus in the yz plane) in the orthotropic inclusions of the Ort_06 and Ort_07 materials, the normalized shear stress, σ xy /σ o xy , was different from the normalized shear stress, σ yz /σ o yz , inside the isotropic prolate spheroidal inclusion under remote shear loading.

A Single Isotropic Oblate Spheroidal Inclusion
In this section, two different oblate spheroidal inclusions are considered: (a) b/a = c/a = 0.5, where a = 6 mm, and (b) b/a = c/a = 0.75, where a = 6 mm (see Figure 5). It should be noted that the length of a (=6 mm) can be arbitrarily chosen. Figures 16 and 17 show a typical discretized model for the single (a) oblate spheroidal inclusion (b/a = c/a = 0.5 where a = 6 mm) and (b) oblate spheroidal inclusion (b/a = c/a = 0.75 where a = 6 mm) used in the VIEM [31], respectively. A total of 7560 standard 20-node quadratic hexahedral elements were used for the single oblate spheroidal inclusion in Figures 16 and 17. The number of elements, 7560, was determined based on a convergence test.
Three different isotropic inclusions (Iso_01, Iso_05 and Iso_06) in Table 2 were used in the numerical calculation. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figures 16 and 17. It should also be noted that the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the isotropic oblate spheroidal inclusions were found to be constant, respectively [1]. Table 22 shows numerical solutions by the volume integral equation method for the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the isotropic oblate spheroidal inclusions. For the inclusions in Iso_01 and Iso_05, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were greater than 1.0, respectively. However, for the inclusion in Iso_06, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were less than 1.0, respectively.   Figure 5). It should be noted that the length of b (= 6 mm) can be arbitrarily chosen. Figures 10 and 11 show a typical discretized model for the single (a) prolate spheroidal inclusion (a/b = c/b = 0.5 where b = 6 mm) and (b) prolate spheroidal inclusion (a/b = c/b = 0.75 where b = 6 mm) used in the VIEM [31], respectively. A total of 7560 standard 20-node quadratic hexahedral elements were used for the single prolate spheroidal inclusion in Figures 10 and 11. The number of elements, 7560, was determined based on a convergence test.
Two different orthotropic inclusions (Ort_06 and Ort_07) in Table 4 were used in the numerical calculation. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figures 10 and 11. It should also be noted that the normalized shear stresses (σxy/σ o xy, σxz/σ o xz and σyz/σ o yz) inside the orthotropic prolate spheroidal inclusions were found to be constant, respectively [1]. Table 21 Figure 26, the cross-section in the xy plane is identical to the cross-section in the xz plane in the oblate spheroidal inclusion, the normalized shear stress, σ xy /σ o xy , was identical to the normalized shear stress, σ xz /σ o xz , inside the isotropic oblate spheroidal inclusion under remote shear loading.
x ≤ 4.5 mm) and (ii) the circumferential direction (0° ≤ θ (see Figure 11) ≤ 360°) of the orthotropic prolate spheroidal inclusions (a/b = c/b = 0.75 where b = 6 mm) under remote shear loading, σ o xy, σ o xz and σ o yz. It was determined in Figures 27 and 28 that the normalized shear stresses (σxy/σ o xy, σxz/σ o xz and σyz/σ o yz) inside the orthotropic prolate spheroidal inclusions are constant in all directions considered. Furthermore, even though, as shown in Figure 26, the cross-section in the xy plane is identical to the cross-section in the yz plane in the prolate spheroidal inclusion, since c55 (shear modulus in the xz plane) is greater than c66 (shear modulus in the xy plane) and c66 is greater than c44 (shear modulus in the yz plane) in the orthotropic inclusions of the Ort_06 and Ort_07 materials, the normalized shear stress, σxy/σ o xy, was different from the normalized shear stress, σyz/σ o yz, inside the isotropic prolate spheroidal inclusion under remote shear loading.

A Single Isotropic Oblate Spheroidal Inclusion
In this section, two different oblate spheroidal inclusions are considered: (a) b/a = c/a = 0.5, where a = 6 mm, and (b) b/a = c/a = 0.75, where a = 6 mm (see Figure 5). It should be noted that the length of a (= 6 mm) can be arbitrarily chosen. Figures 16 and 17 show a typical discretized model for the single (a) oblate spheroidal inclusion (b/a = c/a = 0.5 where a = 6 mm) and (b) oblate spheroidal inclusion (b/a = c/a = 0.75 where a = 6 mm) used in the VIEM [31], respectively. A total of 7560 standard 20-node quadratic hexahedral elements were used for the single oblate spheroidal inclusion in Figures 16 and 17. The number of elements, 7560, was determined based on a convergence test.     Figure 26, the cross-section in the xy plane is identical to the cross-section in the xz plane in the oblate spheroidal inclusion, the normalized shear stress, σxy/σ o xy, was identical to the normalized shear stress, σxz/σ o xz, inside the isotropic oblate spheroidal inclusion under remote shear loading.

A Single Orthotropic Oblate Spheroidal Inclusion
In this section, two different oblate spheroidal inclusions are considered: (a) b/a = c/a = 0.5, where a = 6 mm, and (b) b/a = c/a = 0.75, where a = 6 mm (see Figure 5). It should be noted that the length of a (=6 mm) can be arbitrarily chosen. Figures 16 and 17 show a typical discretized model for the single (a) oblate spheroidal inclusion (b/a = c/a = 0.5 where a = 6 mm) and (b) oblate spheroidal inclusion (b/a = c/a = 0.75 where a = 6 mm) used in the VIEM [31], respectively. A total of 7560 standard 20-node quadratic hexahedral elements were used for the single oblate spheroidal inclusion in Figures 16 and 17. The number of elements, 7560, was determined based on a convergence test.

A Single Orthotropic Oblate Spheroidal Inclusion
In this section, two different oblate spheroidal inclusions are considered: (a) b/a = c/a = 0.5, where a = 6 mm, and (b) b/a = c/a = 0.75, where a = 6 mm (see Figure 5). It should be noted that the length of a (=6 mm) can be arbitrarily chosen. Figures 16 and 17 show a typical discretized model for the single (a) oblate spheroidal inclusion (b/a = c/a = 0.5 where a = 6 mm) and (b) oblate spheroidal inclusion (b/a = c/a = 0.75 where a = 6 mm) used in the VIEM [31], respectively. A total of 7560 standard 20-node quadratic hexahedral elements were used for the single oblate spheroidal inclusion in Figures 16 and 17. The number of elements, 7560, was determined based on a convergence test.
Two different orthotropic inclusions (Ort_06 and Ort_07) in Table 4 were used in the numerical calculation. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figures 16 and 17. It should also be noted that the normalized shear stresses (σxy/σ o xy, σxz/σ o xz and σyz/σ o yz) inside the orthotropic oblate spheroidal inclusions were found to be constant, respectively [1]. Table 23 shows numerical solutions by the volume integral equation method for the normalized shear stresses (σxy/σ o xy, σxz/σ o xz and σyz/σ o yz) inside the orthotropic oblate spheroidal Two different orthotropic inclusions (Ort_06 and Ort_07) in Table 4 were used in the numerical calculation. It should be noted that the VIEM results represent average values of the normalized stresses in all the nodes of the VIEM model in Figures 16 and 17. It should also be noted that the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the orthotropic oblate spheroidal inclusions were found to be constant, respectively [1]. Table 23 shows numerical solutions by the volume integral equation method for the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the orthotropic oblate spheroidal inclusions. For the inclusion in Ort_06, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were greater than 1.0, respectively. However, for the inclusion in Ort_07, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were less than 1.0, respectively. Figure 31 shows numer-ical solutions by the volume integral equation method for the normalized shear stresses (a) σ xy /σ o xy , (b) σ xz /σ o xz and (c) σ yz /σ o yz along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0 • ≤ θ (see Figure 16) ≤ 360 • ) of the orthotropic oblate spheroidal inclusions (b/a = c/a = 0.5 where a = 6 mm) under remote shear loading, σ o xy , σ o xz and σ o yz . Figure 32 shows numerical solutions by the volume integral equation method for the normalized shear stresses (a) σ xy /σ o xy , (b) σ xz /σ o xz and (c) σ yz /σ o yz along (i) the x-axis inside (−6 mm ≤ x ≤ 6 mm) and (ii) the circumferential direction (0 • ≤ θ (see Figure 17) ≤ 360 • ) of the orthotropic oblate spheroidal inclusions (b/a = c/a = 0.75 where a = 6 mm) under remote shear loading, σ o xy , σ o xz and σ o yz . It was determined in Figures 31 and 32 that the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the orthotropic oblate spheroidal inclusions are constant in all directions considered. Furthermore, even though, as shown in Figure 26, the cross-section in the xy plane is identical to the cross-section in the xz plane in the oblate spheroidal inclusion, since c 55 (shear modulus in the xz plane) is greater than c 66 (shear modulus in the xy plane) and c 66 is greater than c 44 (shear modulus in the yz plane) in the orthotropic inclusions of the Ort_06 and Ort_07 materials, the normalized shear stress, σ xy /σ o xy , was different from the normalized shear stress, σ xz /σ o xz , inside the orthotropic oblate spheroidal inclusion under remote shear loading.    Tables 17-22, it was determined that if the inclusion is harder than the matrix, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion are greater than 1.0, respectively. It was also determined that if the inclusion is softer than the matrix, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion are less than 1.0, respectively.
From Figure 26, notable similarities are observed for isotropic inclusions. First, the cross-section in the xy plane of the isotropic prolate spheroidal inclusion is identical to the cross-section in the yz plane and is symmetrical to the cross-sections in the xy and xz planes of the isotropic oblate spheroidal inclusion. Second, the normalized shear stress, σ xy /σ o xy , inside the isotropic prolate spheroidal inclusion is identical to both the normalized shear stress, σ yz /σ o yz , inside the isotropic prolate spheroidal inclusion and the normalized shear stresses, σ xy /σ o xy and σ xz /σ o xz , inside the isotropic oblate spheroidal inclusion under remote shear loading. Third, the cross-section in the xz plane of the isotropic prolate spheroidal inclusion is symmetrical to the cross-section in the yz plane of the isotropic oblate spheroidal inclusion. Fourth, the normalized shear stress, σ xz /σ o xz , inside the isotropic prolate spheroidal inclusion is identical to the normalized shear stress, σ yz /σ o yz , inside the isotropic oblate spheroidal inclusion under remote shear loading.
In contrast, certain differences can be seen for orthotropic inclusions. First, although the cross-section in the xy plane of the orthotropic prolate spheroidal inclusion is still symmetrical to the cross-section in the xy plane of the orthotropic oblate spheroidal inclusion, it is no longer identical to the cross-section in the yz plane of the orthotropic prolate spheroidal inclusion. Second, since the cross-section in the xy plane of the orthotropic prolate spheroidal inclusion is no longer symmetrical to the cross-section in the xz plane of the orthotropic oblate spheroidal inclusion, the normalized shear stress, σ xy /σ o xy , inside the orthotropic prolate spheroidal inclusion is only identical to the normalized shear stress, σ xy /σ o xy , inside the orthotropic oblate spheroidal inclusion under remote shear loading.
Third, since the cross-section in the xz plane of the orthotropic prolate spheroidal inclusion is no longer symmetrical to the cross-section in the yz plane of the orthotropic oblate spheroidal inclusion, the normalized shear stress, σ xz /σ o xz , inside the orthotropic prolate spheroidal inclusion is not identical to the normalized shear stress, σ yz /σ o yz , inside the orthotropic oblate spheroidal inclusion under remote shear loading.      Tables 17-22, it was determined that if the inclusion is harder than the matrix, the normalized shear stresses (σxy/σ o xy, σxz/σ o xz and σyz/σ o yz) inside the inclusion are greater than 1.0, respectively. It was also determined that if the inclusion is softer than the matrix, the normalized shear stresses (σxy/σ o xy, σxz/σ o xz and σyz/σ o yz) inside the inclusion are less than 1.0, respectively.
From Figure 26, notable similarities are observed for isotropic inclusions. First, the cross-section in the xy plane of the isotropic prolate spheroidal inclusion is identical to the cross-section in the yz plane and is symmetrical to the cross-sections in the xy and xz planes of the isotropic oblate spheroidal inclusion. Second, the normalized shear stress, σxy/σ o xy, inside the isotropic prolate spheroidal inclusion is identical to both the normalized shear stress, σyz/σ o yz, inside the isotropic prolate spheroidal inclusion and the normalized shear stresses, σxy/σ o xy and σxz/σ o xz, inside the isotropic oblate spheroidal inclusion under remote shear loading. Third, the cross-section in the xz plane of the isotropic prolate It should be noted that, through numerical analysis using the volume integral equation method, we could quantitatively verify two qualitative predictions: (1) the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the orthotropic spherical inclusions are different from each other, and (2) for orthotropic spheroidal inclusions, there exists only one symmetrical cross-section when the remote loadings are shear (σ o xy , σ o xz and σ o yz ).
It was determined that values of the normalized tensile stress (σ xx /σ o xx ) or the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the isotropic spheroidal inclusions differed significantly from those inside the orthotropic spheroidal inclusions. Therefore, thorough investigation of spheroidal inclusion problems requires stress analysis for both anisotropic spheroidal inclusion problems and isotropic spheroidal inclusion problems.
We also considered multiple isotropic/anisotropic spheroidal inclusions in an infinite isotropic matrix subject to uniform remote tensile loading, σ o xx . In a future paper, the authors will introduce the VIEM solutions of multiple isotropic/orthotropic spheroidal inclusions in an infinite isotropic matrix under arbitrary loading conditions. It is obvious that general characteristics of multiple isotropic/anisotropic inclusion problems cannot be fully analyzed from the basic characteristics of the corresponding single or two isotropic/anisotropic inclusion problems. Therefore, applying multiple inclusion problems to a wide class of real composite materials and structures requires extending the analysis to multiple isotropic/anisotropic inclusions of different shapes.
Both the standard finite element method (FEM) and the boundary element method (BEM) are powerful general-purpose tools in the field of numerical analysis. Since the VIEM is a combination of these two methods, it is also highly beneficial to the field of numerical analysis and can play a very important role in solving "multiple inclusion problems". The authors hope that the results using the VIEM cited in this paper will be used as benchmarked data for verifying the results of similar research using other analytical and numerical methods.

Conclusions
In order to introduce the VIEM as a versatile numerical method for the three-dimensional elastostatic inclusion problem, it was applied to a class of three-dimensional elastostatic inclusion problems. We first considered single isotropic/orthotropic spherical, prolate (with an aspect ratio of 0.5 and 0.75) and oblate (with an aspect ratio of 0.5 and 0.75) spheroidal inclusions in an infinite isotropic matrix under uniform remote tensile loading. Thirteen inclusions with different characteristics were considered in the numerical calculation. Excellent agreement was found between the analytical and numerical solutions using the VIEM for single isotropic spherical inclusion problems. It was determined that the normalized tensile stress (σ xx /σ o xx ) inside the isotropic/orthotropic spherical, prolate and oblate spheroidal inclusions was constant in two different directions (x-axis and circumferential direction). When the inclusion is harder than the matrix, the normalized tensile stress (σ xx /σ o xx ) inside the inclusion can be arranged in ascending order of magnitude: (1) prolate spheroidal inclusion (a/b = c/b = 0.5), (2) prolate spheroidal inclusion (a/b = c/b = 0.75), (3) sphere, (4) oblate spheroidal inclusion (b/a = c/a = 0.75) and (5) oblate spheroidal inclusion (b/a = c/a = 0.5).
We next considered single isotropic/orthotropic spherical, prolate (with an aspect ratio of 0.5 and 0.75) and oblate (with an aspect ratio of 0.5 and 0.75) spheroidal inclusions in an infinite isotropic matrix under remote shear loading. Five inclusions with different characteristics were considered in the numerical calculation. It was determined that the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the isotropic/orthotropic spherical, prolate and oblate spheroidal inclusions were constant in two different directions (x-axis and circumferential direction), respectively. When the inclusion was harder than the matrix, the normalized shear stresses (σ xy /σ o xy , σ xz /σ o xz and σ yz /σ o yz ) inside the inclusion were greater than 1.0, respectively. Furthermore, for isotropic spheroidal inclusions, there existed two identical or symmetrical cross-sections, while for orthotropic spheroidal inclusions, there existed only one symmetric cross-section when the remote loadings were shear (σ o xy , σ o xz and σ o yz ).
It is the authors' hope that the present solutions for various types of inclusions with different material properties under different loading conditions using the parallel volume integral equation method will be established as reference values for verifying the results of other analytical and numerical methods.
It was also determined that applying multiple inclusion problems to a wide class of real composite materials and structures requires extending the analysis to multiple isotropic/anisotropic inclusions of different numbers and shapes. The parallel volume integral equation method (PVIEM) is now generally more applicable and executable than the standard finite element or boundary element methods. Subsequently, the PVIEM can be used to calculate other quantities of practical interest in realistic models of com-posites containing isotropic or anisotropic inclusions of arbitrary shapes under arbitrary loading conditions. It should also be pointed out that, since the VIEM is a combination of the FEM and the BEM, it may have an unknown advantage that neither the FEM model nor the BEM model alone possess. For example, although certain VIEM models are incorrect from the point of view of the standard FEM only, they can be correctly implemented in the VIEM. In a future paper, the authors will attempt to provide more distinct examples to support this new finding. Finally, as a new machine learning-based predictive framework has been proposed for the accurate and efficient evaluation of singular integrals in the boundary element method (BEM) [32], of particular interest to researchers going forward will be the development of a general-purpose machine learning framework for predicting singular integrals [29] in the volume integral equation method.