Finite element modeling of mechanically fastened composite-aluminum joints in aircraft structures

Composite Structures 109 (2014) 198–210

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Finite element modeling of mechanically fastened composite-aluminum joints in aircraft structures Zlatan Kapidzˇic´ a,b,⇑, Larsgunnar Nilsson b, Hans Ansell a a b

Saab AB, SE-581 88 Linköping, Sweden Division of Solid Mechanics, Linköping University, SE-581 83 Linköping, Sweden

a r t i c l e

i n f o

Article history: Available online 11 November 2013 Keywords: Bolted joints Composite-aluminum Finite element modeling Hybrid wing structures

a b s t r a c t A three-dimensional, solid finite element model of a composite-aluminum single-lap bolted joint with a countersunk titanium fastener is developed. The model includes progressive damage behavior of the composite and a plasticity model for the metals. The response to static loading is compared to experimental results from the literature. It is shown that the model predicts the initiation and the development of the damage well, up to failure load. The model is used to evaluate the local force–displacement responses of a number of single-lap joints installed in a hybrid composite-aluminum wing-like structure. A structural model is made where the fasteners are represented by two-node connector elements which are assigned the force–displacement characteristics determined by local models. The behavior of the wing box is simulated for bending and twisting loads applied together with an increased temperature and the distribution of fastener forces and the progressive fastener failure is studied. It is shown that the fastener forces caused by the temperature difference are of significant magnitude and should be taken into account in the design of hybrid aircraft structures. It is concluded that, the account of the non-linear response of the joints results in a less conservative load distribution at ultimate failure load. Ó 2013 Elsevier Ltd. All rights reserved.

1. Introduction The proportion of fiber-reinforced composites in modern aircraft structures is constantly increasing due to their favorable lightweight material properties. At the same time, metals are still required in primary aircraft components e.g. when damage tolerance design is utilized. Consequently, the number of hybrid interfaces where composites and metals are joined, is growing. A common type of such interface is a shear loaded, bolted joint with a titanium bolt and carbon fiber-reinforced polymer (CFRP) and aluminum plates. This study deals with simulation of structures that contain a large number of such joints. Bolted joints are often weak parts of the structure. It is therefore important that they are properly designed from static, fatigue and damage tolerance points of view. From the weight point of view, the joints must not be conservatively over-dimensioned. The need for an effective and accurate analysis method is from these reasons apparent. Traditionally, the development and design of joints have been based on experimental testing of specimens, which has been shown to be both expensive and time-consuming. In recent years, the use of finite element (FE) methods to simulate the behavior of ⇑ Corresponding author at: Saab AB, SE-581 88 Linköping, Sweden. Tel.: +46 (0)13183884. E-mail address: [email protected] (Z. Kapidzˇic´). 0263-8223/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compstruct.2013.10.056

composite joints has increased. In order to replace the experiments, the analysis method should be able to accurately capture a number of issues, such as: three-dimensional state of stress and strain, material behavior, bolt pretension and clamping force, bolt hole clearance, friction, secondary bending effects, contact between the surfaces etc. Another matter, that further complicates the modeling, is that there is no obvious ready-to-use material model for the failure of composites. The complex nature of CFRP failure makes it very hard to derive a universal material model with damage initiation and propagation that works well in all possible situations. Nevertheless, a large number of different criteria for the initiation and development of damage have been proposed over the years. A review with the most common methodologies for modeling of failure can be found in [1]. Another review particularly concerning the failure of bolted composite joints is presented in [2]. Since that review, a vast number of analytical studies of CFRP fastened joints has been published. They range from two-dimensional models with rigid pins representing the bolt [3–5], to three-dimensional models with or without damage development [6–15]. In these studies, bolt pretension, bolt clearance, the effect of countersunk/protruding bolt heads, different failure initiation and damage progression models, by-pass loading and implicit/explicit solution methods have been studied. Also, a large number of experimental investigations has been done, for instance in [16,17]. The level of detail in FE

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210

modeling, in the studies mentioned, requires an extensive computer power in order for the computational time to be reasonable, even for small joint installations. In structures with a very large number of fasteners, the computational and modeling burden tends to become unmanageable if no simplifications are introduced. This has motivated several authors to develop different approaches suited for the incorporation of mechanically fastened joints into FE analysis on a structural level. In most of these methods, the fasteners are represented by structural finite elements. For example, Weyer et al. [18] used connector elements available in the program Abaqus to represent self-piercing rivets in crash analysis. The modeling included elastic, plastic and damage behavior of a single bolted metal sheet specimen. Gray and McCarhty [19] developed a method using beam elements connected to a rigid surface to represent the bolt, and shell elements to model the laminate plates of the joint. Bolt-hole clearance, bolt-torque and friction between the plates were considered, but the damage dissipation in laminates was not accounted for. The functionality of the model was demonstrated on a three-bolt single-lap joint. Another efficient approach was developed by Gray and McCarthy [20], were a user-defined element was implemented in Abaqus. The method is capable of modeling non-linear behavior and failure of composite joints based on semi-empirical approach. A joint with twenty fasteners was used for validation of the method. Ekh and Schön [21] made use of beam elements to represent both bolt and laminates and connector elements to account for bolt-hole clearance and friction. The model was suited for optimization of load transfer, but was limited to single-column joints. A two-dimensional model of a steel–concrete beam with spring elements representing the shear connections was developed by Wang and Chung [22]. A non-linear behavior of the springs, based on experimental data was utilized. The motivation for this work comes from the intrinsic problems of hybrid structures. If materials with different thermal and mechanical properties are mixed in structures, issues that are usually absent in homogeneous structures can arise. Some examples of material properties that can differ for composite laminates and aluminum are: thermal expansion coefficients, failure and fracture mechanisms, response to different types of loading, i.e. tensile versus compressive, fatigue accumulation and scatter, impact resistance, impact residual strength, degree of anisotropy, environmental sensitivity etc. Based on these differences, composite and aluminum materials used in aircraft structures are subject to different design and airworthiness requirements. The issues that arise in the design and certification of hybrid aircraft structures as a result of this mismatch are: thermally induced loads and deformations, multiple failure modes in joints, allowance of buckling and permanent deformations, determination of factors for statistical scatter of material properties, determination of significant load states etc. These issues have to be taken into account in the design and certification of hybrid aircraft structures. The recommended practice for certification of composite assemblies is to begin analyzing and testing at a small specimen level. This way, the risks in technology can be eliminated at an early stage before moving onto testing more complex and expensive structural parts. This is known as the Building Block Approach (BBA) [23]. When it comes to hybrid constitutions this approach may not always be appropriate. As a result of the material mismatch, a hybrid joint may behave in one way when tested as a specimen and in another way when installed in a structure. Based on this argument, the analyzes and tests should be done at a structural level first in order to evaluate the true behavior of a structure. Such testing can be very expensive and a lot of effort can be saved if these tests can be done virtually, i.e. in a computer simulation. In this study, the behavior of a large winglike hybrid box is examined by FE simulations.

