Riveted joints in composites, a practical tool to estimate stresses around the rivet hole

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 Composite Structures, Vol. 263, on May 2021, available at: https://doi.org/10.1016/j.compstruct.2021.113735 Copyright 2021 Elsevier. En idUS Licencia Creative Commons CC BY-NC-ND Riveted joints in composites, a practical tool to estimate stresses around the rivet hole. J.A. Cabrera-González1, G. Vargas2, A. Barroso1* 1 Group of Elasticity and Strength of Materials. School of Engineering. University of Seville, Spain. 2 Department of Engineering. Angelo State University, Texas, USA. *Corresponding autor: [email protected] Abstract Riveted joints have different failure mechanisms, some of them associated to the failure of the rivet itself, and some others associated to the materials to be joined. The failures associated to the material are mainly due to the stress concentrations at the rivet holes. In joints with metallic (ductile) materials, these stress concentrations might not play a relevant role in the failure initiation, but with unidirectional composite (brittle) materials, they do control the failure. The orthotropic behavior of the materials, and the non-linear contact conditions between the rivet and the hole, makes the problem not to have a universal solution. Several Finite Element analyses have been carried out for elementary load cases, different materials and orientations of the fiber, and a twosteps least squares adjustment has been performed to obtain a polynomial expression of the stress distributions for all cases. As a result, a tool is provided, allowing the users to easily estimate the stress distribution (  rr,  ) along the hole contour, and  xy along the tangential lines to the hole circumference, which can be used in failure predictions. All data are available to readers by means of a spreadsheet with all necessary data to estimate the stress distributions. Keywords: Riveted joint, Composite Materials, Failure, Finite Element Analysis. 1. Introduction 1.1. Motivation of the work Materials joined by rivets need, first of all, to be drilled. This fact makes the own drilled material to be a potential location for failure initiation, due to the stress concentrations originated by the hole where the rivets will be installed. While for metallic (isotropic) materials, the shape of the stress field around the hole does not depend on the mechanical properties of the material, for unidirectional composite materials (orthotropic or anisotropic) this stress field shape strongly depends on the anisotropic behavior of the material (the orientation of the fiber). At least, three basic failure mechanism should be taken into account for the materials to be joined, due to the presence of the holes. The first one is the tension failure, typically appearing perpendicularly to the applied load and along the line of the reduced section of the joint where the holes are aligned (Figure 1a). A second failure to be taken into account is the shear failure, which appears when the holes are near the free end of the material (Figure 1b). A third failure is the bearing failure, originated by compression stresses in the contact between the rivet and the wall of the hole in the drilled material (Figure 1c). For sure, there are other failure mechanisms associated to the joining element (the rivet), which are not of interest in the present work. Figure 1. failure mechanisms For each one of these failure mechanisms (tension, shear and bearing) there is a clear stress component responsible for the failure onset (Figure 1), namely  for tension,  xy along the tangential lines of the expected failure path for shear, and  rr (in compression) for bearing. Let us take into account that each rivet row has two kinds of loads, the transfer load (the part of the load which passes to the other material through the rivet) and the by-pass load (which remains in the same material). If the rivet has interference with the hole, there are also self-equilibrated radial loads due to this interference (these loads will not be considered here). Neither load transferred by friction nor thermal loads will be considered in the present work. For isotropic materials, the maximum value of the stress component (normalized with respect to a certain remote nominal stress) and its location along the hole boundary, are typically well known. For example, for the by-pass load component, the maximum circumferential normal stress  is located at 90º with respect to the load direction, and its concentration factor with respect to the remote normal stress depends on the ratio w/d (distance between centers of holes / diameter of the hole) and it is typically slightly higher than 3 (  _max/  nom=3 for an infinite plate width, with a single hole w/d→∞). Nevertheless, this is unfortunately not the case for unidirectional composite laminates, where, both the location and stress concentration factor strongly depend on the fiber orientation and orthotropic mechanical properties.    rr  xy a) b) c) 1.2. Review of previous works To the best author’s knowledge, there are few publications which helps to completely estimate the shape and value of the stresses around the hole for different load cases. In particular, the EUROCOMP Design Code [1] (Chapter 5, pages 172-177) provides, for different elementary load cases: transfer load in tension (and compression), by-pass load in tension (and compression) and shear (with both signs), the values of  and  rr along the hole circumference, and  xy along the tangential lines to the hole, parallel to the load direction. Grüber et al. [2] employed analytical methods in anisotropic multilayered fibre-reinforced composites with pin-loaded holes developed on the basis of layer-related solutions; their analytical results were confirmed in numerical calculations. Savin [3] and Lekhnitskii [4] studied plates with holes of various shapes, mainly in isotropic media, with similar results but different approaches: Savin using an integro-differential method, and Lekhnitskii a series one. Regarding infinite anisotropic plates, other analytic solutions for determining stress and displacement fields of plates with holes were addressed by Hufenbach et al. [5], Daoust and Hoa [6], Ukadgaonker and coworkers [7,8,9], and Sharma [10]. Hufenbach et al. [5] considered several angles of loading and fiber orientations using the conformal mapping and the Schwarz formula to simplify the evaluation of stress functions in the solutions of Savin [3]. Daoust and Hoa [6] used a mapping function only on a unidirectional 0⁰ single layered plate. Ukadgaonker [7,8,9] obtained a solution for a biaxial loading problem based on Gao’s approach [11]. Other authors, like Hirashima [12] obtained the stress and displacement fields using the successive approximation by the point matching technique based on the exact solution of an anisotropic elastic body containing two or more elliptical openings under longitudinal shear applied at infinity. Tan and Gao [13] studied a quadratic isoparametric element formulation of the boundary element method (BEM) for the 2-D analysis of anisotropic bodies with stress concentrations. Hafiani and Dwyer [14] used the edge function method for anisotropic elasticity based on the complex variable formulation and on the superposition of analytical solutions to the field equations. More recently, by Zappalorto [15], obtained the solution for the near and far field in finite orthotropic plates with holes, but the solution for stresses was limited to the mid plane of the plate. Considering orthotropic plates, Wu and Mu [16] studied stress concentration in orthotropic and perforated cylinders comparing a finite element analysis with experimental results. Temiz et al. [17] also conducted finite element analysis of thick composite laminates with a hole under bending loads. Considering that woven fabric composites are orthotropic materials, Pandit et al. [18] and Toubal et al. [19] compared experimental tensile strain fields with approximate theoretical stress solutions on plain-woven fabric composites in the presence of holes; theoretical stresses were determined from existing models such as Lekhnitskii’s model [4]. Karakuzu et al. [20] analyzed experimentally the failure of woven layer composites with holes under pin-joint testing. With respect to composite plates, Ukadgaonker and Rao [7], Kaltakci et al. [22], and Rao et al. [21] studied the stress field around holes in symmetric laminates: the former analyzed triangular holes with blunt corners, the second one plates containing circular holes, and the last one, around rectangular and square cutouts. Other laminate sequences were studied by Rybicki and Schmueser [23] who analyzed stresses around circular holes in laminated plates using a 3-D finite element, with iso-parametric compatible displacement elements; in particular, they investigated the tangential strain distribution around circular holes for a [0/±45/0]s laminate, and the effect of changes in stacking sequence and lay-up angle on the interlaminar normal stress distribution around the free edge of the hole. Lin and Ko [24] used the Airy stress function in terms of full Laurent's series and a complex variable technique, in conjunction with the least-squares "boundary collocation method", to predict laminate strength for various hole sizes, layups, material properties, and loading conditions. Hundman and Horn [25] performed a stress analysis by three approaches (analytical, FEM and experimental) of a composite plate to examine the effects of bonding a plug in a centrally located hole. Kelly and Hallström [26] performed an experimental and numerical study on the effect of laminate stacking sequence and geometry on the bearing strength of fibre reinforced laminates manufactured from non-crimp fabric. Khechai et al. [27] examined cross-and-angle-ply laminated composite plates with single circular holes by the FE method. Similarly, Wu et al. [28] developed an analytical method to calculate stress around interference fit holes on composite pinned plates, and studied the effect of ply properties, ply stacking sequence and load level on stress distribution. In recent years, Pyl et al. [29] experimentally evaluated 3-D printed composite parts reinforced with continuous fibers to reduce stress concentrations around holes. As a summary, stress distributions reported in [1] are quite useful to predict the three potential failure mechanisms appearing in the drilled material. Nevertheless, the information in [1] has the following limitations:  In the case of unidirectional composite materials, results are limited to 0o and 90o in the fiber orientation. When a different orientation is considered, the mid-plane symmetry is lost due to the internal structure of the material and no results are available for these cases.  The shape and values for each stress component is given for composite materials with a volume fraction (vf) of 35 % (had lay-up) and 60 % (pre-pregs), without giving any information of the orthotropic mechanical properties of the materials. Stresses are considered to be equivalent in shape and value when vf is the same, irrespective of the values of the mechanical properties, fact which might be questionable. 1.3. Objective and structure of the work The main objective of the present work is to complete the useful information, as that given in [1], including detailed stress distributions for any fiber orientation (not only 0o and 90o) and for specific material systems (Glass Fiber Reinforced Polymers “GFRP”, Carbon Fiber Reinforced Polymers “CFRP”, and isotropic materials) as well as to provide a practical tool to have a good estimation of the three stress component values, without the need to develop any Finite Element Model (FEM). For that aim, a detailed FEM analysis, considering the nonlinear boundary condition given by the unilateral contact condition between the material and the rivet has been developed and solved for a number of particular cases. In Section 2, the numerical model will be detailed, including some particular discussions about the influence of the presence of the rivet in the obtained stress distributions. In Section 3, the results will be presented and conveniently organized to allow an interpolation procedure to be used for estimating stresses associated to fiber orientations not explicitly included in the analyses. In Section 4 some practical considerations and conclusions will be given. 2. Numerical model 2.1. Description of the model and load cases The FEM numerical model is very simple, but some details will be given in what follows to clearly set the framework of the analysis. Let us remember that, first of all, we need to obtain the reference curves depicted in [1] but without having the information of the material properties, or the explicit boundary conditions of the problems. The geometry is a representative volume of a riveted joint (for one rivet). We have used the same w/d ratio (distance between rivets / rivet diameter) equal to four, used in [1]. In the hole, unilateral frictionless contact has been defined using a circumference representing the boundary of the rivet. Neither clearance nor interference conditions have been considered between the rivet and the hole. As in [1], six elementary load cases have been analyzed, for which the following displacements and loads apply. Figure 2 shows the load cases under study. Figure 2. Elementary load cases considered for the analyses. Transfer load (in tension (Figure 2a) and compression (Figure 2b)): The part of the load passing to the other material in the joint, through the rivet, has one case in tension and another in compression. Due to the nonlinear contact condition between the rivet and the drilled material, both cases are completely different and have to be analyzed separately. The left (right) hand side of the model has a constant normal stress in tension (compression), while the right (left) hand side of the model is stress free. The upper and lower edges have been coupled in the FE model, by pairs of nodes (with the same horizontal, x coordinate), so all degrees of freedom are equal for each pair of coupled nodes. This condition makes the deformed shape of both edges to be equal, and by repetition of this representative volume (in the vertical direction) a multi-row joint can be modelled. No details are given in [1] for the applied boundary conditions at these edges. a) Transfer load (tension) c) By-pass load (tension) e) Shear load (+) b) Transfer load (compression) d) By-pass load (compression) f) Shear load (-) First polynomial adjustment of FEM results. The first step, once we have all the information from the numerical FEM solution, is to make a first least squares fitting of  rr and  along the hole perimeter (at r=R) as a function of the angular polar coordinate (  ) and of  xy as a function of the position along the tangential line to the hole (s) using Equation (1). We have to use the stress vector, the position vector (both obtained for the FEM analysis) and the degree of interpolation in the polyfit function. For each load case, and material, we have carried out this fitting procedure for a wide range of fiber orientations. This is essential to obtain enough information to perform the second interpolation in an accurate way. Figure 5 shows an example of the polynomial fitting for  rr for the transfer load (in tension) case and 0º for the fiber orientation. Both FEM results (in red) and polynomial fitting (in blue) are almost coinciden. The results are symmetrical, as the internal structure of the materials (fiber orientation of the GFRP = 0º) is also symmetric with respect to the horizontal mid-plane of the problem. The contact between the rivet and the composite material hole is lost at approximately ±84º. The stress is normalized by the bearing stress  b, previously introduced. Figure 5. First FEM polynomial fitting (  rr for the transfer load case with the fiber at 0º) Second polynomial adjustment of coefficients of the first polynomial fitting. Having the coefficients for the first least squares step, we must conduct a second one to obtain an expression as a function of the fiber orientation. To accomplish that goal, we have carried out another least squares fitting for each coefficient (of the first polynomial), load case, material, and stress component. For example, we show the following expression for the coefficient 𝑎. 𝑎= 𝑎𝑥+ 𝑎𝑥 + ⋯ + 𝑎𝑥 + 𝑎 (2) where x now represents the fiber orientation (  ). To obtain an accurate representation of each coefficient of the polynomial, we have included a wide range of fiber orientations. For each case, different FEM simulations for different fiber orientations have been carried out to reach a reasonable adjustment. For example, for the radial stress component  rr, the following angles have been evaluated in the FEM simulation: [0, 4, 7, 15, 22, 30, 37, 45, 52, 60, 75, 80, 85, 90]. With the first interpolation for each fiber orientation, we have vectors of coefficients 𝑎 related to each angle. Later, using the polyfit function we can obtain an approximation of these coefficients depending on the fiber orientation. Thanks to that, we have the different coefficients to obtain the stresses for the desired angle. Figure 6 shows (blue circles) the result for coefficient a2 of the polynomial fit of  rr at the previously mentioned 14 fiber orientations, and the red curve is the polynomial fit of these set of points. Figure 6. Polynomial fit for coefficient a2 of the  rr first polynomial fit at different fiber orientations. 3.3. Verification of the model with a fiber orientation not included in the adjustment. One important point is the verification of the adjustments carried out. We must check if the stresses obtained, for a fiber angle not included previously, are accurate. For that reason, we have displayed the stress component  rr at the hole contour (r=R) obtained with FEM and also with the obtained polynomial function for a fiber orientation not used in the least squares fitting (  =56º). As an example, we display in Figure 7 the radial stress (  rr), normalized with the bearing stress (  b), for the Transfer Load (in tension) case for a GFRP oriented at  =56o. Figure 7. Comparison of  rr for FEM results and polynomial fit for a particular case (  =56º). We observe that the approximation is very good for the hole angular range. Loss of symmetry is also visible due to the fiber orientation, information which is not available in [1] which only considers symmetric cases (fiber orientations at 0º and 90º) For the sake of simplicity, we only show this case, but a complete verification campaign has been done for the other load cases, materials and stress components, observing the same excellent agreement between the FEM solution and the polynomial one. 3.4. Presentation of other results. With the numerical simulation carried out, not only the stress profiles have been obtained, a wide variety of information using the calculated curves with the coefficients has also been obtained. First of all, it is possible to represent in the same plot, the stresses for different fiber orientations but for the same material and load case (see Figure 8). We can observe that the stress values are highly influenced by the fiber orientation. For example, for a material oriented at 0o the maximum stress is reached at 0o, while, for a fiber oriented at 90o, the maximum value is near 65o. In Figure 8, the stress  rr in normalized by the bearing stress  b. Figure 8. Radial stress  rr for different fiber orientations (Load Transfer and GFRP). Once we have an accurate polynomial fit approximation, a really useful information can be obtained from the analysis, as for example the value of the maximum stress and the location where it appears for each fiber orientation. Thus, taking into account this two magnitudes, we can estimate which is the most critical fiber orientation. This information is easily obtained through the polynomial approximations of the stresses, which represents one of the main achievements of the work. In the studied case (see Figure 9), the maximum stress appears at  =60o for the carbon fiber. However, for the glass fiber it is reached near 50o. Regarding the location where the maximum stress appears, both materials have a similar behavior. They have a trend equal to the fiber orientation angle except near 90o where it is slightly lower for the CFRP material. Figure 9. Evolution of the maximum  rr stress and its position (  ) with the fiber angle (  ). As mentioned previously, another parameter influenced by the fiber orientation is the loss of contact between the hole and the rivet (see Figure 10), which can also be easily evaluated thanks to the polynomial approximation of the radial stresses (computed at the point where  rr=0). Figure 10. Loss of contact angles vs fiber orientation of the material, for the Load Transfer case. It is important to highlight that all this information is available for each load case, material, stress component and fiber orientation (the latter, for the cases of composite materials) thanks to the obtained coefficients, not being necessary to perform any numerical analyses for that purpose. As mentioned previously, the tool is available in a spreadsheet in which the user has only to select: material, load case, stress component and, in the case of the composite materials, the fiber orientation. Figure 11 shows a general view of the spreadsheet, with the detailed information of all coefficients (in the rows at the bottom), and the stress plot corresponding to the selected problem. Figure 11. Spreadsheet made available to users. Conclusions The contact problem between a rivet and the material being joined has been analyzed. The stresses (  ,  rr and  xy) responsible for the main failure mechanisms have been numerically computed for a wide variety of configurations, including six load cases, three materials, and different fiber orientations. With all these numerical results, a double polynomial fit procedure has been carried out. The first step consisting in obtaining a polynomial approximation of each stress component distribution around the hole (for a particular load case, material and fiber orientation). After computing the numerical solution for different orientations of the fiber angle, a second step is carried out by means of a polynomial fit of each one of the coefficients of the first polynomial approximation. This, allows having a closed stress solution to estimate the stress state around the hole for any fiber orientation. All the data of this two steps procedure have been conveniently organized (and is available) to be used for any reader in a spreadsheet, in which the user needs only to select, the load case, the material, the stress component and the fiber orientation (in the case of a composite material), to have a complete stress estimation around the hole. Other valuable information, such as the loss of contact, maximum stress values, or the location of these maximum stress values are also obtained from the polynomial fitted expressions of the stresses. Acknowledgements The work was supported by Junta de Andalucía (Consejería de Transformación Económica, Industria, Conocimiento y Universidades) and European Regional Development Fund (P18-FR-1928). Data Availability The raw data required to reproduce these findings are available to download from www.germus.es. The processed data required to reproduce these findings are available to download from www.germus.es. References [1] Structural design of polymer composites: Eurocomp design code and handbook. J.L. Clarke (Editor), E & FN Spon, Chapman & Hall, London, ISBN 0419194509. (1996). [2] Grüber B, Hufenbach W, Kroll L, Lepper M, Zhou B. Stress concentration analysis of fibrereinforced multilayered composites with pin-loaded holes. Composites Science and Technology, 67;2007:1439-1450. [3] Savin GN, Stress concentration around holes. New York-Pergamon Press, 1961. [4] Lekhnitskii SG, Theory of elasticity of an anisotropic elastic body. San Francisco- Holden-Day Inc, 1963. [5] Hufenbach W, Schaffer M, Herrmann. Calculation of the stress and displacement field of anisotropic plates with elliptical hole. Ingenieur Arch 1990;60:507–17. [6] Daoust J, Hoa SV, An analytical solution for anisotropic plates containing triangular holes, Compos. Struct. 19;1991:107–130. [7] Ukadgaonker VG, Rao DKN. Stress distribution around triangular holes in anisotropic plates. Composite Structures 45;1999:171-183. [8] Ukadgaonker VG, Rao DKN. A general solution for stresses around holes in symmetric laminates under in-plane loading. Composite Structures 49;2000:339-354. [9] Ukadgaonker VG, Kakhandki V. Stress analysis for an orthotropic plate with an irregular shaped hole for different in-plane loading conditions-part-I. Composite Structures 70;2005:255-274. [10] Sharma DS. Stress concentration around circular/elliptical/triangular cutouts in infinite composite plate. Eng. Lett. 20;2012:1–9. [11] Gao XL. A general solution of an infinite elastic plate with an elliptic hole under biaxial loading. International Journal of Pressure Vessel and Piping 67;1996:95–104. [12] Hirashima KI. Stresses and displacements around openings under longitudinal shear. PROC. OF JSCE 220;1973:131-141. [13] Tan CL, Gao YL. Boundary element analysis of plane anisotropic bodies with stress concentrations and cracks. Composite Structures 20;1992:17-28. [14] Hafiani F, Dwyer JF. Edge function analysis of anisotropic materials with holes and cracks. Computers and Structures 72;1999:779-791. [15] Zappalorto M. Universal equations for the mode I stress distribution in finite size orthotropic plates with blunt notches and holes. Theoretical and Applied Fracture Mechanics 109;2020:102768. [16] Wu HC, Mu B. On stress concentrations for isotropic/orthotropic plates and cylinders with a circular hole. Composites Part B: Engineering 34;2003:127-134. [17] Temiz S, Ozel A, Aydin MD. Applied Composite Materials 10;2003:103-117. [18] Pandit SD, Nishiyabu K, Verpoest I. Strain concentrations in woven fabric composites with holes. Compos Struct 2003;59:361–8. [19] Toubal L, Karama M, Lorrain B. Stress concentration in a circular hole in composite plate. Compos Struct 68;2005:31-36. [20] Karakuzu R, Gulem T, Murat I, Ten B. Failure analysis of woven laminated glass–vinylester composites with pin-loaded hole. Compos Struct 72;2006:27-32. [21] Nageswara Rao DK, Ramesh Babu M, Raja Narender Reddy K, Sunil D. Stress around square and rectangular cutouts in symmetric laminates. Composite Structures 92;2010:2845-2859. [22] Kaltakci MY, Arsalan HM. Stress concentrations of symmetrically laminated composite plates containing circular holes. Iranian Journal of Science and Technology, Transaction B: Engineering 30;2006:447-460. [23] Rybicki Edmund F, Schmueser D. Effect of stacking sequence and lay-up angle on free edge stresses around a hole in a laminated plate under tension. J Compos Mater 12;1978:300–313.