An Efficient Numerical Method for Calculation of Elastic and Thermo-Elastic Fields in a Homogeneous Medium with Several Heterogeneous Inclusions

The work is devoted to calculation of static elastic and thermo-elastic fields in a homogeneous medium with a finite number of isolated heterogeneous inclusions. Firstly, the problem is reduced to the solution of integral equations for strain and stress fields in the medium with inclusions. Then, Gaussian approximating functions are used for discretization of these equations. For such functions, the elements of the matrix of the discretized problem are calculated in explicit analytical forms. The method is mesh free, and only the coordinates of the approximating nodes are the geometrical information required in the method. If such nodes compose a regular grid, the matrix of the discretized problem obtains the Toeplitz properties. By the calculation of matrix-vector products with such matrices, the Fast Fourier Transform technique may be used. The latter accelerates essentially the process of the iterative solution of the disretized problem. The results of calculations of elastic fields in 3D-medium with an isolated spherical heterogeneous inclusion are compared with exact solutions. Examples of the calculation of elastic and thermo-elastic fields in the medium with several inclusions are presented.


Introduction
Calculation of elastic fields in homogeneous materials with isolated heterogeneous inclusions is an important problem of stress analysis of composites and materials with defects.Efficient numerical methods of the solution of this problem are based on the volume integral equations for the fields in heterogeneous media (see, e.g., [1][2][3]).By the use of these equations, the fields inside the inclusions become principal unknowns of the problem.If the fields inside the inclusions are known, the fields in the medium are reconstructed from the original integral equations.Thus, the problem should be solved only in the region occupied by the inclusions.This is the main advantage of the Integral Equation Method (IEM) over the Finite Element Method (FEM), where the fields in the medium as well as in the inclusions are equivalent unknowns, and the solution has to be found in the whole region.The IEM is preferable if inclusions being smaller than the characteristic sizes of the body are situated far from its boundary.
A conventional method of the solution of integral equations is based on the following procedure.The region of integration is divided into a finite number of sub-regions, where the unknown functions are approximated by some standard functions (e.g., polynomial splines, radial functions, wavelets, etc.).After application of the Method of Moments or the Collocation Method the problem is reduced to a finite system of linear algebraic equations with respect to the coefficients of the approximation (the discretized problem) (see, e.g., [3]).The elements of the matrix of this system are integrals over the sub-regions.For conventional approximating functions, a great portion of computer time is spent for numerical calculation of these integrals.The matrices of the discretized problems are usually non-sparse and have large dimensions (if high accuracy of the numerical solutions is required).Computational cost of the solution of linear algebraic sys- tems with such matrices is rather high.In application to the 3D-integral equations of elasticity for a medium with inclusions, this traditional numerical scheme was carried out in [4,5].
In this work, an efficient numerical method for the solution of integral equations of elasticity and thermoelasticity for a 3D-homogeneous medium with a finite number of isolated inclusions is developed.For discretization of these equations, the Gaussian radial approximating functions are used.The theory of approximation by Gaussian and similar functions was developed in [6].For such functions, the elements of the matrix of the discretized problems are calculated in explicit analytical forms.Thus, the time of the construction of this matrix is essentially reduced in comparison with the methods that incorporate conventional approximating functions and require numerical integration by the calculation of the matrix elements.For regular grids of approximating nodes, the matrix of the discretized problem has the Toeplitz structure, and the Fast Fourier Transform (FFT) technique may be used for the calculation of matrix-vector products in the process of iterative solution of the discretized problem.The initial "geometrical" information required in the method is only the coordinates of the approximating nodes but not detailed forms of the mesh cells.Thus, the method is mesh free in fact.
For the numerical solution of 2D-volume integral equations of elasticity, a similar method was developed in [7].In the present work, the method is extended on the 3D-volume integral equations of elasticity and thermo elasticity.
The structure of this paper is as follows.In Section 2, the 3D-volume integral equations of elasticity for strain and stress fields in a homogeneous medium with a set of isolated heterogeneous inclusions are considered.In Section 3, Gaussian approximating functions are used for discretization of these equations.Comparisons of the numerical and exact solutions for a medium with a spherical inclusion which elastic properties vary along the radius are presented in Section 4. A medium with several isolated inclusions is also considered in this section.In Section 5, the problem of thermo-elasticity for the medium with inhomogeneities is considered.In the Conclusions, some details of the proposed method and the area of its application are discussed.

