An analytical solution for multiple inclusions subject to a general applied thermal field

Abstract In this study, we present a general numerical solution to the thermoelastic problem of interacting elliptical inhomogeneities in an infinite matrix material. Our solution provides a framework to study microstructural responses to the external mechanical, thermal, electric, and magnetic fields simultaneously. The proposed semi-analytical method, which enables us to reduce the computational cost and calculate the interactions between many inclusions, is based on conformal mapping, series expansion of the corresponding complex potentials and boundary conditions. The temperature, stress, strain, and displacement functions are derived explicitly in the matrix material and inclusions. The thermal loadings tested in this study include temperature variations and a uniform heat flow. The model is validated against multiple analytical calculations and benchmark numerical problems. The method is ready to apply for simulations with many inclusions.


Introduction
Thermoelasticity of cracks, holes, and second-phase inclusions has been investigated extensively due to the importance of local stress concentrations and their roles in material performance.In this study, we develop a general semi-analytical method that can calculate the effects of thermomechanical loadings on microstructures and can be expanded to add the impact of electric and magnetic fields.
Study of two or multiple inclusions (inhomogeneities) is more challenging.Thermal stresses induced by two inclusions have been calculated by different scientists [26][27][28][29][30][31][32] for almost seven decades.Muramatsu resolved the problem of three holes in a row and an infinite row of holes [33,34].Few studies focus on addressing the thermoelasticity of the multiple inclusions.Dai and Sun derived a series solution to the problem of a plate containing many elliptical inclusions subjected to a steady-state temperature variation field [35].Dai et.al solved the stress and displacement fields of multiple elliptical inclusions with the predefined temperature field using the Faber series expansions for inside the inclusions and Taylor series expansions for the matrix material which requires two truncated numbers.The boundary conditions were implemented with the aid of Fourier series expansions of all terms in the equations to transfer them into the same basis functions.Therefore, this method requires an extra level of approximation as well as the computational cost.Zhang addressed the same issue using the "equivalent inclusion" method and "distributed dislocation" method [36,37].The "equivalent inclusion" method cannot be applied for close-range interactions as each inclusion experiences a nonuniform field.Therefore, this study is only valid for a system with far-distanced inclusions where close interactions can be neglected.
Materials contain multiple inclusions, and understanding their reactions to the external mechanical and thermal fields is required for material design purposes.Increasing the number of impurities elevates the computational costs of modeling, and consequently, predicting the interactions using semi-analytical methods or hybrid modeling becomes more appealing [38][39][40][41][42][43][44].Analytical methods are integrated with the Finite Element Method (FEM) to improve computational costs and scale up simulations.For instance, a family of hybrid-Trefftz elements containing an elliptic hole has been formulated using the complex variable technique [38][39][40][41].Moreover, semi-analytical methods are sufficiently fast and accurate and can be integrated by Machine learning techniques to generate data sets for accelerating efficient microstructural modeling [45][46][47][48].Therefore, calculating interactions between inclusions using CVM not only reduces the computational costs while controlling the accuracy of the calculations, but also it is capable of being incorporated into other analyses to raise the speed of material modeling.
All the reviewed studies, except Zhang's work [36,37], are based on the CVM and selection of various basis functions.The critical point that comes out of these studies is: "what are the basis functions that we can use to be able to implement them to find a general solution for the multiple numbers of inclusions experiencing external fields?"In this study, we aim to build a framework that can fulfill this demand for mechanical and thermal loadings while it can be broadened to include the effects of electric and magnetic fields.The various thermal loadings investigated in this research include temperature variations and a uniform heat flow.The order of this work is as follows.Section 2, the problem is converted to conformal mapping and series expansion of the corresponding complex potentials and boundary conditions.The temperature, stress, strain, and displacement functions are derived explicitly in the matrix material and inclusions.Once the mathematical framework is defined, the algorithm of thermoelasticity is elaborated in Section 4, and the results of numerical test problems are summarized afterwards.

Thermal fields
In quasi-static thermoelasticity problems temperature T and the local heat flux density Q satisfies the following heat conduction equations: where k and S are the thermal conductivity and the heat source respectively.z ¼ x 1 þ ix 2 is the complex variable representing the position vector.The temperature field and the heat flux Q ¼ q 1 þ iq 2 can be expressed by an analytical complex potential function XðzÞ as: Primes denote differentiation with respect to z and a superimposed bar denotes the complex conjugate.Function X depends on thermal loading and geometrical shape of the interfaces (boundaries).Here, we apply the series technique and boundary conditions to determine the function X: Muskhelishvili and Hasebe [7,18] showed the normal flux continuity condition is equivalent to continuity of Im½kXðzÞ at the interface.Thus the continuity equations become: where z 0 is a point on the interface.