199

Fatigue testing of this structure exposed to spectrum loading and an applied temperature, as well as static testing, will be performed in the near future and the results will be compared. The objective of this work is to develop a methodology for simulation of structures that contain hybrid CFRP-aluminum shear loaded bolted joints with damage and failure behavior. The study is divided into two parts. In the first part, the local behavior of a CFRP-aluminum single bolt joint is modeled using the FE code Abaqus [24]. A three-dimensional FE model, using solid elements, is developed. For the composite plate, a progressive damage model (PDM) is implemented and aluminum and titanium parts are modeled using an elasto-plastic material behavior. The resulting joint behavior is compared to results from the experiments performed by Ireman et al. [16] in order to verify the method. In the second part, the method from the first part is used to evaluate the local behavior of several joint configuration that are present in the wing box. These behaviors are then assigned to connector elements available in Abaqus and included into the model of the wing box. Mechanical and thermal loads are applied and the redistribution of fastener loads due to joint failure is examined. 2. Modeling of composite-aluminum bolted joints In this section, the mechanical behavior of a single-lap, shear loaded bolted joint specimen shown in Fig. 1 is assessed. The specimen was previously tested by Ireman et al. [16] where it was denoted AC6III. It contains one plate made of CFRP material HTA/ 6376 with ply thickness of 0.13 mm and a quasi-isotropic stacking sequence ½ð45=0=90Þ4 S . The other plate is made of aluminum AA7475-T76. A single 6 mm countersunk bolt, made of titanium Ti–6Al–4V STA is used with finger-tight pretension. The fitting tolerance between the bolt and hole is according to ISO f7/H10, which gives a clearance of 10–70 lm. The testing with an applied shear load was performed without any lateral support, which means that secondary bending effects were present in the joint. The load–deflection response was monitored at the loading grips and the relative displacement of the plates was obtained by extensiometers attached to the sides of the specimen. The specimen was loaded until final failure occurred in a bearing failure mode. Ireman [12] also performed an FE analysis of the specimen testing but did neither include the failure behavior of the composite nor the plasticity of the metals. 2.1. Finite element modeling In this work, FE modeling of the AC6III specimen is performed using eight-node brick elements with reduced integration, C3D8R in Abaqus. The whole mesh is shown in Fig. 2 and a closer, cutthrough view of the bolt and composite plate is shown in Fig. 3. As can be seen, the bolt, the washer and the nut are considered to be one solid piece. The aluminum plate is meshed with 8 elements in the thickness direction, and in the CFRP plate each ply is represented using one element. In both plates, the mesh is refined close to the contact area between the bolt and the plates. Some mesh convergence studies, based on stresses and strains in this area are performed, but the results are not presented in this paper. For stabilization of the initial loading procedure and to eliminate rigid body modes, weak spring elements are connected to the plates and the bolt in all three directions. Furthermore, a linear spring element is connected in series with the specimen to account for the flexibility of the testing machine, as was done in [25]. In order to eliminate hourglass modes in the reduced integrated elements, Abaqus assigns an hourglass stiffness based on the

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210

200

Fig. 1. Dimensions in mm of the single-bolted joint specimen AC6III tested in [16].

Clamped edge (u=v=w=0)

Composite

Aluminum

z y x

Applied displacement (v=w=0) Fig. 2. The FE model of single-bolted joint from [16].

Table 1 Elasto-plastic properties of AA7475-T76 and Ti–6Al–4V STA. Property

AA7475-T76

Ti–6Al–4V STA

Young’s modulus, E (MPa) Poisson’s ratio m

69,000 0.28

110,000 0.29

Yield stress,

ry (MPa) and

plastic strain,

Fig. 3. Cut-through view of the bolted joint mesh with the aluminum plate removed.

elasticity parameters of the element. However, when the material stiffness is exposed to significant reduction, as in this case, enhanced hourglass control is recommended to stabilize the zero energy modes. The artificial energy introduced with the enhanced hourglass control is checked to ensure that it is negligible and does not affect the accuracy of the solution. The contact between the surfaces in the joint is modeled using the general contact algorithm available in Abaqus. It considers all neighboring surfaces that might come into contact during the analysis and uses finite sliding with a surface-to-surface option. The penalty approach, with default stiffness is used to enforce the contact constraints. A Coulomb friction is assigned to all surfaces with a friction coefficient of 0.2. This value is experimentally estimated by Schön [26] and used by Ireman [12]. Initial clearance is

p

ry

p

ry

p

400 430 490

0 0.002 0.09

950 1034 1103

0 0.002 0.1

eliminated between the bolt and hole edge. The aluminum and titanium are modeled using an elasto-plastic material model with isotropic hardening. The assigned properties are shown in Table 1. For CFRP material, a phenomenological procedure, known as the ply-discount method (PDM), is implemented. A non-linear shear stress–strain relation is assumed. The initial elastic properties and unidirectional strengths are shown in Table 2 and details of the implementation are explained in the next section. As shown in Fig. 2, one free end of the specimen is clamped, while a displacement is incrementally applied to the other free end. The analysis is performed using large displacement formulation and automatic incrementation. 2.2. Progressive damage modeling of CFRP The ply-discount method is a progressive damage procedure that aims to represent the accumulation of damage in a CFRP

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210 Table 2 Lamina elastic and strength properties for HTA7/6376. Property Longitudinal modulus Transverse modulus Out-of-plane modulus In-plane shear modulus Out-of-plane shear modulus Out-of-plane shear modulus In-plane Poisson’s ratio Out-of-plane Poisson’s ratio Out-of-plane Poisson’s ratio Longitudinal tensile strength Longitudinal compressive strength Transverse tensile strength Transverse compressive strength Out-of-plane tensile strength Out-of-plane compressive strength In-plane shear strength Out-of-plane shear strength Out-of-plane shear strength

E11 E22 E33 G12 G13 G23

(GPa) (GPa) (GPa) (GPa) (GPa) (GPa)