Integral Equations of Static Elasticity for a Homogeneous Medium with a Finite Number of Isolated Inclusions
Let an infinite homogeneous medium with the tensor of elastic constants contain a finite number N of isolated inclusions that occupy regions ( ).
Elastic properties inside each region are defined by the tensors , where , , x x x x is a point of the 3D-space.The medium is subjected to an external strain ) field, and the objective is to calculate elastic strain and stress fields in the medium with the inclusions.This problem may be reduced to volume integral equations for the strain   V x be the characteristic function of the region V occupied by the inclusions:  .The strain field   x  in the medium with the inclusions satisfies the following integral equation (see, e.g., [1,8]): Here Summation with respect to repeating tensorial (low) indices is implied.The kernel   K x of the integral operator in this equation is the second derivative of the Green function where   G x is the solution of the following equation: Here ij  is Kronecker's symbol,   x  is Dirac's delta-function.The parentheses in indices mean symmetrization.
A similar equation may be written for the stress field σ(x) in the medium with the inclusions ( [1,8]): Note that the functions   K x and in (2) and (5) behave as and therefore, the integral operators K and S in (1) and (4) are singular.Regularizations of these operators on continuous tensorfunctions with a finite support are indicated in [1,9].Let   ij f x be a smooth tensor-function which Fourier transform is bounded and tends to zero at infinity as 3 k  or faster.In this case, the actions of the operators K and S on such a function are defined by the equations: The integrals in the right-hand sides of these equations exist as ordinary integrals.Here    K k and    S k are the Fourier transforms of the functions If the medium is isotropic, the functions    K k and take the forms: where k is the vector parameter of the Fourier transform, is the scalar product of the vectors k and x, k x  0  and 0  are the Lame constants of the host medium, are the elements of the tensor basis introduced in [1] for presentation of fourth-rank tensors: , , Because the strain   x  and stress

 
x  tensors under the integrals in (1) and ( 4) are multiplied by the function , the elastic fields inside the inclusions (in the region V) are in fact the principal unknowns of the problem.The fields in the medium are reconstructed from the same Equations ( 1) and (4) if the fields inside the inclusions are known.Another important fact that follows from the structure of (1) and ( 4) is that for the numerical solution of this equation, any appropriate region V that includes V may be considered.In particular, by the solution of the problem for an inclusion V of arbitrary shape one can consider a cuboid region V that contains V (see Figure 1).
Unique solutions of ( 1) and ( 4) exist if the tensor of the elastic constants and the inverse tensor do not degenerate inside V [9].

Numerical Solution of the Integral
Equations ( 1) and (4)

Discretization of Equation (1) for Strains by the Gaussian Approximating Functions
We consider Equation (1) for the elastic strains in a me-dium with heterogeneities and following to [6] find its solution in the form of the series: are the nodes of a regular grid that covers a cuboid V that contains the region V occupied by the inclusions, h is the distances between the neighbor nodes, M is the total number of the nodes in V, where , and the tensor Here, the definition (6) of the operator K and Equation (8) for Note that in (14), After introducing a spherical coordinate system in the k-space and integrating firstly over the unite sphere, and then, over the radius k , the integral in the right hand side of (15) is calculated explicitly, and the tensor takes the form: .
Here are the elements of the tensor basis (12).The three scalar functions  in ( 17) and ( 18) have the forms: Here is the probability integral: The system of linear algebraic equations for the unknowns   s  in (13) follows from ( 14) if the latter is satisfied at all the nodes (the Collocation Method).Note that if the nodes of the approximation (13) compose a cubic grid, the coefficients of the approximation   s  coincide with the values of the strain field   (see [6]).As a result, we obtain the system of linear algebraic equations for the coefficients s  in the form: This system may be written in the canonical form as follows: where I is the unit matrix of the dimensions 6 6 M M  , and the vectors of the unknowns X and the right-hand side F are , , , In these equations, , 1,2, , m k M   ; summation from 1 to 3 with respect to repeating indices i, j is implied.These equations follow from ( 23) and ( 27).As it is seen from ( 19)-( 22), the elements of the matrix B in (29)-(34) have simple analytical forms and are calculated fast.

Discretization of Equation (4) for Stresses by the Gaussian Approximating Functions
Let us consider (4) and find its solution in the form similar to (13): Substituting this approximation in (4) and using the Collocation Methods we obtain the following linear algebraic system for the coefficients The last integral is calculated similar to the integral in (15), and the tensor Γ(x) takes the following form: Here the functions Φ 0 , Φ 1 , Φ 2 and 0  , 1  , 2  are defined in ( 18) and ( 19)-( 22).Similar to (23), the discretized Equation (36) may be presented in the canonical form (25).In this case, the vector of unknowns X is composed from the values of the stress tensor ij  at the nodes, the vector of the right

