Recent developments in the evaluation of the 3D fundamental solution and its derivatives for transversely isotropic elastic materials

Recent developments in the evaluation of the 3D fundamental solution and its derivatives for transversely isotropic elastic materials∗ V. Mantiˇc, L. T´avara, J.E. Ortiz, F. Par´ıs Group of Elasticity and Strength of Materials, School of Engineering University of Seville, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain [email protected], [email protected], [email protected], par[email protected] Abstract. Explicit closed-form real-variable expressions of a fundamental solution and its derivatives for three-dimensional problems in transversely linear elastic isotropic solids are presented. The expressions of the fundamental solution in displacements Uik and its derivatives, originated by a unit point force, are valid for any combination of material properties and for any orientation of the radius vector between the source and field points. An expression of Uik in terms of the Stroh eigenvalues on the oblique plane normal to the radius vector is used as starting point. Working from this expression of Uik, a new approach (based on the application of the rotational symmetry of the material) for deducing the first and second order derivative kernels, Uik,j and Uik,jℓ respectively, has been developed. The expressions of the fundamental solution and its derivatives do not suffer from the difficulties of some previous expressions, obtained by other authors in different ways, with complex valued functions appearing for some combinations of material parameters and/or with division by zero for the radius vector at the rotational symmetry axis. The expressions of Uik,Uik,j and Uik,jℓ are presented in a form suitable for an efficient computational implementation in BEM codes. Keywords: transversely isotropic materials, Stroh formalism, fundamental solution, free-space Green’s functions, Somigliana identity, boundary integral equations, boundary element method 1. Introduction A comprehensive review of the history of the development of fundamental solutions and their derivatives for general anisotropic materials, and in particular for the transversely isotropic materials, was recently presented in [1, 2]. Thus, only the most relevant contributions closely related to the present work will be briefly mentioned. ∗Extended version of the Keynote lecture by V. Mantiˇc presented at the International Conference on Boundary Element Techniques IX (BETEQ) in Seville on 10 July, 2008. 1 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) Regarding anisotropic elastic materials, the expression of Uik (i, k = 1,2,3) obtained by Fredholm[3] using a 3D Fourier transformwas used in the pastas the basisfor different contributions aimed to obtain Uik expressions as simple and explicit as possible, e.g. Lifshitz and Rozentsweig [4]. The application of the Stroh formalism to anisotropic elasticity, see Ting [5], to evaluate Uik and its derivatives in 3D has shown to be an important approach that leaded to several substantial contributions starting in the 1970s with Mal´en [6] and Barnett [7], and more recently,Nakamuraand Tanuma [8], Tingand Lee [9], Lee [10, 11]and Tanuma[12]. An implementation in 3-D BEM of expressions of anisotropic fundamental solutions and its derivatives using the Stroh formalism, as introduced by Lee [10], was recently developed by Shiah et. al. [13, 15], Tan et. al. [14] and Buroni et. al. [16]. In the specific case of transversely isotropic materials the expressions deduced by Pan and Chou [17] are usually used in BEM codes [18, 19]. Loloi [20] and Ariza and Dom´ınguez [23] presented the expressions of the strongly singular and hypersingular kernel, respectively, obtained from Pan and Chou’s [17] solution. This solution [17] has several features which make somewhat cumbersome its implementation covering all possible cases: i) expressions depending on the values of △=√C11C33 −C13 −2C44 (positive, negative or zero) and in particular its complex-variable character for △<0; ii) a loss of precision and/or a division by zero for the spherical angle φ=πor φ= 0. In the previous expressions it is assumed that the x3-axis is the rotational symmetry axis of the transversal isotropic material. Although the difficulty with the degeneracy problem at φ=πhas been solved by Loloi [20] by means of an ad hoc approach, the mentioned features may still cause some difficulties in using this expression in further analytic deductions and in BEM development. Other authors as, Ding et. al. [21] also presented a fundamental solution for transversely isotropic media with complex-variable character. While, Yue [22] presented the fundamental solutions in displacement and tractions for different bimaterial systems, that can be used in the displacement boundary integral equation. The aim of the present work is to obtain completely general and closed-form realvariable expressions of Uik(x)and its first and second order derivatives, Uik,j (x)and Uik,jl(x)(i, j, k, l = 1,2,3), valid for any transversely isotropic material and simpler than the available expressions of the fundamental solution and its derivatives, with an arbitrary spatial orientation and allowing for an easy and efficient computational implementation. The use of Uik(x)expression introduced by Ting and Lee [9] and the transformation rule introduced therein, that was corrected and extended in [1], are crucial for the achievement of this aim. 2. Fundamental solution and its derivatives for linear elastic anisotropic materials Consider a cartesian coordinatesystem (x1, x2, x3)with the origin Oand the corresponding spherical coordinate system (r, θ, φ), see Fig. 1. The fundamental solution in displacements Uik originated by the unit concentrated force at O, and its derivatives, can be expressed in terms of the so-called modulation 2 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) functions depending only on the spherical angles θand φin the following form: Uik(r, θ, φ) = 1 4πr b Uik(θ, φ),(1) Uik,j(r, θ, φ) = 1 4πr2b Uik;j(θ, φ),(2) Uik,jℓ(r, θ, φ) = 1 4πr3b Uik;jℓ(θ, φ),(3) where i, j, k, l = 1,2,3and ‘widehat’ symbol denotes the modulation functions of the fundamental solution, and also of its derivatives, with respect to cartesian coordinates. The modulation functions are independent of the radius rin spherical coordinates. While a comma between subscripts denotes a derivative with respect to a cartesian coordinate, a semicolon in a modulation function is used only (mnemotechnically) to represent that this function is associated to the corresponding derivative kernel in (2) or (3). Mal´en [6] and others showed that b Uik(θ, φ) = Hik(θ, φ),(4) where Hik is the Barnett-Lothe tensor of the Stroh formalism of the anisotropic elasticity. A novel approach for obtaining an analytic expression of Hik was introduced by Ting and Lee [9]. By applying standard rules for differentiation with respect to cartesian coordinates in (2) and (3) we obtain b Uik;j(θ, φ) =Hik,j (r, θ, φ)r−Hik(φ, θ)r,j,(5) b Uik;jℓ(θ, φ) =Hik,jℓ(r, θ, φ)r2−Hik,j(r, θ, φ)r,ℓr−Hik,ℓ(r, θ, φ)r,j r −Hik(θ, φ) (δjℓ −3r,jr,ℓ),(6) where r,j(θ, φ) = rj(r, θ, φ)/r. By conveniently introducing new modulation functions obtained from the derivatives of Hik(θ, φ)as Hik,j (r, θ, φ) = r−1b Hik;j(θ, φ),(7) Hik,jℓ(r, θ, φ) = r−2b Hik;jℓ(θ, φ),(8) the modulation functions of the derivatives of Uik can be expressed in terms of the modulation functions Hik,b Hik;jand b Hik;jℓ, b Uik;j(θ, φ) = b Hik;j(θ, φ)−Hik(θ, φ)r,j ,(9) b Uik;jℓ(θ, φ) = b Hik;jℓ(θ, φ)−b Hik;j(θ, φ)r,ℓ −b Hik;ℓ(θ, φ)r,j −Hik(θ, φ)(δjℓ −3r,jr,ℓ).(10) Thus, in order to achieve expressions of the modulations functions b Uik,b Uik;jand b Uik;jℓ suitable for implementationin BEM codes we need to deduce general and as simple and compact as possible expressions of modular functions Hik,b Hik;jand b Hik;jℓ. As will be shown, for transversely isotropic materials simple and compact expressions of these modulation functions are only needed for a particular coordinate plane including the x3- axis, the plane x2= 0 being used in the present work. 3 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) 3. Modulation function of the fundamental solution and of its derivatives for linear elastic transversally isotropic materials 3.1. Transformation laws The key idea in the present approach for the evaluation of the modulation functions b Uik;j and b Uik;jℓ is to generalize the procedure proposed by Ting and Lee [9] for the evaluation of b Uik =Hik. This procedure will be sketched briefly in what follows. Let us define a vector b xcorresponding to a vector x= (x1, x2, x3)as b x= (r12,0, x3),where r12 =qx2 1+x2 2.(11) Let c= cos φ=x3/r and s= sin φ=r12/r, the angle φ,0≤φ≤π, being shown in Fig. 1. Then, defining n= (c, 0,−s)and m= (0,1,0),[n,m,b x/r]forms a righthanded triad. Figure 1: Points xand b xin spherical coordinates associated to a transversely isotropic material. Assuming that the x3-axis is the rotational symmetry axis of the material considered here, general expressions of the modulation functions (tensors) b Uik,b Uik;jand b Uik;jℓ for any x, in terms of cos and sin functionsof spherical angles φand θof x, see Fig. 1, can be obtained from the expressions of these modulation functions for the corresponding point b xby the following transformation of the components of these functions: b Uik(x) = ΩiaΩkb b Uab(b x),(12) b Uik;j(x) = ΩiaΩkbΩjc b Uab;c(b x),(13) b Uik;jℓ(x) = ΩiaΩkbΩjcΩℓd b Uab;cd(b x),(14) 4 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) where the rotation matrix Ωis defined as Ω(θ) =   cos θ−sin θ0 sin θcos θ0 0 0 1 ,(15) the angle θ,0≤θ < 2π, being shown in Fig. 1. Hence, if expressions of the modulation functions for a point b x(situated on the symmetry plane x2= 0) suitable for an efficient computer implementation are available, then it can be expected that the above transformation rules will lead to an efficient computer implementation of these modulation functions for any general point xas well. As will be shown in what follows this is the case of transversely isotropic materials with the x3-axis being rotational symmetry axis. It should be stressed that the transformation rules in (12-14) do not hold for general anisotropic materials (or any other orientation of transversely isotropic materials) as these rules are exclusively associated to the above mentioned rotational symmetry of the material with respect to the x3-axis. With reference to the plane x2= 0 (i.e. θ= 0) it will be useful to evaluate the particular expressions of Ω(θ)and its derivatives with respect to θat θ= 0, Ω(0) =   1 0 0 0 1 0 0 0 1 ,Ω′(0) =   0−1 0 100 000 ,(16) Ω′′ (0) =  −1 0 0 0−1 0 0 0 0 . The expressions of the derivatives with respect to cartesian coordinates of a modulation function f(θ, φ)for θ= 0 take the simple form: f,1=r−1c∂φf, f,2=r−1s−1∂θf, f,3=−r−1s∂φf, (17) and f,11 =r−2c2∂φφf−2cs∂φf, f,22 =r−2s−2∂θθf+cs−1∂φf, f,33 =r−2s2∂φφf+ 2cs∂φf, f,12 =r−2cs−1∂φθf−s−2∂θf,(18) f,13 =r−2−cs∂φφf−(c2−s2)∂φf, f,23 =−r−2∂φθf. The components of the unit vector b xare r,1=s, r,2= 0, r,3=c. (19) 3.2. Modulation function b Uik Explicit expressions for the non-zero terms of the modulation function b Uik(θ= 0, φ) = Hik(θ= 0, φ)in the plane x2= 0, can be obtained as shown in [9], see also [1], H(b x) =   H11 0H13 0H22 0 H13 0H33 ,(20) 5 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) where Hik =Hik(0, φ)are expressed as H11 =1 C66β3 +C44c2+C33s2 C11C44gh −f ξ, H22 =1 C11g+f ξ, H33 =1 gh h+c2 C44 +s2 C11 , H13 =˜ H13s, (21) and ˜ H13 =(C13 +C44)c C11C44gh , β3=C66c2+C44s2 C66 1/2 , h=c4+ηs2c2 C11C44 +C33s4 C11 1/2 , g=2(h+c2) + ηs2 C11C44 1/2 , (22) ξ=g(h+gβ3+β2 3), f=h+c2 C66 +gh C66β3 +C33s2 C11C44 , η=C11C33 −C2 13 −2C13C44, (23) β3,h,gand ξbeing bounded and positive dimensionless functions given in terms of c and s[9]. This is a key issue for the expressions of modulation functions presented in the following sections, as these dimensionless functions are the only functions appearing in the denominators. This fact will guarantee the non-singular character of the modulation function expressions presented. Finally, bringing together Equations (1), (4), (12), (15) and (20-23), an explicit and completelygeneralexpressionfor the fundamentalsolutionU(x)in atransverselyisotropic material is obtained. The form of this expression suitable for a computational implementation obtained by performing explicitly the transforms indicated in (12) is given in [1]. The presence of several zero components in H(b x)and Ωleads to simple and short expressions of the components of H(x). 3.3. Modulation function b Uik;j In order to evaluate b Uik;j(θ= 0, φ), we first evaluate the derivatives indicated in (5). By applying formulas in (17) to (4), the following expressions of b Hik;j(θ= 0, φ)defined 6 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) in (7) are easily obtained in terms of derivatives of Hik(θ, φ): b Hik;1 =c∂φHik(0, φ), b Hik;2 =s−1∂θHik(0, φ),(24) b Hik;3 =−s∂φHik(0, φ). The required expressions of the first order partial derivatives of Hik(θ, φ)with respect to angles θand φcan be evaluated from (12) in view of (4) as follows: ∂φHik(θ, φ) =Ωim(θ)Ωjn(θ)∂φHmn(0, φ)(25) ∂θHik(θ, φ) = Ω′ im(θ)Ωjn(θ) + Ωim(θ)Ω′ jn(θ)Hmn(0, φ).(26) For the particular case of θ= 0, and in view of (16), derivatives (25) can be expressed as ∂φHik(0, φ) =Ωim(0)Ωjn(0)∂φHmn(0, φ) =  H′ 11 0H′ 13 0H′ 22 0 H′ 13 0H′ 33  ,(27) where H′ mn =∂φHmn(θ= 0, φ)represent the derivativesof the correspondingfunctions defined in (21) and will be expressed explicitly later on, and derivatives (26) as ∂θHik(0, φ) = Ω′ im(0)Ωjn(0) + Ωim(0)Ω′ jn(0)Hmn(0, φ) =  0H11 −H22 0 H11 −H22 0H13 0H13 0 ,(28) where Hmn =Hmn(0, φ). The completely general and closed-form expressions of b Uik;j(b x)can now be easily obtained by substitution of (24) with (27-28) and (19-20) into (9) considered for a point b x(i.e. θ= 0). The expressions obtained are presented in a compact form suitable for computer implementation: b U11;1 =H′ 11c−H11s, b U12;2 =˜ H12s, b U11;3 =−H′ 11s−H11c, b U22;1 =H′ 22c−H22s, b U23;2 =˜ H13,b U22;3 =−H′ 22s−H22c, b U33;1 =H′ 33c−H33s, b U33;3 =−H′ 33s−H33c, b U13;1 =H′ 13c−H13s, b U13;3 =−H′ 13s−H13c, (29) where H′ 11 =˜ H11s, H′ 22 =˜ H22s, H′ 33 =˜ H33s, (30) 7 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) ˜ H11 =−˜ β3 C66β2 3−C44c2+C33s2 C11C44gh ˜ h h+˜g g!+2 (C33 −C44)c C11C44gh −1 ξ ˜ f−˜ ξf ξ!, ˜ H22 =−˜g C11g2+1 ξ ˜ f−˜ ξf ξ!, ˜ H33 =−1 gh 2c C11 +˜ h−2c C44 !−H33 ˜ h h+˜g g!,(31) ˜ H12 =C33 C11C44gh −ηc2+C33C44s2 C11C44 (h+c2)1 C11gh +g C66β3ξ, +1 ξη−2C33C66 C11C44C66 +C44g C2 66β3, H′ 13 =C13 +C44 C11C44gh c2−s2−γcs, γ=g′ g+h′ h, and ˜ β3=(C44 −C66)c C66β3 , ˜ h=1 hηc C11C44 c2−s2+2C33cs2 C11 −2c3, ˜g=1 gηc C11C44 +˜ h−2c,(32) ˜ ξ=g2 (C44 −C66)c C66 +˜ h+β3˜g+˜ β3g+˜gξ g, ˜ f=˜ h−2c C66 +1 C66β3 h˜g+˜ hg −˜ β3gh β3!+2C33c C11C44 . In (32), ‘tilde’ functions are related to β′ 3, h′, g′, ξ′, and f′, that represent the first order derivatives with respect to the angle φ(at θ= 0) of the corresponding functions defined in (22-23), in the following way: β′ 3=˜ β3s, h′=˜ hs, g′= ˜gs, ξ′=˜ ξs, f′=˜ fs. (33) ‘Tilde’ functions are introduced in order to avoid the presence of divisions by s, that will lead eventually in divisions by zero when φ= 0. The remaining components of b Uik;j(b x)vanish, b U12;1 =b U23;1 =b U11;2 =b U13;2 =b U22;2 =b U33;2 =b U12;3 =b U23;3 = 0.(34) It should be mentioned that in the original expressionof b U12,2(b x), directly obtainedby differentiating (24), the term (H11 −H22)/s appeared, which would lead to zero divided by zero for points at the x3-axis. This term, which has a finite limit value for φ→0or π, has been rewritten in the form ˜ H12s, which is well defined for any point with r > 0. 8 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) 3.4. Modulation function b Uik;jℓ To determine b Uik;jℓ(θ= 0, φ), we first evaluate the second order derivatives appearing in (6). Using formulasin (18) applied to (4), the following expressions of b Hik;jℓ(θ= 0, φ) defined in (8) are easily obtained in terms of derivatives of Hik(θ, φ): b Hik;11 =c2∂2 φφHik(0, φ)−2cs∂φHik(0, φ), b Hik;22 =s−2∂2 θθHik(0, φ) + cs−1∂φHik(0, φ), b Hik;33 =s2∂2 φφHik(0, φ) + 2cs∂φHik(0, φ), b Hik;12 =cs−1∂2 φθHik(0, φ)−s−2∂θHik(0, φ),(35) b Hik;13 =−cs∂2 φφHik(0, φ)−(c2−s2)∂φHik(0, φ), b Hik;23 =−∂2 φθHik(0, φ). Recall the following symmetry relations, due to the interchanging the order of taking partial derivatives, b Hik;jℓ =b Hik;ℓj .(36) The required expressions of the second order partial derivatives of Hik(θ, φ)with respect to angles θand φcan be evaluated from (12) in view of (4) as follows: ∂2 φφHik(θ, φ) =Ωim(θ)Ωjn(θ)∂φφHmn(0, φ),(37) ∂2 φθHik(θ, φ) = Ω′ im(θ)Ωjn(θ) + Ωim(θ)Ω′ jn(θ)∂φHmn(0, φ),(38) ∂2 θθHik(θ, φ) = Ω′′ im(θ)Ωjn(θ) + 2Ω′ im(θ)Ω′ jn(θ) + Ωim(θ)Ω′′ jn(θ)Hmn(0, φ). (39) For the particular case of θ= 0, and in view of (16), derivatives (37-39)can be expressed as: ∂2 φφHik(0, φ) =Ωim(0)Ωjn(0)∂φφHmn(0, φ) =  H′′ 11 0H′′ 13 0H′′ 22 0 H′′ 13 0H′′ 33  ,(40) ∂2 φθHik(0, φ) = Ω′ im(0)Ωjn(0) + Ωim(0)Ω′ jn(0)∂φHmn(0, φ) =  0H′ 11 −H′ 22 0 H′ 11 −H′ 22 0H′ 13 0H′ 13 0 ,(41) ∂2 θθHik(0, φ) = Ω′′ im(0)Ωjn(0) + 2Ω′ im(0)Ω′ jn(0) + Ωim(0)Ω′′ jn(0)Hmn(0, φ) = −2(H11 −H22) 0 −H13 0 2(H11 −H22) 0 −H13 0 0  ,(42) 9 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) Figure 8: U11,11 in plane x2= 0, for △<0. Figure 9: U11,11 in plane x2= 0, for △= 0. 16 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) Figure 10: U11,11 in plane x2= 0, for △>0. Figure 11: U13 in plane x2= 0, for △<0. 17 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) Figure 12: U13 in plane x2= 0, for △= 0. Figure 13: U13 in plane x2= 0, for △>0. 18 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) Figure 14: U13,1in plane x2= 0, for △<0. Figure 15: U13,1in plane x2= 0, for △= 0. 19 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) Figure 16: U13,1in plane x2= 0, for △>0. Figure 17: U13,11 in plane x2= 0, for △<0. 20 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) Figure 18: U13,11 in plane x2= 0, for △= 0. Figure 19: U13,11 in plane x2= 0, for △>0. 21 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) Figure 20: U23,2in plane x2= 0, for △<0. Figure 21: U23,2in plane x2= 0, for △= 0. 22 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) Figure 22: U23,2in plane x2= 0, for △>0. Figure 23: U33 in plane x2= 0, for △<0. 23 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) Figure 24: U33 in plane x2= 0, for △= 0. Figure 25: U33 in plane x2= 0, for △>0. 24 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) Figure 26: U33,1in plane x2= 0, for △<0. Figure 27: U33,1in plane x2= 0, for △= 0. 25 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) implicit none c−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− c Model parameters c−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− c Ela stic s t i f f n e s s constants real ∗8 c11 , c12 , c13 , c33 , c44 , c66 c H tensor and Displacement fundamental solutio n real ∗8 Htheta (3 ,3) , gu (3 ,3) c First and second order d er ivat i v e of the c displacement fundamental soluti on real ∗8 upf (3 ,3 ,3) , uppf (3 ,3 ,3 ,3) c First and second order d e r ivat i ves (when x2=0) real ∗8 upg (3 ,3 ,3) , uppg (3 ,3 ,3 ,3) c First and second order d e riv a t iv es modulation c functi ons ( indepedent of r ) real ∗8 up (3 ,3 ,3) , upp (3 ,3 ,3 ,3) c cos [ phi ] and sin [ phi ]; cos [ theta ] and sin [ theta ] real ∗8 c , s , ct , st c Aux . variables real ∗8 eta , h , g , beta3 , xi , f , x1 , x2 , x3 ,gamma c Derivatives of aux . variables real ∗8 hp , gp , beta3p , xip , fp real ∗8 hg , gg , beta3g , xig , fg real ∗8 hpp , gpp , beta3pp , xipp , fpp c extra variables real ∗8 pi , rt , r12 c New aux . variables real ∗8 H11 , H22 , H33 , H13 real ∗8 H11p , H22p , H33p , H13p , H12g , H13g real ∗8 H11pp , H22pp , H33pp , H13pp , H11g , H22g , H33g c Radius vector between the source point and f i e l d point real ∗8 x (3) c−−−−−−−−−−−−−−−−−−− Local variables −−−−−−−−−−−−−−−−− integer i , j , k , l ,m, n , flagsd c ====================================================== pi =4.d0∗datan ( 1. d0 ) rt =dsqrt (x(1)∗∗2+x(2)∗∗2+x (3)∗ ∗ 2) x1=x (1 ) ; x2=x (2 ); x3=x (3) c−−− Calculate aux variables −−− r12= dsqrt ( x(1)∗∗2+x (2)∗ ∗ 2) c=x ( 3 ) / r t s=r12 / r t eta=c11∗c33−c13∗∗2−2.∗c13∗c44 beta3 = ds qrt ( ( c66∗c∗∗2+c44∗s∗ ∗ 2)/ c66 ) h=dsqr t ( c∗∗4+( eta∗s∗s∗c∗c / ( c11∗c44 )) +( c33 ∗( s ∗ ∗ 4)/ c11 )) g=dsqrt (2∗( h+c∗∗2)+ eta ∗s∗ ∗ 2/( c11∗c44 )) xi=g∗(h+g∗beta3 +beta3 ∗beta3 ) c−− Calculate Hx matrix components −−− H11 =(1/( c66∗beta3 ) ) + ( ( c44∗c∗c+c33∗s∗s ) / ( c11∗c44∗g∗h )) 32 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) &−(1/ xi )∗(( h+c∗c ) / c66+g∗h / ( c66∗beta3 )+c33∗s∗s / ( c11∗c44 )) H22 =(1/( c11∗g ) )+( 1/ xi )∗(( h+c∗c )/ c66+g∗h / ( c66∗beta3 )+c33∗s∗s / &(c11∗c44 ) ) H33 =(1/( h∗g ))∗( ( h+c∗c ) / c44+s∗s / c11 ) H13=c∗s∗( c13+c44 ) / ( c11∗c44∗g∗h) c−− Calculate Cos ( thet a ) and Sin ( theta ) rotatio n matrix −−− ct =x ( 1) / r12 s t=x ( 2) / r12 c−−− Calculate Htheta tensor components −−− c E x p li c it r ot ati on ( f a st e r than i m p l ic i t ) Htheta (1 ,1)=H11∗ct ∗∗2+H22∗st ∗∗2 Htheta (1 ,2)=(H11−H22)∗ct ∗st Htheta (1 ,3)=H13∗ct Htheta (2 ,2)=H22∗ct ∗∗2+H11∗st ∗∗2 Htheta (2 ,3)=H13∗st Htheta (3 ,3)=H33 Htheta (2 ,1)= Htheta (1 ,2) Htheta (3 ,1)= Htheta (1 ,3) Htheta (3 ,2)= Htheta (2 ,3) gu=Htheta / ( 4 . d0∗pi∗r t ) c−−− Calculate aux variables −−− f =(h+c∗c ) / C66+g∗h / ( C66∗beta3 )+C33∗s∗s / ( C11∗C44) H12g=C33 / ( C11∗C44∗g∗h)−( et a ∗c∗∗2+C33∗C44∗s∗ ∗ 2)∗(1/( C11∗g∗h)+ &g / ( C66∗beta3 ∗xi ) ) / ( C11∗C44∗(h+c ∗∗2))+(( eta −2∗C33∗C66 ) / &(C11∗C44∗C66 )+(C44∗g ) / ( C66∗∗2∗beta3 ) ) / xi beta3g =(C44−C66)∗c / ( C66∗beta3 ) hg =(( ( eta ∗c )∗( c∗∗2−s∗ ∗ 2) )/( C11∗C44 )+(2.∗C33∗c∗s∗∗2)/C11− &2.∗c∗∗3)/ h gg =(( eta ∗c ) / ( C11∗C44)+hg−2.∗c )/ g xig=g∗(2∗c∗(C44−C66 ) / C66+hg+beta3 ∗gg+beta3g∗g)+gg∗xi / g fg=2∗c∗C33 / ( C11∗C44)+( hg−2∗c ) / C66+ &(h∗gg+hg∗g−beta3g ∗g∗h / beta3 ) / ( C66∗beta3 ) H33g=(2∗c /C11+(hg−2∗c ) / C44 ) / ( h∗g)−H33∗(gg / g+hg / h) H22g=−gg / ( C11∗g∗∗2)+( fg−xig∗f / xi ) / xi H11g=−beta3g / ( C66∗beta3 ∗∗2)−(C44∗c∗∗2+C33∗s∗∗2)∗( gg / g+hg / h ) / &(C11∗C44∗h∗g)+2∗(C33−C44)∗c / ( C11∗C44∗h∗g)−(fg−xig∗f / xi ) / xi beta3p=beta3g∗s 33 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) hp=hg∗s gp=gg∗s xip=xig∗s fp=fg∗s H11p=H11g∗s H22p=H22g∗s H33p=H33g∗s gamma=gp / g+hp / h H13g=c ∗(C13+C44 ) / ( C11∗C44∗g∗h ) H13p=(C13+C44 )∗( c∗∗2−s∗∗2−c∗s∗gp / g−c∗s∗hp / h ) / ( C11∗C44∗h∗g ) c−− Fir st Order Derivatives of GU −− (d2U / dxidxj ) when x2=0 upg =0. d0 upg (1 ,1 ,1)=H11p∗c−H11∗s upg (2 ,2 ,1)=H22p∗c−H22∗s upg (3 ,3 ,1)=H33p∗c−H33∗s upg (1 ,3 ,1)=H13p∗c−H13∗s upg (3 ,1 ,1)= upg (1 ,3 ,1) upg (2 ,3 ,2)=H13g upg (1 ,2 ,2)=H12g∗s upg (2 ,1 ,2)= upg (1 ,2 ,2) upg (3 ,2 ,2)= upg (2 ,3 ,2) upg(1 ,1 ,3)= −H11p∗s−H11∗c upg(2 ,2 ,3)= −H22p∗s−H22∗c upg(3 ,3 ,3)= −H33p∗s−H33∗c upg(1 ,3 ,3)= −H13p∗s−H13∗c upg (3 ,1 ,3)= upg (1 ,3 ,3) c−−− Calculate up tensor components from upg −−− c E x p li c it r ot ati on ( f a st e r than i m p l ic i t ) up (1 ,1 ,1)= ct ∗( upg (1 ,1 ,1)∗ct ∗∗2+(2∗upg (1 ,2 ,2)+ upg (2 ,2 ,1))∗st ∗∗2) up (1 ,1 ,2)= st ∗(( upg (1 ,1 ,1) −2∗upg (1 ,2 ,2))∗ct ∗∗2+upg (2 ,2 ,1)∗st ∗∗2) up (1 ,1 ,3)= upg (1 ,1 ,3)∗ct ∗∗2+upg (2 ,2 ,3)∗st ∗∗2 up (1 ,2 ,1)= st ∗( upg (1 ,2 ,2)∗st ∗∗2+( upg (1 ,1 ,1) −upg (1 ,2 ,2) &−upg (2 ,2 ,1))∗ct ∗∗2) up (1 ,2 ,2)= ct ∗( upg (1 ,2 ,2)∗ct ∗∗2+( upg(1 ,1 ,1) −upg (1 ,2 ,2) &−upg (2 ,2 ,1))∗st ∗∗2) up (1 ,2 ,3)= ct ∗st ∗( upg (1 ,1 ,3) −upg (2 ,2 ,3)) up (1 ,3 ,1)= upg (1 ,3 ,1)∗ct ∗∗2+upg (2 ,3 ,2)∗st ∗∗2 up (1 ,3 ,2)= ct ∗st ∗( upg (1 ,3 ,1) −upg (2 ,3 ,2)) up (1 ,3 ,3)= upg (1 ,3 ,3)∗ct up (2 ,2 ,1)= ct ∗(( upg(1 ,1 ,1) −2∗upg (1 ,2 ,2))∗st ∗∗2+upg (2 ,2 ,1)∗ct ∗∗2) up (2 ,2 ,2)= st ∗( upg (1 ,1 ,1)∗st ∗∗2+(2∗upg (1 ,2 ,2)+ upg (2 ,2 ,1))∗ct ∗∗2) up (2 ,2 ,3)= upg (1 ,1 ,3)∗st ∗∗2+upg (2 ,2 ,3)∗ct ∗∗2 up (2 ,3 ,1)= ct ∗st ∗( upg (1 ,3 ,1) −upg (2 ,3 ,2)) up (2 ,3 ,2)= upg (1 ,3 ,1)∗st ∗∗2+upg (2 ,3 ,2)∗ct ∗∗2 up (2 ,3 ,3)= upg (1 ,3 ,3)∗st up (3 ,3 ,1)= upg (3 ,3 ,1)∗ct 34 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) up (3 ,3 ,2)= upg (3 ,3 ,1)∗st up (3 ,3 ,3)= upg (3 ,3 ,3) up (2 ,1 ,1)=up (1 ,2 ,1) up (2 ,1 ,2)=up (1 ,2 ,2) up (2 ,1 ,3)=up (1 ,2 ,3) up (3 ,1 ,1)=up (1 ,3 ,1) up (3 ,1 ,2)=up (1 ,3 ,2) up (3 ,1 ,3)=up (1 ,3 ,3) up (3 ,2 ,1)=up (2 ,3 ,1) up (3 ,2 ,2)=up (2 ,3 ,2) up (3 ,2 ,3)=up (2 ,3 ,3) upf=up / ( 4 . d0∗pi ∗r t ∗∗2) c−−− Calculate NEW aux variables −−− beta3pp=−beta3p ∗∗2/ beta3 +(C44−C66 )∗( c∗∗2−s∗ ∗ 2)/( C66∗beta3 ) hpp=(−hp∗∗2−2∗h∗∗2+2∗c∗∗2∗s∗∗2∗(3+3∗C33 /C11−eta / ( C11∗C44))+ &eta ∗(( c∗∗2−s∗ ∗ 2)∗ ∗ 2)/( C11∗C44 ) ) / h gpp=(−gp∗∗2+( c∗∗2−s∗∗2)∗( eta / ( C11∗C44)−2)+hpp ) / g xipp =2∗gp ∗( xip /g−gp∗xi / ( g∗∗2))+ g ∗(2∗(C44−C66 )∗( c∗∗2−s∗ ∗ 2)/C66 &+hpp+2∗beta3p ∗gp+beta3 ∗gpp+beta3pp∗g)+gpp∗xi / g fpp=2∗C33∗( c∗∗2−s∗ ∗ 2)/( C11∗C44 )+( hpp −2∗(c∗∗2−s∗ ∗ 2))/ C66+ &(2∗hp∗gp+h∗gpp+hpp∗g−(beta3pp∗g∗h+beta3p ∗gp∗h+beta3p ∗hp∗g ) / beta3 &+2∗g∗h∗beta3p ∗∗2/( beta3 ∗ ∗ 2 )) /( C66∗beta3 ) &−(h∗gp+hp∗g)∗beta3p / ( C66∗beta3 ∗∗2) H13pp=−H13p∗gamma+(C13+C44)∗(−4∗c∗s−((c∗∗2−s∗∗2)∗gp+c∗s∗gpp ) / g− &((c∗∗2−s∗∗2)∗hp+c∗s∗hpp )/ h+c∗s∗(( gp / g )∗∗2+( hp / h )∗ ∗ 2)) &/(C11∗C44∗h∗g) H33pp=−(gp / ( h∗g∗∗2)+ hp / ( g∗h∗∗2))∗(2∗c∗s / C11+(hp−2∗c∗s ) / C44)+ &(2∗(c∗∗2−s∗ ∗ 2)/ C11+(hpp −2∗(c∗∗2−s∗ ∗ 2))/ C44 ) / ( g∗h)−H33p∗gamma− &H33 ∗(( gpp∗g−gp ∗ ∗ 2)/( g∗∗2)+( hpp∗h−hp ∗ ∗ 2)/( h∗∗2)) H22pp =(2.∗( gp∗∗2)−gpp∗g ) / ( C11∗g∗∗3)+ xip ∗( f∗xip / xi−fp ) / ( xi ∗∗2)+ &(fpp −(( xipp∗f+xip∗fp )∗xi−f∗xip ∗ ∗ 2)/( xi ∗ ∗ 2))/ xi H11pp=(2∗beta3p∗∗2−beta3pp∗beta3 ) / ( C66∗beta3 ∗∗3)−4∗c∗s∗(C33−C44) &∗gamma / ( C11∗C44∗g∗h )+(C44∗c∗∗2+C33∗s∗∗2)∗(gamma∗∗2− &(gpp∗g−gp ∗ ∗ 2)/( g∗∗2)−(hpp∗h−hp ∗ ∗ 2)/( h ∗ ∗ 2)) /( C11∗C44∗g∗h)+ &2∗(C33−C44 )∗( c∗∗2−s∗ ∗ 2)/( C11∗C44∗g∗h)−xip ∗( f∗xip / xi−fp ) / ( xi∗∗2)− &(fpp −(( xipp∗f+xip∗fp )∗xi−f∗xip ∗ ∗ 2)/( xi ∗ ∗ 2))/ xi c−− Second Order Derivatives of GU −− (d2U / dxidxj ) when x2=0 35 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) uppg =0. d0 uppg (1 ,1 ,1 ,1)= H11pp∗c∗∗2−4∗H11p∗c∗s−H11∗c∗∗2+2∗H11∗s∗∗2 uppg (1 ,1 ,1 ,3)= −H11pp∗c∗s−2∗H11p∗( c∗∗2−s∗∗2)+3∗H11∗c∗s uppg (1 ,1 ,2 ,2)=H11g∗c−2∗H12g−H11 uppg (1 ,1 ,3 ,3)= H11pp∗s∗∗2+4∗H11p∗c∗s+2∗H11∗c∗∗2−H11∗s∗∗2 uppg (1 ,2 ,1 ,2)=( H11g−H22g)∗c−H12g−H11+H22 uppg (1 ,2 ,2 ,3)= −H11p+H22p−H12g∗c∗s uppg (1 ,3 ,1 ,1)= H13pp∗c∗∗2−4∗H13p∗c∗s−H13∗c∗∗2+2∗H13∗s∗∗2 uppg (1 ,3 ,1 ,3)= −H13pp∗c∗s−2∗H13p∗( c∗∗2−s∗∗2)+3∗H13∗c∗s uppg (1 ,3 ,2 ,2)= −H13g∗(2∗s+gamma∗c)−H13 uppg (1 ,3 ,3 ,3)= H13pp∗s∗∗2+4∗H13p∗c∗s+2∗H13∗c∗∗2−H13∗s∗∗2 uppg (2 ,2 ,1 ,1)= H22pp∗c∗∗2−4∗H22p∗c∗s−H22∗c∗∗2+2∗H22∗s∗∗2 uppg (2 ,2 ,1 ,3)= −H22pp∗c∗s−2∗H22p∗( c∗∗2−s∗∗2)+3∗H22∗c∗s uppg (2 ,2 ,2 ,2)=H22g∗c+2∗H12g−H22 uppg (2 ,2 ,3 ,3)= H22pp∗s∗∗2+4∗H22p∗c∗s+2∗H22∗c∗∗2−H22∗s∗∗2 uppg (2 ,3 ,2 ,3)= −H13p−H13g∗c uppg (3 ,3 ,1 ,1)= H33pp∗c∗∗2−4∗H33p∗c∗s−H33∗c∗∗2+2∗H33∗s∗∗2 uppg (3 ,3 ,1 ,3)= −H33pp∗c∗s−2∗H33p∗( c∗∗2−s∗∗2)+3∗H33∗c∗s uppg (3 ,3 ,2 ,2)= −H33+H33g∗c uppg (3 ,3 ,3 ,3)= H33pp∗s∗∗2+4∗H33p∗c∗s+2∗H33∗c∗∗2−H33∗s∗∗2 uppg (2 ,3 ,1 ,2)= uppg (1 ,3 ,2 ,2) c E x p li c it r ot ati on ( f a st e r than i m p l ic i t ) upp (1 ,1 ,1 ,1)= uppg (1 ,1 ,1 ,1)∗ct ∗∗4+uppg (1 ,1 ,2 ,2)∗ct ∗∗2∗st ∗∗2+4 &∗uppg (1 ,2 ,1 ,2)∗ct ∗∗2∗st ∗∗2+uppg (2 ,2 ,1 ,1)∗ct ∗∗2∗st ∗∗2+uppg (2 ,2 ,2 ,2) &∗st ∗∗4 upp (1 ,1 ,1 ,2)= uppg (1 ,1 ,1 ,1)∗st ∗ct ∗∗3−uppg (1 ,1 ,2 ,2)∗st ∗ct ∗∗3 &−2∗uppg (1 ,2 ,1 ,2)∗s t ∗ct ∗∗3+2∗uppg (1 ,2 ,1 ,2)∗ct ∗st ∗∗3+uppg (2 ,2 ,1 ,1) &∗ct ∗st ∗∗3−uppg (2 ,2 ,2 ,2)∗ct ∗st ∗∗3 upp (1 ,1 ,1 ,3)= uppg (1 ,1 ,1 ,3)∗ct ∗∗3+2∗uppg (1 ,2 ,2 ,3)∗ct ∗st ∗∗2 &+uppg (2 ,2 ,1 ,3)∗ct ∗st ∗∗2 upp (1 ,1 ,2 ,2)= uppg (1 ,1 ,1 ,1)∗ct ∗∗2∗st ∗∗2+uppg (1 ,1 ,2 ,2)∗ct ∗∗4 &−4∗uppg (1 ,2 ,1 ,2)∗ct ∗∗2∗st ∗∗2+uppg (2 ,2 ,1 ,1)∗st ∗∗4+uppg (2 ,2 ,2 ,2) &∗ct ∗∗2∗st ∗∗2 upp (1 ,1 ,2 ,3)= uppg (1 ,1 ,1 ,3)∗st ∗ct ∗∗2−2∗uppg (1 ,2 ,2 ,3)∗st ∗ct ∗∗2 &+uppg (2 ,2 ,1 ,3)∗st ∗∗3 upp (1 ,1 ,3 ,3)= uppg (1 ,1 ,3 ,3)∗ct ∗∗2+uppg (2 ,2 ,3 ,3)∗s t ∗∗2 upp (1 ,2 ,1 ,1)= uppg (1 ,1 ,1 ,1)∗st ∗ct ∗∗3+uppg (1 ,1 ,2 ,2)∗ct ∗s t ∗∗3 &−2∗uppg (1 ,2 ,1 ,2)∗s t ∗ct ∗∗3+2∗uppg (1 ,2 ,1 ,2)∗ct ∗st ∗∗3−uppg (2 ,2 ,1 ,1) &∗st ∗ct∗∗3−uppg (2 ,2 ,2 ,2)∗ct ∗st ∗∗3 upp (1 ,2 ,1 ,2)= uppg (1 ,1 ,1 ,1)∗ct ∗∗2∗st ∗∗2−uppg (1 ,1 ,2 ,2)∗ct ∗∗2∗st ∗∗2 36 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) &+uppg (1 ,2 ,1 ,2)∗ct ∗∗4−2∗uppg (1 ,2 ,1 ,2)∗ct ∗∗2∗st ∗∗2+uppg (1 ,2 ,1 ,2) &∗st ∗∗4−uppg (2 ,2 ,1 ,1)∗ct ∗∗2∗st ∗∗2+uppg (2 ,2 ,2 ,2)∗ct ∗∗2∗st ∗∗2 upp (1 ,2 ,1 ,3)= uppg (1 ,1 ,1 ,3)∗st ∗ct ∗∗2−uppg (1 ,2 ,2 ,3)∗st ∗ct ∗∗2 &+uppg (1 ,2 ,2 ,3)∗st ∗∗3−uppg (2 ,2 ,1 ,3)∗s t ∗ct ∗∗2 upp (1 ,2 ,2 ,2)= uppg (1 ,1 ,1 ,1)∗ct ∗s t ∗∗3+uppg (1 ,1 ,2 ,2)∗st ∗ct ∗∗3 &+2∗uppg (1 ,2 ,1 ,2)∗s t ∗ct ∗∗3−2∗uppg (1 ,2 ,1 ,2)∗ct ∗st ∗∗3−uppg (2 ,2 ,1 ,1) &∗ct ∗st ∗∗3−uppg (2 ,2 ,2 ,2)∗s t ∗ct ∗∗3 upp (1 ,2 ,2 ,3)= uppg (1 ,1 ,1 ,3)∗ct ∗s t ∗∗2+uppg (1 ,2 ,2 ,3)∗ct ∗∗3 &−uppg (1 ,2 ,2 ,3)∗ct ∗st ∗∗2−uppg (2 ,2 ,1 ,3)∗ct ∗st ∗∗2 upp (1 ,2 ,3 ,3)= uppg (1 ,1 ,3 ,3)∗ct ∗st−uppg (2 ,2 ,3 ,3)∗ct ∗st upp (1 ,3 ,1 ,1)= uppg (1 ,3 ,1 ,1)∗ct ∗∗3+3∗uppg (1 ,3 ,2 ,2)∗ct ∗st ∗∗2 upp (1 ,3 ,1 ,2)= uppg (1 ,3 ,1 ,1)∗st ∗ct ∗∗2−2∗uppg (1 ,3 ,2 ,2)∗st ∗ct ∗∗2 &+uppg (1 ,3 ,2 ,2)∗st ∗∗3 upp (1 ,3 ,1 ,3)= uppg (1 ,3 ,1 ,3)∗ct ∗∗2+uppg (2 ,3 ,2 ,3)∗s t ∗∗2 upp (1 ,3 ,2 ,2)= uppg (1 ,3 ,1 ,1)∗ct ∗s t ∗∗2+uppg (1 ,3 ,2 ,2)∗ct ∗∗3 &−2∗uppg (1 ,3 ,2 ,2)∗ct ∗st ∗∗2 upp (1 ,3 ,2 ,3)= uppg (1 ,3 ,1 ,3)∗ct ∗st−uppg (2 ,3 ,2 ,3)∗ct ∗st upp (1 ,3 ,3 ,3)= uppg (1 ,3 ,3 ,3)∗ct upp (2 ,2 ,1 ,1)= uppg (1 ,1 ,1 ,1)∗ct ∗∗2∗st ∗∗2+uppg (1 ,1 ,2 ,2)∗s t ∗∗4 &−4∗uppg (1 ,2 ,1 ,2)∗ct ∗∗2∗st ∗∗2+uppg (2 ,2 ,1 ,1)∗ct ∗∗4+uppg (2 ,2 ,2 ,2) &∗ct ∗∗2∗st ∗∗2 upp (2 ,2 ,1 ,2)= uppg (1 ,1 ,1 ,1)∗ct ∗s t ∗∗3−uppg (1 ,1 ,2 ,2)∗ct ∗s t ∗∗3 &+2∗uppg (1 ,2 ,1 ,2)∗s t ∗ct ∗∗3−2∗uppg (1 ,2 ,1 ,2)∗ct ∗st ∗∗3+uppg (2 ,2 ,1 ,1) &∗st ∗ct∗∗3−uppg (2 ,2 ,2 ,2)∗s t ∗ct ∗∗3 upp (2 ,2 ,1 ,3)= uppg (1 ,1 ,1 ,3)∗ct ∗s t ∗∗2−2∗uppg (1 ,2 ,2 ,3)∗ct ∗st ∗∗2 &+uppg (2 ,2 ,1 ,3)∗ct ∗∗3 upp (2 ,2 ,2 ,2)= uppg (1 ,1 ,1 ,1)∗st ∗∗4+uppg (1 ,1 ,2 ,2)∗ct ∗∗2∗st ∗∗2 &+4∗uppg (1 ,2 ,1 ,2)∗ct ∗∗2∗st ∗∗2+uppg (2 ,2 ,1 ,1)∗ct ∗∗2∗st ∗∗2 &+uppg (2 ,2 ,2 ,2)∗ct ∗∗4 upp (2 ,2 ,2 ,3)= uppg (1 ,1 ,1 ,3)∗st ∗∗3+2∗uppg (1 ,2 ,2 ,3)∗st ∗ct ∗∗2 &+uppg (2 ,2 ,1 ,3)∗st ∗ct ∗∗2 upp (2 ,2 ,3 ,3)= uppg (1 ,1 ,3 ,3)∗st ∗∗2+uppg (2 ,2 ,3 ,3)∗ct ∗∗2 37 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) upp (2 ,3 ,1 ,1)= uppg (1 ,3 ,1 ,1)∗st ∗ct ∗∗2−2∗uppg (1 ,3 ,2 ,2)∗st ∗ct ∗∗2 &+uppg (1 ,3 ,2 ,2)∗st ∗∗3 upp (2 ,3 ,1 ,2)= uppg (1 ,3 ,1 ,1)∗ct ∗s t ∗∗2+uppg (1 ,3 ,2 ,2)∗ct ∗∗3 &−2∗uppg (1 ,3 ,2 ,2)∗ct ∗st ∗∗2 upp (2 ,3 ,1 ,3)= uppg (1 ,3 ,1 ,3)∗ct ∗st−uppg (2 ,3 ,2 ,3)∗ct ∗st upp (2 ,3 ,2 ,2)= uppg (1 ,3 ,1 ,1)∗st ∗∗3+3∗uppg (1 ,3 ,2 ,2)∗st ∗ct ∗∗2 upp (2 ,3 ,2 ,3)= uppg (1 ,3 ,1 ,3)∗st ∗∗2+uppg (2 ,3 ,2 ,3)∗ct ∗∗2 upp (2 ,3 ,3 ,3)= uppg (1 ,3 ,3 ,3)∗st upp (3 ,3 ,1 ,1)= uppg (3 ,3 ,1 ,1)∗ct ∗∗2+uppg (3 ,3 ,2 ,2)∗s t ∗∗2 upp (3 ,3 ,1 ,2)= uppg (3 ,3 ,1 ,1)∗ct ∗st−uppg (3 ,3 ,2 ,2)∗ct ∗st upp (3 ,3 ,1 ,3)= uppg (3 ,3 ,1 ,3)∗ct upp (3 ,3 ,2 ,2)= uppg (3 ,3 ,1 ,1)∗st ∗∗2+uppg (3 ,3 ,2 ,2)∗ct ∗∗2 upp (3 ,3 ,2 ,3)= uppg (3 ,3 ,1 ,3)∗st upp (3 ,3 ,3 ,3)= uppg (3 ,3 ,3 ,3) c Symmetry upp (1 ,1 ,2 ,1)= upp (1 ,1 ,1 ,2); upp (1 ,1 ,3 ,1)= upp (1 ,1 ,1 ,3) upp (1 ,1 ,3 ,2)= upp (1 ,1 ,2 ,3); upp (1 ,2 ,2 ,1)= upp (1 ,2 ,1 ,2) upp (1 ,2 ,3 ,1)= upp (1 ,2 ,1 ,3); upp (1 ,2 ,3 ,2)= upp (1 ,2 ,2 ,3) upp (1 ,3 ,2 ,1)= upp (1 ,3 ,1 ,2); upp (1 ,3 ,3 ,1)= upp (1 ,3 ,1 ,3) upp (1 ,3 ,3 ,2)= upp (1 ,3 ,2 ,3); upp (2 ,2 ,2 ,1)= upp (2 ,2 ,1 ,2) upp (2 ,2 ,3 ,1)= upp (2 ,2 ,1 ,3); upp (2 ,2 ,3 ,2)= upp (2 ,2 ,2 ,3) upp (2 ,3 ,2 ,1)= upp (2 ,3 ,1 ,2); upp (2 ,3 ,3 ,1)= upp (2 ,3 ,1 ,3) upp (2 ,3 ,3 ,2)= upp (2 ,3 ,2 ,3); upp (3 ,3 ,2 ,1)= upp (3 ,3 ,1 ,2) upp (3 ,3 ,3 ,1)= upp (3 ,3 ,1 ,3); upp (3 ,3 ,3 ,2)= upp (3 ,3 ,2 ,3) c Reciprocity upp (3 ,1 ,1 ,1)= upp (1 ,3 ,1 ,1); upp (3 ,1 ,1 ,2)= upp (1 ,3 ,1 ,2) upp (3 ,1 ,1 ,3)= upp (1 ,3 ,1 ,3); upp (3 ,1 ,2 ,1)= upp (1 ,3 ,2 ,1) upp (3 ,1 ,2 ,2)= upp (1 ,3 ,2 ,2); upp (3 ,1 ,2 ,3)= upp (1 ,3 ,2 ,3) upp (3 ,1 ,3 ,1)= upp (1 ,3 ,3 ,1); upp (3 ,1 ,3 ,2)= upp (1 ,3 ,3 ,2) upp (3 ,1 ,3 ,3)= upp (1 ,3 ,3 ,3) upp (3 ,2 ,1 ,1)= upp (2 ,3 ,1 ,1); upp (3 ,2 ,1 ,2)= upp (2 ,3 ,1 ,2) upp (3 ,2 ,1 ,3)= upp (2 ,3 ,1 ,3); upp (3 ,2 ,2 ,1)= upp (2 ,3 ,2 ,1) upp (3 ,2 ,2 ,2)= upp (2 ,3 ,2 ,2); upp (3 ,2 ,2 ,3)= upp (2 ,3 ,2 ,3) upp (3 ,2 ,3 ,1)= upp (2 ,3 ,3 ,1); upp (3 ,2 ,3 ,2)= upp (2 ,3 ,3 ,2) upp (3 ,2 ,3 ,3)= upp (2 ,3 ,3 ,3) 38 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) upp (2 ,1 ,1 ,1)= upp (1 ,2 ,1 ,1); upp (2 ,1 ,1 ,2)= upp (1 ,2 ,1 ,2) upp (2 ,1 ,1 ,3)= upp (1 ,2 ,1 ,3); upp (2 ,1 ,2 ,1)= upp (1 ,2 ,2 ,1) upp (2 ,1 ,2 ,2)= upp (1 ,2 ,2 ,2); upp (2 ,1 ,2 ,3)= upp (1 ,2 ,2 ,3) upp (2 ,1 ,3 ,1)= upp (1 ,2 ,3 ,1); upp (2 ,1 ,3 ,2)= upp (1 ,2 ,3 ,2) upp (2 ,1 ,3 ,3)= upp (1 ,2 ,3 ,3) uppf =upp / ( 4 . d0∗pi∗r t ∗∗3) return end Acknowledgments This work was supported by the Consejer´ ıa de Innovaci´ on, Ciencia y Empresa de la Junta de Andaluc´ ıa, Projects of Excellence TEP-1207, TEP-2045 and TEP-4051, and by the Ministerio de Ciencia e Innovaci´ on, Spain, Research Projects TRA2006-08077 and MAT2009-14022. References [1] T´avara L, Ortiz J E, Mantiˇc V, Par´ıs F. Unique real-variable expressions of displacement and traction fundamentalsolutions coveringall transversely isotropic materials for 3D BEM. International Journal for Numerical Methods in Engineering;74:776– 798, 2008. [2] T´avara L, Mantiˇc V, Ortiz J E, Par´ıs F. Unique real-variable expressions of the integral kernels in the Somigliana stress identity covering all transversely isotropic elastic materials for 3D BEM. Computer Methods in Applied Mechanics and Engineering (submitted). [3] Fredholm I. Sur les ´equations de l’´equilibre d’um corps solide ´elastique. Acta Mathematica;23:1–42, 1900. [4] Lifshitz IM, Rozentsweig LN. Construction of the Green tensor for the fundamental equation of elasticity theory in the case of unbounded elastically anisotropic medium. Zhurnal Eksper. i Teoreticheskoi Fiziki;17:783–791, 1947. [5] Ting TCT. Anisotropic Elasticity Theory and Applications. Oxford University Press: Oxford, 1996. [6] Mal´en K. A unified six-dimensional treatment of elastic Green’s functions and dislocations. Physica Status Solidi B;44: 661–672, 1971. [7] Barnett DM. The precise evaluation of derivatives of the anisotropic elastic Green’s functions. Phys. Stat. Sol. (b);49:741748, 1972. 39 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) [8] Nakamura G, Tanuma K. A formula for the fundamental solution of anisotropic elasticity. Quarterly Journal of Mechanics and Applied Mathematics;50:179–194, 1997. [9] Ting TCT, Lee VG. The three-dimensional elastostatic Green’s function for general anisotropic linear elastic solids. Quarterly Journal of Mechanics and Applied Mathematics;50:407–426,1997. [10] Lee VG. Explicit expression of derivatives of elastic Green’s functions for general anisotropic materials. Mechanics Research Communications;30: 241-249, 2003. [11] Lee VG. Derivatives of the three-dimensional Greens functions for anisotropic materials. International Journal of Solids and Structuress;46: 3471-3479, 2009. [12] Tanuma K. Stroh Formalism and Rayleigh Waves Journal of Elasticity;89:1–3, 2007. [13] Shiah YC, Tan CL, Lee VG. Evaluation of explicit-form fundamental solutions for displacements and stresses in 3D anisotropic elastic solids. Computer Modeling in Engineering and Sciences;34:205–226,2008. [14] Tan CL, Shiah Y C, Lin CW. Stress analysis of 3D generally anisotropic elastic solids using the boundary element method. Computer Modeling in Engineering and Sciences;41:195–214,2009. [15] Shiah YC, Tan CL, Lee VG. Internal point solutions for displacements and stresses in 3D anisotropicelastic solids using the boundaryelement method.Computer Modeling in Engineering and Sciences;69:167–197, 2010. [16] Buroni FC, Ortiz JE, S´aez A. Multiple pole residue approach for 3D BEM analysis of mathematical degenerate and non-degenerate materials. International Journal for Numerical Methods in Engineering;86:1125–1143,2011. [17] Pan Y, Chou T. Point force solution for an infinite transversely isotropic solid. Journal of Applied Mechanics;98:608–612, 1976. [18] S´aez A, Ariza MP, Dom´ınguez J. Three-dimensional fracture analysis in tranversely isotropic solids. Engineering Analysis with Boundary Elements;20:287–298,1997. [19] Pan E, Yuan GF. Boundary element analysis of three-dimensional cracks in anisotropic solids. International Journal for Numerical Methods in Engineering; 48:211–237, 2000. [20] Loloi M. Boundary integral equation solution of three-dimensionalelastostatic problems in transversely isotropic solids using closed-form displacement fundamental solutions. International Journal for Numerical Methods in Engineering;48:823– 842, 2000. 40 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012) [21] Ding HJ, Liang J, Chen B. The unit point force solution for both isotropic and transversely isotropic media. Communications in Numerical Methods in Engineering; 13:95–102, 1997. [22] Yue ZQ. Elastic fields in two joined transversely isotropic solids due to concentrated forces. International Journal of Engineering Science;33:351–369,1995. [23] Ariza MP, Dom´ınguez J. Boundary element formulation for 3D transversely isotropic cracked bodies. International Journal for Numerical Methods in Engineering;60:719–753,2004. [24] Par´ıs, F, Ca˜nas, J., Boundary Element Method, Fundamentals and Applications. OUP: Oxford, 1997. [25] Aliabadi MH. The Boundary Element Method, Volume 2 Applications in Solids and Structures. John Wiley & Sons: Chichester, 2002. [26] Sutradhar A, Paulino G, Gray L. Symmetric Galerkin Boundary Element Method. Springer: Berlin, 2008. 41 V. Mantic et al./ Electronic Journal of Boundary Elements, Vol. 10, No. 1, pp. 1-41 (2012)