m12 m13 m23 XT XC YT YC ZT ZC S12 S13 S23

(MPa) (MPa) (MPa) (MPa) (MPa) (MPa) (MPa) (MPa) (MPa)

141 10 11 5.2 5.2 3.9 0.3 0.5 0.5 2250 1400 65 300 30 344 120 80 80

laminate. In a typical procedure, the load is increased until a certain failure criterion is fulfilled at a material point. Depending on the activated failure mode, specific stiffness parameters are reduced at that material point. In this way, a decrease of load-bearing capacity is achieved in a direction consistent with the triggered failure mode. Although very simple and effective, the ply-discount method suffers from a major drawback. Even if the failure criteria employed may be physically based, the choice of degraded material parameters and the level of degradation is completely arbitrary. A good understanding of the failure conditions and damage development is required to produce accurate results. Also, the results predicted should preferably be correlated to experimental evidence to ensure their validity. Despite its shortcomings, the ply-discount method has been successfully implemented by numerous authors to different joint configurations. The most commonly used failure criteria in three-dimensional models e.g. [10] are the criteria proposed by Hashin [27]. In two-dimensional models, e.g. [3–5] the applied criteria used are the ones proposed by Chang and Chang [28] and Chang and Lessard [29]. In the current study, the failure criteria proposed by Olmedo and Santiuste [6] are utilized. These criteria are based on previously proposed criteria by Chang and Chang [28] and Chang and Lessard [29], which were extended in [6] to include the out-ofplane stresses. Both the extended and the original criteria include the following non-linear relation between the in-plane shear strain and the shear stress, as derived by Hahn and Tsai [30]

c12 ¼

1 s12 þ as312 G12

ð1Þ

where G12 is the in-plane shear modulus and a is an experimentally determined constant. If a ¼ 0, the non-linear strain–stress dependency is eliminated. It is assumed that the non-linearity mainly is caused by micro-cracking in the matrix. Therefore, similar relations can be assumed to hold for the out-of-plane shear strains, as was done by Olmedo and Santiuste [6]. In order to implement Eq. (1) into the FE-model, the non-linear term is rewritten such that a numerically stable algorithm is obtained, as described in [24]. The shear stress at the end of the load increment i is then given as a linear function of the shear strain ðiþ1Þ sðiþ1Þ ¼ ð1  dÞG12 c12 12

ð2Þ

where ðiÞ 2

d¼

ðiÞ 3

ðiÞ

3aG12 ðs12 Þ  2aðs12 Þ =c12 1 þ 3aG12 ðs

ðiÞ 2 12 Þ

ð3Þ

201

The parameter d can be seen as a damage parameter that degrades the initial shear modulus G12 to its current value ð1  dÞG12 . In [6], the proposed criteria include fiber tensile and compressive failure, matrix in-plane tensile and compressive failure, matrix out-of-plane tensile and compressive failure and fiber–matrix shearing failure. Fiber tensile failure is assumed to have occurred in a material point if the following condition is satisfied for r11 P 0

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u  s212 s2 u 2 þ 34 as412 2G1313 þ 34 as413 u r11 2G12 t ¼1 þ 2 þ 2 S12 S XT þ 34 aS412 2G1313 þ 34 aS413 2G12

ð4Þ

The third term under the square root is due to the contribution from the out-of-plane shear stress to the fiber failure. If a ¼ 0, and the shear strength S12 ¼ S13 , the criterion reduces to the tensile fiber mode criterion proposed by Hashin [27]. Originally, Chang and Chang [28] assumed that S12 used should be the in situ ply shear strength measured from a cross-ply laminate, since that value may differ significantly from the strength measured from unidirectional plies. Fiber compressive failure is assumed to take place if the compressive stress (r11  0) is equal to the fiber buckling strength

jr11 j ¼1 XC

ð5Þ

This is the simplest and the most commonly used assumption, where obviously, the influence of the other stress components is neglected. The matrix in-plane failure criterion is met if

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u s223 u 2 s212 þ 3 as4 þ 34 as423 u r22 12 2G12 2G23 4 t ¼1 þ 2 þ 2 S12 S Y þ 34 aS412 2G2323 þ 34 aS423 2G12

ð6Þ

where Y ¼ Y T for r22 P 0 and Y ¼ Y C for r22  0. The same comments following the fiber tensile failure criterion are valid here, except that the criterion reduces to the Hashin and Rotem criterion [31] for a ¼ 0. Matrix out-of-plane failure is considered by expanding the Hashin [27] criterion for delamination initiation with the non-linear shear stress/strain behavior, which gives

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u s223 u 2 s213 þ 3 as4 þ 34 as423 u r33 13 4 2G13 2G23 t ¼1 þ 2 þ 2 S13 S Z þ 34 aS413 2G2323 þ 34 aS423 2G13

ð7Þ

where Z ¼ Z T for r33 P 0 and Z ¼ Z C for r33  0. This criterion is not present in the original set of criteria from [28,29]. It is included in [6] to account for the damage caused by effects of tightening torque in pretightened joints. Fiber–matrix shearing failure is initiated for r11  0 when

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u  s212 s2 u 2 þ 34 as412 2G1313 þ 34 as413 u r11 2G t ¼1 þ 212 þ 2 S12 S XC þ 34 aS412 2G1313 þ 34 aS413 2G12

ð8Þ

As in Eq. (4) and (6) the out-of-plane shear stress is added to the original criterion. It is obvious that the criteria (8) and (5) are identical if the shear stresses vanish. Consequently, the latter criterion should also enclose the fiber compression. However, fiber buckling may follow subsequent to fiber–matrix shearing since only the shear stiffnesses degrade after fiber–matrix shearing failure. If a failure initiation criterion is met at a material point, certain stiffness parameters of that material point are reduced. The simplest alternative would be to reduce all parameters to zero, which

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210

