Singularity parameter determination in adhesively bonded lap joints for use in failure criteria

ACCEPTED MANUSCRIPT 1 Singularity parameter determination in adhesively bonded lap joints for use in failure criteria. A. Barroso, V. Manti, F. París. ([email protected], [email protected], [email protected]) Group of Elasticity and Strength of Materials, University of Seville, Spain. ABSTRACT A procedure for the singularity characterization of anisotropic multimaterial corners which typically appear in adhesively bonded lap joints between metals and composites is presented and implemented in the present work. The characterization in terms of characteristic exponents (stress singularities), characteristic functions and generalized stress intensity factors allows the definition of an oriented test program in order to analyze the suitability of a singularity based failure criterion for structures of this type. Keywords: A) Adhesive joints, C) Interfacial strength, C) Failure criterion, C) Stress singularity. 1. INTRODUCTION Adhesively bonded joints between metallic sheets and composite laminas generate critical points, where the Linear Theory of Elasticity predicts unbounded stresses. These singularity stress fields are due to abrupt changes in geometry and/or abrupt changes in the elastic properties between the different materials, Sinclair [1]. In the particular case of adhesively bonded single (or double) lap joints between an aluminium sheet and a CFRP laminate, these critical points can be easily identified in Fig. 1a. Fig 1 With a polar coordinate system (r, θ ) located at the multimaterial corner tip (Fig. 1b), where θ i-1 and θ i are the angles of the i-th material wedge (N wedges), the asymptotic singularity stress field in the vicinity of the corner can be assumed to have the following series expansion structure with separation of variables (some particular corner configurations, which will not be considered here, do not admit the separation of variables, see for example Joseph and Zhang [2]). ACCEPTED MANUSCRIPT 2 ,,,),( )2( ),( ,,,),( )2( ),( 11 11 θθ π θ θθ π θσ λ λ λ rjig rK ru rjif r K r ik n k k i ijk n k k ij k k k == ==   = − = − (1) where n is the number of terms considered in the series expansion, λ k are the characteristic exponents, leading to singular terms when 0< λ k<1 and regular terms if λ k1, fijk( θ ) and gik( θ ) are the characteristic functions: smooth functions inside each material wedge which contain the angular dependency of the stress and displacement components respectively associated to the singularity mode λ k and Kk are the Generalized Stress Intensity Factors (GSIFs) which represent the weight of each term in the asymptotic series expansion. The characterization of the singular stress field is fundamental if the initiation of failure (an initial crack or localized damage whose propagation implies the complete failure of the joint) can be assumed to occur due to that singular stress field. In the same way as Linear Elastic Fracture Mechanics (LEFM) accurately predicts the propagation of a crack when the SIF (Kk) reaches a certain allowable value (the fracture toughness of the material), it seems reasonable to consider the possibility of a GSIF based failure criterion (see Gradin & Groth [3], Groth [4], Hattori [5], Reedy [6], Reedy & Guess [7], Akisanya [8], Yang & Munz [9], Lefebvre & Dillard [10], Ishii et al [11], Wang & Rose [12], Dunn et al. [13], Penado [14], Lazzarin et al [15], Leguillon & Yosibash [16] and Quaresimin & Ricotta [17] among others). To analyze the suitability of such an approach using the concepts of the LEFM, it is necessary to have the tools to evaluate λ k, fijk( θ ), gik( θ ) and Kk. Only after this careful characterization can experimental tests be carried out in an oriented way, changing in a controlled manner the values of the singularity parameters, in order to observe their influence (if it occurs) in a failure initiation. In successive sections we will introduce the analytical procedure to determine the characteristic exponents and characteristic functions (Section 2), and the numerical models (in BEM) as well as the implemented procedure based on a least squares determination of the GSIFs (Section 3). The accuracy of the method is checked with some benchmark problems, and some examples regarding the characterization of adhesively bonded lap joints between metal and composites are presented (Section 4). With these tools, a preliminary experimental test program can be specifically designed to investigate the influence of the singularity parameters on failure initiation. 2. CHARACTERIZATION OF THE SINGULARITY MODES The characteristic exponents λ k and the characteristic functions fijk( θ ) and gik( θ ) only depend on material elastic properties, local geometry and local boundary conditions. ACCEPTED MANUSCRIPT 3 Considering homogeneous linear elastic materials, perfect adhesion between the material wedges, homogeneous boundary conditions (at external faces) and plane states (plane stress, plane strain or generalized plane strain), at the present moment there is enough information and efficient analytical and computational tools for the evaluation of accurate values of λ k, see Dempsey and Sinclair [18,19], Pageau et al. [20], Ting [21], Manti et al. [22] and Poonsawat et al. [23,24] among others. In this work a semianalytical tool, developed by the authors, Barroso et al. [25], is employed. Some of the main characteristics and improvements of the present tool in comparison with the previous works are explained briefly in the following. All linear elastic material behaviour laws (from isotropic to fully anisotropic) can be included simultaneously in the analysis. Depending on the approach, it has often been observed in previous works that only material behaviour laws of one kind (only isotropic, or only orthotropic) have been considered in the analysis. When considering anisotropic materials and using the powerful and elegant Lekhnitskii-Eshelby-Stroh formalism for anisotropic elasticity: Lekhnitskii [26], Eshelby et al. [27] and Stroh [28,29], or simply, Stroh formalism in what follows, these limitations were usually associated to the fact that isotropic materials are, in the framework of the Stroh formalism, mathematically degenerate materials and the formalism becomes much more cumbersome, see Ting and Hwu [30] and Wang and Ting [31]. A particular effort has been focused on obtaining explicit expressions of the Stroh matrices (A and B) for transversely isotropic materials, a material behaviour which can be assumed to represent the unidirectional long fibre composite laminas. These expressions have also been properly orthogonalized and normalized for the powerful orthogonality and closure relationships of the formalism to be fulfilled. These expressions, Barroso et al. [32], will help to develop further application of the Stroh formalism to transversely isotropic materials. It is of major importance to notice that mathematically degenerate cases can appear when considering transversely isotropic materials (see Tanuma [33]) irrespective of the value of the elastic constants of the material and only depending on the relative orientation of the material, a fact that has also been considered in the program implementation. All typical orthogonal homogeneous boundary conditions can be imposed at the external faces of the multimaterial corner (see Table 1), and closed corners (called also cross points, all materials bonded, no external faces) can also be analyzed. The matrices in Table 1 are defined in terms of the outward normal (n) and tangential vectors (sr,s3) (see Fig. 2) at the external faces ( θ 0, θ N) with 0 and I being respectively the null and identity matrices. Fig 2 Table 1 ACCEPTED MANUSCRIPT 4 Following an original idea by Ting [21], the concept of transfer matrix has been further applied to degenerate and extraordinary degenerate materials, so that independently of the number and type of materials in the corner, the characteristic equation from which the characteristic exponents are obtained is derived in a simple manner, always given by the vanishing condition of a 3×3 determinant in the case of open corners (or 6×6 determinant for closed corners). The characteristic equation written as: 0K =)( ˆ)2( λ N, (2) with ),()()()( )()()()( ˆ 0 )4( 0 )3( 0 )2( 0 )1()2( θθθθ θθθθ ϕϕϕ ϕ T uNN T NN T uNNu T NNuN DKDDKD DKDDKDK ++ ++= (3) where )( e au θ D and )( e a θ ϕ D are defined in Table 1 and )(i N Kare the submatrices of N K which is the transfer matrix of the entire corner and is obtained as the product of the transfer matrices of the material wedges EN (N being the number of material wedges in the corner): 121 ··...·· EEEEK − = NNN . (4) The structure of a transfer matrix Ei is detailed in Ting [21] for non-degenerate cases and in Barroso et al. [25] for all cases of mathematical degeneracy. The characteristic functions fijk( θ ), gik( θ ) associated to λ k can also be numerically evaluated, always bearing in mind that GSIFs have to be defined together with these characteristic functions for the stress field to be completely and coherently defined, see Yosibash & Szabo [34] and also Pageau et al. [35]. In short, the tool presented in this section permits an analytical knowledge of the characteristic exponents and functions, which will be considered to be known in the next section. 3. GENERALIZED STRES INTENSITY FACTORS In the evaluation of the GSIFs, the global geometry and far field loading are now necessary, thus numerical models or experimental tests (e.g. using photoelasticity) have to be used. Additionally, in most cases, a post processing of the data is necessary, the accuracy of the results being finally dependent on the “quality” of both processes. For example, due to the asymptotic character of the elastic representation in (1), GSIFs obtained from the solution far from the corner tip can not be representative of the real elastic solution, while GSIFs obtained extremely near the corner tip can be affected by numerical errors associated to the model discretization. ACCEPTED MANUSCRIPT 5 There exist techniques which overcome these difficulties by means of path-independent integrals evaluated in the neighbourhood of the corners, see for example Qian and Akisanya [36], Wu [37], Banks-Sills and Sherer [38] and Cisilino & Ortiz [39] among others. Further information on these different techniques can be found in Helsing and Jonsson [40] and Sinclair [1]. Unlike the typically high accuracy found in the literature regarding the evaluation of characteristic exponents, the accuracy in the evaluation of GSIFs is substantially worse. More reliable results for benchmark problems are needed which could be used as reference values in new evaluation methods. Some interesting comments regarding the validity of published numerical results in the literature can be found in Helsing and Jonsson [41]. In the present work, with the previous knowledge of accurate values of the characteristic exponents and characteristic functions, a simple procedure has been implemented based on a least squares method in terms of displacements and/or stresses. The numerical models have been performed using a BEM code developed by Graciani [42]. BEM models are more suitable (when compared with FEM models) due to their higher accuracy in solving problems involving singular stresses. Extending the approach developed by Munz and Yang [43] using FEM, the present work uses a least squares method to evaluate the GSIF values. An error function J (5a) is defined as summing squares of differences between the displacements, Ju (5b), and/or tractions, Jt (5c), obtained from the numerical model (uBEM and TBEM) and the corresponding values obtained from the asymptotic series expansion in (1) (useries and Tseries). The function J depends on the values of 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. ( ) )0,(,,..., 01 ≥ + = babJaJKKKJ tuqn (5a) [ ] 2 1 0 1 01 ),,...,,,(),(  = = = −= A B j M n qnjn series jn BEM uKKKruruJ α αα θθ (5b) [ ] 2 1 0 1 1),...,,,(),(  = = = −= A B j M n njn series jn BEM tKKrTrTJ α αα θθ (5c) The only unknowns in J are the GSIFs, Kk (k=1,...,n) and the rigid body translations K0q (if necessary). In the J expression three summatories are considered: Α represents the number of displacement components (A=1 for ur, u θ or uz only, A=2 for ur and u θ , ur and uz or u θ and uz and A=3 for ur, u θ and uz) B defines the number of edges used in the evaluation of the error (with N material wedges, there exist N+1 edges where the BEM model gives direct results, without internal points) and finally M indicates the number of nodes along an edge used in the evaluation of J. An approximation for Kk (and K0q) is obtained from the solution of the following linear equation system: ACCEPTED MANUSCRIPT 6 0 ),...,( 1= ∂ ∂ k n K KKJ , (k=1,...,n) 0 ),,...,( 0 01 = ∂ ∂ q qn K KKKJ , (q=x,y,z) (6) The authors have successfully checked the procedure, defined by (5) and (6), using displacements and stresses as well. In the next section the benchmark problems present corners with free-free boundary conditions, where only displacements have been used in the evaluation of the GSIFs, and, for comparison purposes, the rest of the problems analyzed have been evaluated using only displacements (a,b)=(1,0). 4. BENCHMARK PROBLEMS AND APPLICATIONS To analyze the accuracy of the implemented method, two simple problems have been selected, a 90º notch with tensile loadings and an isotropic free-free bimaterial corner. Although the method is independent of the complexity of the problem (number and type of materials, boundary conditions, etc.), it seemed reasonable to check it at least in these two problems. Nevertheless, the lack of benchmark problems in the literature to check with is evident. As mentioned before, the values of the GSIFs depend on of the definition of the characteristic functions (Yosibash and Szabo [34]), and many authors do not present both of them, the difficulty of such comparisons also being due to the lack of uniformity and information in the normalization/standardization of the presented results as indicated by Pageau et al [35]. 4.1 Single-edge 90º notched specimen subjected to tensile loading. This simple problem (Fig. 3) has been analyzed with a/w=0.5 and h/w=1, with a=5cm. The BEM model has 709 nodes and 709 linear elements, with a size along the outer edges of 0.1 cm. Along the edges converging at the corner the element adjacent to the tip has a size of 10-7 cm, growing with a factor of 1.5 until a size of 0.1 cm is reached, with a total of 100 elements along these faces. Fig 3 This symmetric load configuration, which generates a mode I (opening) in the neighbourhood of the corner, has been analyzed in detail by Helsing and Jonsson [40] from where the reference value of KI=4.295886967699 has been taken. Using ur and u θ (A=2) and both edges (N=2) for the evaluation of the error J in (5), Fig. 4 represents the absolute value of the relative error (up to 5%) between the obtained solution and the reference considering: a) three and b) seven terms in the asymptotic series expansion in (1). It is important to notice that horizontal axes do not represent distances to the notch tip but only node numbers, and although element length far from the tip is constant, element lengths near to the notch tip change gradually. ACCEPTED MANUSCRIPT 7 All possible combinations resulting from a consecutive group of nodes have been considered, the right hand side axis being the initial node number and the left hand side axis the final node number (the bottom half part of plots in Fig.4 then being meaningless). It is clear that major errors are associated to: • A small group of nodes far away from the notch tip (the peak at the left hand side of Fig. 4) where the more regular terms of the asymptotic expansion and also possibly the external boundary obviously affect the solution. • A small group of nodes very close to the notch tip (the peak at the right hand side of Fig. 4), where numerical errors associated to the discretization affect the solution. Fig 4 The more terms of the asymptotic expansion are considered in the analysis, the lower the errors associated to the far group of nodes are achieved. The previous consideration is clear when comparing Fig. 4a) and 4b). In both cases, it can be observed that, from an engineering point of view, for most of the “reasonable” group of nodes (avoiding small groups of nodes far away and very close to the notch tip) results are accurate enough (<3%). For the same configuration, and using a path-independent technique, Ortiz et al. [44], the errors obtained are always below 1%. 4.2 Isotropic bimaterial corner. The tools developed in Sections 2 and 3 have also been used to analyze the bimaterial corner in Fig. 5, (with h=1) from Qian and Akisanya [36]. Fig 5 The reference solution reported by Qian and Akisanya was obtained for Dundurs parameters given by α =0.8 and β =0.2 which indicate that the material in the top part of Fig. 5 is about ten times stiffer than the material in the bottom part. Using a FEM model and path-independent integrals, finally transformed into domain integrals, Qian and Akisanya obtained for the first three terms (only the first one being singular) the following values for λ k and ak (ak being a dimensionless Kk): λ 1=0.6747 (a1=0.6301), λ 2=1.1637 (a2=-0.3671) and λ 3=1.5938 (a2=0.5443). The BEM model used in the present work has linear elements with a size of 0.025h far away from the corner and 10-8h at the corner tip. The characteristic exponents obtained using the semi-analytical tool developed in Barroso et al. [25] were λ 1=0.673473, λ 2=1.167477, λ 3=1.589147, very similar to those obtained by Qian and Akisanya. The ACCEPTED MANUSCRIPT 8 evaluation of the GSIFs was also performed using (5) and (6) with ur and u θ (A=2), the three edges which converge at the corner (N=3) and an average value from the flattest part of reasonable groups of nodes (the central top part of the corresponding representation, similar to Fig 4b), avoiding far and close nodes to the corner tip. Considering the first three terms of the asymptotic representation, the obtained value was a1=0.673688, which differs 6.92% from the result by Qian and Akisanya. Ortiz et al. [44] for the same problem and a path-independent technique reported an average value of a1=0.67829, which differs 7.65% from Qian and Akisanya and 0.68% from the present results. The good agreement reported by Qian and Akisanya (using FEM) in the range 10-4h and 10-1h is improved by the present results, see Fig. 6 (for internal points at r=0.1h) with a BEM model. The range used with the least squares method in the present work is 5.818585·10-6h<r<0.6224868h, which leads to a closer agreement between numerical and analytical results. Fig 6 Once the tool has been demonstrated to be accurate in these particular problems, it will be used to characterize the multimaterial corners which typically appear in adhesively bonded joints between metals and composites. Note that only stresses (a,b)=(0,1) have also been used in a least squares expression (5), instead of only displacements (a,b)=(1,0), obtaining similar results with differences below 2%. In this particular case, as the corner faces are stress free boundaries, only the common interface (at 0º) could be used to evaluate Kk by the least squares using σ θθ . 4.3 Adhesively bonded joint between a metal and a composite. Different test configurations can be specifically designed to analyze the suitability of a singularity based failure initiation criterion by means of previous knowledge of the singularity parameters. For example, as GSIF dimensions depend on the associated characteristic exponents λ k, only GSIFs associated to the same value of λ k can be strictly compared. Accordingly, the following four different joint configurations (Fig. 7, where Lo is the overlap length and the numbers between brackets at the caption are the thicknesses in mm) with the same local corner configuration can be analyzed (only the overlap zone is detailed) in which the inner adherent is the metallic sheet and the first outer material of the composite adherent is a 0º layer: To show the capabilities of the previously introduced tools, the numerical analysis of configuration in Fig. 7a) is presented. The bonded double-lap joint is loaded in shear by means of a tensile force (see Fig. 8, where the units of the axes are the dimensions of the model in mm). The symmetry of the double-lap joint configuration is used to model ACCEPTED MANUSCRIPT 9 half of the geometry. Both the undeformed and the deformed geometry (x20) are shown in Fig. 8. Linear elements have been used and a progressive refinement towards the inner corner has been applied up to a final size of 10-8 mm at the corner tip, a total of 1484 nodes having been used in the BEM model. Fig 7 The effect of curing residual stresses due to the polymerization temperature of the adhesive (115ºC) and the different thermal expansion coefficients (CFRP, adhesive and aluminium) has also been taken into account in the BEM model. The thermomechanical properties of the three materials used in the analysis are respectively: CFRP (x1 being the fiber direction, E11=141.3 GPa, E22=E33=9.58 GPa ν 12= ν 13=0.3, ν 23=0.32, G12=5.0 GPa, α 1=-1·10-6 ºC-1, α 2= α 3=26·10-6 ºC-1), aluminium (E=68.67 GPa, ν =0.33, α =24.5·10-6 ºC-1), epoxy adhesive (3.0 GPa, ν =0.35, α =45·10-6 ºC-1). The applied tensile load is 400 N (200 N at each laminate at the right hand side). Fig 8 Using the semi-analytical tool previously introduced (Barroso et al. [25]) the characteristic exponents ( λ k) and functions (f( θ ) and g( θ )) have been obtained. The GSIFs associated to both mechanical and thermal loads have been evaluated using expressions (5) and (6) once the BEM model has been solved. Results for the inner closed corner (corner "b" in Fig. 8 with 90º wedge of [0º] CFRP and 270º wedge of adhesive) and the outer corner at the left hand side of the adhesive spew fillet (corner "a" in Fig. 8 with 180º aluminium wedge and 23º adhesive wedge) are shown in Table 2 (Kg associated to the characteristic exponent λ k=1 representing the rigid body rotation). Table 2 The values of Kk are standardized following Pageau et al. [35] in such a way that σ θθ | θ =0º=K/(2 π r)1-λ. In the particular case of the closed corner, and due to the low thickness of the adhesive layer (≈0.1 mm), the least squares adjustment has been performed between r=10-6 and r=0.03 mm (measured from the corner tip) to avoid the nearest nodes to the corner tip (between 10-8 and 10-6 mm) and those which could be far away from the asymptotic representation (>1/3 of the adhesive thickness). The results for r=0.0194 mm are plotted in Fig. 9, where circles are the BEM solution (at internal points) and the angle is measured from the horizontal interface between the CFRP laminate and the adhesive. Separate contributions of each term (K1, K2 and K3) are plotted separately, together with results considering two and three terms. Fig 9 ACCEPTED MANUSCRIPT 16 -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. 6- Comparison of stresses in the bimaterial problem. a) b) c) d) Fig. 7- a) Al (3.2)-[0º]8 (1.6) Lo=12.5mm, b) Al (3.2)-[0º]16 (2.9) Lo=12.5mm, c) Al (3.2)-[0º]12 (2.2) Lo=20mm, d) Al (3.2)-[0º2,90º2]s (2.2) Lo=12.5mm. Aluminium adhesive 23º Aluminium adhesive 23º a λ 0 o adhesive 0 o adhesive c λ =0.763236 Corner “a” Corner “b” aluminium (3.2 mm thickness) [0º] 8 CPFR laminate (1.6 mm thickness) adhesive (0.1 mm thickness) 50 mm 50 mm12.5 mm Fig. 8- Geometry (top) and detail of the BEM model (bottom). ACCEPTED MANUSCRIPT 17 0o adhesive 0o adhesive c λ =0.763236 Aluminium adhesive 23º Aluminium adhesive 23º a λ λ 1=0.763236 K1=-0.00427067 λ 1=0.887773 K1=0.00134598 λ 2=0.889389 K2=0.0387508 λ 2=1.680930 K2=0.00118731 λ 3=1.10698 K3=-0.0171819 λ =1 Kg=-0.00394706 λ =1 Kg=-0.000905341 Table 2.- Values of λ λλ λ and K for the analyzed corners. σ σσ σ θ θθ θ -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0 30 60 90 120 150 180 210 240 270 300 330 360 angle (degrees) MPa Mode 1 Mode 2 Mode 3 Mode 1 + Mode 2 Mode 1 + Mode 2 + Mode 3 BEM results Term 1 + Term 2 Term 1 + Term 2 + Term 3 Term 1 Term 2 Term 3 Fig. 9- σ σσ σ θθ θθθθ θθ in corner "b" (90º wedge [0º] CFRP - 270º wedge of adhesive). -15 -10 -5 0 5 10 15 20 -15 -10 -5 0 5 10 15 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1 0 1 2 -4 -3 -2 -1 0 1 2 3 4 -4 -2 0 2 4 Mode 1 σ θθ >0 σ θθ <0 Mode 2 σ θθ >0 Mode 3 σ θθ <0 -15 -10 -5 0 5 10 15 20 -15 -10 -5 0 5 10 15 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -2 -1 0 1 2 -4 -3 -2 -1 0 1 2 3 4 -4 -2 0 2 4 Mode 1 σ θθ >0 σ θθ <0 Mode 2 σ θθ >0 Mode 3 σ θθ <0 term 2 term 3 term 1 Fig. 10- σ σσ σ θθ θθθθ θθ for each term in the “90º [0º] CFRP - 270º adhesive” corner.