Numerical Solution of the System (25)
It follows from ( 17)-( 22) and ( 29)-(34) that B in (25) is a non-sparse matrix which dimensions may be very large if high accuracy of the solution is required.For the solution of linear algebraic systems with such matrices, only iterative methods are efficient.For instance, if the Minimal Residue Method (see, e.g., [10]) is used, the n-th iteration of the solution of ( 25) is calculated as follows: with the initial values X (0) , Y (0) of the vectors X and Thus, the vector  

Y
is to be multiplied by the matrix B at every step of the iteration process.For nonsparse matrices of large dimensions, calculation of such a product is an expensive computational operation.If, however, a regular grid of approximating nodes is used, the volume of calculations is reduced substantially.Let us consider the product BY in detail.For the matrix B that corresponds to (28)-(34), this product is a combination of the following sums: where the tensor   ijkl x ∏ is defined in (17).For a regular node grid with a step h, the coordinates of every node   r x can be presented in the form: , , Here , , 1,2, m n p   are integers, 1 , 2 , 3 are minimal values of the node coordinates in the cuboid V which sides are parallel to the axes 1

L L L
x , 2 x , 3 x .Thus, the position of every node may be defined by 3 integers (m, n, p).Connection between the one and three-index numerations of the nodes may be introduced by the equation: .
Here N 1 , N 2 , N 3 is the number of the nodes along the corresponding sides of V, In the threeindex numeration, the sum (45) is presented as follows: , , , , , It is seen from the last equation that the object has the Toeplitz structure: it depends on the differences of the indices: , and .As a result, the Fourier Transform Technique can be used for the calculation of the triple sums in (49), and therefore, of the matrix-vector products.Application of the FFT algorithms for calculation of these sums essentially accelerates the iterative process (42)-(44).In addition, one has to keep in the computer memory not all the matrices ijkl of the dimensions but only one column and one row of every such a matrix: the object that has the dimensions and is calculated once for all the iteration process.The details of the FFT algorithm for the calculation of the matrix-vector products are described in [7,11].

An Isolated Heterogeneous Inclusion in a Homogeneous Elastic Medium
Let us apply the method to the calculation of elastic fields inside a spherical isotropic inclusion of a radius a with radially varying elastic properties.The elastic fields will be considered in a Cartesian coordinate system ( 1 2 3 , , x x x ) with the origin at the center of the inclusion.The medium is isotropic with the Young module E 0 and the Poisson ratio 0 .First, we consider an inclusion with parabolic change of the Young module along the radius: .The regular grids of approximating nodes with the steps 0.1 h a  (the total number of the nodes is 9261 M  ), and 0.0317 h a  (M = 262 144) were considered.The graphs in Figure 4 describe the distribution of the shear stress 12  along the x 1 and x 3 axes if the medium is subjected to the external shear stress tensor The bold lines in Figures 2-4 are exact distributions of the components of the elastic stress field obtained by the method presented in [8].It is seen that for the grid step 0.0317 h a  , the numerical solutions are practically coincide with the exact ones.The influence of the parameter H in (35) on the numerical results is not essential if 0.7 2 H   , and 1 H  is taken in the calculations.The same problem was solved with the help of Equation (1) for strains.The medium was subjected to the    of the strain tensor inside the inclusion are in Figures 5-6. Figure 7 shows the distribution of the shear strain 12  by application of the external shear strain ten- The solid lines in these figures are exact distributions of the corresponding components of the strain tensor.It is seen that similar to the equation for stresses, the numerical solution coincides practically with the exact one for the step of the node grid 0.0317 h a  . But the number of iterations in the process (42)-( 44) turns out to be almost two times more than by the use of Equation ( 4) for stresses.This fact reflects general situation: the iteration process based on Equation (4) for stresses converges faster than the same process based on Equation (1) for strains if the inclusion is more rigid than the medium.In the opposite case, when the inclusion is softer than the medium, the iteration process based on the equation for strains converges faster than the same process based on the equation for stresses.
Note that if elastic moduli are not continuous on surfaces inside the inclusion or on the inclusion boundary, some components of the stress and strain tensors have jumps on these surfaces.
In this case, for accurate description of the elastic fields, the grid of the approximating nodes should be sufficiently fine.In the next example, we consider elastic fields in a spherical inclusion of the radius a that consists of a central kernel in the region  (M = 8 120 601) were used for the construction of the graphs in Figures 8-10.The bold lines in these figures are exact solutions obtained by the method presented in [8].It is seen from these figures that the numerical solutions tend to the exact one when the step h decreases.
The stress distributions in a layered spherical inclusion which Young moduli are more than the Young module of the medium are presented in Figures 11-13.In this case, and in Figure 13 to the external shear stress field For the numerical solutions, Equation ( 4) for stresses was used with the node grid steps 0.1 h a  , 0.02 h a  , and 0.01 h a  .Figures 2-13 demonstrate convergence of the numerical solution to the exact distribution of the elastic fields inside the inclusions when step h decreases.The iterative scheme (42)-( 44) converges for any finite values of the elastic constants of the inclusion and the matrix.The number of the iterations grows together with the contrast in the elastic properties of the matrix and the inclusion as well as with the number M of the approximating nodes.Note that for not very large contrast in the elastic properties of the medium and the inclusion, the well-known Conjugate Gradient Method (see, e.g., [10]) applied to the solution of (23) or (36) converges faster than the Minimal Residue Method (42)-(44).But for large contrasts in the properties, the Conjugate Gradient Method may di-verge, meanwhile the Minimal Residue Method keeps converging.

