Next Article in Journal
Solutions of the Two-Wave Interactions in Quadratic Nonlinear Media
Next Article in Special Issue
Numerical Analysis of Convection–Diffusion Using a Modified Upwind Approach in the Finite Volume Method
Previous Article in Journal
A Class of Sixth Order Viscous Cahn-Hilliard Equation with Willmore Regularization in ℝ3
Previous Article in Special Issue
On Semi-Analytical Solutions for Linearized Dispersive KdV Equations
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Department of Mechanical and Design Engineering, Hongik University, Sejong City 30016, Korea
*
Author to whom correspondence should be addressed.
Mathematics 2020, 8(11), 1866; https://doi.org/10.3390/math8111866
Submission received: 18 September 2020 / Revised: 20 October 2020 / Accepted: 22 October 2020 / Published: 26 October 2020
(This article belongs to the Special Issue Numerical Modeling and Analysis)

Abstract

:
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.

1. 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.

2. 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 umo(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:
u m ( x ) = u m o ( x ) + R [ δ ρ ω 2 g i m ( ξ , x ) u i ( ξ ) δ c i j k l g i , j m ( ξ , x ) u k , l ( ξ ) ] d ξ
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 gim(ξ,x) denotes the elastodynamic Green’s function for the unbounded homogeneous matrix material. The value gim(ξ,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,jm(ξ,x) will be in a later section.
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:
u m ( x ) = u m o ( x ) R δ c i j k l g i , j m ( ξ , x ) u k , l ( ξ ) d ξ
where the integral is over entire space R and δcijkl = cijkl(1) − cijkl(2). The value gim(ξ,x) is the elastostatic Green’s function (or Kelvin’s solution) for the unbounded homogeneous matrix material. The value gim(ξ,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 gim(ξ,x) and gi,jm(ξ,x) while the elastostatic volume integral Equation (2) requires only gi,jm(ξ,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].
It should be noted that Lee and his co-workers e.g., [4,5,7,10,11,14] have been developing the VIEM based on numerical modeling and analysis, while Buryachenko e.g., [12,15,16] has been performing research more mathematically.
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.

3. Numerical Analysis Based On the Three-Dimensional Volume Integral Equation Method

3.1. 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, σoxx, as shown in Figure 4 (see also Figure 2b).
We first consider a single isotropic cubic inhomogeneity (−4 mm ≤ x, y, z ≤ 4 mm) in the unbounded isotropic matrix. Four models with a different number of quadratic hexahedral elements for the isotropic inhomogeneity are used for the convergence test; 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). Figure 4 shows typical discretized models used in the VIEM [17].
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 Table 1 and Table 2. In order to distinguish different material properties easily, we assign a distinct material name (mat_01, mat_02, mat_03, ---) to each material.

3.1.1. Volume Integral Equation for an Unbounded Isotropic Matrix Containing a Single Isotropic Cubic Inhomogeneity

The coordinate axes are x1, x2, x3. The stress–strain relationships (σα = cαβεβ (α, β = 1, 6)) for the isotropic inhomogeneities can be expressed in the form:
{ σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 } = [ λ + 2 μ λ λ 0 0 0 λ λ + 2 μ λ 0 0 0 λ λ λ + 2 μ 0 0 0 0 0 0 μ 0 0 0 0 0 0 μ 0 0 0 0 0 0 μ ] { ε 1 ε 2 ε 3 ε 4 ε 5 ε 6 } .
It should be noted that the stress–strain relationships are also referred to as constitutive equations.
In Equation (3), the stiffness constants cij are related to the engineering elastic constants through c11 = c22 = c33 = (λ + 2μ), c12 = c13 = c23 = λ, and c44 = c55 = c66 = μ.
The displacement components in the volume integral Equation (2) for multiple isotropic inhomogeneities can be expressed in the form:
u 1 ( x ) = u 1 o ( x ) V { δ ( λ + 2 μ ) g 1 , 1 1 u 1 , 1 + δ λ ( g 1 , 1 1 u 2 , 2 + g 2 , 2 1 u 1 , 1 ) + δ λ ( g 1 , 1 1 u 3 , 3 + g 3 , 3 1 u 1 , 1 ) + δ ( λ + 2 μ ) g 2 , 2 1 u 2 , 2 + δ λ ( g 2 , 2 1 u 3 , 3 + g 3 , 3 1 u 2 , 2 ) + δ ( λ + 2 μ ) g 3 , 3 1 u 3 , 3 + δ μ [ g 2 , 3 1 ( u 2 , 3 + u 3 , 2 ) + g 3 , 2 1 ( u 2 , 3 + u 3 , 2 ) ] + δ μ [ g 1 , 3 1 ( u 1 , 3 + u 3 , 1 ) + g 3 , 1 1 ( u 1 , 3 + u 3 , 1 ) ] + δ μ [ g 1 , 2 1 ( u 1 , 2 + u 2 , 1 ) + g 2 , 1 1 ( u 1 , 2 + u 2 , 1 ) ] } d ξ 1 d ξ 2 d ξ 3
u 2 ( x ) = u 2 o ( x ) V { δ ( λ + 2 μ ) g 1 , 1 2 u 1 , 1 + δ λ ( g 1 , 1 2 u 2 , 2 + g 2 , 2 2 u 1 , 1 ) + δ λ ( g 1 , 1 2 u 3 , 3 + g 3 , 3 2 u 1 , 1 ) + δ ( λ + 2 μ ) g 2 , 2 2 u 2 , 2 + δ λ ( g 2 , 2 2 u 3 , 3 + g 3 , 3 2 u 2 , 2 ) + δ ( λ + 2 μ ) g 3 , 3 2 u 3 , 3 + δ μ [ g 2 , 3 2 ( u 2 , 3 + u 3 , 2 ) + g 3 , 2 2 ( u 2 , 3 + u 3 , 2 ) ] + δ μ [ g 1 , 3 2 ( u 1 , 3 + u 3 , 1 ) + g 3 , 1 2 ( u 1 , 3 + u 3 , 1 ) ] + δ μ [ g 1 , 2 2 ( u 1 , 2 + u 2 , 1 ) + g 2 , 1 2 ( u 1 , 2 + u 2 , 1 ) ] } d ξ 1 d ξ 2 d ξ 3
u 3 ( x ) = u 3 o ( x ) V { δ ( λ + 2 μ ) g 1 , 1 3 u 1 , 1 + δ λ ( g 1 , 1 3 u 2 , 2 + g 2 , 2 3 u 1 , 1 ) + δ λ ( g 1 , 1 3 u 3 , 3 + g 3 , 3 3 u 1 , 1 ) + δ ( λ + 2 μ ) g 2 , 2 3 u 2 , 2 + δ λ ( g 2 , 2 3 u 3 , 3 + g 3 , 3 3 u 2 , 2 ) + δ ( λ + 2 μ ) g 3 , 3 3 u 3 , 3 + δ μ [ g 2 , 3 3 ( u 2 , 3 + u 3 , 2 ) + g 3 , 2 3 ( u 2 , 3 + u 3 , 2 ) ] + δ μ [ g 1 , 3 3 ( u 1 , 3 + u 3 , 1 ) + g 3 , 1 3 ( u 1 , 3 + u 3 , 1 ) ] + δ μ [ g 1 , 2 3 ( u 1 , 2 + u 2 , 1 ) + g 2 , 1 3 ( u 1 , 2 + u 2 , 1 ) ] } d ξ 1 d ξ 2 d ξ 3
where u1(x), u2(x) and u3(x) are the three-dimensional displacement components, δcαβ = cαβ(1) − cαβ(2) (α, β = 1, 6), where cαβ(1) represent the elastic stiffness constants of the isotropic inhomogeneities while cαβ(2) denote those for the isotropic matrix material; δc11 = δc22 = δc33 = (λ1 + 2μ1) − (λ2 + 2μ2), δc12 = δc13 = δc23 = λ1 − λ2, and δc44 = δc55 = δc66 = μ1 − μ2.
In Equations (4)–(6), gim(ξ,x) shows the Green’s function for the unbounded isotropic matrix material and is presented by Banerjee [18] and Pao and Varatharajulu [19] as:
g 1 1 = 1 16 π ( 1 ν ) μ r [ ( x 1 ξ 1 ) 2 r 2 + ( 3 4 ν ) ]
g 2 1 = g 1 2 = 1 16 π ( 1 ν ) μ r [ ( x 1 ξ 1 ) ( x 2 ξ 2 ) r 2 ]
g 3 1 = g 1 3 = 1 16 π ( 1 ν ) μ r [ ( x 1 ξ 1 ) ( x 3 ξ 3 ) r 2 ]
g 2 2 = 1 16 π ( 1 ν ) μ r [ ( x 2 ξ 2 ) 2 r 2 + ( 3 4 ν ) ]
g 3 2 = g 2 3 = 1 16 π ( 1 ν ) μ r [ ( x 2 ξ 2 ) ( x 3 ξ 3 ) r 2 ]
g 3 3 = 1 16 π ( 1 ν ) μ r [ ( x 3 ξ 3 ) 2 r 2 + ( 3 4 ν ) ]
where r = | xξ | = ( x 1 ξ 1 ) 2 + ( x 2 ξ 2 ) 2 + ( x 3 ξ 3 ) 2 , ν is Poisson’s ratio, and μ is the shear modulus.
Also, gi,jm(ξ,x) is presented as:
g 1.1 1 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) 3 r 3 + ( 1 4 ν ) ( x 1 ξ 1 ) r }
g 1.2 1 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) 2 ( x 2 ξ 2 ) r 3 + ( 3 4 ν ) ( x 2 ξ 2 ) r }
g 1.3 1 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) 2 ( x 3 ξ 3 ) r 3 + ( 3 4 ν ) ( x 3 ξ 3 ) r }
g 2.1 1 = g 1 , 1 2 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) 2 ( x 2 ξ 2 ) r 3 ( x 2 ξ 2 ) r }
g 2.2 1 = g 1 , 2 2 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) ( x 2 ξ 2 ) 2 r 3 ( x 1 ξ 1 ) r }
g 2.3 1 = g 1 , 3 2 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) ( x 2 ξ 2 ) ( x 3 ξ 3 ) r 3 }
g 3.1 1 = g 1 , 1 3 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) 2 ( x 3 ξ 3 ) r 3 ( x 3 ξ 3 ) r }
g 3.2 1 = g 1 , 2 3 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) ( x 2 ξ 2 ) ( x 3 ξ 3 ) r 3 }
g 3.3 1 = g 1 , 3 3 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) ( x 3 ξ 3 ) 2 r 3 ( x 1 ξ 1 ) r }
g 2.1 2 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) ( x 2 ξ 2 ) 2 r 3 + ( 3 4 ν ) ( x 1 ξ 1 ) r }
g 2.2 2 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 2 ξ 2 ) 3 r 3 + ( 1 4 ν ) ( x 2 ξ 2 ) r }
g 2.3 2 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 2 ξ 2 ) 2 ( x 3 ξ 3 ) r 3 + ( 3 4 ν ) ( x 3 ξ 3 ) r }
g 3 , 1 2 = g 2 , 1 3 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) ( x 2 ξ 2 ) ( x 3 ξ 3 ) r 3 }
g 3 , 2 2 = g 2 , 2 3 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 2 ξ 2 ) 2 ( x 3 ξ 3 ) r 3 ( x 3 ξ 3 ) r }
g 3 , 3 2 = g 2 , 3 3 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 2 ξ 2 ) ( x 3 ξ 3 ) 2 r 3 ( x 2 ξ 2 ) r }
g 3.1 3 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 1 ξ 1 ) ( x 3 ξ 3 ) 2 r 3 + ( 3 4 ν ) ( x 1 ξ 1 ) r }
g 3.2 3 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 2 ξ 2 ) ( x 3 ξ 3 ) 2 r 3 + ( 3 4 ν ) ( x 2 ξ 2 ) r }
g 3.3 3 = 1 16 π ( 1 ν ) μ r 2 { 3 ( x 3 ξ 3 ) 3 r 3 + ( 1 4 ν ) ( x 3 ξ 3 ) r }
It should be noted that for the applied uniformly remote stress, the displacement vector u° is of the form:
u 1 o = c 1 x , u 2 o = c 2 y , u 3 o = c 3 z
where the constants C1–C3 are related to the tensile and shear components of the applied uniformly remote stress.