Stress and displacement fields
Muskhelishvili [7] derived the general solution for the problems of plane elasticity in terms of the two complex potentials functions uðzÞ and wðzÞ: Following Kushch [49], we modify fields using an equation wðzÞ ¼ wðzÞ þ zu 0 ðzÞ: The stress tensor components r ij , displacement vector u ¼ ðu 1 , u 2 Þ, considering the effect of thermal expansion, and traction forces can be expressed as follows: where l is the shear modulus, is the Poisson's ratio, and j is Kolosov's constant equal to j ¼ 3 À 4 and b l ¼ 2a l ð1 þ Þ for plane strain, j ¼ ð3 À Þ=ð1 þ Þ and b l ¼ 2a l for plane stress; a l is the linear coefficient of thermal expansion.ÀY þ iX are traction forces and b l Ð XðzÞdz is the thermal displacement.

Single inclusion
Consider an infinite elastic matrix containing an elliptical elastic inclusion with different thermomechanical properties.Without loss of generality, we assume that inclusion centered in the point z ¼ 0 and Coordinate axes are directed along the major (l 1 ) and minor (l 2 ) axes of the ellipse.The aspect ratio and the inter-foci distance of the ellipse are e ¼ l 2 =l 1 < 1 and respectively.Here, we use an elliptic coordinate frame with radial f and angular g coordinates.Elliptic coordinates n ¼ f þ ig are translated into Cartesian coordinates by z ¼ dcoshðnÞ: f ¼ f 0 coincides with the boundary of elliptic inclusion and equals to The ambient temperature field T far ¼ Re½X far ðzÞ is the analytical function.The elliptical inclusion produces the local, vanishing at infinity thermal stress distribution due to the thermal incompatibility of the two phases.Both phases are assumed to be thermally and mechanically homogeneous and isotropic.The two materials are assumed to be perfectly bonded along the interface and have perfect thermal contact which implies temperature and normal heat flux continuity through the interfaces.We solve the problem by deriving the temperature field and calculating the corresponding displacement fields using Eq. ( 6).Indexes "0" and "1" represent the matrix material and inclusion related quantities.Following the Kushch [50,51] solution, we expand the temperature function T 1 in the form of complex function X 1 as Knowing T 1 is finite inside the inclusion [7] results in G n ¼ G Àn and Likewise, by assuming finite far-field temperature T far in a vicinity of inhomogeneity, the complex field becomes: The total temperature in the matrix is: where X s is the disturbance temperature field induced by the inclusion and it is vanishing at infinity.G n , F n and S n are the unknown complex coefficients and v ¼ exp ðnÞ ¼ z=d þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðz=dÞ 2 À 1 q : Unknown coefficients are derived by fulfilling the boundary conditions given in Eq. ( 3). S 0 ¼ 0 is required to have the vanishing disturbance fields at infinity.Consequently, are only function of angular coordinate.

Resolving temperature equations
i) The temperature and heat flux continuity equations are simplified using the orthogonality of Fourier harmonics and become: Where k1 ¼ k 1 =k 0 : Eliminating G n from Eqs. ( 12) and ( 13) results in: Thus, the problem of one inclusion subjected to the external temperature field can be calculated.

Single inclusion: Thermal displacement
Thermal fields displacement inside the inclusion is: For fields in the matrix we have: where S n ¼ 0 for n 0: The singular term in displacement field, named thermal dislocation [3], was eliminated by a specific selection of the complex potentials [3,6,32].Then for the multiinclusion system, the new complex functions with singular terms have to be expanded in another inclusion coordinate system.Here, we expand the displacement singular term and integrate it with our previous method [52] to derive the Thermal and mechanical external field effects in one system.It also provides the framework to expand the method to integrate it with external electroand magneto-static effects.According to Eq. ( 6), ln v is the only term which is not a factor of v 6n : Fourier series expansion will simplify it on the boundary of the ellipse where v ¼ v 0 e ig ¼ v 0 t: The Fourier series expansions of ln can be given by means of the boundary value of Faber series.
Using the potentials and the boundary conditions, the expansion coefficients can be derived.First, displacement continuity across the interface results to: For more detail on derivation of the equation please see Eq. ( 10) in [53] and Eq. ( 4).An infinite sets of linear equations obtained by equating the coefficients of equal Fourier harmonics t n ¼ e ing : Displacement and traction continuity yield to: Where A n and B n are coefficients of the disturbance field in the matrix and C n and D n are fields coefficients in the inclusion region, a n and b n are fields coefficients in the matrix resulting from the far-field load and finally, S n , F n and G n are the coefficient of thermal fields in the matrix material and inclusion consequently.The two Eqs.(19) and (20) together with two equations from the continuity of traction at the boundary (21) provide a linear set of equations to calculate the expansion coefficients.From here, solving the problem is equivalent to deriving the expansion coefficients of Equations.Once the expansion coefficients are obtained, the potential functions and consequently displacement and stress fields can be calculated.
Another set of linear equations can be derived using the below scalar potentials in such a way that cancels dislocation terms and avoid Fourier series expansions of the logarithmic term of the thermal displacement.Thus the S 1 ðÀ1Þ n v n 0 term will be removed from the linear sets of equations and A 0 ¼ À b 0l S 1 d 2ðj 0 þ1Þ :