Let a homogeneous medium contain three isolated spherical inclusions of equal radii
We assume that the centers of the inclusions are at the vertices of an equilateral triangle with the coordinates: The numerical solution is constructed in the cuboid V that contains all these inclusions ( 1 : 1 ) (Figure 14).The Young moduli of the medium and the inclusions are E 0 and E, and 0 0.001 E E  .x  that go through the centers of the inclusions are in Figures 16-18, correspondingly.For the solution, Equation (4) for stresses was used.

The Problem of Thermo-Elasticity
We consider a homogeneous elastic medium with several Here   If the medium is subjected to an external stress field and a temperature field ,2   T x , the stress tensor in the medium satisfies the following system of differential equations: Here e  is the tensor of elastic deformation that defines the stress tensor σ according to Hooke's law, is the temperature deformation of the medium, and   0 0, , 0, The system (56) may be interpreted as a system of differential equations for internal stresses in a homogeneous medium with the constant stiffness tensor in the presence of dislocation moments of the density 0 C   m x .The solution of (56) may be presented in the integral form (see [8,12]):   stresses in a medium with heterogeneities.
Let be a constant temperature field, and for , the medium be free of the temperature stresses.In this case, the value of the last integral in (59) depends on the conditions at infinity.If deformation is not restricted at infinity, this integral is equal to zero (see [8]): Meanwhile, for the condition that the total deformation of the medium is equal to zero, we have The right-hand side in this equation is the stress field in an infinite homogeneous medium with the properties ,  This system may be presented in the canonical form (25) which left-hand side coincides with the discretised form of Equation ( 4) for stresses, and the components of the vector F in its right-hand side are defined as follows: Because the objects in (64) for have the Teoplitz properties, the FFT technique may be used for the calculation of the vector F.     Note that the integral equation for the stress field in an elasto-plastic heterogeneous medium follows from (57) and (58) if the temperature deformation T is changed for the plastic deformation p  .Because the plastic deformation is a non-linear function of stresses, the final equation will be non-linear with respect to the stress tensor.The conventional procedure of linearization of elasto-plastic problems is well known.The process of loading of the medium is divided into small intervals, and the problem is considered as linear in every such an interval.The same algorithm may be applied to the solution of the in-tegral equation for the stresses in elasto-platic heterogeneous media.