3.1.2. 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:
{ σ 1 σ 2 σ 3 σ 4 σ 5 σ 6 } = [ c 11 c 12 c 13 0 0 0 c 12 c 22 c 23 0 0 0 c 13 c 23 c 33 0 0 0 0 0 0 c 44 0 0 0 0 0 0 c 55 0 0 0 0 0 0 c 66 ] { ε 1 ε 2 ε 3 ε 4 ε 5 ε 6 } .
The stiffness constants cij are related to the engineering elastic constants through:
c 11 = 1 ν 23 ν 32 E 2 E 3 Δ , c 12 = ν 21 + ν 31 ν 23 E 2 E 3 Δ = c 21 , c 13 = ν 31 + ν 21 ν 32 E 2 E 3 Δ = c 31 , c 22 = 1 ν 31 ν 13 E 3 E 1 Δ , c 23 = ν 32 + ν 31 ν 12 E 3 E 1 Δ = c 32 , c 33 = 1 ν 12 ν 12 E 1 E 2 Δ , c 44 = μ 23 , c 55 = μ 31 , c 66 = μ 12 ,
where,
Δ = 1 ν 12 ν 21 ν 23 ν 32 ν 31 ν 13 2 ν 12 ν 23 ν 31 E 1 E 2 E 3 .
In the above equations, E1, E2, E3 are Young’s moduli in 1, 2 and 3 directions, respectively, μ23, μ31 and μ12 are the shear moduli in the 2-3, 3-1 and 1-2 planes, respectively. νij is Poisson’s ratio for transverse strain in the j-direction when stressed in the i-direction. It should be noted that the elastic moduli satisfy the reciprocal νij/Eiji/Ej.
The displacement components in the volume integral Equation (2) for multiple orthotropic inhomogeneities can be expressed in the form:
u 1 ( x ) = u 1 o ( x ) V { δ c 11 g 1 , 1 1 u 1 , 1 + δ c 12 ( g 1 , 1 1 u 2 , 2 + g 2 , 2 1 u 1 , 1 ) + δ c 13 ( g 1 , 1 1 u 3 , 3 + g 3 , 3 1 u 1 , 1 ) + δ c 22 g 2 , 2 1 u 2 , 2 + δ c 23 ( g 2 , 2 1 u 3 , 3 + g 3 , 3 1 u 2 , 2 ) + δ c 33 g 3 , 3 1 u 3 , 3 + δ c 44 [ g 2 , 3 1 ( u 2 , 3 + u 3 , 2 ) + g 3 , 2 1 ( u 2 , 3 + u 3 , 2 ) ] + δ c 55 [ g 1 , 3 1 ( u 1 , 3 + u 3 , 1 ) + g 3 , 1 1 ( u 1 , 3 + u 3 , 1 ) ] + δ c 66 [ g 1 , 2 1 ( u 1 , 2 + u 2 , 1 ) + g 2 , 1 1 ( u 1 , 2 + u 2 , 1 ) ] } d ξ 1 d ξ 2 d ξ 3
u 2 ( x ) = u 2 o ( x ) V { δ c 11 g 1 , 1 2 u 1 , 1 + δ c 12 ( g 1 , 1 2 u 2 , 2 + g 2 , 2 2 u 1 , 1 ) + δ c 13 ( g 1 , 1 2 u 3 , 3 + g 3 , 3 2 u 1 , 1 ) + δ c 22 g 2 , 2 2 u 2 , 2 + δ c 23 ( g 2 , 2 2 u 3 , 3 + g 3 , 3 2 u 2 , 2 ) + δ c 33 g 3 , 3 2 u 3 , 3 + δ c 44 [ g 2 , 3 2 ( u 2 , 3 + u 3 , 2 ) + g 3 , 2 2 ( u 2 , 3 + u 3 , 2 ) ] + δ c 55 [ g 1 , 3 2 ( u 1 , 3 + u 3 , 1 ) + g 3 , 1 2 ( u 1 , 3 + u 3 , 1 ) ] + δ c 66 [ g 1 , 2 2 ( u 1 , 2 + u 2 , 1 ) + g 2 , 1 2 ( u 1 , 2 + u 2 , 1 ) ] } d ξ 1 d ξ 2 d ξ 3
and,
u 3 ( x ) = u 3 o ( x ) V { δ c 11 g 1 , 1 3 u 1 , 1 + δ c 12 ( g 1 , 1 3 u 2 , 2 + g 2 , 2 3 u 1 , 1 ) + δ c 13 ( g 1 , 1 3 u 3 , 3 + g 3 , 3 3 u 1 , 1 ) + δ c 22 g 2 , 2 3 u 2 , 2 + δ c 23 ( g 2 , 2 3 u 3 , 3 + g 3 , 3 3 u 2 , 2 ) + δ c 33 g 3 , 3 3 u 3 , 3 + δ c 44 [ g 2 , 3 3 ( u 2 , 3 + u 3 , 2 ) + g 3 , 2 3 ( u 2 , 3 + u 3 , 2 ) ] + δ c 55 [ g 1 , 3 3 ( u 1 , 3 + u 3 , 1 ) + g 3 , 1 3 ( u 1 , 3 + u 3 , 1 ) ] + δ c 66 [ g 1 , 2 3 ( u 1 , 2 + u 2 , 1 ) + g 2 , 1 3 ( u 1 , 2 + u 2 , 1 ) ] } d ξ 1 d ξ 2 d ξ 3
where u1(x), u2(x) and u3(x) are the three-dimensional displacement components, δcαβ = cαβ(1) − cαβ(2) (α, β = 1, 6), where cαβ(1) represent the elastic stiffness constants of the orthotropic inhomogeneities while cαβ(2) denote those for the isotropic matrix material; δc11 = c11 − (λ2 + 2μ2), δc22 = c22 − (λ2 + 2μ2), δc33 = c33 − (λ2 + 2μ2), δc12 = c12 − λ2, δc13 = c13 − λ2, δc23 = c23 λ2, and δc44 = c44 − μ2, δc55 = c55 − μ2, δc66 = c66 − μ2.
In Equations (13)–(15), gim(ξ,x) is Green’s function for the unbounded isotropic matrix material. Thus, the VIEM does not need to use Green’s function for the orthotropic material of the inhomogeneity.