corresponds to a sudden brittle failure with no energy absorption in the damage process. This is not a very realistic assumption, since it has been confirmed in many experimental studies that the stress degrades gradually after the onset of failure. In their work, Zang and Rowland [11] and Tserpes et al. [25] concluded that the reduction of the parameters to zero gives a too severe degradation when compared to experimental results. It also causes convergence problems in the FE-solution procedure because the elements in the area of large contact pressure tend to become highly distorted. A more realistic set of degradation rules were introduced by Tan [32] and extended by Camanho [10] to three dimensions. Similar rules were used in [4,6,11,25] all with good correlation to experimental results. The same degradation factors are used in this work and are shown in Table 3. The progressive damage model, with the failure criteria described in this section, is implemented using the Abaqus subroutine USDFLD. This subroutine allows the user to assign the material properties as functions of user defined field variables, which in turn can be functions of field quantities at material points. In this case, these variables are defined as the failure initiation criteria presented in the equations above and are functions of the stress tensor components at material points. The analysis is driven by incremental increase of the applied load and the strain and stress fields are evaluated at the end of each increment. When the stresses fulfill a failure criterion at a material point, the material properties are degraded. At this point, iterations should be performed at the same load level in order to capture the redistribution of the stresses to other material points which then also may fail. However, if the load increment size is very small, these additional iterations may be omitted [8]. The analyzes preformed in this study utilize the automatic time stepping control with maximal load increment of one per cent of the total applied load and are terminated when severe convergence difficulties are encountered. 2.3. Model validation The results from the FE-analyzes are compared to experimental results performed by Ireman et al. [16]. In the experimental set up, the behavior of the joint was evaluated by measuring the applied force–displacement response and the relative displacements of the plates. The latter measurements were performed at low load levels and are not used here for comparison with FE-results. In Fig. 4, the applied force–displacement response from experiments is compared with three different versions of FE-results. In the first version, only the linear elastic material model is utilized in all components, i.e. no composite damage or metal plasticity is involved. The second version includes the elasto-plastic material model for metal parts and in the third version both the composite damage model and the metal plasticity are present. Initially, the response is linear, both in the experiment and in the analyzes, indicating the absence of pretension and bolt hole clearance effects. At approximately 7 kN, the experimental curve and the curve with composite damage and metal plasticity start to deflect from the linear behavior. The other two curves continue linearly until the curve with plasticity deflects slightly from the linear one. This indicates that the effects of damage accumulation in

Table 3 Degradation rules in different failure modes. Failure mode

E11

E22

E33

G12

G13

G23

m12

m13

m23

Matrix failure Fiber failue Fiber–matrix shear Shear non-linearity

– 0.14 – –

0.4 0.4 – –

0.4 0.4 – –

– 0.25 0.25 1d

– 0.25 0.25 –

0.2 0.2 – –

0 0 0 –

0 0 0 –

0 0 – –

18

Bolt shear failure

16

Experiment AC6III FE, no damage, no plasticity

14

FE, no damage and plasticity FE, damage and plasticity

12

F [kN]

202

10 8 6 4 2 0

0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

u [mm] Fig. 4. Comparison between experimental results and FE-results.

the composite are initiated at at about 7 kN. The experimental curve reaches a maximal load just below 14 kN, where ultimate bearing failure in the composite plate occurs and the load carrying capacity starts to diminish. At the same point, the FE-analysis starts having convergence difficulties and the run is terminated shortly afterwards. Fig. 4 also shows the load level at which the ultimate shear failure of the bolt is expected to occur. This level is simply the material shear strength times the area of the bolt cross-section. It is evident that the load in the experiment never quite reaches this level, which confirms that the joint failed in the bearing failure mode. However, it is reasonable to expect that the failure mode would be bolt shear failure for the same bolt diameter and a significantly thicker laminate. Fig. 5 shows the cross-section of the composite plate at 14 kN and the elements where the failure criteria are met and material properties degraded. It is concluded that the PDM material model used is able to capture the failure initiation point and the damage development in the composite up to ultimate failure load. The progression of failure beyond this point is not modeled. The degradation parameters used here were successfully employed in other studies, e.g. in [10], for similar quasi-isotropic laminates but for different types of bolt installations. However, it is uncertain whether the same parameters would be appropriate for other types of lay-ups and laminate materials. Therefore, experimental verification of the chosen parameters should always be performed for new laminate configurations. Alternative approaches to modeling damage and subsequent loss of stiffness in composite materials have been developed in the context of continuum damage mechanics (CDM). After failure initiation, the damage development due to each failure mechanism is typically controlled by a damage variable. The material stiffness properties are successively reduced depending on the current damage state, which in turn depends on a chosen damage evolution law. This law is thus a key parameter in the framework of CDM. The choice of law is often based on some kind of physical or energy-conserving ground. Usually, once the damage variable has reached its maximal allowable value the material stiffness is exhausted. A large number of CDM models are published in the literature. The implemented material models are often demonstrated on regular coupons or specimens with open hole tension or compression and they seem to work well. Bearing failure of bolted joints, on the other hand, is problematic in the sense that the area of damage concentration is also the area of contact between the

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210

FV1 (Avg: 75%) +1.0e+00 +5.0e−01 +0.0e+00

FV2 (Avg: 75%) +1.0e+00 +5.0e−01 +0.0e+00

(a) Matrix damage

203

FV3 (Avg: 75%) +1.0e+00 +5.0e−01 +0.0e+00

(b) Fiber-matrix shear damage

(c) Fiber damage

Fig. 5. Cross-section of the composite plate at 14 kN and the damage field variables.

bolt and the hole edge. A complete degradation of stiffness and the subsequent element deletion in the area of contact creates convergence problems and yields poor results. From a physical point of view and based on the experimental observations, material removal is not very realistic. The material that has failed due to the bearing failure mode seem to have some residual stiffness left but there is no obvious way to model this kind of behavior and further research is required. 3. Structural modeling This section deals with the FE-modeling of a wing-like CFRPaluminum structure, that is to be fatigue tested at elevated temperature. The purpose of designing, analyzing and testing such a structure is to address the issues of hybrid structures described in the introducing part of this article. The focus is directed on modeling and studying the failure behavior of structural joints that contain a large number of bolts and on the redistribution of loads due to joint failure. The structure is a box of dimensions 3300x630x150 mm and consists of four CFRP skins, eight C-shaped aluminum spars, six aluminum rib sections and one aluminum splice section, see Fig. 6. The material properties are the same as in Table 1 and Table 2. All parts are held together by bolted joints with titanium bolts. All four skins have the same lay-ups and thicknesses and each skin has three areas with different thickness: the skin area overlapping the splice has the thickness of 8.32 mm, the area between the splice and rib Section 6.24 mm and the area outside the rib Section 7.28 mm. Each composite skin area is