Conclusions
An efficient numerical method for the solution of 3Dproblems of elasticity and thermo-elasticity for a medium with isolated heterogeneous inclusions is proposed.The method may be applied to the solution of other problems of mathematical physics that can be reduced to volume integral equations.In particular, in [13,14], the method was used for the solution of the problems of electromagnetic wave scattering on perfectly conducting screens and 3D-dielectric bodies.
The main difficulty in the numerical solution of 3Dintegral equations of mathematical physics is a large number of approximation functions (approximating nodes) that should be taken in order to achieve acceptable accuracy.In the present method, this difficulty is overcome by the following means.First, the Gaussian approximating functions allow us to obtain the elements of the matrix of the discretized problems in closed analytical forms and thus, to calculate them fast.Second, the use of regular node grids provides Toeplitz properties to the matrix of the discretized problem.As a result, the number of independent elements of this matrix is essentially reduced, and in addition, the FFT algorithms may be applied to calculating matrix-vector products with such matrices.
For some integral equations of mathematical physics or other types of approximating functions, the elements of the matrix of the discretized problem are not calculated explicitly but expressed via a finite number of standard one-dimensional integrals that can be tabulated (see [6,13]).So, in these cases, the matrices of the discretized problems are also calculated fast.
For regular node grids and any type of identical approximating functions centered at the nodes, the matrices of the discretized problems obtain the Toeplitz properties.Thus, in these cases, the possibility of exploiting FFT algorithms for the calculation of the matrix-vector products always exists.
There are various ways of improvement of the method.As it is shown in [6], the Gaussian functions multiplied with special polynomials increase the precision of the approximation.If such functions are used in the framework of the method, all the algorithms presented in this work will be the same, and only the form of the functions 0  , 1  , and 2  in ( 19)-(22) will change.
Note that in the framework of the method, many parts of the problem may be calculated in parallel.For instance, the elements of the block-matrices in (29)-(34) and products of these matrices with vectors may be carried out independently, and at the same time.This is another advantage of the method and a source of accelerating of the computation process.
into the integral in (1) leads to the following equation: 27) Here T  is the transposition operation.The matrix B in (25) has the dimensions 6 6 M M  and consists of 36 sub-matrices pq b of the dimensions M hand side F consists of the components of the external stress field 0 ij  at the nodes similar to (26), (27).The matrices pq in (28) are defined in (29)-(34), where b of the x 1 -axis, 0  is a scalar.In this case, the Young module is a continuous function together with the components right-hand sides in these figures, the functions and 22 2 on their left-hand sides.The graphs in Figures 2-3 correspond to the numerical solution of Equation (4) for stresses inside the cube V:

Figure 2 .
Figure 2. The dependence of the component σ 11 of the stress field in a spherical inclusion with parabolic distribution of the elastic properties along the radius.The medium is subjected to a uniaxial stress field in the direction of the x 1 -axis.

Figure 3 .
Figure 3.The same as in Figure 2 for the component 22  .

Figure 4 .
Figure 4.The same as in Figure 2 for the component 12  , but the medium is subjected to a constant shear stress 0 12  .

Figure 5 .
Figure 5.The dependence of the component ε 11 of the strain field in a spherical inclusion with parabolic distribution of the elastic properties along the radius.The medium is subjected to a uniaxial strain field in the direction of the x 1 -axis.

.
The graphs in Figures 11-12 correspond to the uniaxial external stress field

Figure 8 .
Figure 8.The dependence of the component ε 11 of the elastic strain field inside a spherical inclusion with step-wise change of the elastic properties along the radius.The medium is subjected to a uniaxial strain field in the direction of the x 1 -axis.

Figure 16 .
Figure 16.The distribution of the component σ 33 of the stress field along the lines x 2 = 0.6, x 3 = 0 (left-hand side) and x 2 = -0.3,x 3 = 0 (right-hand side) that pass through the centers of the inclusions; the medium is subjected to a uniaxial stress field in the direction of the x 3 -axis.

Figure 17 .
Figure 17.The same as in Figure 16 for the stress component σ 11 .

V x is the
characteristic function of the region occupied by the inclusions, the tensors and 0 the operator of incompatibility ( ijk  is the Levi-Civita tensor).If we introduce an auxiliary tensor ' may be rewritten in the following form:


is an external stress field applied to the medium.The kernel   S x of the integral operator in this equation coincides with the generalized function   S x in (5), integration in (58) is spread over the entire medium.Substitution of the tensors and and (53) into (58) yields 59) If 0 T  , this equation coincides with Equation (4) for . 0 .

0 C 0 
subjected to a constant temperature field T by the condition that the total deformation of the medium is zero.In the absence of external stresses 0 ij  and restrictions at infinity, (59) takes the form similar to (35) and application of the Collocation Method lead to the following discretized form of this equation:


Distributions of temperature stresses in the medium with three spherical inclusions considered in Subsection 4.2 are presented in Figures19-24.The Young's module E 0 , Poisson ratio ν 0 , and the coefficient of thermo expansion

Figure 19 .
Figure 19.The distribution of the component σ 11 of the thermo-stress tensor along the axes that go through of the centers of the inclusions shown in Figure 14.

Figure 20 .
Figure 20.The distribution of the component σ 11 of the thermo-stress tensor in the plane x 3 = 0 for the medium with three inclusions shown in Figure 14.

Figure 21 .
Figure 21.The same as in Figure 19 for the thermo-stress component σ 22 .

Figure 22 .
Figure 22.The same as in Figure 20 for the thermo-stress component σ 22 .