A least squares procedure for the evaluation of multiple generalized stress intensity factors at 2D multimaterial corners by BEM

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, Vol. 36, Issue 3, on March 2012, available at: https://doi.org/10.1016/j.enganabound.2011.09.011 Copyright 2011 Elsevier. En idUS Licencia Creative Commons CC BY-NC-ND 1 A least squares procedure for the evaluation of multiple generalized stress intensity factors at 2D multimaterial corners by BEM. A. Barroso, E. Graciani, V. Mantič, F. París Group of Elasticity and Strength of Materials, School of Engineering, University of Seville, Spain. Abstract. Detailed characterization of linear elastic stress states at corners and crack tips requires knowledge of the stress singularity orders, the characteristic angular functions and the generalized stress intensity factors (GSIF). Typically a high accuracy is found in the literature for the evaluation of the stress singularity orders and characteristic angular functions (numerically computed from analytical expressions in most cases). Nevertheless, GSIF values, evaluated by means of a numerical model using FEM or BEM and usually by postprocessing the results, are often reported with a lower level of confidence. A robust procedure is presented in this work for the evaluation of the GSIF at multimaterial corners. The procedure is based on a simple least squares technique involving stresses and/or displacements, computed by BEM, at the neighbourhood of the corner tip. A careful verification of the robustness and accuracy of the procedure using a few benchmark problems in the literature has been carried out. Applications of the procedure developed to the evaluation of GSIFs appearing at corners in metal-composite adhesive lap joints are presented. Keywords: Stress intensity factor, linear elastic anisotropic material, multimaterial wedge, singularity analysis. 1.- Introduction Problems having abrupt changes in geometry and/or material properties at some singular points solved under the assumptions of linear elasticity present unbounded stresses (referred to as stress singularities), see Wieghardt (1907), Williams (1952), Vasilopoulos (1988), Dempsey and Sinclair (1979, 1981), Sinclair (2004) and Paggi and Carpinteri (2008). Neighbourhood of such a point is usually referred to as corner, multimaterial corner if several materials meet at this point or also cross points if the singular point is located inside the domain. If the conditions for small scale yielding apply, the obtained singular elastic solution at a certain distance to the corner tip is representative of the real stress state. In this case, the linear elastic solution using a polar coordinate system (r,  ) centred at the corner tip, see Figure 1, admits, with the exception of some degenerate cases, a representation for the stresses  ij and displacements ui by the following asymptotic series expansion, see Kondratev (1967), Costabel and Dauge (1993), Knees and Sändig (2006) and Barsuk (2010) for a rigorous mathematical justification: 2 1 i N  0  N  i  i-1 x 1 x 2 1 i N  0  N  i  i-1 x 1 x 2  r Fig. 1.- Multimaterial corner. )(),( )(),( 1 11     m i M m mi m ij M m m ij grKru f r K r m m       , (i,j=r,  z) (1) in which  m are the singularity exponents, )(  m ij f and )(  m i g are respectively the characteristic angular functions for stresses and displacements, and Km, the weights of the terms, are the so-called Generalized Stress Intensity Factors (GSIF). Thus, the set of Km (m=1,...,M) defines the local elastic state at a corner. Notice that each term in (1) represents a solution of the governing partial differential equations in the corner domain. We will refer to each term as a mode similar to that used in the case of cracks. It is assumed that  m are naturally ordered fulfilling Re  m  Re  m+1. Logarithmic terms have not been considered in this work for the sake of brevity; see Sinclair (1999) for further information. However, the procedure for the evaluation of GSIFs presented here may be easily generalized to include these logarithmic terms. When a  is a complex number, as in the case of interface cracks,  =  R+i  I (where  R and  I are real numbers), the associated GSIF is also a complex number K=KR+iKI (where KR and KI are real numbers). In such a case, two terms can be included in (1). In the representation of the stresses, one term would be equal to   )(Re 1   ijR frK  and the other would be   )(Im 1   ijI frK , while in the representation of the displacements, one term would be equal to   )(Re   iR grK and the other   )(Im   iI grK . The characteristic angular functions )(  m ij f and )(  m i g together with the singularity exponents  m depend only on the local geometry, local material properties and the type of local boundary conditions. The GSIFs, m K, in addition depend on the global geometry and prescribed boundary conditions, their values being proportional to the magnitude of boundary conditions for linear elastic solutions. The rigid body motions are included in (1) for  m=0 (translations) and  m=1 (rotations) with the appropriate definition of m i g and the corresponding m ij f=0. Note that Re  m>0 for the other modes in (1). 3 The terms in (1) with Re  m>0, where  m =1-  m, thus 0<Re  m<1, are called singular, giving rise to unbounded stresses as r0+. In this case,  m are referred to as the stress singularity orders. Failure predictions at these corner-points admit several approaches, most of them being based either on allowable values of these GSIFs or on allowable values of field variables evaluated by means of the local stress or strain fields governed by the GSIFs. The evaluation of these GSIFs is then of crucial importance. For the evaluation of the GSIFs, the global geometry and far field loading must be considered. Thus, in general, numerical models by means of FEM or BEM or experimental tests (e.g. using photoelasticity) have to be used. Techniques for the evaluation of GSIFs can be roughly divided into four basic groups, according to the local/global character of the procedure and to the necessity, or otherwise, of post-processing of the FEM or BEM analysis results. Local techniques are sensitive to the accuracy of the numerical solution values for stresses and/or displacements close to the corner tip, while global techniques, working also, or only, with the elastic solution far from the corner tip, typically by making use of conservative integrals, are thus less sensitive to the solution accuracy at the corner tip. Regarding the second criterion for the classification, those techniques using postprocessing of basic field variables do not need to be incorporated into the numerical codes (FEM, BEM), but do not typically have as good accuracy as the methods which directly incorporate the singularity shape functions into the discretization, usually requiring a modification of the numerical code. Some references belonging to these groups are included in Table 1, which do not aim to be an exhaustive review. Further information can be found in Helsing and Jonsson (2002a), Sinclair (2004) and Paggi and Carpinteri (2008). Particular examples of GSIF evaluation involving materials having non-isotropic constitutive laws can be found in Quaresimin and Ricotta (2006), Nomura et al. (2009) using the H-integral and including thermal stresses, or the interaction integral by Cisilino and Ortiz (2005). Local techniques Global techniques Evaluation of GSIFs using postprocessing of numerical (FEM, BEM) results Based, e.g., on least squares adjustment (Chan et al. 1970, Munz and Yang, 1993), Barroso et al (2004), Liu et al (2008). Based, e.g., on conservative integrals, Banks-Sills (1997), Sinclair et al (1984), Omer & Yosibash (2005). Incorporation of singularity shape functions in the problem discretization Quarter point elements (Henshell & Shaw 1975, Barsoum 1976, Gray et al. 2003) and other singularity elements (Wait, 1978), Functions in the whole domain (Helsing and Jonsson, 2002a) and (Marin, Lesnic and Mantič, 2004) Table 1.- Classification of procedures for GSIF evaluation. Typically, high accuracy can be found in the literature regarding the evaluation of singularity exponents  m and characteristic functions )(  m ij f and )(  m i g in (1), see Wieghardt (1907), Williams (1952) and Vasilopoulos (1988) for single isotropic corners, Dempsey and Sinclair (1979, 1981) for multimaterial isotropic corners, Mantič 4 (1997) for single orthotropic corners, Pageau et al (1995) for single anisotropic corners, Mantič et al (2003) for multimaterial anisotropic antiplane corners, Ting (1997), Barroso et al (2003), Hwu et al. (2003) and Yin (2003) for multimaterial anisotropic corners, which can be computed by finding the roots of an analytical function. However, the accuracy in the evaluation of GSIFs is substantially worse as numerical models of the global elastic problem and postprocessing of the results are needed. Reliable and accurate results for benchmark problems which could be used as reference values in new evaluation methods are needed. Some interesting comments regarding the validity of the published numerical results in the literature can be found in Helsing and Jonsson (2002b). The approach proposed and explored in this paper is aimed at being an accurate and easy-to-use procedure for the evaluation of GSIFs. It can be located in the above classification in the first row (see Table 1), with no need to modify the FEM or BEM code applied, and between the two columns, as both near and far field data of the analysis results can be used. It should be noted that no special need for an accurate solution from the FEM or BEM analysis step in the neighbourhood of the corner tip is required, and it has a certain robustness with the use of the far-field data (far from the corner tip). The paper is divided into five main sections. Section 2 deals with the description of the basic features of the Boundary Element code (BEM) applied, París and Cañas (1997) and Graciani et al. (2005), used for the numerical models in this work. Section 3 describes the postprocessing procedure for the evaluation of GSIFs, which is based on a least squares method in terms of the displacements and/or stresses using the results obtained along the boundary edges of the corner and also the common edges between material wedges in the case of a multimaterial corner. The proposed postprocessing procedure for evaluation of GSIFs Km (m=1,...,M) is particularly well suited for numerical solutions arising from BEM models as the displacements and stress vector are the basic field variables used in the procedure. Nevertheless, no limitation appears in using the numerical solutions obtained by FEM models. Difficulties with the possible ill-conditioning of the resulting linear system are also analyzed in this section. Section 4 is devoted to the evaluation of the implemented procedure by means of the well known benchmark problems from the literature, performing parametric analyses to check the accuracy and robustness of the procedure. Section 5 applies the procedure for the evaluation of GSIFs in configurations of multimaterial corners appearing in adhesively bonded double-lap joints between aluminium and carbon fibre laminates. A final section summarizes the main features of the work. 2.- Main features of the multidomain BEM code The present study involves the numerical analysis of 2D plane elasticity problems including one or multiple materials, with isotropic or orthotropic behaviours, in which singular stresses are present. Two-dimensional BEM is employed in its traditional collocation formulation using continuous linear elements to mesh the boundaries of the problem in the displacement Somigliana identity. 5 The so-called displacement BIE can be obtained from the Somigliana identity by a wellknown asymptotic analysis, París and Cañas (1997), Mantič and París (1995), Mantič and París (2004), yielding:   KK KK j K ij KK j K ij K j K ij dtUduTuC   )()(),()()(),()()( yyyxyyyxxx (2) where y x j i , ,  ; x and y are, respectively, the collocation and integration points, both placed at the boundary K  of the solid K under study; K j u and K j t are, respectively, the displacement and traction components; K ij C is the coefficient tensor of the free term; finally K ij U and K ij T are, respectively, the components of the fundamental solution in displacements and tractions. Superscript K is employed to denote the possibility of analyzing problems having multiple domains with possibly different linear elastic material behaviours. The well-known analytic expressions of K ij C, K ij U and K ij T for isotropic materials can be found, for example, as in París and Cañas (1997). Simple explicit complex representations of K ij U and K ij T were introduced by Mantič and París (1995), along with a new analytical formula of the corresponding K ij C. Finally, explicit expressions of K ij U and K ij T for anisotropic mathematically non-degenerate and degenerate materials (i.e., with distinct and repeated roots of the sextic characteristic equation of the material) were presented in Mantič and París (2004). In problems involving multiple materials, perfect bonding along the interfaces is assumed in the present work. The equilibrium and compatibility conditions are imposed using a weak formulation, see Blázquez et al (1998), Graciani et al (2005) and Blázquez et al (2006). This formulation permits the use of non-conforming boundary meshes along bonded boundaries. In such a case, although equilibrium and compatibility are not strictly imposed at the nodes of the boundary mesh, a smooth numerical solution is obtained which accurately fulfils the equilibrium and compatibility along the boundary in an average sense. Conversely, if conforming boundary meshes are employed (as in the examples presented in this paper), the use of the weak formulation of the equilibrium and compatibility conditions yields a system of equations which is a linear combination of the strong formulation of the equilibrium ( 0 A i B itt ) and compatibility (A i B iuu ) conditions at the nodes of the mesh of the bonded boundaries. In some cases, depending on the local geometry and the boundary conditions, more than two unknowns may appear at certain nodes. In that case additional collocation equations are employed inside the elements adjacent to these nodes to ensure that a square matrix is obtained in the final system of equations, see París and Cañas (1997). If the boundary conditions permit the appearance of rigid body motions, the system of equations is appropriately modified by means of the Method F2 described in Blázquez et al (1996), to obtain a non-singular square matrix in the final system of equations. 3.- Evaluation of the GSIFs 6 In the present work, the only unknowns which will be considered in the local asymptotic series expansion (1) are the GSIFs Km (m=1,...,M). In the general procedure proposed and implemented by Barroso et al (2003), the analytic expressions of the characteristic angular functions are given in terms of the singularity exponents  m, which are computed as real or complex roots of an analytic function. A procedure based on a least squares method in terms of displacements and/or stresses will be described and studied in what follows to determine Km. The least squares method in displacements has been proved, Munz and Yang (1993), to be a simple and reliable method for the simultaneous evaluation of multiple GSIFs. In the present work, this approach has been extended to consider both displacements and/or stresses in order to evaluate the GSIF values. The present procedure, introduced by Barroso et al (2004) can also be considered a generalization of the work by Liu et al. (2008), as in the present proposal there is no limit to the number of stress singularities to be considered and these stress singularities, together with the characteristic angular functions, are obtained semi-analytically. The proposed least squares fitting procedure is based on the minimization of the difference between the numerical solution (obtained by BEM in this work) and the analytical one (1), given in terms of unknown values of Km (m=1,...,M). Hereinafter the BEM solution will be denoted by superscript "BEM" and the analytical one by superscript "series". The sum of squares of differences in values of the displacements, BEMseries uu  , and tractions, BEMseries tt  , at nodes placed along the radial edges of the corner are summed together. N  and Nr (see Figure 2) denote respectively the number of radial edges to be taken into account and the number of nodes taken at different distances along each radial edge. The so-called error function J, defined as       )0,(,...,,...,,..., 111  baKKbJKKaJKKJ MMuM  (3) is built up with the contribution of the following sums of squares of differences in displacements Ju and the stress vector J  :     )0(),(),...,,,(,..., 2 , 1 1 11           aruKKruaKKJ r N i N j ji BEM Mji series Mu r (4)     )0(),(),...,,,(,..., 2 , 1 1 11           brtKKrtbKKJ r N i N j ji BEM Mji series M r (5) The expression of J depends on the values of weights a and b, the usual choices being (a,b)=(0,1), (a,b)=(1,0) or (a,b)=(l2,  2) where l and  are some characteristic length and stress values which make stresses and displacements have values of the same order of magnitude. The dimensionless weights a  and b  allow the isolated components of the nodal displacements and tractions to be used. 7 j=1 j=2 j=3 j=4 j=5 j=6 N r j=1 j=2 j=3 j=4 j=5 j=6=N  N r N r N r N r N r Fig. 2.- Definition of parameters used in the least squares procedure. The use of BEM is particularly suitable for the evaluation of the error function J as the numerical results can be computed for the nodes defining the mesh. Nevertheless, no additional difficulty appears for the generalization of the error function J when working with an arbitrary set of nodes at any location inside a corner neighbourhood, a situation which will be useful in the case of having FEM results. A set of GSIFs Km (m=1,...,M) which minimizes the error function J in (3) is obtained by solving the following linear system of equations: 0 ),...,( 1   j M K KKJ (j=1,...,M) (6) The procedure admits complex values of the GSIFs, as detailed in Section 1. In such a case, two terms have to be included in (1) and correspondingly, two terms, associated to the real and imaginary parts of the GSIFs, have also to be taken into account in (6). Note that additional data could be obtained from the BEM results at some internal points in the domain. However, the idea herein is to carry out an additional check of the procedure by trying to represent the stresses and displacements inside the corner neighbourhood domain using only BEM results at the boundary mesh nodes in the fitting process. This will give confidence in the overall process, as no domain data have been used to fit domain variables. The resulting linear system arising from (6) will be calculated explicitly in the following. Only the displacement component ur will be considered in the evaluation of the error function J (i.e. ar=1, a  =0, a=1, b=0), for the sake of simplicity. The displacement ur at a point (usually a node) n defined by the radius ri and polar angle  j using M terms in the asymptotic series expansion representation (1) can be expressed by: m M m nmj m r M m imji series rKagrKru m   11 )(),(   (7) 8 where )( j m rinm gra m    is the coefficient associated to the Km-term of the series expansion of ur in (1) evaluated at the point n with the polar coordinates (ri,  j). The derivative with respect to Kj of the error function J built with the displacement component ur evaluated at N nodes (n=1,...,N) is a linear function of Kj (j=1,...M) which can be expressed as follows: 0)(2)(2 1 11 1                     N n M m BEM rnjnjmnm N n nj M m BEM rmnm j nuaaKaanuKa K J (8) where ),()( ji BEM r BEM rrunu   denotes the value of displacement ur at the point n, usually coinciding (as mentioned previously) with a boundary element node. Expression (8) can be written in matrix notation as: bAKAA ··· TT  (9) where A is a (N×M) matrix, N and M being respectively the number of points (nodes) taken for building the error function J and M the number of terms in the series expansion (the number of GSIF values) (1). As long as the number of nodes N is greater than the number of terms considered for the displacement representation M, A is a rectangular matrix with more rows than columns and thus expected to have full rank M. K (M×1) is the vector of unknowns (GSIF values) and b (N×1) is the vector of numerical results for ur at the chosen points.            NMN M aa aa    1 111 A,            M K K  1 K,            N b b  1 b (10) It is clear that the square (M×M) matrix AT·A can have rank=M if and only if the number of nodes N is equal to or greater than the number M of GSIFs to be calculated. Only in that case can the inverse of AT·A exist and be computed. Equation (9) is the typical matrix expression that appears when solving an overdetermined linear system bKA  · (11) using the 2-norm minimization 2 ·min bKA K, see Golub and van Loan (1996) or Trefethen and Bau (1997). Notice that the fulfilment of (11) corresponds to the vanishing differences in the error function J, i.e. J=0. The solution b of the full rank least squares problem is theoretically unique. Nevertheless, due to the nature of the matrix components )( j m rinm gra m   , the evaluation of the matrix AT·A using only nodes very close to the corner tip has been shown to give rise to ill-conditioned matrices with numerically computed rank(AT·A)<M. This ill-conditioning includes the cases in which the number of nodes N 15 c) d) 0 1 2 3 4 5 0 1 2 3 4 5 e) f) Fig. 5.- Errors in KI considering: a-b) M=3, c-d) M=5 and e-f) M=7 terms in the series expansion (1), a), c) and d) pictures including node distances and b), d) and f) including node numbers. With reference to the influence of the terms of the series expansion used, in general, the higher the number of terms (M), the lower the error found. Acceptable results (errors <1%) are found in most of the reasonable groups of nodes, avoiding small groups (NM) far away from the notch tip. In view of the error distributions in Figures 5a, c and e it seems that the variation of error for groups of nodes with approximately the same position of the central point (the arithmetic mean of the largest and smallest distances of nodes in the group, respectively, rmax and rmin, to the corner point), rc=(rmax+rmin)/2, is relatively small. 16 0 5 10 15 20 25 30 35 40 45 50 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 r c /L error (%) 3 terms 7 terms 5 terms 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 r c /L error (%) 3 terms 7 terms 5 terms Fig. 6.- Errors in KI considering M=3, 5 and 7 terms as a function of rc/L. This observation has motivated the representation in Figure 6, where the error is plotted for all groups of nodes considered in Figures 5 as a function of the central point distance to the corner tip rc. These plots confirm that rc/L is a very relevant parameter and that for small M, a suitable rc/L should be relatively small, while for larger M values rc/L can be larger. In particular, in the present problem the optimal values of rc/L are around 0.1, 0.2 and 0.30.4, respectively, for M=3, 5 and 7. In general, groups of nodes with rc/L>0.5 are not suitable for KI evaluation in this problem. 4.3.- The isotropic free-free bimaterial corner. The problem, depicted in Fig. 7, of an isotropic free-free bimaterial corner subjected to tensile loading has been taken from Qian and Akisanya (1999), a configuration which was analyzed later by Ortiz et al (2005) and Alcalá (2008). This problem is still simple enough to have a relatively simple analytic solution of the asymptotic field of stresses and displacements in terms of  m and Km.  h h h  1=90º  2=120º h  r Fig. 7.- Bimaterial free-free isotropic corner (Qian and Akisanya, 1999). 17 Qian and Akisanya reported results for the following Dundurs parameters:  =0.8 and  =0.2, which correspond to a system in which the upper material is about ten times stiffer than the lower material. The expressions for displacements and stresses by Qian and Akisanya (using their nomenclature) are:    N k m ijkk m ij frHr k 1 1)(),(   ,   N k m ikk m igrHru k 1 )(),(   (12) where ),,,,( 21 1   kkk ahH k  . (13) In (12-13) m is the material, Hk is the GSIF associated to the corresponding  k, ak is the dimensionless GSIFs (equal to Hk when h=1 mm and  =1 MPa), )(  m ijk f and )(  m ik g are known characteristic angular functions, see Qian and Akisanya (1999),  and  are the Dundurs parameters and  1 and  2 are the solid angles for both solids (90º and 120º respectively in the present case). The corner in Figure 7 has only one singular mode (0<  1<1), Qian and Akisanya reporting the first three singularity exponents:  1=0.6747,  2=1.1637 and  3=1.5938. These values are in relatively good agreement with those obtained using the analytical procedure (including a numerical computation of roots of an analytic function) by Barroso et al. (2003):  1=0.673473 (0.18% difference),  2=1.167477 (0.32% difference) and  3=1.589147 (0.29% difference). The GSIF values obtained by Qian and Akisanya using FEM and by means of a pathindependent contour integral (finally transformed into a domain integral) effectively show a reasonable path independency when evaluated at different distances from the corner tip: 0.0053h-0.0063h, 0.0217h-0.0255h and 0.0869h-0.11h as shown in Table 4. The average values and variation coefficients reported for the three domains are respectively a1=0.6301 (VC=0.009%), a2=-0.3671 (VC=0.129%) and a3=0.5443 (VC=0.798%). r=0.0053h-0.0063h r=0.0217h-0.0255h r=0.0869h-0.11h average a1 0.6301 0.6301 0.6300 0.6301 a2 -0.3666 -0.3673 -0.3675 -0.3671 a3 0.5407 0.5430 0.5491 0.5443 Table. 4.- Results for a1, a2 and a3, by Qian and Akisanya (1999). This problem has also been analyzed, using other techniques, by Ortiz et al (2005) and Alcalá (2008). A summary of their results is presented prior to the present results. The idea is to compare all of them, as each one uses a different approach to evaluate the GSIFs. Ortiz et al. (2005) reported for the same configuration, and using a domain-independent integral in conjunction with a 3D BEM model considering plane strain boundary 18 conditions, the values for a1 and a2 shown in Table 5 (where z/t is the relative depth and r/h the radial distance to the notch). The average values and variation coefficients are 1 a=0.67829 (VC=0.16%) and 2 a=-1.54431 (VC=1.44%). z/t=0.0 z/t=0.5 z/t=1.0 average a1 r/h=0.2 0.67700 0.67697 0.67669 0.67829 r/h=0.3 0.67877 0.67876 0.67853 r/h=0.4 0.67938 0.67938 0.67918 a2 r/h=0.2 -1.52424 -1.51860 -1.50673 -1.54431 r/h=0.3 -1.55476 -1.55223 -1.54716 r/h=0.4 -1.56654 -1.56558 -1.56294 Table. 5.- Results for a1 and a2 by Ortiz et al. (2005). Alcalá (2008) also analyzed the same problem, using the H-integral, see Bueckner (1973), with a 2D BEM model, the values for ai (i=1,2,3) being evaluated at the same distances from the corner tip as used by Qian and Akisanya. The results, in Table 6, show average values for a1=0.6788 (VC=0.075%), a2=-1.538 (VC=2.62%) and a3=0.5630 (VC=60.7%) or a3=0.7732 (VC=3.6%), eliminating the first two results for a3 (at r=0.0053h and r=0.0063h), which do not show the uniformity of the other four results. r=0.0053h r=0.0063h r=0.0217h r=0.0255h r=0.0869h r=0.11h average a1 0.6784 0.6781 0.6788 0.6788 0.6795 0.6792 0.6788 a2 -1.4746 -1.5006 -1.5585 -1.5581 -1.5686 -1.5683 -1.538 a3 -0.0174 0.3026 0.7553 0.7433 0.7942 0.7999 0.7732 Table. 6.- Results for a1, a2 and a3, by Alcalá (2008). Results from Qian and Akisanya, Ortiz and Alcalá will be compared with those obtained in the present work. With reference to the present method, the main characteristics of the BEM model used to evaluate the GSIF values are the following: linear continuous boundary elements, element size far from the corner tip is 0.025h, element size at the corner tip is 10-8h and a ratio of 1.5 between the length of consecutive elements is used for the progressive refinement along the radial edges in a neighbourhood of the corner tip. The edges converging to the corner have, respectively, 66 (  =90º), 66 (  =0º) and 72 (  =-120º) elements. The deformed and undeformed shapes are shown in Fig. 8, together with a detail of the discretization at the corner tip. 19 -1,00 -0,50 0,00 0,50 1,00 1,50 -0,75 -0,25 0,25 0,75 1,25 -0,20 -0,10 0,00 0,10 0,20 0,30 0,40 0,50 -0,20 -0,10 0,00 0,10 0,20 0,30 Deformada (x10000) Indeformada Deformed shape (x10 4 ) Undeformed Fig. 8.- Deformed and undeformed shapes of the bimaterial problem. Table 7 shows average values for (a1, a2, a3) and their corresponding variation coefficients (in %) using in the least squares procedure the displacement components ur and u  at the three edges emerging from the corner (  = -120º, 0º and 90º), thus N  =3, N=66 and M=3 in equation (3). The results in each row in Table 7 correspond to all groups of nodes including the range of distances indicated. Six different ranges have been chosen to observe the influence of nodes close to the corner tip. Unlike the values of a1, values of a2 and a3 show significant variations when using nodes close to the corner tip. As soon as these nodes are excluded from the fitting procedure, the results show a clear decrease in the value of the variation coefficient, which corresponds to the flattest part of the ai (i=1,2,3) plots. range (mm) a1 (VC%) a2 (VC%) a3 (VC%) 8.74·10 - 6 - 4.91·10 - 1 0.678908 (1.25%) -2.55 (>100%) -2.64 (>100%) 6.65·10 - 5 - 4.91·10 - 1 0.677977 (1.08%) -2.09 (54%) 1.50 (>100%) 5.05·10 - 4 - 4.91·10 - 1 0.676982 (0.92%) -1.86 (23%) 1.14 (>100%) 3.83·10 - 3 - 4.91·10 - 1 0.675973 (0.77%) -1.75979 (8.6%) 0.90 (38%) 2.91·10 - 2 - 3.28·10 - 1 0.675493 (0.58%) -1.72064 (4.2%) 0.82861 (8.5%) 1.64·10 - 1 - 2.91·10 - 1 0.675046 (0.43%) -1.71727 (3.9%) 0.81757 (4.5%) selected values 0.675046 -1.71727 0.81757 Table. 7.- Results for a1, a2 and a3. Results in Table 7 show the importance of choosing an appropriate averaging area in order to correctly evaluate the values of ai. As observed in Table 7, the number of nodes to avoid, located very close to the corner, is higher for the evaluation of a3 than for a2 and for a1 respectively, as the stress singularity order is lower. 20 The "selected values" row appearing in Table 7 has been chosen based on results from Section 4.2, Figure 6. It shows that the results with error below 1% can be found, irrespective of the number of terms included in the asymptotic series expansion, for rc/L ratios between 0.1 and 0.2, which is basically the choice in Table 7 (0.164-0.291). As a first observation, see Table 8, the results from Ortiz, Alcalá and those obtained in this work are much closer (at least in a2 and a3) than those obtained by Qian and Akisanya. In fact, Alcalá (2008), using the explicit expressions by Qian and Akisanya, obtains results which are much more similar to those obtained using the present method than to those by Qian and Akisanya (with great differences in a2 and a3). a1 a2 a3 Qian & Akisanya (1999) 0.6301 -0.3671 0.5443 Ortiz et al (2005) 0.67829 -1.54431 --- Alcalá (2008) 0.6788 -1.538 0.7732 Present work 0.675046 -1.71727 0.81757 Table. 8.- Summary of average results for a1, a2 and a3. To see the influence of using stresses and/or displacements in the procedure, data from the interface edge have been used (where both displacements and stresses are not null). A group of 25 consecutive nodes, Nr=25 in (4-5), at  =0º in the range (10-6h<r<0.02h) (nodes between number 10 and number 35, from a total of 66 along the edge) and M=3 has been used. Table 9 shows the results for a1 with different combinations of displacement and stress components, making use of the appropriate values for a, b in the evaluation of (3-5). components a1 ur 0.656855 (a=1, b=0) u  0.661111 (a=1, b=0) ur and u  0.662099 (a=1, b=0)  0.686260 (a=0, b=1)  r  0.634967 (a=0, b=1)   and  r  0.685272 (a=0, b=1) ur, u  ,    and  r  0.674215 (a=10 12 , b=10 - 8 ) ur, u  and  r  0.648634 (a=10 12 , b=10 - 4 ) ur, u  and   0.674305 (a=10 12 , b=10 - 8 ) Table. 9.- Results for a1 using different displacement and stress components. Taking the results of a1 from Table 7 (between 0.675046 and 0.678908) for all averaging areas, results from Table 8 show that the closest result is obtained when using all displacement and stress components (ur, u  ,  and  r  ) giving a value of a1=0.674215. To explore the capability of obtaining accurate results at a certain location using data from other locations, three different choices for the least squares fitting have been selected (Fig. 9) for the stress and displacement representation of the Qian and Akisania problem (Fig. 7) at points along a circumference with r = 0.1 h from the corner tip. The three completely different choices of nodes used in the procedure are: a) nodes along the 21 interface (Fig. 9a) in the range r ≤ 0.02 h, which represents 1/5 of the distance at which  ij and ui will be represented in Figures 10-12, b) nodes along the interface (Fig. 9b) in the range r ≤ 0.5 h, five times the distance at which  ij and ui will be represented, and c) nodes along the three edges of the problem (Fig. 9c), in the range r ≤ 0.02 h. All these three options use exactly the same BEM mesh, the only difference consisting in the segment of radial edges from which the BEM solution is used in the evaluation of ai (i=1,2,3). data for least squares: three edges at r<0.02 h at  =-120º, 0º and 90º representation at r=0.1 h data for least squares: interface r<0.02 h at  =0º representation at r=0.1 h a) 0.5 h data for least squares: interface at r<0.5 h at  =0º representation at r=0.1 h b) c) Fig. 9.- Data ranges and representations used for results in the isotropic bimaterial corner. The displacements (ur, u  ) and stresses (  r  ,  ) obtained by the present procedure using the three choices defined in Fig. 9 are shown, respectively, in Figs. 10-12. In these figures, the dots represent the BEM solution at some internal points uniformly distributed along a circumference of radius r=0.1h from the corner tip, while the continuous lines are the series expressions for displacements and stresses (1) using the first three terms in addition to the rigid body rotation term. Approximations of displacements and stresses using data from Fig. 9a and b, and shown in Fig. 10 and 11, are very good and similar in accuracy. When data from the three radial edges are included in the procedure (Fig. 9c) the comparison with the analytical solution is much better (Fig. 12) than in the two cases in which only data from one radial edge were used with data too close to and too far from the notch tip (rc0.01h and rc0.25h respectively, 0.1 and 2.5 times the distance from the radius r=0.1h at which  ij and ui are represented). The procedure has shown a very robust behaviour. u r u  mm angle MPa  r    angle   Fig. 10.- Displacements (ur, u  ) and stresses (  ,  r  ) at r=0.1h using displacement data from Fig.9a 22 u r u  mm angle MPa   angle  r   Fig. 11.- Displacements (ur, u  ) and stresses (  ,  r  ) at r=0.1h using displacement data from Fig.9b u r u  mm angle MPa   angle  r    Fig. 12.- Displacements (ur, u  ) and stresses (  ,  r  ) at r=0.1h using displacement data from Fig.9c The angular variation of stresses at r=0.1 h is shown in Fig.13, the dots being the BEM solution ),(  r BEM at internal points uniformly distributed along the circumference with radius r=0.1 h, the dashed line given by the ),(  r series solution for ai (i=1,2,3) by Qian and Akisanya (1999) and the continuous line given by the present results for ai (i=1,2,3), see Table 8. It can be seen that the stress component  rr is not continuous at the common interface (  =0º). -1.000 -0.500 0.000 0.500 1.000 1.500 2.000 2.500 3.000 -120 -90 -60 -30 0 30 60 90  MPa BEM Qian & Akisanya Present results rr      r Fig. 13.- Comparison of numerical and serial stresses in the bimaterial joint at r=0.1h. 23 The approximation ),(  r series obtained using the present approach for ai evaluation fits very well with the BEM solution ),(  r BEM , while larger differences can be observed with the solution by Qian & Akisanya. A difference of around 7% exists in the value of a1 between results by Qian and Akisanya and the ones obtained in the present work, while huge differences are found in a2 and a3 (see Table 8). These differences must nevertheless be carefully interpreted, as the characteristic angular functions for displacements (ur, u  ) and stresses (  ,  r  ) associated to terms 2 and 3 are quite similar, as can be observed in Fig. 14 (the modes have been conveniently multiplied for comparison purposes). -8 -6 -4 -2 0 2 4 6 8 -120 -90 -60 -30 0 30 60 90 ur (II) ur (III) u  (II) u  (III) -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 -120 -90 -60 -30 0 30 60 90   (II)   (III)  r  (II)  r  (III) (2)(2) (2)(2) (3)(3) (3)(3) (2)(2) (2)(2) (3) (3) a) b) Fig. 14.- Comparison of characteristic angular functions 2 and 3 a) displacements and b) stresses. Although, strictly speaking, the characteristic angular functions are linearly independent for different terms, if two of them are qualitatively similar in shape, as seems to be the case here (see Fig. 14), the comparison of particular values of GSIF for each term may have significant differences but the global contribution of both terms together can be quite similar. The values of a2 and a3 by Qian and Akisanya and those obtained in the present work are very different (Table 8), which may lead to the conclusion that one of the two results is wrong. Nevertheless, the complete stress representation in Fig. 13 does not differ as much as is indicated by the great difference in the obtained values of ai (i=1,2,3), giving some support to the above mentioned idea. This fact can make it difficult to separate these modes and evaluate correctly the values of a2 and a3. In any case, the values obtained with the present procedure, for the three ai (i=1,2,3), are in close agreement with the results obtained by Ortiz et al (2005) using a different BEM code and Alcalá (2008) using explicit expressions from Qian and Akisanya, giving the impression that results obtained in the present work with the proposed procedure are reliable and can be used as a reference in this benchmark problem. 5.- Application: Singularity analysis of composite-metal adhesive lap joints Once the present postprocessing procedure for evaluation of GSIFs has been verified in Section 4, it will be applied to a corner configuration which typically appears in structural adhesively bonded lap joints between metals and composites. In Fig. 15 a double lap joint [Aluminium (Al) (3.2 mm thickness) - Carbon Fiber Reinforced Plastic 24 (CFRP) [0º]12 (2.2 mm thickness)] is depicted. For the sake of clarity, the bimaterial corners have been conveniently magnified. metal adhesive metal adhesive 0 o adhesive 0 o adhesive adhesive 0 o adhesive 0 o metal adhesive metal adhesive CFRP [0] 12 (2.2 mm thickness) Aluminium (3.2 mm thickness) 12.5 mm ABCD Adhesive layer thickness = 0.1 mm Fig. 15. Geometry of the AL-CFRP double-lap joint. The mechanical properties of the materials are: Al (E=68.67 GPa,  =0.33), epoxy adhesive (E=3.0,  =0.35) and 0º CFRP (E11=141.3 GPa, E11=E22=9.58 GPa,  12=  13=0.3,  23=0.32, G12=G13= 5.0 GPa, G23=3.5 GPa). Taking, for example, the inner closed corner, marked in Fig. 15 with the letter B (90º wedge of the CFRP laminate and 270º wedge of the adhesive layer and fillet), where failure typically initiates, see Barroso et al. (2009a), and using the tool developed by the same authors in Barroso et al. (2003), the following singularity exponents were obtained:  1=0.763236,  2=0.889389,  3=1.10698. This last one does not correspond to a singular mode as  >1, see the stress representation in (1). Small scale yielding assumption has been verified to be valid for this case in Barroso et al (2009b). The search for the GSIFs has been attempted using a BEM model and employing linear elements with a progressive refinement of the mesh towards the corner tip. The size of the smallest element adjacent to the corner is 10-8 mm. A uniform tensile normal stress of 125 MPa is applied at the left boundary of the aluminium plate. A detail of the joint overlap zone of the BEM model and its deformed shape are shown in Fig. 16. Only a half-part of the problem has been modelled due to the existence of symmetry. 0.00 0.50 1.00 1.50 2.00 2.50 3.00 3.50 4.00 45.0 47.0 49.0 51.0 53.0 55.0 57.0 59.0 61.0 63.0 65.0 laminado de carbono [0º]12 Aluminio Deformada (x20) Deformed shape (x20) CFRP laminate [0] 12 Aluminium Fig. 16. BEM model (detail of the overlap zone).