quasi-isotropic and symmetrically stacked using a layup sequence ½ð45=0=90Þn S , where n is the number of sub-sequences. Notice that this is the same layup as in the bolted joint in the previous section. Between each skin area there is a ply drop-off region. Table 4 shows the thicknesses of the aluminum parts. All bolts have a 6 mm diameter except for the splice-skin bolts which are 8 mm. The area between the ribs is the area of main interest in the experiment. In this area, all the bolts in CFRP-aluminum joints are countersunk. In the experimental set up, the wing box will be placed on a steel support base, simply supported at the ends on both upper and lower side. The mechanical loads will be applied via two rectangular steel frames, clamped around the section of the box at rib section positions. Four load actuators, two on each frame, will be used to apply the load in the transverse direction. The diagonally placed actuators will be assigned the same load, which will result in a four point bending/twisting situation, see Fig. 7. Moreover, the testing set up will at the same time be placed in a furnace that generates temperatures up to 90 °C during the fatigue and static testing. 3.1. The finite element model of the structure The FE-model of the experimental set up is shown in Fig. 7. All aluminum parts are represented as their middle surfaces and meshed with shell elements S4R. The CFRP skins are modeled with continuum shell elements, SC8R, that discretize the entire threedimensional geometry. One element with stacked lay-up formulation is used in the thickness direction. The steel frames are

Fig. 6. Studied wing box structure with dimensions 3300  630  150 [mm]. One of the upper skins is not shown and the other one is partially hidden for easier visibility of the inner parts. The bolts connecting the skins with the inner structure are indicated in the figure.

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210

204

Table 4 Aluminum component thicknesses. The flanges of each component are the parts connected to the skins, webs are the longitudinal parts and the stiffeners are the parts that are perpendicular to the webs, see Fig. 6.

Spar Rib Splice

Flange

Web

Stiffener

4.5 4.5 6

3 3 3

2 3 3

modeled using solid brick elements and for the supports, rigid elements are used. All parts are assigned linear elastic material behavior. Contact conditions are enforced between all parts in the model with the general contact algorithm. The solution of the FE-model is obtained in several consecutive general steps in Abaqus/Standard. In the first step, a pretension load is applied to the vertical braces of the steel frames, while the rigid supports are constrained in all directions. The purpose of the pretension is to hold the frame in place while the actuator loads are applied, so that the forces can be transferred to the wing box via contact. For this purpose, a sufficiently high pretension load is chosen to ensure that the contact is maintained between the frame and the box during loading. No estimation of the pretension load is made in correlation with the experimental set up. In the second step, a constant temperature difference DT = 50 °C is applied to all parts. Vertical point loads representing the actuator forces are applied directly onto the steel frame in the third step, see Fig. 7.

from the connector element nodes to the nodes of the parts that are to be joined. In this way the connector elements representing the bolts become independent of the surrounding mesh, which makes the modeling work easier. There are two kinds of element interfaces in the structural model: shell-to-shell (S–S) interface and continuum shell-to-shell (C–S) interface. In the S–S interface, the shell elements placed at the middle surfaces of the plates are linked by the connector elements. In the C–S interface, the shell elements placed at the middle of the aluminum plate are linked to the lower side of the continuum shell element representing the composite. This set-up is illustrated in Fig. 8. Table 5 summarizes the different bolted joint configurations and interfaces present in the FE-model of the wing box. The local behavior for each connector element, in terms of force–displacement characteristics, is determined from a local solid FE-model of a single bolt joint. This model is principally the same as the specimen model described previously. However, now the plates are square with 4D long sides, where D is the bolt diameter, as shown in Fig. 9. The small in-plane dimensions of the plates are chosen in order to eliminate the secondary bending effects. It is assumed that these effects also are negligible in the wing box, due to the out-of-plane deflection restrictions imposed on the joints by lateral support of the surrounding structure. However, all singlelap joints have an inherent load path eccentricity, even if the secondary bending is eliminated. This causes the bolt to tilt and bend,

3.2. Bolted joints Each bolt in the structure is represented by Abaqus connector element CONN3D2 of type BUSHING. This is a two-node line element, that has six degrees of freedom (DOF) in each node and connects two parts at their interface. It is implemented by a Lagrange multiplier technique, where the fastener forces are solved as the constraint forces in the interconnections. The nodes of the connector element are not necessarily placed coincidentally with the nodes at the interface of the bolted parts. Instead, the Abaqus feature *FASTENER is used to deploy distributing coupling elements

Fig. 8. Shell-to-shell (S–S) and continuum shell-to-shell (C–S) interfaces.

Fig. 7. FE-model of the wing box structure. One of the skins is hidden for easier visibility of the inner structure.

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210 Table 5 The configurations of the bolted joints in the wing box.

205

35

Plate materials

Plate thicknesses (mm)

Bolt diameter (mm)

Element interface

Skin-Spar Skin-Splice Skin-Rib Spar-Splice Spar-Rib

CFRP-Al CFRP-Al CFRP-Al Al-Al Al-Al

6.24–4.5 8.32–6 7.28–4.5 3–3 3–3

6 8 6 6 6

C–S C–S C–S S–S S–S

Bolt shear failure

30 25

Pretension No pretension

F [kN]

Bolted articles

20

Semi−empirical

15 10

F

5 0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

u [mm]

F

C B A

Fig. 9. Local solid model of the single bolt joint.

which in turn gives an uneven contact pressure distribution over the thickness which affects the damage distribution. Another effect that is overlooked is the influence of the by-pass loads on the force–displacement characteristics of a joint. By-pass loads certainly affect the damage propagation in the composite plate [13] but it is unclear to what extent they also influence the flexibility. Indeed, the by-pass loads are present in the wing box joints. Moreover, the proportion of by-pass to transferred loads vary for different bolts and depend on the applied load level. This means that the connector element characteristics would have to be bolt-dependent and solution-dependent in order to accommodate for any by-pass load effects on the flexibility. The bolts in the wing box will be torque tightened in the experiment. Therefore, the pretension effects on the joint flexibility in the local model is briefly studied. Fig. 10 shows the applied force–displacement response of a local solid bolt model from the skin-splice joint, with and without a pretension load on the bolt. The applied pretension force is 15 kN and is included in the model by an Abaqus built-in feature for pretension of sections. The figure also shows the linear semi-empirically established response proposed by Huth [33], which is often used in industrial applications in load distribution analyzes. The load level, where the bolt shear failure is expected, is also shown in the figure. The effects of the tightening force on the force–displacement behavior are apparent. In the initial loading, the pretension force results in a higher load transfer for the same displacement due to the friction forces. With further loading the two curves have the same, almost linear, slope until the curve with no pretension starts to deflect. This is due to the initiation of the damage in the composite plate. The curve with pretension departs from linear at a higher load level. At the maximal load level, the two curves converge to one, due to the loss of frictional contact. These effects have also been confirmed experimentally by Ireman et al. [16]. It is also interesting to notice that the load level, to which both curves are converging, is the predicted bolt shear failure load. As anticipated, for thicker laminates the failure mode is bolt shear failure. Since the ultimate failure level of the joint seems to be unaffected