Single inclusion: Verification
Here we test and verify the method with couple of benchmark problems.

i. Uniform temperature change
Temperature variation causes displacement fields in the system and can be calculated as: where X far ¼ DT is a constant.We test the method on a circular inclusion located at the origin of the coordinate system.Indexes "1" and "0" denote inclusion and matrix material.Figure 1 illustrates generated stresses in the matrix material when the system is subjected to a uniform temperature variation.Radial and hoop stresses in the matrix material differ less than 4.4% compared with the analytical analysis [54,55] where r rr ¼ pR 2 =r 2 and The method is verified with a crack (a collapsed elliptical cavity with aspect ratio e ¼ 0:01=0:5) under uniform heat flux CzÞ and F n ¼ d n1 d C=2; d ij is the Kronecker delta.The stress fields (Figure 2) are compared with an adiabatic crack under uniform heat flux in [56] and [57].Finite-sized effects are not included in this study, and discrepancies between the results in the vicinity of the sample's edges are expected to be seen.The thermo-mechanical properties of the material used in this work are summarized in Table 1.
We assume that cracks and holes are traction-free surfaces without deformation.This condition imposes a singular term in stress fields of matrix materials.The predicted SIF becomes:   Ref [58] Ref [56,57] Poisson ratio which is the same formula as the analytical SIF derived by Koizumi [6].In Table 2 the normalized SIFs are compared with those of other simulations of a rectangular plate with a slant crack [56,58].Results of simulations show the same trend for SIF variation with crack orientation however, the discrepancy in the magnitude is associated with the different applied boundary conditions.

Finite array of inclusion
When there are N elliptical inclusions in the material, we define N parallel local coordinate systems (x 1p , x 2p ), each with its origin at the centroid of the ellipse; another parallel coordinate system is designated as the global coordinate system (x 1 , x 2 ).Defining the coordinates of a point expressed in the global coordinate system ðx 1 , x 2 Þ relates to the coordinates of the point expressed in the local system ðx 1p , x 2p Þ by z ¼ z p þ Z p , where Z p is the coordinate of the centroid of the pth inclusion in the global coordinate system.
To describe the differences in the relative orientations of the ellipses another set of N local rotated coordinate systems (y 1p , y 2p ) is required, aligned along the major and minor axes of each ellipse.If the coordinate axes (y 1p , y 2p ) are related by an anti-clockwise rotation of h p relative to the local coordinate axes (x 1p , x 2p ) then y p ¼ y 1p þ iy 2p ¼ e Àih p z p : The conformal transformation y p ¼ D p coshn p , where 2D p is the distance between the foci of the pth ellipse, introduces the local elliptic coordinates n p ¼ f p þ ig p : Then z p ¼ D p e ih p coshn p ¼ d p coshn p (Figure 3).

Temperature fields
Using the linear superposition principle, the general solution T 0 ðzÞ for the multiply connected matrix in the global, and local z q coordinate system, is written as: where X far , and X dis are complex temperature fields caused by an external temperature field, and a total disturbance field induced by all inclusions respectively.By applying Eq. (A.1) to X s we have: where S q n 0forn 0 and h pq ¼ h q À h p , and This is the regular part of the solution, so The same formula is derived in Appendix A and Appendix B using the displacement continuity at the ellipses interfaces.By substituting Eq. ( 27) into Eq.( 14) for qth inclusion's interface, we come to the following compact formula:

Stress and displacement fields
The stress and displacement fields in and around the inclusions are induced by the far-field stresses and temperature.In the previous section, we found the temperature fields of a finite number of inclusions.
Using the potentials stress and displacement fields are predictable by generalizing Eqs. ( 19)-( 21) to finite number of inclusions as below: where index q refers to the qth inclusion.From this point onwards the problem reduces to solving a set of linear equations where the unknown coefficients can be obtained.Once all coefficients are derived, the other scalar potentials and consequently the stress and displacement fields can be calculated.Another set of linear equations can be derived using the below scalar potentials to avoid Fourier series expansions of the logarithmic term of the thermal displacement A numerical solution with an arbitrary accuracy can be obtained using the truncation method retaining a finite number of n max harmonics and hence a finite number of equations.The convergence of the solutions can be addressed by first the convergence of the series solutions for temperature calculations and second the convergence of stress fields solutions.
The convergence rate of Eq. ( 7) is substantially affected by the problem parameters such as the inclusions aspect-ratio, and the minimum distance between inclusions.Convergence accelerates for e ! 1, and a higher number of harmonics requires for collapsed inclusions and cracks.According to Kushch's study [50], for a finite number of circular inclusions with a minimum distance of r min 0:1R or volume content of 0.7, the truncated number of n max ¼ 20 provides the results with at least 5-digit accuracy.The form of Eqs. ( 30)- (33) are remarkably simple if it is not simplest, and they are based on the Multipole method framework with modified far-field loading [53].Using the MM framework is rational as the validity and the convergence of the series solution are already tested and verified for non-overlapping and neighboring inclusions or holes approaching each other.For two close inclusions with r min ¼ 0:1R, as many as 25 harmonics are required to get a practically convergent solution.Therefore, the value of n max ¼ 25 is adapted for the subsequent calculations.

Numerical results
The problem illustrated in Figure 4 with the thermo-elastic constants of the matrix and inclusions listed in Table 3 are selected to assess and compare the method.The system is subjected to a uniform temperature change DT ¼ 20: Figure 4 shows the stress fields obtained by the presented method outside the inclusions are in good agreement with the results of the FEA and Zhang method [37].Although close interactions were not considered in [37] work, it can predict far interactions between inclusions with fixed boundary conditions accurately.Our method calculated the close interactions which cause the non-uniform field inside the centered inclusion.However, because inclusions are positioned far from each other, stress fields inside other inclusions are nearly uniform and derived external fields become similar to the Zhang predictions.We excluded the fields around the singular points and branchcut in the inclusions.
Figure 5 shows the stress fields of the same problem with the traction-free boundary conditions.Our simulation indicates that the stress fields are non-uniform inside the inclusions, and therefore, the "equivalent inclusion method" cannot be used.However, our method can predict the accurate stress fields which agree well with FEM simulations.
We verify the method with another mathematical model [32] and FEM simulations for two circular holes with radius R 1 =R 2 ¼ 1 separated with d ¼ 0:4R 1 and 0:1R 1 : The heat conduction and thermal stresses are calculated for material with thermal conductivity k ¼ 10 W=mK,   6 and 7. Distributions of temperature derived by the three methods agree well.Comparing the results of the two theoretical methods shows that, the magnitude of the stress fields is nearly identical, however, stress distributions seem dissimilar.The differences can be associated with the plotting method or the selections of the singularity term in the complex potential functions as Song selected the potentials in a way that they cancel the thermal dislocation terms [32].
External thermal and mechanical loading generate contradictory stress field distributions.For example, in the case of two holes, the external mechanical field imposes the stress amplifications between holes [52,59] while Figure 6 indicates stress shielding caused by the thermal field.The ultimate stress fields are determined from the entire interactions and this must be incorporated in designing materials to optimize their functionality such as brittleness, strength or toughness depending on the desired operation.
This study aims to further elevate our knowledge in understanding the localized thermal stresses around the holes and inclusion and their effects on crack initiation and propagation.Previously, we developed a method to study fracture in the vicinity of the circular inclusions [52,60].In our future work, we will integrate the mathematical framework derived from this work with our crack propagation method to study systems containing cracks, holes and inclusions.

