A fast and non-degenerate scheme for the evaluation of the 3D fundamental solution and its derivatives for fully anisotropic magneto-electro-elastic materials

Depósito de Investigación de la Universidad de Sevilla https://idus.us.es/ This is an Accepted Manuscript of an article published by Elsevier in Engineering Analysis with Boundary Elements, Volume 105, on August 2019, available at: https://doi.org/10.1016/j.enganabound.2019.04.010 Copyright 2019 Elsevier. En idUS Licencia Creative Commons CC BY-NC-ND Afastandnon-degeneratescheme for the evaluation of the 3D fundamental solution and its derivatives for fully anisotropic magneto-electro-elastic materials✩ Federico C. Buronia,∗,CristianoUbessi b,GabrielHattori c, Rogério J. Marczakb,AndrésSáez d aDepartment of Mechanical Engineering and Manufacturing, Universidad de Sevilla. Camino de los Descubrimientos s/n , Seville E-41092, Spain bDepartment of Mechanical Engineering - DEMEC/PROMEC, Universidade Federal do Rio Grande do Sul. Sarmento Leite 425, 90050-160, Porto Alegre,Brazil cDepartment of Engineering, University of Cambridge, CB2 1PZ, Cambridge, UK dDepartment of Continuum Mechanics and Structural Analysis, Universidad de Sevilla. Camino de los Descubrimientos s/n , Seville E-41092, Spain Abstract Anewexpressionforthefundamentalsolutionisintroduced,presenting three relevant characteristics: (i)itisexplicitintermsoftheStroh’seigenvalues, (ii)itremainswell-definedwhensomeStroh’seigenvaluesarerepeated, and (iii)itisexact. Afastandrobustnumericalschemefortheevaluation of the fundamental solution and its derivatives developed fromdoubleFourier series representations is presented. The Fourier series representation is possible due to the periodic nature of the solution. The attractiveness of this series solution is that the information of the material properties is contained only in the Fourier coefficients, while the information of the dependence of the evaluation point is contained in simple trigonometric functions. This implies that any order derivatives can be determined by spatial differentiation of the trigonometric functions. Moreover, Fourier coefficients need to be obtained only once for a given material, leading to an efficientmethodology. ✩Special issue on the Tribute of the Brazilian Computational Mechanics Community to Carlos Brebbia: Edited by Ney Augusto Dumont and José Telles ∗corresponding author Email address: [email protected] (Federico C. Buroni) Preprint submitted to EABE April 1, 2019 The robustness of the scheme arises from the properties (i)and(ii)ofthe new expression for the fundamental solution, which is used tocomputethe Fourier coefficients. The proposed approach combines the clean structure of the Stroh formalism with the simplicity of Fourier expansions, addressing the old drawbacks of anisotropic fundamental solutions. Key words: explicit expressions, mathematical degenerate materials,Stroh formalism, three-dimensional magnetoelectroelasticity,Green’sfunctions, Fourier series representation 1. Introduction In the last few years, the Magneto-Electro-Elastic (MEE) coupling presented in composites consisting of Piezoelectric (PE) and Piezomagnetic (PM) phases has been the focus of intensive research due to thenumerous emerging possibilities at multiple material scales. These coupling phenomena have also used for modern applications in various devices: sensors, actuators and smart structures are just some of the current applications. New possibilities in fabrication and design of advanced materials of general anisotropy and novel micro/nano devices and structures require robust and efficient numerical methods for the analysis of anisotropic multicoupling materials. Historic developments, challenges and perspectives of MEE materials and its applications can be found in recent reviews [1, 2]. Fundamental solutions can be used to solve sophisticated boundary value problems. They are a central subject in modelling the physical and mechanical behaviour of solids. For instance, they are required to obtain stresses due to internal defects in materials. They are also required for the analysis of general problems by some numerical methods, such as theBoundary Element Method (BEM) [3–5], to which Professor Carlos Brebbia has made extensive contributions [6]. In all these methodologies, fundamental solutions are typically evaluated thousands or millions of times. Advances in the use of these numerical methods are limited by the availability ofefficient and robust numerical schemes for the corresponding fundamentalsolutions.Additionally, these solutions require extensive computations, which can render the numerical method infeasible for large problems. All multifield materials have some degree of anisotropy. Hence, the development of anisotropic fundamental solutions with coupled behaviour has been built upon previous research for purely elastic anisotropic elasticity his2 torically. A recent comprehensive review of this development can be found in Muñoz-Reja et al. [7], so we focus our attention in the main drawbacks of existing three-dimensional (3D) anisotropic MEE solutions. The main interest in finding fundamental solutions in a suitable form for their numerical implementation has been the computational cost. For comparison purposes, the general anisotropic elastic case contains 21 different elastic constants, however the amount of calculations for the fundamental solutionisalreadylarge. For the MEE case, with 75 different material constants, the computational cost becomes a key issue. In general, fundamental solutions are available in the literature either in integral form or in terms of explicitexpressions. On one hand, fundamental solution and its derivatives may be presented in integral form (see the work by Han [8] and references therein). This approach is called implicit and involves numerical integration of regular functions. When the integral part of the fundamental solution is arranged in order to be dependent only on the direction of the evaluation vector but not on its modulus (see for instance deduction by Wang & Achenbach [9] for the elastic case), the smoothness of the kernels depends on the degree of anisotropy. In alternative integral forms, corresponding kernels may become less smooth and more difficult to integrate for small modulus oftheposition vector [10]. At any event, regarding the BEM implementation of the fundamental solutions, the numerical computation of the involvedintegralsmay lead in many cases to inefficient and computationally costly schemes. On the other hand, explicit expressions are more desirable, as provided by Pan [11] for the fundamental solution of MEE materials and by Buroni and Sáez [12] for its derivatives. In this approach, one needs to solve an eigenproblem for each evaluation point, and one deals with the typical problem of mathematical degeneracy which can arise in explicit formulations. Buroni and Sáez [12] have presented a set of explicit solutions, bothforthefundamental solution and its derivatives, which consider all possible kind of degeneracies – however one needs to know which solution to apply, and this is not a trivial task. In addition to the above-mentioned contributions, additional effort has been invested in recent years for constructing adequate explicit solutions. For instance, see the works of Xie et al. [13] and [14] for PE and MEE materials, respectively, where both solutions degenerate mathematically. It is remarkable the explicit approach recentlyproposedbyXieet al. [15] for PE materials, which is valid for degenerate cases. However, to the authors’ best knowledge no non-degenerate explicit fundamental solution has been derived for MEE materials. AsolutionconstructedbyRadon-Strohfor3 malism [10, 17] has been proposed and two different evaluationapproaches, complex form and real form, have been suggested and discussedinHsuet al. [16]. Between these two approaches, the one based upon thecomplex form also has the feature of single time evaluation, a conceptsimilartothe presented in this study. Moreover, some approximate approaches such as the data base and interpolation scheme proposed by Wilson & Cruse [18] have also beenproposed. See for instance the work by Muraishi [19] for PE materials. Another interesting approach based on the relationship between the Rayleigh expansion and Fourier representation, which leads to a series representation in terms of spherical harmonics has been proposed recently [20]. In the context of contact and crack problems see the work by Fabrikant [21] and for thermo-MEE fundamental solutions the work by Pasternak et al. [22]. In summary, existing fundamental solutions exhibit at leastoneofthefollowing inconveniences: (i)theyareimplicit,or(ii)theyaremathematically degenerated, preventing the use of such solutions in the general case, and/or (iii)theyareevaluatedbyanapproximatedprocedure.Inthiscontext, the aim of this work is to provide a numerical scheme suitable for numerical implementation of the 3D fundamental solution and its firstandsecond-order derivatives for PE, PM and MEE materials with general anisotropy. The remainder of the paper is organized as follows: In Section 2, basic equations and notation for magnetoelectroelasticity are introduced.Inordertoovercome the above mentioned limitations, in Section 3 a new expression for the fundamental solution in MEE materials is presented which hasthreerelevant characteristics: (i)itisexplicitintermsoftheStroh’seigenvalues,(ii)it remains well-defined when some Stroh’s eigenvalues are equal(mathematical degeneracy) or nearly equal (quasi-mathematical degeneracy), and (iii)itis exact. Next, by using this new expression and a representation of the solution based on double Fourier series, a very fast and robust numerical scheme for the evaluation of the 3D fundamental solution and its derivatives is presented in Section 4. Some numerical results validate the proposed approach while showcasing the attained accuracy. We discuss the conclusions in Section 5. 2. Basic equations of linear magnetoelectroelasticity Let xi(i=1,2,3)beaCartesiancoordinatesysteminthree-dimensions. The extended notation introduced by Barnett & Lothe [23] for PE materials is very convenient for the purpose of this work. In this way, the linear MEE 4 problem can be formulated in an elastic-like fashion by extending the elastic displacement field vector uiwith the addition of the electric potential ϕand the magnetic potential ϑas [24] uJ=⎧ ⎨ ⎩ ujJ⩽3 ϕJ=4 ϑJ=5, (1) and by defining an extended elasticity tensor with the following components [24] CiJKm = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ cijkm J, K ⩽3 emij J⩽3; K=4 eikm J=4; K⩽3 qmij J⩽3; K=5 qikm J=5; K⩽3 −λim J=4; K=5or J=5; K=4 −ϵim J, K =4 −µim J, K =5. (2) where cijkl,ϵil and µil denote the components of the elastic stiffness tensor at constant electric and magnetic fields, the dielectric permittivity tensor at constant strains and magnetic fields, and the magnetic permeabilities tensor at constant stresses and electric displacements, respectively; eijk,qijk and λil are the PE coupling coefficients at constant magnetic fields, PMcoupling coefficients at constant electric fields and ME coupling coefficients at constant strains and electric fields, respectively. We assume isothermal conditions. The material constants tensors show the following symmetry conditions cijkl =cjikl =cijlk =cklij,e kij =ekji,q kij =qkji,(3) ϵkl =ϵlk,λkl =λlk,µ kl =µlk. Due to these symmetries, CiJKm =CmKJi is satisfied. Moreover, the elastic constant, dielectric permittivity, magnetic permeabilityandMEcoupling tensors are positive definite, i.e. cijkmγijγkm >0,ϵijEiEj>0,µ ijHiHj>0,λimEiHj>0(4) ∀γkm,E j,H j∈R;γkm =γmk =0,E j=0,H j=0. 5 and no constraint is imposed on the PE and PM coupling tensors.Equation (4) is known as the strong convexity condition and it is equivalent to positive definiteness of the internal energy function. In the definitions above, the lowercase (elastic) and uppercase (extended) subscripts take values 1, 2, 3 and 1, 2, 3 (elastic), 4 (electric), 5 (magnetic), respectively. As pointed out by Fan [25] for piezoelectricity, these matrix representations are not tensors. So one has to be careful whenchangingcoordinates systems (as the ones used in Section 4). Then, usingtheintroduced matrix representation, elliptic equilibrium equations fortheelastic,electric and magnetic problems in terms of the extended displacementscanberecast in a similar way to Navier’s equation of elasticity as CiJKmuK,mi +fJ=0,(5) where fJis the extended body force vector, defined as fJ=⎧ ⎨ ⎩ fjJ⩽3 −feJ=4 −fmJ=5, (6) being fi,feand fmthe three components of body forces, the electric charge density and the electric current density, respectively. As usual, comma denotes differentiation. Note that uncoupled problems, i.e.,purelyelastic, electric and/or magnetic, can be considered by setting the corresponding coefficients eijk,qijk and/or λil to zero. 3. Non-degenerate fundamental solution The fundamental solution is defined as a two-point second-order tensor in a five-dimension space with components UKP such that satisfies the elliptic partial differential equations (5) where the generalisedbodyforcevector corresponds to a point load fJ=δJPδ(x−x′)being δ(x−x′)the Dirac delta function located at the source point x′and δJK the five-dimension Kronecker delta. In homogeneous media the fundamental solutions depends on the relative vector x−x′so, for simplicity it is considered that the Cartesian coordinate system has the origin at the source point x′,thusthefundamental solution is just a function of the evaluation point x.Foraphysicalinterpretation see [12]. 6 The fundamental solution can be expressed as a singular term by a modulation function Has UJK(x)= 1 4πrHJK(x)(7) where x=rˆe with r=|x|=0.ThemodulationfunctionHJK(x)depends on the direction of xbut not on its modulus, so HJK(x)=HJK(ˆe ).This function can be put in the context of the Stroh formalism [26] being known as one of the three extended Barnett-Lothe tensors, which is symmetric and H(ˆe )=H(−ˆe ).Hence,U(x)is also symmetric and even, i.e.: UJK(x)=UJK(−x).(8) Therefore, the following parity relationships for the derivatives of U(x) UJK,m(x)=−UJK,m(−x)(9) and UJK,mn(x)=UJK,mn(−x).(10) are satisfied. The tensor HJK can be evaluated as [12, 27] HJK(ˆe )=1 π%+∞ −∞ Γ−1 JK(p)dp, (11) with ΓJK(p)=QJK +(RJK +RKJ)p+TJKp2,(12) and QJK =CiJKmninm,R JK =CiJKmnimm,T JK =CiJKmmimm,(13) where niand miare the components of any two mutually orthogonal unit vectors such that (n,m,ˆe )is a right-handed triad. Note that QJK and TJK are symmetric like their elastic counterparts [26], but the MEE coupling cause the loss of positive definiteness of these matrices. However, as Ting [26] shows for piezoelectric materials, it can be proven thatQJK,TJK and ΓJK are non-singular in this case, so their inverses are guaranteed. Moreover, Hand Uare independent of the choice of the unit vectors mand non the oblique plane. 7 The kernel in equation (11) is a single-valued holomorphic function in the upper complex half-plane except at the five complex poles with positive imaginary part and their conjugates that corresponds to the roots of the ten-order polynomial equation |Γ(p)|=0.(14) The determinant in (14) can be factorised as |Γ(p)|=|T| 5 & ξ=1 (p−pξ)(p−¯pξ),(15) where pξare known as the Stroh’s eigenvalues and the bar over pξdenotes the complex conjugate and Tis defined in equation (13). Stroh’s eigenvalues can be obtained as the roots of tenth-order characteristic equation (14) as well as solving the eigen-problem, as described in Appendix A. Then, assuming that all Stroh’s eigenvalues are different, the integration in equation (11) can be done by the Cauchy’s residue theory to yield HJK(ˆe )= 2i |T| 5 ' α=1 ˆ ΓJK(pα) (pα−¯pα) 5 ( ξ=1 ξ=α (pα−pξ)(pα−¯pξ) ,(16) where ˆ ΓJK is the adjugate of ΓJK defined as ΓPJ(p)ˆ ΓJK(p)=|Γ(p)|δPK and i=√−1.Clearly,thisexpressionisnotvalidfordegeneratecaseswhen there are repeated Stroh’s eigenvalues. Following an idea of[27],equation (16) can be algebraically modified in order to obtain a well-defined solution, valid even for repeated Stroh’s eigenvalues of any multiplicity, as explained below. The adjugate matrix ˆ ΓJK(p)is a polynomial in pof degree eight. Let ˆ ΓJK(p)= 8 ' n=0 pnˆ Γ(n) JK,(17) where ˆ Γ(n) JK (n=0,...,8) are real symmetric matrices which only depend on the material properties and the position vector ˆe .These matrices can be computed in a straightforward way in terms of the Stroh matrices (13) 8 Then, points on the x3-axis such as ˆe (3) =(0,0,1)Tor ˆe (3) =(0,0,−1)T,in the new coordinate system x∗ i(i=1,2,3)theybecomeˆe ∗ (3) =(−1,0,0)Tor ˆe ∗ (3) =(1,0,0)T,respectively. Thisartificecanbeusedsincetheexpressions in equations (28) and (30) are well-defined for vectors like ˆe ∗ (3).Modulation functions ˜ UIKj and ˜ UIKjl are not tensors. They can be transformed according to the following rules ˜ UIKj(ˆe (3))=Ω(2)|π 2 AI Ω(2)|π 2 BK Ω(2)|π 2 cj ˜ U∗ ABc(ˆe ∗ (3)),(32) and ˜ UIKjl(ˆe (3))=Ω(2)|π 2 AI Ω(2)|π 2 BK Ω(2)|π 2 cj Ω(2)|π 2 dl ˜ U∗ ABcd(ˆe ∗ (3)).(33) Naturally, since Ω(2)|π 2is not necessarily a material symmetry operation, the MEE properties must be expressed in the x∗ i(i=1,2,3)systemas c∗ ijkl =Ω(2)|π 2 ia Ω(2)|π 2 jb Ω(2)|π 2 kc Ω(2)|π 2 ld cabcd,(34) e∗ ijk =Ω(2)|π 2 ia Ω(2)|π 2 jb Ω(2)|π 2 kc eabc,(35) q∗ ijk =Ω(2)|π 2 ia Ω(2)|π 2 jb Ω(2)|π 2 kc qabc,(36) ϵ∗ ij =Ω(2)|π 2 ia Ω(2)|π 2 jb ϵab,(37) λ∗ ij =Ω(2)|π 2 ia Ω(2)|π 2 jb λab,(38) µ∗ ij =Ω(2)|π 2 ia Ω(2)|π 2 jb µab,(39) in order to compute the solutions in the x∗ i(i=1,2,3)coordinatesystem. Table 2 compares the fundamental solution for fully MEE Material C (properties listed in table B.7) at point x=(1,1,1)Tobtained with the Fourier series approach, the non-degenerate approach and the solution by Buroni & Sáez [12] as reference. Table 3 presents the results obtained with the proposed formulation and its comparison with the results by Buroni & Sáez and finite difference approach for the UJK,3derivatives at point x=(1,1,1)T [12]. Table 4 presents the compararison for the UJK,12 derivatives at the same point with the finite difference approach by using the Buroni & Sáez solution [12]. For all cases, 256 Gauss points have been used for integration of the coefficients (21), and the series take 20 terms. Very goodagreement is observed amongst the solutions for all cases. 15 (JK)Buroni & Sáez (2010) present work (non degenerate) present work (Fourier Series ρ=2) 11 8.7225398002742391 ×10−48.7225310879278819 ×10−48.7225302692569269 ×10−4 12 6.4154620180000892 ×10−56.4154582776641815 ×10−56.4154583909996026 ×10−5 13 3.1446825949164037 ×10−43.1446793068922219 ×10−43.1446789927188727 ×10−4 14 −1.2680519416303164 ×10−3−1.2680506415118376 ×10−3−1.2680505211279109 ×10−3 15 −5.1231568254934725 ×10−7−5.1231513273905938 ×10−7−5.1231508865968785 ×10−7 22 7.6090426614139801 ×10−47.6090350460375387 ×10−47.6090344129206722 ×10−4 23 4.9374127465890965 ×10−54.9374095746782780 ×10−54.9374098591871851 ×10−5 24 4.5315107542192896 ×10−44.5315058718214098 ×10−44.5315053684158277 ×10−4 25 2.8343787079494569 ×10−72.8343751170138947 ×10−72.8343748212786792 ×10−7 33 1.0440331060742933 ×10−31.0440321207110419 ×10−31.0440320278586277 ×10−3 34 −1.8426611715379502 ×10−3−1.8426594309862083 ×10−3−1.8426592695127788 ×10−3 35 −6.5916283018091271 ×10−7−6.5916214947368479 ×10−7−6.5916209154751710 ×10−7 44 1.2157849166703620 ×10−21.2157837419870272 ×10−21.2157836354733913 ×10−2 45 1.8466901349619211 ×10−61.8466881829979327 ×10−61.8466880294716113 ×10−6 55 2.5182299673182738 ×10−72.5182273648575246 ×10−72.5182271523844513 ×10−7 Table 2: Green’s function UJK at evaluation point x=(1,1,1)Tfor MEE Material C. (JK)Finite difference Buroni & Sáez (2010) present work (Fourier Series ρ=2) 11 −3.675418036705983 ×10−13 −3.675418038147529 ×10−13 −3.6754151848859640 ×10−13 12 −1.022473809428751 ×10−13 −1.022473811350401 ×10−13 −1.0224716811402825 ×10−13 13 −1.130914219850401 ×10−14 −1.130914522684852 ×10−14 −1.1309163083534299 ×10−14 14 3.207066889615237 ×10−43.207067004117469 ×10−43.2070743668963542 ×10−4 15 2.240137639332812 ×10−72.240137681772090 ×10−72.2403482815477269 ×10−7 22 −3.482702552773726 ×10−13 −3.482702562831942 ×10−13 −3.4827117098878988 ×10−13 23 2.457210719924869 ×10−15 2.457208785075626 ×10−15 2.4578986477356111 ×10−15 24 −1.386199675011690 ×10−4−1.386199583437589 ×10−4−1.3862246603200998 ×10−4 25 −1.123708543004768 ×10−7−1.123708508292927 ×10−7−1.1236385421649195 ×10−7 33 −8.404303395022170 ×10−14 −8.404303328191928 ×10−14 −8.4043404516964220 ×10−14 34 1.966253560794704 ×10−41.966253482908411 ×10−41.9662681422586747 ×10−4 35 1.850871822783363 ×10−71.850871866296891 ×10−71.8506905085345953 ×10−7 44 −3.346553299576044 ×106−3.34655331275473 ×106−3.3465596382245759 ×106 45 −7.696119549791547 ×102−7.69611965386665 ×102−7.6960757074656776 ×102 55 −1.323349050039724 ×102−1.32334904512777 ×102−1.3233341260687896 ×102 Table 3: Derivative of Green’s function UJK,3at evaluation point x=(1,1,1)Tfor MEE Material C. In order to evaluate the performance of the double Fourier series approach, the following errors schemes are devised: e(K,α) int (λIJ) := α / m=−α α / n=−α000λ(m,n) IJ −λR(m,n) IJ 000 2α2,(40) where λR(m,n) IJ is a reference matrix coefficient evaluated with a high number of Gauss abscissas, and λ(m,n) IJ is the same coefficient being evaluated with a 16 (JK)Finite difference with Buroni & Sáez (2010) present work (Fourier Series ρ=2) 11 1.7612993919238762 ×10−41.7613069813436755 ×10−4 12 1.4372187047786456 ×10−41.4372186057878247 ×10−4 13 −2.4864503755635603 ×10−4−2.4864148494725331 ×10−4 14 4.3863676517483414 ×10−44.3861250202513018 ×10−4 15 −5.3289261534288785 ×10−4−5.3288733367096404 ×10−4 22 −3.8721191709814931 ×10−4−3.8719435901746122 ×10−4 23 −2.6527046054561354 ×10−4−2.6527176860501949 ×10−4 24 −6.1505992611293401 ×10−4−6.1506203927697429 ×10−4 25 3.6085826796324048 ×10−43.6084569437355527 ×10−4 33 −5.5435380481296344 ×10−8−5.5444351257591942 ×10−8 34 −1.9217495103716011 ×10−7−1.9217332626512755 ×10−7 35 −2.8879775065048592 ×10−7−2.8879064153501508 ×10−7 44 1.3028359742552897 ×10−71.3028391615667287 ×10−7 45 2.1695644783532599 ×10−82.1693234709859794 ×10−8 55 −1.8092540754710937 ×10−7−1.8092827923235988 ×10−7 Table 4: Second derivative of Green’s function UJK,12 at evaluation point x=(1,1,1)T for MEE Material C. Korder rule; and eS(UIJ) := 1 1UIJ −UR IJ1 1S ∥UR IJ∥S ,(41) where ∥·∥S=2S2|·|dω,S2denotes a unit sphere in R3. The error e(K,α) int defined in equation (40) is used in order to measure the accumulated error in the integration of the Fourier coefficients for a given number of terms α.Figures2-7showtheevolutionofthemeanerrorofall N×Ncomponents (N=3,4or5forelastic,PEorMEE)oftheFourier coefficients as 'e(K,α) int := N / I=1 N / J=1 e(K,α) int (λIJ) N2.(42) The reference solution is computed with 20 terms and 256 Gausspoints are used in order to obtain Fourier coefficients (21). Figure 2 refers to PE Material A with ρ=1.Infigure2(a)weshowinlogscalethemeanerror vs. number of Gauss points used to compute each of the Fourier coefficients (21). Each curve corresponds to a αnumber of terms included. The same data are shown in figure 2 (b), this time as a family of curves of constant number of KGauss points as a function of the number of terms α.It can be observed that if 64 Gauss points are used the error remains below 10−10 for any number of terms included in the series. Figure 3 shows the same evolution of the error for the same material but using the formulation with 17 Figure 2: Mean e(K,α) int for 16 components the fundamental solution for Material A. Formulation for ρ=1.(a)Errorvs.numberofGausspoints.(b)Errorvs.numberofterms in the series ρ=2.Comparisonwithpreviousfigure2leadtotheconclusionthatless Gauss points are needed for the same accuracy if the present formulation with ρ=2for the Fourier series approach is used. Figures 4-7 presents similar results for materials B and C. These numerical tests suggest that –for practical applications– 64 Gauss points are sufficient to evaluate accurately the Fourier coefficients (21), in agreement with [29]. In order to show the convergence behaviour of the Fourier series approach for ρ=1and ρ=2,figures8-10illustrateinlogscalethemeanvaluesof the error /eSfor Materials A, B and C, respectively. In this case the error is defined as 'eS:= N / J=1 J / I=1 eS(UIJ) 1 2N(N+1) .(43) For all cases, Fourier coefficients (21) have been computed with 128 Gauss points. As a reference for computing eSthe non-degenerate solution given by (7), (18) and (19) has been used. The integration on the unitsphereS2 has been performed with standard double Gaussian quadraturewith64×64 points. The three figures illustrate the fast convergence that shows the Fourier expansion by taking ρ=2when compared with convergence for ρ=1.Itisshownthathigheraccuracyisobtainedforagivencut-offof the series, or equivalently, for a given degree of accuracy less terms need to be included into the series by taking into account this simpledetailonthe periodicity of the Barnett-Lothe tensor. 18 Figure 3: Mean e(K,α) int for 16 components the fundamental solution for Material A. Formulation for ρ=2.(a)Errorvs.numberofGausspoints.(b)Errorvs.numberofterms in the series Figure 4: Mean e(K,α) int for 25 components the fundamental solution for Material B. Formulation for ρ=1.(a)Errorvs.numberofGausspoints.(b)Errorvs.numberofterms in the series 19 Figure 5: Mean e(K,α) int for 25 components the fundamental solution for Material B. Formulation for ρ=2.(a)Errorvs.numberofGausspoints.(b)Errorvs.numberofterms in the series Figure 6: Mean e(K,α) int for 25 components the fundamental solution for Material C. Formulation for ρ=1.(a)Errorvs.numberofGausspoints.(b)Errorvs.numberofterms in the series 20 Figure 7: Mean e(K,α) int for 25 components the fundamental solution for Material C. Formulation for ρ=2.(a)Errorvs.numberofGausspoints.(b)Errorvs.numberofterms in the series Figure 8: Mean values /eSof the 10 different components the fundamental solution for Material A. Formulation for ρ=1,2. 21 Figure 9: Mean values /eSof the 15 different components of the fundamental solution for Material B. Formulation for ρ=1,2. Figure 10: Mean values /eSof the 15 different components of the fundamental solution for Material C. Formulation for ρ=1,2. 22 5. Conclusions The 3D extended displacement fundamental solution and its firstand second-order derivatives for PE, PM and MEE materials have been obtained and its effective implementation further discussed in this paper. The new expression for the fundamental solution is (i)explicitintermsoftheStroh’s eigenvalues, (ii)itremainswell-definedwhensomeStroh’seigenvaluesare equal (mathematical degeneracy) or nearly equal (quasi-mathematical degeneracy), and (iii)itisexact. Werealisethatthemathematicaldegeneracy presented in previous explicit formulations can be removed by factorization of the denominator in eq. (16), since this comes from the mathematical structure of the solution and not from physical arguments. This solution is used as building block for the development of an alternative efficient approach based on double Fourier series representations. The Fourierseriesrepresentation is possible due to the periodic nature of the solution.Themainbenefit from this series solution is that the information of the material properties is contained only in the Fourier coefficients, while the information of the dependence of the evaluation point position is contained in simpletrigonometric functions. This results in two advantages: first, any order derivatives can be determined by simple spatial differentiation of the trigonometric functions. We present results for firstand second-order but higher-order derivatives can be built in a straightforward way using this methodology if required. Second, the Fourier coefficients need to be obtained only once for a given material, leading to a very efficient methodology for numerical implementations. We have shown that, exploiting the π-periodicity of variable φ,betteraccuracyis obtained for a given number of terms, i.e. convergence is improved (in some cases by several orders of magnitude). Fourier expansion representation is real-valued, which is an important feature for numerical applications. The robustness of the scheme arise from the fact that, due to the properties (i) and (ii)ofthenewexpressionforthedisplacementfundamentalsolution, the Fourier coefficients can be computed for any general anisotropic coupled material with any kind of mathematical degeneracy in the Stroh context. In summary, fundamental solutions for anisotropic materials deal with two main drawbacks, the mathematical degeneracy and its overly complex and computationally expensive structure. In this work we developed a scheme for the evaluation of 3D fully anisotropic fundamental solution for MEE materials merging the best of two worlds: the clean structure oftheStrohformalism along with the simplicity of Fourier expansions. These developments 23 are expected to help to mitigate the mentioned old drawback offundamental solutions. Acknowledgements The authors would like to dedicate this work to the memory of Prof. Carlos A. Brebbia (1938-2018), whose relevant research contributions helped to forge the Boundary Element Method and were key in its further development. This work was supported by the Ministerio de Economía, Industria y Competividad of Spain under project DPI2017-89162-R. R.J.M. would like to acknowledge CNPq, project no.310649, for the support. Appendix A. Computation of the Stroh’s eigenvalues In this work, the Stroh’s eigenvalues are obtained numerically by solving the linear eigen-problem [26]: )N1N2 N3NT 1*)a b*=p)a b*,(A.1) where N1=−T−1RT,N2=T−1,N3=RT−1RT−Q,(A.2) with Q,Rand Tdefined in (13) and the superscript Tdenoting transpose. In this implementation the GEEVX subroutine of LAPACK library has been used in order to compute the corresponding eigenvalues. Then, the five complex eigenvalues with positive imaginary part are the so-called Stroh’s eigenvalues. The remainder fives are their complex conjugates. Appendix B. Materials In this appendix we summarize the material properties used inthiswork. Materials A and B are transversely isotropic and therefore satisfied the following relations: c1212 =c1111 −c1122 2,c 1313 =c2323,c 2222 =c1111,c 2233 =c1133 e322 =e311,e 223 =e113,q 322 =q311,q 223 =q113 (B.1) λ22 =λ11,ϵ22 =ϵ11,µ 22 =µ11. Non-vanishing components for Materials A, B and C are presented in Tables B.5, B.6 and B.7, respectively. 24