Fig. 10. Response of the local solid bolt model in the skin-splice joint as the applied force vs. the displacement of the loaded edge, see Table 5 for dimensions.

by pretension, the pretension is neglected in connector element characterization. The displacement that is used to characterize the behavior of the connector element is extracted form the local solid model in different ways for the two element interface types. For the S–S type the connector element displacement is defined as the relative displacement of the points A and C in Fig. 9, and for the C–S interface type the relative displacement of points A and B is extracted. In this way the relative displacement of the connector element nodes is related to the points where the element is attached to the plates. The force–displacement curve extracted from the solid model is given to the connector element as elastoplastic behavior assigned to the element shear resultant force as function of the resultant in-plane displacement. The resultants are formed as vector sums of the in-plane force and displacement components. Remaining DOFs are assigned rigid behavior except for the drilling DOF, i.e. rotation around the element line, which is released. Once the load level of the shear resultant reaches the bolt shear failure load, the connector damage functionality in Abaqus is activated and the load is linearly decreased to zero over a very short displacement range. This range is set to 1% of the maximal plastic displacement. The sudden loss of load carrying capacity signifies the nature of fastener shear failure. When the resultant force has dropped to zero the connector element is removed from the model. Fig. 11 shows the force–displacement curves for the three joint configurations present in the area between the rib sections. The connector elements outside the rib sections are expected to be exposed to forces within the linear elastic regime and are therefore assigned linear elastic behavior according to Huth [33]. Damage behavior in the connector elements is characterized by softening response that causes local instabilities and convergence difficulties in the implicit code. To overcome this issue two features available in Abaqus/Standard are activated: viscous regularization and stabilization. The first one is applied on connector element level as a regularization of the elements constitutive response by introducing viscosity into the damage evolution law. The second feature provides stabilization of the solving process by adding viscous forces to the global equilibrium equations. The amount of viscosity introduced into the solution is in both cases controlled by user-defined parameters, which have to be chosen carefully. If the viscosity effect is too small it might not eliminate the convergence problems, and too much of it might distort the

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210

206

Skin−Splice

35

Spar−Splice Bolt shear failure, 8mm bolt

30

F [kN]

25 20 Bolt shear failure, 6mm bolt

15 10 5 0

3.3. Structural analysis results

Skin−Spar

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Δu [mm] Fig. 11. Elasto-plastic and damage behavior of the connector elements representing the bolts in the area between the rib sections. For configuration of the joints, see Table 5 and for the location, see Fig. 6.

solution. Anyway, the energy dissipated by viscous regularization and stabilization should be checked a posteriori to make sure that it is negligible compared to the total energy.

This section presents the responses of three different bolt groups in the wing-box, see Fig. 12(a). The first group is denoted ’’Skin-Spar bolt group’’ and connects two of the skins to one of the outer spar pairs and to the splice. The bolts are placed along x-direction in a single row containing 130 bolts. The second bolt group is the Skin-Splice bolt group and it consists of four bolt rows in y-direction, with 16 bolts in each row. This group connects two of the composite skins to the splice flange. There are four bolts in this group that also belong to the first bolt group and these are the ones placed at the edge of the splice flange, between the outer spars. Both of the mentioned bolt groups are situated on the tensile loaded side of the wing box, which corresponds to the lower side of a real wing. The third bolt group is the Spar-Splice bolt group and it connects each of the eight spars to the splice stiffeners with 8 joints, each containing 5 bolts placed in a row in the z-direction. Only 5 of these bolt rows are fully visible in Fig. 12(a) and the remaining 3 are indicated by arrows. All fasteners and bolt rows are numbered within their own bolt group and Fig. 12(b) shows the numbering in the splice area. Fig. 13 shows the distribution of resultant shear bolt forces in the three bolt groups due to an applied temperature of DT = 50 °C. Each bolt row in all three bolt groups displays a U-shaped force distribution, which seems to be symmetric with respect to the coordinate axis. In Fig. 13(a) disturbances in the two U-shapes can be noticed. This is due to the expansion of the

Fig. 12. Studied bolt groups with row and fastener numbering.

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210

Skin−Splice bolt group

4000

4000

3500

3500

3000

3000

2500

2500

F [N]

F [N]

Skin−Spar bolt group

2000

2000

1500

1500

1000

1000

500

500

0

0

20

40

60

80

207

100

0

120

row 1 0

row 2 16

row 3 32

Bolt nr

row 4 48

64

Bolt nr

(b) Skin-Splice bolt group, 4 bolt rows in y-direction.

(a) Skin-Spar bolt group, 1 bolt row in x-direction.

Spar−Splice bolt group 4000 3500 3000

F [N]

2500 2000 1500 1000 500 0

row 1 row 2 row 3 row 4 row 5 row 6 row 7 row 8 0

5

10

15

20

25

30

35

40

Bolt nr

(c) Spar-Splice bolt group, 8 bolt rows in z-direction. Fig. 13. Distribution of resultant shear bolt force due to applied temperature of DT = 50 °C in all three bolt groups studied.

aluminum ribs in the y-direction, whereby a force in the same direction is transferred via the skin-spar bolts at rib positions. From the overall maximal force levels and Fig. 11, the conclusion is drawn that all bolts still are in the linear elastic regime at the currently applied temperature. Fig. 14 presents the distribution of bolt forces due to applied temperature and mechanical loads. Each graph contains two curves. The curve, denoted ’’Failure onset’’, shows the force distribution just before the first fastener fails, when the elements in the area between the ribs are assigned the behavior from Fig. 11. The other curve shows the force distribution at the same applied load level, but with linear elastic behavior assigned to all connector elements. The applied actuator load levels at the onset of first bolt failure are P 1 ¼ 253:3 kN and P 2 ¼ 130:6 kN, which gives the bending and twisting moments M b ¼ 134:9 kN m and M t ¼ 72:9 kN m in the cross-section of the area between the ribs. These loads are 1.43 times higher than the dimensioning loads for ultimate failure when Saab in-house methods for dimensioning are used. It is evident that the effects of thermally induced loads cannot be neglected. It can be seen from Fig. 14 that a linear elastic behavior assigned to the fasteners leads to an overestimation of the maximal bolt loads at the onset of the first fastener failure which therefore is conservative. It is also noticed that several bolts in bolt groups Skin-Splice and Spar-Splice approach their corresponding ultimate failure force at 30 kN and 17 kN, respectively. Most of the bolts in bolt group Skin-Spar are still in the linear elastic regime, as anticipated.