Conclusion and further work
In this paper, we developed a numerical solution to simulate the thermal interactions between a finite array of arbitrarily oriented cracks and inclusions with varying aspect ratios and material properties subjected to external stress and temperature fields.In this method, we expanded the complex potentials of stress and temperature fields into a truncated Faber series in the local coordinate system of the inclusions using the same basis functions.Firstly, we derive expansion coefficients of the temperature field by applying thermal boundary conditions.Secondly, we apply the stress and displacement boundary conditions in assistance with the Fourier series expansions of the logarithmic term of the thermal displacement to calculate the expansion coefficients of the complex potentials using the system of 4n m axn linear equations.Here, n is the number of inclusions, and n m ax is a truncated number of Faber series expansions.The method provides an accurate semi-analytical tool to derive stress and displacement fields, as well as stress intensity factors.A set of benchmark problems with known analytical and computational results was selected to assess and verify the method.This study provides a basis for further investigation of microstructure response to the external fields, including electric and magnetic fields.The method has applications to study polycrystalline and composite materials, particularly in heat conduction and thermal stress of the fiber-reinforced composite and porous materials.Forthcoming studies using hypergeometric functions to expand v Àn p ¼ ½z p =d p 6 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðz p =d p Þ 2 À 1 q Àn , then Taylor expansion of the z p in the vicinity of the z q , and finally re-expansion of z n q to hypergeometric functions.The obtained series of expansion coefficients is and, Converting a "singular" elliptical harmonic v Àn p to a series of hypergeometric functions and applying an additional Taylor expansion provides an additional geometric restriction on the convergence of g pq nm series.Therefore, for the nearest neighbors another computationally expensive formula [66] shall be used.This integral form of coefficient converged for any two non-overlapping ellipses.
Àncosh À1 dq cos g q þzpq dp cos ðmg q Þdg q (A.5) Another useful formula that is used in this context is a differentiation of the re-expansion formula with respect to the z p : @v Àn p @z p ¼ X k l pq nk ðv q Þ Àk (A.6) The above mentioned re-expansion method, is used to evaluate fields in a system containing a finite array of aligned inclusions.u s p ðy p Þ ¼ ju s p ðy p Þ À ðy p À y p Þu s p ðy p Þ 0 À w s p ðy p Þ þ b lp e Àih p Ð X p s ðz p Þdz p u r pq ðy q Þ ¼ ju r pq ðy q Þ À ðy q À y q Þu r pq ðy q Þ 0 À w s q ðy q Þ þ b lp e Àih q Ð X pq r ðz q Þdz q : (B. 3) The scalar potentials that solve the problem in the p and q coordinate systems are: u r pq ðy q Þ ¼ P m6 a mpq v Àm q w r pq ðy q Þ ¼ P m6 b mpq À ma mpq sinhf 0q sinhn q v q v 0q À v 0q v q " # v Àm q (B.5) By equating similar terms in u s p ðz p Þ and u r pq ðz q Þ and transforming the basic functions, the expansion coefficients a mpq and b mpq in the q coordinate system can be derived in terms of A np and B np : Also, b lp Ð X p s ðz p Þdz p b lp Ð X pq r ðz q Þdz q X p s ðz p Þ X pq r ðz q Þ (B.6) Which is expected as T p ðz p Þ ¼ T pq ðz q Þ:

Figure 1 .
Figure 1.Depiction of normalized radial stresses r rr =pR 2 in a) 1D (black and red lines are derived from this study and analytical analysis), and b) 2D.

Figure 2 .
Figure 2. Distribution of stress fields r yy , r xy , and r xx around the a) crack Ref [56] and b) crack like cavity.

Figure 4 .
Figure 4. a) A square plane containing 5 inclusions subjected to DT: Von Mises stress obtained by b) presented method, c) FEM,and d) equivalent inclusion method methods[37].

Figure 5 .
Figure 5.A square plane containing 5 inclusions with traction free boundary conditions subjected to DT ¼ 20K temperature variation.Von Mises stress obtained by a) COMSOL FEM b) presented method.

Figure 6 .
Figure 6.Distributions of temperature and von Mises stress around two circular holes a) [32], b) our method, and c) COMSOL FEM.

Table 1 .
Material Properties; simulations are independent of thermal conductivity.

Table 2 .
SIFs of a slant crack calculated in this study compared with other simulations of a rectangular plate of width 2 W, length 4 W and a central slant crack of length 2a ¼ 0:6W: Figure 3. Local and global coordinate systems.

Table 3 .
Thermo-elastic constants of the matrix and inclusions.Young's Module E ¼ 200 GPa, Poisson ratio ¼ 0:3, and expansion coefficient a ¼ 10 À5 1=K parameters.The remote heat flux with C 2 ¼ 10 4 K=m is applied in the x 2 direction.Distributions of temperature and von Mises stress are illustrated in Figures By re-expanding the function X p s ðz p Þ in qth coordinate system we have: Daphne Jackson Trust.Data reported in this paper can be obtained upon request by emailing the corresponding authors or d.balint@imperial.ac.uk. the