3.1.3. 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, gim(ξ,x) behaves as 1/r while its derivatives behave as 1/r2 as r → 0. It ought to be noted that the VIEM involves only gim(ξ,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, Ω1f and Ω2f, as shown in Figure 5b. If the singular point is located at a midside node, then Ωf is divided into three triangular prisms, Ω1f, Ω2f and Ω3f, 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 Ω11f maps onto Ω11f. Subsequently, Ω11f maps onto Ω11f in Figure 5d. The last mapping stands for a linear transformation carrying out numerical integrations in a standard way.
Figure 6 shows a typical hexahedral element used in the VIEM. Standard 20-node quadratic hexahedral elements were used in finite element discretization of the VIEM models. Figure 7 shows a detailed explanation for a subdivision of each triangular prism into three tetrahedral subcells after a subdivision of Ωf′ into two (Ω1f′ and Ω2f′) or three triangular prisms (Ω1f′, Ω2f′ and Ω3f′) (see Figure 5b) for the hexahedral element used in the VIEM (Figure 6).
Table 3 and Table 4 show a detailed explanation for a subdivision of each triangular prism into three tetrahedral subcells after a subdivision of Ωf′ into two or three triangular prisms (see Figure 5b,c), respectively, for the hexahedral element used in the VIEM (Figure 6). Table 3 shows a subdivision of each triangular prism into three tetrahedral subcells after a subdivision of Ωf′ into two triangular prisms (see Figure 5b,c). This subdivision technique is applied to all the corner nodes, 1, 3, 5, 7, 13, 15, 17 and 19 of the hexahedral element shown in Figure 6. Thus, the total number of cases in Table 3 increases to 48 (8 nodes × 6 cases). Table 4 shows a subdivision of each triangular prism into three tetrahedral subcells after a subdivision of Ωf′ into three triangular prisms (see Figure 5b,c). The same subdivision technique is applied to all the midside nodes, 2, 4, 6, 8, 9, 10, 11, 12, 14, 16, 18 and 20 of the hexahedral element shown in Figure 6. Thus, the total number of cases in Table 4 rises to 108 (12 nodes × 9 cases). It should be noted that, in Table 3 and Table 4, 1, 2 or 1, 2 and 3 in the triangular prism column and 1, 2 and 3 in the tetrahedron column are arbitrarily selected. There is no particular order in the numbers, 1, 2 or 1, 2 and 3. Therefore, in order to compute singular integrals accurately, it is necessary to investigate 156 total cases (48 cases in Table 3 and 108 cases in Table 4) per hexahedral element.

3.1.4. 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 (σxxoxx) 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 σoxx = 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 n3 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 (=143) elements.
Figure 9a shows the normalized tensile stress component (σxxoxx) 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 (σxxoxx) 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.
Figure 10, Figure 11 and Figure 12 show the normalized tensile stress component (σxxoxx) along the x-axis and the XX-axis (−8 mm ≤ x, XX ≤ 8 mm) for Model_8 × 8 × 8, Model_10 × 10 × 10, Model_14 × 14 × 14 and Model_16 × 16 × 16 under uniform remote tensile loading for the three different orthotropic inhomogeneities (mat_01, mat_02 and mat_03), respectively. The remote applied load was assumed to be σoxx = 143.1 GPa.
Figure 10 a shows the normalized tensile stress component (σxxoxx) along the x-axis inside the orthotropic cubic inhomogeneity #1 (mat_02 in 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 (σxxoxx) 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.
Figure 11a and Figure 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. Figure 11b and Figure 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 (σxxoxx) for the isotropic cubic inhomogeneity in Figure 9 appears to be considerably different from those of the orthotropic cubic inhomogeneities in Figure 10, Figure 11 and Figure 12.

3.2. A single Spherical Inhomogeneity Problem

3.2.1. 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, σxxo, 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 (σxxoxx) 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 σoxx = 150.0 GPa.
It should be noted that the normalized tensile stress component (σxxoxx) 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:
σ i j i = c i j k l ( 2 ) ( ε k l c ε k l * ) = c i j k l ( 2 ) ( S k l p q ε p q * ε k l * )
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):
S 1111 = S 2222 = S 3333 = ( 7 5 v ) 15 ( 1 v ) , S 1122 = S 2233 = S 3311 = S 1133 = S 3322 = S 2211 = ( 5 v 1 ) 15 ( 1 v ) , S 1212 = S 2323 = S 3131 = ( 4 5 v ) 15 ( 1 v ) .
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.