A slight increase in the applied load leads to the first ultimate bolt failure, at bolt nr 16 and 32 in rows 1 and 2 at the edge of the Skin-Splice bolt group, see Fig. 12(b). Fig. 15 shows the force redistribution in the Skin-Splice and the Spar-Splice joints. As a result of loss of load carrying capacity in the failed bolts, the forces are redistributed to the neighboring fasteners, i.e. the fasteners close to the end of each row of Skin-Splice bolt group and in the Spar-Splice joint, especially in row 4. The forces in the fasteners at the very end of rows 3 and 4 in Skin-Splice bolt group however, tend to decrease. The reason for this is that the loss of a load path at the end of rows, 1 and 2 causes the load to redistribute towards the middle of the splice, thus decreasing the transferred load through the end of the rows 3 and 4. Further progression of the fastener failure in the Skin-Splice bolt group is shown in Fig. 16. It can be seen from Fig. 16(a) that at this point, three fasteners in each of the rows 1 and 2 have completely failed and that two more are on their way to fail. The loads in the fasteners at the end of rows 3 and 4 have decreased further. The effects of load redistribution to other fasteners is more pronounced. It is expected that the amount of total fastener shear load transferred over the splice joint is the same before and after the progressive fastener failure. However, there are several factors that can disturb this equilibrium. One is that the stabilization option used to damp the local instabilities dissipates some of the total energy and consequently can lower the transferred forces. This effect is however negligible in this case, since the dissipated energy by stabilization is confirmed to be very low, under 1% of the total

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210

208

x 10

Skin−Spar bolt group

4

x 10

3.5

3.5

Linear elastic connector Failure onset

3

3 2.5

F [N]

F [N]

2.5 2

2

1.5

1.5

1

1

0.5

0.5

0

Skin−Splice bolt group

4

0

20

40

60

80

100

row 1

0 0

120

row 2 16

Bolt nr

row 4 48

64

Bolt nr

(a) Skin-Spar bolt group.

(b) Skin-Splice bolt group. Spar−Splice bolt group

4

2.5

row 3 32

x 10

F [N]

2

1.5

1

0.5

0

0

row 1 row 2 row 3 row 4 row 5 row 6 row 7 row 8 5 10 15 20 25 30 35 40

Bolt nr

(c) Spar-Splice bolt group. Fig. 14. Distribution of resultant shear bolt force due to applied temperature of DT = 50 °C and mechanical bending/twisting loads at first bolt failure onset.

Skin−Splice bolt group

4

x 10

4

2.5

x 10

Spar−Splice bolt group

Failure onset 2 bolts failed

3.5

2

3

F [N]

F [N]

2.5 2

1.5 1

1.5 1

0.5 0.5 0 0

row 1

row 2 16

row 3 32

row 4 48

64

Bolt nr

(a) Skin-Splice bolt group.

0 0

row 1 row 2 row 3 row 4 row 5 row 6 row 7 row 8 5 10 15 20 25 30 35 40

Bolt nr

(b) Spar-Splice bolt group.

Fig. 15. Distribution of resultant shear bolt force due to applied temperature of DT = 50 °C and mechanical bending/twisting loads shortly after the first bolt failure.

strain energy. Another effect that can lower the total fastener load is that the thermally induced forces decrease with fewer fasteners. In other words, when ten of the outer fasteners are broken, some of the energy stored in the joint by thermal expansion is released and the remaining fastener forces decrease. Due to the loss of the fasteners at the end of rows 1 and 2, the third effect is that the outer Spar-Splice starts twisting and joint row 4 starts transferring load via the splice stiffener by contact and by axial bolt forces. The axial

force is not one of the components that contributes to the shear resultant force and this effect is therefore not evident from Fig. 16(b). However, the influence of shifted load transfer on load redistribution to the fasteners outside the splice area can be seen in Fig. 17. Due to the loss of Skin-Splice fasteners the load is redistributed to Spar-Splice joint row 4, which leads to increased shear forces in the part of Skin-Spar joint, bolts 50–63.

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210

Skin−Splice bolt group

4

x 10

4

2.5

Failure onset 10 bolts failed

3.5

x 10

209

Spar−Splice bolt group

2

3

1.5

F [N]

F [N]

2.5 2 1.5

1

1 0.5 0.5 0

row 1 0

row 2 16

row 3 32

row 4 48

64

0

0

row 1 row 2 row 3 row 4 row 5 row 6 row 7 row 8 5 10 15 20 25 30 35 40

Bolt nr

Bolt nr

(a) Skin-Splice bolt group.

(b) Spar-Splice bolt group.

Fig. 16. Distribution of resultant shear bolt force due to applied temperature of DT = 50 °C and mechanical bending/twisting loads after failure of 10 fasteners in Skin-Splice bolt group.

4

Skin−Spar bolt group

x 10 3.5

Failure onset 10 bolts failed

3

F [N]

2.5 2 1.5 1 0.5 0

0

20

40

60

80

100

120

Bolt nr Fig. 17. Distribution of resultant shear bolt force in Skin-Spar bolt group due to applied temperature of DT = 50 °C and mechanical bending/twisting loads after failure of 10 fasteners in Skin-Splice bolt group.

4. Conclusions Local behavior of a bolted composite-aluminum shear lap joint is modeled using PDM in the composite plate. The resulting force– displacement characteristics correlate reasonably well with the experimental results found in the literature. It is concluded that PDM, despite its shortcomings, can be used to determine the behavior of composite failure in bearing failure mode. The initiation of composite damage is observed locally at low load levels, but the effect of the damage on the force–displacement curve is not present until approximately half the maximal load. This behavior should be carefully considered in the design situations, when load levels for permanent deformations and ultimate failure are determined. Mixing of aluminum and composite materials in structures gives rise to effects that might be absent otherwise. One of the effects is that thermally induced loads arise in hybrid assemblages as a result of the difference in thermal expansion properties of the materials. It is shown that the fastener forces caused by temperature difference are of significant magnitude and should be considered in the design of hybrid structures. It is also noticed that different material failure behaviors affect the local stiffness of bolted joints and consequently the load distribution in the structure. Because of the size of the numerical problem, solid

