Boundary Element Analysis of Three-Dimensional Exponentially Graded Isotropic Elastic Solids

Copyright c 2007 Tech Science Press CMES, vol.22, no.2, pp.151-164, 2007 Boundary Element Analysis of Three-Dimensional Exponentially Graded Isotropic Elastic Solids R. Criado1,J.E.Ortiz 1,V.Manti ˇ c1,L.J.Gray 1,2and F. París1 Abstract: A numerical implementation of the Somigliana identity in displacements for the solution of 3D elastic problems in exponentially graded isotropic solids is presented. An expression for the fundamental solution in displacements, Uj, was deduced by Martin et al. (Proc. R. Soc. Lond. A,458, pp. 1931–1947, 2002). This expression was recently corrected and implemented in a Galerkin indirect 3D BEM code by Criado et al. (Int. J. Numer. Meth. Engng., 2008). Starting from this expression of Uj,anew expression for the fundamental solution in tractions Tjhas been deduced in the present work. These quite complex expressions of the integral kernels Ujand Tjhave been implemented in a collocational direct 3D BEM code. The numerical results obtained for 3D problems with known analytic solutions verify that the new expression for Tjis correct. Excellent accuracy is obtained with very coarse boundary element meshes, even for a relatively high grading of elastic properties considered. Keyword: functionally graded materials, boundary element method, three-dimensional elasticity, Somigliana identity, fundamental solution in tractions. 1 Introduction Functionally Graded Materials (FGMs) [Suresh and Mortensen (1998)] represent a new generation of composites, having a continuous variation of apparent material properties obtained through a progressive variation of their microstructural composition. Stress concentrations appearing at 1School of Engineering, University of Seville, Camino de los Descubrimientos s/n, Sevilla, E-41092, Spain 2Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831-6367, USA material discontinuities in various applications (for example, thermal barrier coatings) can be avoided or diminished using FGMs. The first numerical studies of FGMs have been carried out using the Finite Element Method (FEM) [Lee and Erdogan (1998); Santare and Lambros (2000); Anlas et al. (2000); Kim and Paulino (2002a,b); Paulino and Kim (2004)] due to its capability to include, relatively easily, variation of material properties. The Meshless Local Petrov-Galerkin (MLPG) Method [Atluri and Shen (2002); Sládek et al. (2004, 2005, 2007); Ching and Chen (2006)] has also recently shown to be a versatile approach to solve different problems including non-homogeneous media. The Boundary Element Method (BEM) [París and Cañas (1997); Aliabadi (2002)] is another technique for elastic analysis, capable of solving problems with material and geometrical discontinuities, e.g., crack growth and contact, and also very suitable for flaw detection and shape optimization. Nevertheless, an adaption of BEM to nonhomogeneous media is a hard task, as fundamental solutions or Green’s functions (corresponding to concentrated loads or sources) for such media are difficult to obtain. Fundamental solutions for heat transfer problems in non-homogeneous media have been presented in Shaw and Makris (1992); Clements (1998); Gray et al. (2003); Berger et al. (2005) and Kuo and Chen (2005) and implemented in BEM codes by Gray et al. (2003); Clements and Budhi (1999); Núñez et al. (2002); Paulino et al. (2002) and Sutradhar and Paulino (2004). Green’s functions due to surface loads in non-homogeneous half-spaces have been deduced in Han et al. (2006) and Seyrafian et al. (2006). Fundamental solutions for 2D and 3D elastic 152 Copyright c  2007 Tech Science Press CMES, vol.22, no.2, pp.151-164, 2007 problems in exponentially graded isotropic materials have been deduced only recently in Chan et al. (2004) and Martin et al. (2002). These solutions have not as yet been checked computationally, to the knowledge of the present authors, which can be due to the fact that implementing them in a BEM code is far from straightforward. In the present work the displacement fundamental solution Ujl corresponding to a point force in a 3D exponentially graded elastic isotropic media, developed originally in Martin et al. (2002) and corrected by Criado and co-workers (2005, 2008), is employed in the form presented in Criado and co-workers (2005, 2008). Moreover, a new expression of the corresponding traction fundamental solution Tjl is presented herein, and both functions have been implemented in a 3D collocational BEM code. To check the correctness of the kernel function expressions and to prove their suitability to be implemented in a BEM code, and also to check the overall BEM implementation, two 3D problems with known analytic solutions for exponentially graded materials have been analysed by this BEM code. 2 Properties of Elastic Exponentially Graded Isotropic Materials The fourth rank tensor of elastic stiffnesses cijk (i,j,k,=1,2,3)for an exponentially graded material varies according to the following law: cijk(x)=Cijkexp(2 β β β ·x),(1) where xis a point in the material and the vector β β β defines the direction and exponential variation of grading, β = β β β . According to (1), points situated in a plane perpendicular to β β β have the same stiffnesses, Cijkgiving the stiffnesses in the plane including the origin of coordinates. In the case of isotropic materials, the Lamé constants λ and μ satisfy cijk(x)= λ (x) δ ij δ k+ μ (x) δ ik δ j+ δ i δ jk, (2) where δ ij is Kronecker delta, and hence for exponential grading λ (x)= λ 0exp(2 β β β ·x)and μ (x)= μ 0exp(2 β β β ·x).(3) Here λ 0and μ 0are the Lamé constants on the plane that includes the origin of coordinates. It is easy to check, that λ (x)/ μ (x)= λ 0/ μ 0= 2 ν /(1−2 ν ), ν being the (constant) Poisson ratio defined as ν = λ 0/2( λ 0+ μ 0). Note that the assumption of a constant Poisson ratio appears to be reasonable for many real graded materials, and is a commonly used approximation in an analysis of graded materials. 3 Elastic Fundamental Solution in 3D Exponentially Graded Isotropic Materials 3.1 Displacement fundamental solution According to Martin et al. (2002), the displacement fundamental solution can be written as U(x,x)=exp{− β β β ·(x+x)} ·U0(x−x)+Ug(x−x),(4) where Uj(x,x)gives the j-th displacement component at xdue to a unit point force acting in the -direction at point x,andU0is the weakly singular Kelvin fundamental solution associated to a homogenous isotropic material defined by λ 0and μ 0(see París and Cañas (1997); Aliabadi (2002)). The so-called grading term Ug j(x−x)=−1 4 πμ 0r1−e− β r δ j +Aj(x−x)(5) is bounded for r→0 and vanishes for β =0, r= rwhere r=x−x. Let an orthogonal system of coordinates (˜x1,˜x2,˜x3), whose origin is placed at x,be defined by the orthonormal right-handed triad {n,m,ˆ β β β },whereˆ β β β =(ˆ β 1,ˆ β 2,ˆ β 3)= β β β / β ,and nand mare orthonormal vectors in the plane perpendicular to β β β . Let the following spherical coordinate system (r,Θ,Φ)be associated to this Analysis of Isotropic Elastic Solids 153 coordinate system: r·n=rsinΘcosΦ, r·m=rsinΘsinΦ r·ˆ β β β =rcosΘ, (6) where 0 ≤Θ≤ π and 0 ≤Φ≤2 π . According to Criado and co-workers (2005, 2008) the term Ajl is composed of the following five integrals: Ajl =− β 4 π (1− ν ) μ 0 I1 − β 2 π 2(1− ν ) μ 0 (I2−I3+I4−I5,),(7) where I1= 2 ∑ s=0 2 ∑ n=0 π /2 0 R(n) se−|k|ysIn(Kys)sin θ d θ , (8) I2= 2 ∑ s=0 π /2 θ m R(0) ssin θ  π /2 η m sinhΨsd η d θ ,(9) I3= 2 ∑ s=0 π /2 θ m R(2) ssin θ  π /2 η m sinhΨscos2 η ·d η d θ , (10) I4= 2 ∑ s=1 π /2 θ m M(1) ssin θ  π /2 η m coshΨssin η ·d η d θ , (11) I5= 2 ∑ s=1 π /2 θ m ˜ M(1) ssgn(k)sin θ  π /2 η m sinhΨs ·sin η d η d θ , (12) the extensive notation introduced in this equation being now defined. First, In(x)denotes the modified first kind Bessel function of order n, I1(Kys)= 2 π  π /2 0 sinh(Kyssin η )sin η d η ,(13) In(Kys)= 2 π in π /2 0 cosh(Kyssin η )cosn η d η , n=0,2.(14) The integration limits θ mand η m(0≤ θ m, η m≤ π 2)are defined by θ m(Θ)=1 2 π −Θ, |k(r,Θ, θ )|=K(r,Θ, θ )sin η m(Θ, θ ),(15) where k(r,Θ, θ )= β rcos θ cosΘand K(r,Θ, θ )= β rsin θ sinΘ, and the range of θ guarantees that η mis well defined. The argument of the hyperbolic functions is Ψs(r,Θ, θ , η ) =K(r,Θ, θ )ys( θ )(sin η m(Θ, θ )−sin η ),(16) where the functions ysare given by y0=1, y1( θ )=q( θ )+q2( θ )−1, y2( θ )=q( θ )−q2( θ )−1, (17) with q( θ )≥1definedas q( θ )=1+2 ν 1− ν sin2( θ ).(18) The functions R(n) sand M(n) sare given by R(0) s=M(0) s,R(2) s=−M(2) s,s=0,1,2, (19) R(1) s=−M(1) s+ M(1) ssgn(k),s=1,2, (20) M(n) 0=fn(1) 2D(1),M(n) s=fn(ys) (1−y2 s)D(ys), (21) n=0,2ands=1,2, M(1) s=f1(ys) D(ys), M(1) s=˜ f1(ys) D(ys),s=1,2, (22) while the functions fiare defined by f0(x)=1 28 ν x4−(−x2+1)(−2x2q+1) ·(njn+mjm)sin2 θ +8 ν x4sin2 θ +(−x2+1)[−x2−(−2x2q+1)cos2 θ ]ˆ β jˆ β , (23) 154 Copyright c  2007 Tech Science Press CMES, vol.22, no.2, pp.151-164, 2007 f1(x)=x3(4 ν −1)(sjˆ β −ˆ β js)sin θ ,(24) ˜ f1(x)=−1 2(sjˆ β +ˆ β js)(−2x2q+1)sin2 θ ,(25) f2(x)=−1 2[8 ν x4−(−x2+1)(−2x2q+1)] ·nj(ncos2Φ+msin2Φ) +mj(nsin2Φ−mcos2Φ)sin2 θ ,(26) sj(Φ)=njcosΦ+mjsinΦ,(27) and the polynomials D(x)and D(x)by D(x)=x4−2x2q+1andD(x)=−4x3+4xq. (28) Notice that D(x)is not the derivative of D(x). A discussion of the properties of the fundamental solution Ujl and some aspects of the above expression, together with recommendations for its numerical evaluation can be found in Criado and co-workers (2005, 2008). 3.2 Traction fundamental solution The direct boundary integral equation for surface displacement requires the displacement fundamental solution, and the corresponding traction fundamental solution. The starting point in the evaluation of tractions in an exponentially graded material due to a unit point force is the differentiation of the fundamental solution in displacements Ujl. These derivatives are used to determine the corresponding strains, and then employing the constitutive law with the tensor of elastic stiffnesses given in (2-3), the corresponding stresses can be obtained. Differentiation of (4) yields ∂ Uj ∂ xk (x,x)=exp − β β β ·(x+x) · ∂ U0 j ∂ xk (x−x)+ ∂ Ug j ∂ xk (x−x) − β kUj(x,x).(29) Although the derivative of U0 jis strongly singular, this term eventually produces the Kelvin traction kernel for a homogeneous material; the expressions can be found in París and Cañas (1997) and Aliabadi (2002). The derivative of Ug jis weakly singular and can be expressed, in view of (5), as ∂ Ug j ∂ xk (x−x)= − δ j 4 πμ 0e− β r( β r,k) r−(1−e− β r)r,k r2 + ∂ Aj ∂ xk (x−x),(30) where the derivative of Ajl is, according to (7), decomposed into the sum of the derivatives of the integrals Ii ∂ Ajl ∂ xk (x−x)=− β 4 π (1− ν ) μ 0 ∂ I1 ∂ xk − β 2 π 2(1− ν ) μ 0 ∂ I2 ∂ xk − ∂ I3 ∂ xk + ∂ I4 ∂ xk − ∂ I5 ∂ xk. (31) Note that the weakly singular character of ∂ Ug j/ ∂ xkdirectly follows from the boundedness of Ug jand the Gauss divergence theorem. When differentiating Ii(i=2,...,5), involving double integrals with respect to η and θ , it should be taken into account that while their superior limits are constant, their inferior limits are varying with the positions of the field and source points, x and x, as follows: •Inferior limit of the integral in θ : θ m= θ m(x,x), •Inferior limit of the integral in η : η m= η m(x,x, θ ). Thus, derivatives of these double integrals are evaluated by applying the following rule twice: d dx B A(x) f(x,t)dt=B A(x) ∂ f(x,t) ∂ xdt−f(x,A(x))dA dx. (32) By also taking into account that η m( θ = θ m)= π /2andΨs( η = η m)=0, see Criado (2005), the following expressions are obtained after some al- Analysis of Isotropic Elastic Solids 155 gebraic manipulations: ∂ I1 ∂ xk = 2 ∑ s=0 2 ∑ n=0 π /2 0 e−|k|yssin θ  ∂ R(n) s ∂ xk In(Kys) +R(n) s−ys ∂ |k| ∂ xk In(Kys)+ ∂ In(Kys) ∂ xkd θ (33) ∂ I2 ∂ xk = 2 ∑ s=0 π /2 θ m sin θ  ∂ R(0) s ∂ xk π /2 η m sinhΨsd η +R(0) s π /2 η m coshΨs ∂ Ψs ∂ xk d η d θ (34) ∂ I3 ∂ xk = 2 ∑ s=0  π /2 θ m sin θ  ∂ R(2) s ∂ xk π /2 η m sinhΨscos2 η d η +R(2) s π /2 η m coshΨscos2 η∂ Ψs ∂ xk d η d θ (35) ∂ I4 ∂ xk = 2 ∑ s=1 π /2 θ m sin θ  ∂ M(1) s ∂ xk π /2 η m coshΨssin η d η +M(1) s π /2 η m sinhΨssin η∂ Ψs ∂ xk d η − ∂η m ∂ xk sin η md θ (36) ∂ I5 ∂ xk = 2 ∑ s=1 π /2 θ m sin θ · ∂ ˜ M(1) s ∂ xk sgn(k) π /2 η m sinhΨssin η d η +˜ M(1) ssgn(k) π /2 η m coshΨssin η∂ Ψs ∂ xk d η d θ (37) where ∂ R(n) s ∂ xk =0,n=0andn=1 with s=0,(38) =− ∂ M(1) s ∂ xk − ∂  M(1) s ∂ xk sgn(k),(39) n=1 with s=1,2, =− ∂ M(2) s ∂ xk ,n=2,(40) ∂ Ψs ∂ xk = ∂ K ∂ xk ys(sin η m−sin η )+Ky scos η m ∂η m ∂ xk , (41) ∂ In ∂ xk =2 π (−1)n/2 π /2 0 cos(n η )sinh(Ky ssin η ) (42) · ∂ K ∂ xk yssin η d η ,n=0,2, =2 π  π /2 0 sin( η )cosh(Ky ssin η ) ∂ K ∂ xk ys (43) ·sin η d η ,n=1. The derivative of η mis expressed as ∂η m ∂ xk =1 Kcos η m ∂ |k| ∂ xk −sin η m ∂ K ∂ xk,(44) where ∂ k ∂ xk = β  ∂ r ∂ xk cos θ cosΘ+rcos θ∂ cosΘ ∂ xk, (45) ∂ K ∂ xk = β  ∂ r ∂ xk sin θ sinΘ+rsin θ∂ sinΘ ∂ xk. (46) The derivatives of M(n) sand  M(n) sappearing in 156 Copyright c  2007 Tech Science Press CMES, vol.22, no.2, pp.151-164, 2007 the above expressions are given by: ∂ M(n) s ∂ xk =0,n=0,(47) =1 D(ys) ∂ f1 ∂ xk (ys),n=1,(48) =1 D(1) ∂ f2 ∂ xk (1),n=2 with s=0, (49) =1 (1−ys2)D(1) ∂ f2 ∂ xk (1),(50) n=2 with s=1,2, ∂  M(1) s ∂ xk =1 D(ys) ∂ ˜ f1 ∂ xk (ys),(51) where ∂ f1 ∂ xk (x)=x3(4 ν −1) ∂ sj ∂ xk ˆ β l− ∂ sl ∂ xk ˆ β jsin θ (52) ∂ f2 ∂ xk (x)=−0.5[8 ν x4−(−x2+1)(−2x2q+1)] (53) ·njnl ∂ cos2Φ ∂ xk +ml ∂ sin2Φ ∂ xk +mjnl ∂ sin2Φ ∂ xk −ml ∂ cos2Φ ∂ xksin2 θ (54) ∂ ˜ f1 ∂ xk (x)=−0.5 ∂ sj ∂ xk ˆ β l+ ∂ sl ∂ xk ˆ β j ·(−2x2q+1)sin2 θ .(55) Finally, ∂ sj ∂ xk =nj ∂ cosΦ ∂ xk +mj ∂ sinΦ ∂ xk (56) ∂ ∂ xk = ∂ ∂ ˜xj ∂ ˜xj ∂ xk =Ljk ∂ ∂ ˜xj ,(57) where L1k=nk,L2k=mk,L3k=ˆ β k,and ∂ cosΘ ∂ ˜xj = δ j3 r−r3 r2 ∂ r ∂ ˜xj ,(58) ∂ sinΘ ∂ ˜xj =1 rr2−r2 3r ∂ r ∂ ˜xj −r3 δ j3(59) −r2−r2 3 r2 ∂ r ∂ ˜xj , ∂ cosΦ ∂ ˜xj = δ j1 r2−r2 3 (60) −r1 (r2−r2 3)3/2r ∂ r ∂ ˜xj −r3 δ j3, ∂ sinΦ ∂ ˜xj = δ j2 r2−r2 3 (61) −r1 (r2−r2 3)3/2r ∂ r ∂ ˜xj −r3 δ j3. Note that the finite integrals in (33-37) can be computed by the standard Gaussian quadrature, except for the integral in (36) where the function ∂η m ∂ xkappears. This function is weakly singular for θ → θ m(including a square root singularity) and accurate values of this integral can be achieved, e.g., by the singularity subtraction technique or a special Gauss-Jacobi quadrature. For further details and a comparison of both techniques, which showed that the special quadrature is slightly more accurate, see Criado (2005). The strains Eij(x,x)associated with the fundamental solution in displacements Uj(x,x)are given by Eij(x,x)=1 2 ∂ Ui ∂ xj (x,x)+ ∂ Uj ∂ xi (x,x), (62) and incorporating them into the constitutive law defining the elastic stiffnesses (2-3), yields the corresponding stresses Σij=2 μ (x)Eij(x,x)+ λ (x)Ekk(x,x) δ ij. (63) Then, substituting (62) into (63) and using (29) Analysis of Isotropic Elastic Solids 157 yields Σij(x−x)= exp  β β β ·(x−x)Σ0 ij(x−x)+Σg ij(x−x), (64) where the strongly singular term Σ0 ij(x−x)represents the stress tensor σ ij at xoriginated by a unit point force in direction at xin the homogeneous elastic isotropic material having Lamé constants μ 0and λ 0(see París and Cañas (1997); Aliabadi (2002)). The weakly singular grading term Σg ij(x−x)is expressed as: Σg ij(x−x)= μ 0 ∂ Ug i ∂ xj + ∂ Ug j ∂ xi − β iU0 j+Ug j− β jU0 i+Ug i + λ 0 ∂ Ug k ∂ xk − β kU0 k+Ug k δ ij.(65) It is important to notice the differences between the exponential prefactors in the expressions of the kernels Ujand Σijin (4) and (64) and their implications for the behaviour of these kernels. First, these prefactors imply that while Uj(x,x) depends on the position of both the source and field points, xand x,Σij(x−x)depends only on the difference of these points. This difference can be explained by the fact that the equilibrium equation for the displacement fundamental solution Ujincludes the graded material properties [Martin et al. (2002)], whereas the equilibrium equation for the stress fundamental solution Σij: ∂ Σij ∂ xj + δ i δ (x−x)=0,(66) with δ (x)being the three-dimensional Dirac delta, is naturally independent of these material properties. Second, taking the source point xfixed (e.g., at the origin of coordinates), the exponential prefactor in (4) diminishes the values of the displacements Ujin stiffer zones and increases them in more compliant zones, whereas the exponential prefactor in (64) increases the values of the stresses Σijin stiffer zones and diminishes them in more compliant zones. The traction fundamental solution Ti(x,x),associated with the unit outward normal vector n(x), is obtained from Σij(x−x)by the Cauchy lemma: Ti(x,x) =Σij(x−x)nj(x)(67) =exp  β β β ·(x−x)T0 i(x,x)+Tg i(x,x), (68) where, as for the stress, T0 i(x,x)represents the well-known strongly singular fundamental solution in tractions for a homogeneous material (with parameters μ 0and λ 0) [París and Cañas (1997); Aliabadi (2002)], and Tg i(x,x)is the weakly singular grading term obtained from Σg ij(x−x), Tg i(x,x)=Σg ij(x−x)nj(x). The behaviour of the stresses originated by a unit point force is illustrated in Figure 1, where three-dimensional plots of some components of the stress fundamental solution Σij(x−x)are shown, the origin of coordinates being coincident with the force application point, i.e. x=0. All the components presented show some singular behaviour at x=0. Nevertheless, as can be guessed looking at the values and/or variations of the stresses plotted, these singularities are not always of the same kind, which will be explained in the following. The differences between these plots in Figure 1 can be more easily understood if the symmetry and skew-symmetry planes in the corresponding elastic problems for exponentially graded and homogeneous materials are considered. The stronger singular behaviour of Σ0 ijat x=0 in comparison with Σg ijis dominant in Figures 1 a) and f). Slight differences between behaviour of stresses in these plots can be attributed to the effect of the exponential prefactor, not present in Figure 1 a), and to different contributions due to Σg ij. The components Σ0 ijcorresponding to plots in Figures 1 b) to e) vanish (which can be shown by applying standard symmetry arguments), hence the stress values shown in these plots are exclu- 158 Copyright c  2007 Tech Science Press CMES, vol.22, no.2, pp.151-164, 2007 (a) Σ111(x1,x2,0)(b) Σ231(x1,x2,0) (c) Σ113(x1,x2,0)(d) Σ123(x1,x2,0) (e) Σ333(x1,x2,0)(f) Σ333(0,x2,x3) Figure 1: Components of the Σijk kernel for μ 0=4, ν =0.35 and grading parameter β β β =(0,0,0.1). Analysis of Isotropic Elastic Solids 159 sively given by the components of Σg ijtimes the exponential prefactor. Finally, in order to show in a very clear way the differences between the stress fundamental solutions in homogeneous and exponentially graded materials, components Σ333(0,0,x3)and Σ0 333(0,0,x3)are plotted in Figure 2. It can be seen how the skew-symmetric character of stresses in the homogeneous material is changed in presence of a grading, in agreement with the above explained effect of the exponential prefactor. Figure 2: Component Σ0 333(0,0,x3)(corresponding to the homogeneous material) and Σ333(0,0,x3)for μ 0=4, ν =0.35 and grading parameter β β β =(0,0,5). 4 Boundary Element Method The boundary integral formulation for an isotropic, exponentially graded body Ω⊂R3 with a bounded (Lipschitz and piecewise smooth) boundary ∂ Ω=Γwill be briefly discussed in this section. The derivation follows the standard procedures for a homogeneous material [París and Cañas (1997); Aliabadi (2002)]. Starting from the 2nd Betti Theorem of reciprocity of work for a graded material, one can derive the corresponding Somigliana identity: Ci(x)ui(x)+− Γ Ti(x,x)ui(x)dS(x) =Γ Ui(x,x)ti(x)dS(x),(69) expressing the displacements ui(x)at a domain or boundary point x∈Ω∪Γin terms of the boundary displacements ui(x)and tractions ti(x), x∈Γ. Notice that zero body forces have been considered in (69). The strongly singular traction kernel integral is evaluated in the Cauchy principal value sense, and Ci(x)= lim ε →0+S ε (x)∩Ω Ti(x,x)dS(x)(70) is the coefficient tensor of the free term, S ε (x)being a spherical surface of radius ε centered at x and oriented by the unit normal vectors pointing to the center. It is important to note that, despite the complexity of the Tikernel expression, this evaluation is not a problem. The weakly singular grading term and the exponential prefactor in (68) will play no role in the limit procedure in (70). Thus, the value of Ciin (70) coincides with the value of Cifor the homogeneous isotropic material whose properties are defined by the Lamé constants λ 0and μ 0,i.e., Ci(x)= lim ε →0+S ε (x)∩Ω T0 i(x,x)dS(x).(71) Hence, Ci(x)= δ ifor x∈Ω,Ci(x)=1 2 δ ifor x∈Γsituated at a smooth part of Γ,andforan edge or corner point of Γ,Ci(x)is given by the size, shape and spatial orientation of the interior solid angle at x. A general explicit expression of the symmetric tensor Ci(x)in terms of the unit vectors tangential to the boundary edges and the unit outward normal vectors to the boundary surfaces at xcan be found in Mantiˇc (1993). For an unbounded domain Ω, the form of the Somigliana identity in (69) holds if the following radiation condition [Costabel and Dauge (1997)] for the displacement solution uiand the corresponding traction solution tiis fulfilled for any fixed x∈Ω: lim ρ →∞Γ ρ Ui(x,x)ti(x)−Ti(x,x)ui(x)dS(x) =0,(72) where Γ ρ is the spherical surface of radius ρ centered at the origin of coordinates. In fact, (72) implies that the integral in (72) vanishes for any sufficiently large Γ ρ including Γin its interior and