3.2.2. A Single Orthotropic Spherical Inhomogeneity

We consider a single orthotropic spherical inhomogeneity in an unbounded isotropic matrix under uniform remote tensile loading, σxxo, 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 (σxxoxx) 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 σoxx = 143.1 GPa. It should be noted that the normalized tensile stress component (σxxoxx) inside the orthotropic spherical inhomogeneities was found to be constant [23].

4. 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 (σxxoxx) inside the spherical inhomogeneities was also found to be constant for the single orthotropic spherical inhomogeneity problems.
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.

Author Contributions

Conceptualization, J.L.; VIEM Analysis, J.L.; Investigation, J.L. and M.H.; Validation, M.H. and J.L.; Writing-original draft, J.L. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the Ministry of Trade, Industry and Energy (MOTIE), KOREA, through the Education Program for Creative and Industrial Convergence (No. N0000717) and 2019 Hongik University Research Fund. This work was also supported by the National Supercomputing Center with supercomputing resources including technical support (No. KSC-2020-CRE-0107) and the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2019R1H1A2039647).

Acknowledgments

The authors would like to thank all reviewers and journal editors for their constructive and motivating comments. The authors also thank Jason P. Buschman at Hongik University (Sejong Campus) for his thorough proofreading of the revised manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bose, S.K.; Mal, A.K. Elastic waves in a fiber-reinforced composite. J. Mech. Phys. Solids 1974, 22, 217–229. [Google Scholar] [CrossRef]
  2. Maslov, B.P. Stress concentration in an isotropic matrix bonded by anisotropic fibers. Sov. Appl. Mech. 1987, 23, 963–969. [Google Scholar] [CrossRef]
  3. Yang, R.B.; Mal, A.K. The effective transverse moduli of a composite with degraded fiber-matrix interfaces. Int. J. Eng. Sci. 1995, 33, 1623–1632. [Google Scholar] [CrossRef]
  4. Lee, J.K.; Mal, A.K. A volume integral equation technique for multiple inclusion and crack interaction problems. J. Appl. Mech. 1997, 64, 23–31. [Google Scholar] [CrossRef]
  5. Lee, J.K.; Mal, A.K. A volume integral equation technique for multiple scattering problems in elastodynamics. Appl. Math. Comput. 1995, 67, 135–159. [Google Scholar] [CrossRef]
  6. Nimmer, R.P.; Bankert, R.J.; Russel, E.S.; Smith, G.A.; Wright, P.K. Micromechanical modeling of fiber/matrix interface effects in transversely loaded SiC/Ti-6-4 metal matrix composites. J. Compos. Tech. Res. 1991, 13, 3–13. [Google Scholar]
  7. Lee, J.; Mal, A. Characterization of matrix damage in metal matrix composites under transverse loads. Comput. Mech. 1998, 21, 339–346. [Google Scholar] [CrossRef]
  8. Aghdam, M.M.; Falahatgar, S.R. Micromechanical modeling of interface damage of metal matrix composites subjected to transverse loading. Compos. Struct. 2004, 66, 415–420. [Google Scholar] [CrossRef]
  9. Mal, A.K.; Knopoff, L. Elastic wave velocities in two component systems. J. Inst. Math. Appl. 1967, 3, 376–387. [Google Scholar] [CrossRef]
  10. Lee, J.K.; Lee, H.C.; Jeong, H.G. Multiple scattering using parallel volume integral equation method: Interaction of SH waves with multiple multilayered anisotropic elliptical inclusions. Math. Probl. Eng. 2015. [Google Scholar] [CrossRef] [Green Version]
  11. Lee, J.K.; Lee, H.M.; Mal, A. A mixed volume and boundary integral equation technique for elastic wave field calculations in heterogeneous materials. Wave Motion 2004, 39, 1–19. [Google Scholar] [CrossRef]
  12. Buryachenko, V.A. Micromechanics of Heterogeneous Materials; Springer: New York, NY, USA, 2007. [Google Scholar]
  13. Advances in Computers and Information in Engineering Research; Michopoulos, J.G.; Rosen, D.W.; Paredis, C.J.; Vance, J.M. (Eds.) ASME Press: New York, NY, USA, 2020; Volume 2, in press. [Google Scholar]
  14. Lee, J.K.; Oh, S.M.; Mal, A. Calculation of interfacial stresses in composites containing elliptical inclusions of various types. Eur. J. Mech. A Solid. 2014, 44, 17–40. [Google Scholar] [CrossRef]
  15. Buryachenko, V.A.; Bechel, V.T. A series solution of the volume integral equation for multiple-inclusion interaction problems. Compos. Sci. Technol. 2000, 60, 2465–2469. [Google Scholar] [CrossRef]
  16. Buryachenko, V.A. Solution of general integral equations of micromechanics of heterogeneous materials. Int. J. Solids Struct. 2014, 51, 3823–3843. [Google Scholar] [CrossRef] [Green Version]
  17. PATRAN User’s Manual, version 2008 r2; MSC Sotware: Newport Beach, CA, USA, 2008.
  18. Banerjee, P.K. The Boundary Element Methods in Engineering; McGraw-Hill: London, UK, 1993. [Google Scholar]
  19. Pao, Y.-H.; Varatharajulu, V. Huygens’ principle, radiation conditions, and integral formulas for the scattering of elastic waves. J. Acoust. Soc. Am. 1976, 59, 1361–1369. [Google Scholar] [CrossRef]
  20. Li, H.B.; Han, G.M.; Mang, H.A. A new method for evaluating singular integrals in stress analysis of solids by the direct boundary element method. Int. J. Numer. Methods Eng. 1985, 21, 2071–2098. [Google Scholar] [CrossRef]
  21. Shibata, M.; Ono, K. Stress concentration due to a prolate spheroidal inclusion. J. Compos. Mater. 1978, 12, 132–138. [Google Scholar] [CrossRef]
  22. Eshelby, J.D. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proc. Royal Soc. A 1957, 241, 376–396. [Google Scholar]
  23. Ma, H. Solutions of Eshelby-Type Inclusion Problems and a Related Homogenization Method Based on a Simplified Strain Gradient Elasticity Theory. Ph.D. Thesis, Department of Mechanical Engineering, Texas A&M University, College Station, TX, USA, 2010. [Google Scholar]