modeling including failure behavior is not suitable for modeling larger structures including many joints. Some simplifications are necessary and in this work connector elements are used to represent the fasteners and the behavior that they exhibit. The force– displacement response of connector fasteners is determined from a local solid fastener model including PDM for the composite plate. In this study, the important effects of bolt hole clearance were not included. However, connector elements do have a capability to include these kinds of behavior and this should be further looked into. Also, the effects of by-pass loads on the local force–displacement characteristics of hybrid joints need to be further studied. The fastener load distribution due to bending/twisting and thermally induced loads in a large hybrid wing box is studied. It is shown that considering the non-linear response of the joints results in less conservative load distributions in typical wing structures on ultimate failure load level. On operational load levels, the composite damage induced by repeated loading, might lead to bolt hole widening and consequently, to load redistribution. This phenomenon should be studied further. Acknowledgments The authors would like to acknowledge support from the HYBRIS-project, Swedish Defense Material Administration, and Saab AB who funded the work presented in this article. References [1] Orifici AC, Herszberg I, Thomson RS. Review of methodologies for composite material modelling incorporating failure. Compos Struct 2008;86(1-3): 194–210. [2] Camanho PP, Mattews FL. Stress analysis and strength prediction of mechanically fastened joints in FRP: a review. Compos A: Appl Sci Manuf 1997;28(6):529–47. [3] Dano ML, Gendron G, Picard A. Stress and failure of mechanically fastened joints in composite laminates. Compos Struct 2000;50(3):287–96. [4] Dano ML, Kamal E, Gendron G. Analysis of bolted joints in composite laminates: Strains and bearing stiffness predictions. Compos Struct 2007;79(4):562–70. [5] Xiao Y, Ishikawa T. Bearing strength and failure behavior of bolted composite joints (part II: modeling and simulation). Compos Sci Technol 2005;65(7– 8):1032–43. [6] Olmedo A, Santiuste C. On the prediction of bolted single-lap composite joints. Compos Struct 2012;94(6):2110–7. [7] Hühne C, Zerbst A-K, Kuhlmann G, Steenbock C, Rolfes R. Progressive damage analysis of composite bolted joints with liquid shim layers using constant and continuous degradation models. Compos Struct 2010;92(2):189–200. [8] McCarthy CT, McCarthy MA, Lawlor VP. Progressive damage analysis of multibolt composite joints with variable bolt-hole clearances. Composites: Part B 2005;36(4):290–305. [9] Linde P, Pleitner J, De Boer H, Carmone C. Modelling and simulation of fiber metal laminates. In: ABAQUS users conference; 2004. [10] Camanho PP, Mattews FL. A progressive damage model for mechanically fastened joints in composite laminates. J Compos Mater 1999;33:2248–80.

210

Z. Kapidzˇic´ et al. / Composite Structures 109 (2014) 198–210

[11] Zhang J, Rowland J. Damage modeling of carbon-fiber reinforced polymer composite pin-joints at extreme temperatures. Compos Struct 2012;94(8):2314–25. [12] Ireman T. Three-dimensional stress analysis of bolted single-lap composite joints. Compos Struct 1998;43(3):195–216. [13] Rosales-Iriarte F, Fellows NA, Durodola JF. Failure prediction in carbon composites subjected to bearing versus bypass loading. J Compos Mater 2012;46(15):1859–78. [14] Egan B, McCarthy CT, McCarthy MA, Gray PJ, Frizzell RM. Modelling a singlebolt countersunk composite joint using implicit and explicit finite element analysis. Comput Mater Sci 2012;64:203–8. [15] Chishti M, Wang CH, Thomson RS, Orifici AC. Numerical analysis of damage progression and strength of countersunk composite joints. Compos Struct 2012;94(2):643–53. [16] Ireman T, Ranvik T, Eriksson I. On damage development in mechanically fastened composite laminates. Compos Struct 2000;49(2):151–71. [17] Chishti M, Wang CH, Thomson RS, Orifici AC. Experimental investigation of damage progression and strength of countersunk composite joints. Compos Struct 2012;94(3):865–73. [18] Weyer S, Hooputra H, Zhou F. Modeling of self-piercing rivets using fasteners in crash analysis. In: ABAQUS users conference; 2006. [19] Gray PJ, McCarthy CT. A global bolted joint model for finite element analysis of load distributions in multi-bolt composite joints. Composites: Part B 2010;41(4):317–25. [20] Gray PJ, McCarthy CT. A highly efficient user-defined finite element for load distribution analysis of large-scale bolted composite structures. Compos Sci Technol 2011;71(12):1517–27. [21] Ekh J, Schön J. Finite element modeling and optimization of load transfer in multi-fastener joints using structural elements. Compos Struct 2008;82(2):245–56.

[22] Wang AJ, Chung KF. Advanced finite element modelling of perforated composite beams with flexible shear connectors. Eng Struct 2008;30(10):2724–38. [23] MIL-HDBK-17-3F. Military handbook, polymer matrix composites. U.S. Department of Defense; 2000. [24] Abaqus Analysis User’s Manual, Simulia; 2011. . [25] Tserpes KI, Labeas G, Papanikos P, Kermanidis Th. Strength prediction of bolted joints in graphite/epoxy composite laminates. Composites: Part B 2002;33:521–9. [26] Schön J. Coefficient of friction for aluminum in contact with carbon fiber epoxy composite. Tribol Int 2004;37(5):395–404. [27] Hashin Z. Failure criteria for unidirectional fibre composites. J Appl Mech 1980;47:329–34. [28] Chang FK, Chang KY. A progressive damage model for laminated composites containing stress concentrations. J Compos Mater 1987;21(9):834–55. [29] Chang FK, Lessard LB. Damage tolerance of laminated composites containing an open hole and subject to compressive loadings: part I – Analysis. J Compos Mater 1991;25(1):2–43. [30] Hahn HT, Tsai SW. Nonlinear elastic behavior of unidirectional composite laminates. J Compos Mater 1973;7(1):102–18. [31] Hashin Z, Rotem A. A fatigue failure criterion for fiber reinforced materials. J Compos Mater 1973;7:448–64. [32] Tan SC. A progressive failure model for composite laminates containing openings. J Compos Mater 1991;25(5):556–77. [33] Huth H. Influence of fastener flexibility on the prediction of load transfer and fatigue life for multiple-row joints. Fatigue in mechanically fastened composite and metallic joints. ASTM STP 927, Philadelphia; 1986. p. 221–50.