Figure 1. Transverse cross-section of a SiC/Ti-15-3 composite.
Figure 1. Transverse cross-section of a SiC/Ti-15-3 composite.
Mathematics 08 01866 g001
Figure 2. Geometry of the general (a) elastodynamic and (b) elastostatic problem.
Figure 2. Geometry of the general (a) elastodynamic and (b) elastostatic problem.
Mathematics 08 01866 g002
Figure 3. The field point and the interior point (including the boundary point).
Figure 3. The field point and the interior point (including the boundary point).
Mathematics 08 01866 g003
Figure 4. Typical discretized cubic models (−4 mm ≤ x, y, z ≤ 4 mm) in the volume integral equation method (VIEM). (a) A single cube (Model_10 × 10 × 10); (b) A single cube (Model_16 × 16 × 16).
Figure 4. Typical discretized cubic models (−4 mm ≤ x, y, z ≤ 4 mm) in the volume integral equation method (VIEM). (a) A single cube (Model_10 × 10 × 10); (b) A single cube (Model_16 × 16 × 16).
Mathematics 08 01866 g004
Figure 5. Application of the tetrahedron polar co-ordinates to quadratic, hexahedral, isoparametric singular elements [20]. (a) Mapping of Ωf onto Ωf; (b) subdivision of Ωf′ at the singular point 1; (c) subdivision of triangular prism Ω1f′ into three tetrahedral; (d) mapping of Ω11f′ onto Ω11f″ and of Ω11f″ onto Ω11f‴.
Figure 5. Application of the tetrahedron polar co-ordinates to quadratic, hexahedral, isoparametric singular elements [20]. (a) Mapping of Ωf onto Ωf; (b) subdivision of Ωf′ at the singular point 1; (c) subdivision of triangular prism Ω1f′ into three tetrahedral; (d) mapping of Ω11f′ onto Ω11f″ and of Ω11f″ onto Ω11f‴.
Mathematics 08 01866 g005
Figure 6. A typical quadratic hexahedral element used in the VIEM.
Figure 6. A typical quadratic hexahedral element used in the VIEM.
Mathematics 08 01866 g006
Figure 7. A subdivision of each triangular prism into three tetrahedral subcells after a subdivision of Ωf′ into two or three triangular prisms (see Figure 5b,c). (a) A subdivision of each triangular prism into three tetrahedral subcells at the singular point 15 after a subdivision of Ωf′ into two triangular prisms (Ω1f′ and Ω2f′) (see Figure 5b,c); (b) a subdivision of each triangular prism into three tetrahedral subcells at the singular point 14 after a subdivision of Ωf′ into three triangular prisms (Ω1f′, Ω2f′ and Ω3f′) (see Figure 5b,c).
Figure 7. A subdivision of each triangular prism into three tetrahedral subcells after a subdivision of Ωf′ into two or three triangular prisms (see Figure 5b,c). (a) A subdivision of each triangular prism into three tetrahedral subcells at the singular point 15 after a subdivision of Ωf′ into two triangular prisms (Ω1f′ and Ω2f′) (see Figure 5b,c); (b) a subdivision of each triangular prism into three tetrahedral subcells at the singular point 14 after a subdivision of Ωf′ into three triangular prisms (Ω1f′, Ω2f′ and Ω3f′) (see Figure 5b,c).
Mathematics 08 01866 g007
Figure 8. A typical single cubic model (−4 mm ≤ x, XX ≤ 4 mm) (between two green lines) including the outside of the cube (−8 mm ≤ x ≤ −4 mm and 4 mm ≤ x ≤ 8 mm) in the VIEM. (a) A single cubic model (−4 mm ≤ x, XX ≤ 4 mm) with 1000 elements (Model_10 × 10 × 10) (between two green lines) including the outside of the cube (−8 mm ≤ x ≤ −4 mm and 4 mm ≤ x ≤ 8 mm). 1000 elements are used for the discretization of the outside of the cube; (b) a single cubic model (−4 mm ≤ x, XX ≤ 4 mm) with 4096 elements (Model_16 × 16 × 16) (between two green lines) including the outside of the cube (−8 mm ≤ x ≤ −4 mm and 4 mm ≤ x ≤ 8 mm). 4096 elements are used for the discretization of the outside of the cube.
Figure 8. A typical single cubic model (−4 mm ≤ x, XX ≤ 4 mm) (between two green lines) including the outside of the cube (−8 mm ≤ x ≤ −4 mm and 4 mm ≤ x ≤ 8 mm) in the VIEM. (a) A single cubic model (−4 mm ≤ x, XX ≤ 4 mm) with 1000 elements (Model_10 × 10 × 10) (between two green lines) including the outside of the cube (−8 mm ≤ x ≤ −4 mm and 4 mm ≤ x ≤ 8 mm). 1000 elements are used for the discretization of the outside of the cube; (b) a single cubic model (−4 mm ≤ x, XX ≤ 4 mm) with 4096 elements (Model_16 × 16 × 16) (between two green lines) including the outside of the cube (−8 mm ≤ x ≤ −4 mm and 4 mm ≤ x ≤ 8 mm). 4096 elements are used for the discretization of the outside of the cube.
Mathematics 08 01866 g008
Figure 9. Numerical results using the volume integral equation method (VIEM) for the normalized tensile stress component (σxxoxx) along the x-axis and the XX-axis (−8 mm ≤ x, XX ≤ 8 mm) for Model_8 × 8 × 8, Model_10 × 10 × 10, Model_14 × 14 × 14 and Model_16 × 16 × 16 under uniform remote tensile loading. (a) x-axis; (b) XX-axis.
Figure 9. Numerical results using the volume integral equation method (VIEM) for the normalized tensile stress component (σxxoxx) along the x-axis and the XX-axis (−8 mm ≤ x, XX ≤ 8 mm) for Model_8 × 8 × 8, Model_10 × 10 × 10, Model_14 × 14 × 14 and Model_16 × 16 × 16 under uniform remote tensile loading. (a) x-axis; (b) XX-axis.
Mathematics 08 01866 g009aMathematics 08 01866 g009b
Figure 10. Numerical results using the volume integral equation method (VIEM) for the normalized tensile stress component (σxxoxx) along the x-axis and the XX-axis (−8 mm ≤ x, XX ≤ 8 mm) for Model_8 × 8 × 8, Model_10 × 10 × 10, Model_14 × 14 × 14 and Model_16 × 16 × 16 under uniform remote tensile loading. (a) x-axis; (b) XX-axis.
Figure 10. Numerical results using the volume integral equation method (VIEM) for the normalized tensile stress component (σxxoxx) along the x-axis and the XX-axis (−8 mm ≤ x, XX ≤ 8 mm) for Model_8 × 8 × 8, Model_10 × 10 × 10, Model_14 × 14 × 14 and Model_16 × 16 × 16 under uniform remote tensile loading. (a) x-axis; (b) XX-axis.
Mathematics 08 01866 g010
Figure 11. Numerical results using the volume integral equation method (VIEM) for the normalized tensile stress component (σxxoxx) along the x-axis and the XX-axis (−8 mm ≤ x, XX ≤ 8 mm) for Model_8 × 8 × 8, Model_10 × 10 × 10, Model_14 × 14 × 14 and Model_16 × 16 × 16 under uniform remote tensile loading. (a) x-axis; (b) XX-axis.
Figure 11. Numerical results using the volume integral equation method (VIEM) for the normalized tensile stress component (σxxoxx) along the x-axis and the XX-axis (−8 mm ≤ x, XX ≤ 8 mm) for Model_8 × 8 × 8, Model_10 × 10 × 10, Model_14 × 14 × 14 and Model_16 × 16 × 16 under uniform remote tensile loading. (a) x-axis; (b) XX-axis.
Mathematics 08 01866 g011
Figure 12. Numerical results using the volume integral equation method (VIEM) for the normalized tensile stress component (σxxoxx) along the x-axis and the XX-axis (−8 mm ≤ x, XX ≤ 8 mm) for Model_8 × 8 × 8, Model_10 × 10 × 10, Model_14 × 14 × 14 and Model_16 × 16 × 16 under uniform remote tensile loading. (a) x-axis; (b) XX-axis.
Figure 12. Numerical results using the volume integral equation method (VIEM) for the normalized tensile stress component (σxxoxx) along the x-axis and the XX-axis (−8 mm ≤ x, XX ≤ 8 mm) for Model_8 × 8 × 8, Model_10 × 10 × 10, Model_14 × 14 × 14 and Model_16 × 16 × 16 under uniform remote tensile loading. (a) x-axis; (b) XX-axis.
Mathematics 08 01866 g012
Figure 13. A typical discretized spherical model in the volume integral equation method (VIEM). (a) A spherical model; (b) an inside view of a spherical model.
Figure 13. A typical discretized spherical model in the volume integral equation method (VIEM). (a) A spherical model; (b) an inside view of a spherical model.
Mathematics 08 01866 g013
Figure 14. Comparison of the numerical results using the volume integral equation method (VIEM) and the analytical solutions for the normalized tensile stress component (σxxoxx) along the x-axis inside the isotropic spherical inclusions with a radius of 7 mm under uniform remote tensile loading.
Figure 14. Comparison of the numerical results using the volume integral equation method (VIEM) and the analytical solutions for the normalized tensile stress component (σxxoxx) along the x-axis inside the isotropic spherical inclusions with a radius of 7 mm under uniform remote tensile loading.
Mathematics 08 01866 g014
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.
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.
Mathematics 08 01866 g015
Table 1. Material properties of the isotropic matrix and the isotropic inhomogeneity.
Table 1. Material properties of the isotropic matrix and the isotropic inhomogeneity.
(Unit: GPa)Isotropic Matrix (mat_01)Isotropic Inhomogeneity (mat_01)
λ67.34176.06
μ37.88176.06
Table 2. Material properties of the isotropic matrix and the isotropic and orthotropic cubic inhomogeneity.
Table 2. Material properties of the isotropic matrix and the isotropic and orthotropic cubic inhomogeneity.
Isotropic Matrix
(mat_01, 02, 03, 04)
Isotropic Inhomogeneity
(mat_01)
Orthotropic Inhomogeneity #1
(mat_02)
Orthotropic Inhomogeneity #2
(mat_03)
Orthotropic Inhomogeneity #3
(mat_04)
c11 [GPa]143.10528.18139.54279.08418.61
c12 [GPa]67.34176.063.907.8011.70
c13 [GPa]67.34176.063.907.8011.70
c22 [GPa]143.10528.1815.2830.5645.83
c23 [GPa]67.34176.063.296.599.88
c33 [GPa]143.10528.1815.2830.5645.83
c44 [GPa]37.88176.065.9011.8017.70
c55 [GPa]37.88176.065.9011.8017.70
c66 [GPa]37.88176.065.9011.8017.70
Table 3. A subdivision of each triangular prism into three tetrahedral subcells at each singular point after a subdivision of Ωf′ into two triangular prisms (see Figure 5b,c). The singular point is located at the corner point of the hexahedral element shown in Figure 6.
Table 3. A subdivision of each triangular prism into three tetrahedral subcells at each singular point after a subdivision of Ωf′ into two triangular prisms (see Figure 5b,c). The singular point is located at the corner point of the hexahedral element shown in Figure 6.
Singular PointTriangular PrismTetrahedronNode 1Node 2Node 3Node 4
11115719
2117519
31131719
2113517
2115317
31131517
31137113
2319713
33151913
2135719
2317519
33151719
51157119
2511319
35131719
2151313
2531513
35151713
71171313
2731513
37151913
2173515
2751715
37171915
13111317153
2135173
313153
211319175
2137195
313175
15111519175
2157195
315375
211513197
2151137
315317
17111715133
2171313
317153
211713191
2171971
317751
19111917155
2191535
319375
211915133
2191313
319173
Table 4. A subdivision of each triangular prism into three tetrahedral subcells at each singular point after a subdivision of Ωf′ into three triangular prisms (see Figure 5b,c). The singular point is located at the midside of the hexahedral element shown in Figure 6.
Table 4. A subdivision of each triangular prism into three tetrahedral subcells at each singular point after a subdivision of Ωf′ into three triangular prisms (see Figure 5b,c). The singular point is located at the midside of the hexahedral element shown in Figure 6.
Singular PointTriangular PrismTetrahedronNode 1Node 2Node 3Node 4
21123515
2251715
32171415
2125719
2217519
32141719
3127113
2219713
32141913
41141313
2431513
34151613
2147113
2419713
34161913
3145719
2417519
34161719
61163515
2651715
36171815
2161315
2613115
36181315
3167113
2619713
36181913
81181315
2813115
38201315
2183515
2851715
38172015
3185717
2871917
38192017
9119131517
29191317
39121917
21915317
293517
3951217
319315
29175
397125
101110131519
210151719
310171119
211011319
2107119
31011719
3110317
210537
3101157
111111171913
211151713
311101513
211119713
2117113
31111013
3111751
211531
3113101
121112171915
212191315
31213915
211251715
2123515
3129315
3112753
212173
312913
14111417155
2141535
314325
211419175
2147195
314275
311413197
2141137
314217
16111619177
2161757
316547
211613197
2161137
316417
311615131
2163151
316431
18111817153
2185173
318653
211815133
2181313
318163
311813191
2181971
318761
20112015133
2201313
320183
212017153
2205173
320853
312019175
2207195
320875
Table 5. Material properties of the isotropic matrix and the isotropic spherical inhomogeneity.
Table 5. Material properties of the isotropic matrix and the isotropic spherical inhomogeneity.
Isotropic Matrix
(mat_05)
Isotropic Inhomogeneity
(mat_05)
Isotropic Matrix
(mat_06)
Isotropic Inhomogeneity
(mat_06)
λ [GPa]75.0150.075.0375.0
μ [GPa]37.575.037.5187.5
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Lee, J.; Han, M. Three-Dimensional Volume Integral Equation Method for Solving Isotropic/Anisotropic Inhomogeneity Problems. Mathematics 2020, 8, 1866. https://doi.org/10.3390/math8111866

AMA Style

Lee J, Han M. Three-Dimensional Volume Integral Equation Method for Solving Isotropic/Anisotropic Inhomogeneity Problems. Mathematics. 2020; 8(11):1866. https://doi.org/10.3390/math8111866

Chicago/Turabian Style

Lee, Jungki, and Mingu Han. 2020. "Three-Dimensional Volume Integral Equation Method for Solving Isotropic/Anisotropic Inhomogeneity Problems" Mathematics 8, no. 11: 1866. https://doi.org/10.3390/math8111866

APA Style

Lee, J., & Han, M. (2020). Three-Dimensional Volume Integral Equation Method for Solving Isotropic/Anisotropic Inhomogeneity Problems. Mathematics, 8(11), 1866. https://doi.org/10.3390/math8111866

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop