The Effects of Fatigue Cracks on Fastener Loads during Cyclic Loading and on the Stresses Used for Crack Growth Analysis in Classical Linear Elastic Fracture Mechanics Approaches

High strength threaded fasteners are widely used in the aircraft industry, and service experience shows that for structures where shear loading of the joints is significant, like skin splices, fuselage joints or spar caps-web attachments, more cracks are initiated and grow from the edges of the fastener holes than from features like fillets radii and corners or from large access holes. The main causes of this cracking are the stress concentrations introduced by the fastener holes and by the threaded fasteners themselves, with the most common damage site being at the edge of the fastener holes. Intuitively, it is easy to visualize that after the crack initiation, during the growth stages, some of the load transferred initially by the fastener at the cracked hole will decrease, and it will be shed to the adjacent fasteners that will carry higher loads than in uncracked condition. Using currently available computer software, the method presented in this paper provides a relatively quick and quantitatively defined solution to account for the effects of crack length on the fastener loads transfer, and on the far field and bypass loads at each fastener adjacent to the crack. At each location, these variations are determined from the 3-dimensional distribution of stresses in the joint, and accounting for secondary bending effects and fastener tilt. Two cases of a typical skins lap splice with eight fasteners in a two rows configuration loaded in tension are presented and discussed, one representative for wing or fuselage skins configurations, and the second case representative for cost effective laboratory testing. Each case presents five cracking scenarios, with the cracks growing from approx. 0.03 inch to either the free edge, next hole or both simultaneously.


Introduction
With the rise of fracture mechanics, the management and certification of safety critical or economically significant structures has switched towards crack propagation models using linear elastic fracture mechanics (LEFM) which correlates parameters like crack tip stress intensity factor to its growth rate. At the same time, new stronger materials have led to higher working stresses and, in many situations, having higher crack growth rates, so once initiated, the duration to a critical size is reduced [1].
High strength threaded fasteners are widely used in the aircraft industry, and service experience shows that for structures where shear loading of the joints is significant, like skin splices, fuselage joints or spar caps-web attachments, more cracks are initiated and grow from the edges of the fastener holes than from features like fillets radii and corners or from large access holes. The main causes of this cracking, particularly in the case of single shear joints, are the stress concentrations introduced by the fastener holes and by the threaded fasteners themselves with the most common damage site being at the edge of the fastener holes.
These types of joints have been addressed extensively by researchers and engineers using experimental and analytical techniques (e.g. [2] [3] [4] [5]), and two major steps, with their desired methods and data, have been identified for the analysis of the fatigue performance of fastened joints. These steps are: • The calculations of the force distributions in the joint, with methods and data to account for fastener flexibility, eccentricities in the load paths and the support from surrounding structure, and • The calculation of the fatigue quality, which requires the influence from the bypass loads, the loads transferred by the fasteners, the moments due to eccentricities and the effects of bearing and friction interplay to be known for all the materials, hole quality and fastener types in the joint.
A systematic review of some approaches used for fatigue analysis is presented by E. Jarfall in [6]. In reference [7], the same author defines a stress severity factor concept, and shows that both the far field and bearing contact loads have a significant influence on the fatigue strength of a shear joint. Methods and test data for shear and tension joints have been developed by many major aircraft manufacturers and/or research organizations. Some are published in the technical literature [8]- [13]. In all these methods, a clear distinction is made between hole filling rivets and threaded fasteners.
In a typical LEFM crack growth analysis, the effects of crack length and the in- known as "the beta functions". The load transferred at the fastener is assumed to remain constant as one crack grows to the next nearby feature and, after that, it is often assumed that there is no more load transfer through that fastener.
Intuitively it is easy to visualize that after the crack initiation, during the growth stages, some of the load transferred initially by the fastener at the cracked hole will decrease, and it will be shed to the adjacent fasteners that will carry higher loads than in the uncracked condition. Concurrently with the changes in the load transfers at the fasteners, the bypass and even the far field loads will have to change to maintain the loads balance.
Using currently available computer software, the method presented in this paper provides a relatively quick and quantitatively defined solution to account for the effects of crack length on the fastener loads transfer at each fastener adjacent to the crack in a single shear loaded joint.
Two cases of a typical skins lap splice with eight fasteners in a two rows configuration loaded in tension are presented and discussed. The first case is representative for wing or fuselage skins and the second case is representative for cost effective laboratory testing. Both models present five cracking scenarios, with the cracks growing from approximately 0.03 inch to either the free edge, the next hole or both simultaneously.
The approach enables development of simple tabular/look up parametric functions for updating the reference stresses used with the "beta functions".

Brief Overview of the Fastener Modelling Approach and Typical Results
Contemporary finite element packages like Nastran or MSC Marc, Abaqus or Ansys have incorporated dedicated methods and/or elements to enable the modelling and evaluation of the local stress/strain fields at and around the pre-tensioned fasteners. The methods available involve splitting of the fastener shank in two parts and the "shortening" of the shank with a set of unidirectional multipoint constraint (MPC) equations until the desired pretension force is reached. After that the two parts are "locked" with the shortening remaining fixed while the other loads are being applied. A brief summary of this approach as used in MSC Nastran [14] is presented in Figure 1.
For all subsequent load cases, the control grid is fixed in all three directions (1,2,3). This approach is very effective when the fasteners are subject to separating tensile forces and when the components remain relatively fixed in operation, with no lateral loads requiring the fasteners to sustain shear forces. It is not compatible with large displacements nor with fasteners loaded in shear.
Another fastener modelling technique, introduced in reference [15], maintains the continuity of the pin without changing its installation clearance (does  not use thermal contraction 1 ) and it is compatible with both large parts displacements and with the transfer of shear forces. For this approach, the bolt pin is kept continuous and a separate mesh is used for the nut. After the pretension step, the bolt and the nut are glued together along the common thread area and this condition is kept for all the subsequent loading steps. The initial pretension is introduced in the first step through two RBE3 elements, one acting on the pin shank and the other on the nut. A very short overview of this technique with typical results relevant to pretension and bearing loads and to the separation threshold [8] is shown in Figure 2. 1 The use of thermal contraction for pretension would require a stiffness based analysis, similar to that presented in [8] for each joint configuration to establish the temperature differential ΔT and it is not suitable when operational thermal considerations are of importance to structural or fatigue strength.

Geometry and FEA Representations
The fastener modelling technique of reference [15] was used for a typical lap splice with eight fasteners in a two rows configuration, loaded in tension only.
The analysis was done for two distinct cases, one representative for wing or fuselage skins marked as model C and the second case representative for cost effective laboratory testing and marked as model D.
Both models are identical in terms of geometry, materials, fixed end supports and applied loads. The difference between models is in the boundary conditions applied along the lateral edges, and they are shown in Figure 3.
For the aircraft configuration, the long edges are fixed in direction 1 (Ox) to approximate the constraint imposed by the adjacent structure, whereas model D This simulates either the rest of the skin on model C or the jaws of the loading machine in model D.
The external load was applied in four increments of 0.5%, 10%, 90% and 100%, is 278.8 lbf (1240 N), corresponding to 63 lbf/in (11.03 N/mm) running load or to a uniform gross area tension stress of 1.0 ksi (6.895 MPa) in the top plate.
The FEA analysis was conducted using shell elements located at the mid planes of the two plates in Nastran solution 400 with non-linear material and  The geometry, materials and the specification for fasteners (including the close tolerance values 2 ) are defined in Figure 4. In the same figure is shown the numbering used for the fastener holes and/or the fastener pins.
In Figure 5 is shown the meshing pattern used. The crack length increment used was defined by the sizes of the edges of the shell elements along the Ox axis.
All units for the input/output were in the imperial system (lbf, in, psi, lbf/in etc.).
For each fastener, four contact bodies were defined (Body_Pin1_Head, Body_Pin1_Thread, Body_Nut1_Head, Body_Nut1_Thread, etc.) and one contact body for each of the plates. Both models contained 21 contact zones defined for the initial set-up and setting in position, and 41 contact zones for the loading stage. With the exception of the contacts between the threads of fasteners and collars (glued), the friction was not included. This is because for conservatism, the friction is neglected in the structural analysis of aircraft structures. For the same reason the data in this paper does not include the initial pretension. Although the effects of pretension are significant for the crack initiation and growth, the aim of this paper is to address and quantify the effects of crack length on the forces (and stresses) used in popular current crack growth analysis programs.
The contacts at the threads of pins and collars were "glued" at the initial set-up (with the card BCTABL1 set to zero).

Cracking Configurations Evaluated and Software Used for Pre and
Post Processing Both models were evaluated for five cracking scenarios, with the cracks growing from approx. 0.03 inch to either the free edge, to the next hole or both simultaneously. The five configurations considered are shown in Table 1 The generation of Patran databases, the first stage of post processing and the extraction of results for each of the runs were conducted using dedicated PCL session files. The second stage of post processing contained the analysis and systematization of the results followed by data reduction for interpretation and generation of several levels of summary tables and graphs. It was conducted in batch mode for both models in Excel using VBA subroutines.
The following results were extracted and post processed for each crack length in all configurations of Table 1: 2 In close tolerance the hole sizes are defined by the standardized Transition Fit classification for steel bolts in aluminium structure or by the Close Ream classification for steel bolts in steel or Ti structures. In the cases presented the fit between the pin and the hole in the FEA models was a clearance of 0.002" in diameter.   • Contact forces normal to the pins: these forces were derived from the Partran's outputs for "Resultants Contact Forces", this time at the pin nodes.
There are 12 equally spaced layers of elements on each pin. The purpose of extracting the pin loads was to evaluate the distribution of the bearing forces at each hole at various depth levels of the plates.  • Far field and bypass element centroidal stresses were extracted for the upper and lower surfaces of P1 at each hole, as σ x , σ y and τ xy strictly in CID 0, as max-min principal invariants and as Tresca/max shear. These stresses were "averaged" for the band widths as described in section 3.3.
To keep uniformity with the bolt constant 4 /fastener flexibility 5 concepts of Tate and Rosenfeld/L. Jarfall, the traditional form of the gross area stress 3 The axial running loads were defined as the distributed axial force acting per unit of length (analogue to the classical definition of shear flow) given by the ratio between the sum of the axial forces of the ten elements and their total length in direction Ox. Similar for lateral running loads. 4 Bolt constant was defined as a linear approximation of the total displacement caused by fastener deformation, fastener tilt and the deformation of the fastener holes (all three deformations are nonlinearly dependent on the applied load). There are several semi-empirical formulae used in the industry for its calculations. It is commonly known as the "fastener flexibility constant" and is used to approximate the fastener loads in a 1D beam idealization.

Effect of Crack Configuration and Crack Length on the Forces Transferred by the Fasteners
For splices in tension, there is a distinct difference in the patterns of fastener load redistributions as functions of crack length between the single tip and double tip cracks. Also, the lateral restraints play a role in the way the loads from the fastener at the cracked hole are shed to adjacent fasteners. However, regardless of configuration and crack length, the variations of fastener loads are ± 10%.
Generally the fastener at the cracked hole becomes less loaded and the adjacent fasteners in same row experience load increases. The load shedding extends also to the fasteners in the second row, but the absolute changes in values are smaller.
An envelope of percent variations from the baseline values is presented in Table   2.
The most significant observations from an engineering perspective are: • For a single sided crack initiating from an edge hole and growing towards the free edge (Config. 1), the fastener load is at its minimum value just before the crack breaks though the end ligament. When the crack breaks through, the load at that fastener returns close to its value in baseline configuration. Noting that the fastener head and the collar continue hold the two plates together, this behavior was reconciled by considering the end supports provided by the plate for this fastener which, during the growth stage, are analogous to a beam with both edges restrained, and after breakthrough, the support configuration is single sided ("the hook effect" 6 ). As the ligament at the free edge breaks, the bypass loads in the plate at the other side of the cracked hole increase and the load transferred by fastener #1 increases slightly over its baseline value.
• If, after the break through, a new crack initiates at the other side of the hole and grows inwards as in Config 2, the variations in fastener loads are similar to those in Config 1 but, after the crack breaking into the second hole, the loads at both fasteners remain below their baseline values.
• For one tip cracks initiating from holes remote from the free edges, the variations in the fastener loads are of similar ratios and, after the fracture of the ligament, the hook effects are smaller with both fasteners carrying between 90% and 95% of their initial load.
• The largest reductions in the fastener loads are for double sided cracks. After braking one or two ligaments the fastener at the cracked hole carries 20% less load than in baseline configuration. This load remains unchanged until a new crack starts to grow. 6 Another visualization analogy may be made with a lug, when the crack breaks through one side, the loads are reacted in bending and shear like in a hook.   For quick reference the final baseline fastener loads are shown in Table 3. All values are in CID 0 (the coordinate system shown in Figure 6(c)).

Effect of Crack Length on the Through Thickness Distribution of Bearing Forces between the Fastener Pins and the Hole Surfaces
The results for the bearing distributions at the pin-hole interfaces account for the fastener load transfer, for the eccentricity of the joint (secondary bending), bolt tilting and for the effects of the contact forces between the pins and the surfaces of the holes in the plates and between the fastener heads and the plates, albeit the fact that 2D elements located at midplane were used for the plates. They are however preliminary results because the pretension clamping forces, friction and the factors influencing the tractions forces between the faying surfaces (like cladding, sealants, etc) have not been included.      For this reason the data reduction was somehow limited, with the results being post-processed only on a pin -by-pin basis for each crack configuration at 12 equidistant levels along the pin. Typical results for the Oy components of the bearing forces are shown in Figure 12 below. They are for two initial crack lengths, for the last increment before ligament fracture and after its fracture.
Because the plates were idealized using 2D elements located at the midplane of the plates and because, as shown in Figure 6, the faying surface between the top and bottom plates is between levels 4 and 5 of the pin, the bearing loads at the bottom of the hole in plate 1 needed to be multiplied by a factor of approx. 1.5.
This was necessary in order to account for fact that the contact area in the plate is thinner than five pin levels 7 . A similar approach was applied for the distribution of the bearing loads in lower plate.
After the above correlation between the bearing loads in the pins and those at the hole surface, the results showed that: • the bearing loads at the hole increase continuously from the fastener head (or from the collar) towards the separation surface. The maximum values are at the faying surface (after applying the factor for the transfer of the pin loads to the 2D plate) regardless of the cracked or uncracked condition. This aspect of the contact loads at pin-hole interfaces is described in more detail in publications like [7] [9] [16]. It is attributed in large portion to secondary bending and, to a lesser degree, to bolt tilting; • the presence of the cracks at the hole or at full length have very little influence on the through thickness distribution of bearing loads (or stresses). At the faying surface they were less than 3% different from the baseline configuration, and followed the same pattern as the fastener shear load (including the "hook effect" after fracture), • at the fastener hole next to the cracked one, because the total shear force increases with the crack length, the maximum bearing loads at the faying surface increase also by around 7% from their values in baseline configuration shown in Table 3, • for the laboratory testing configuration, the distribution of the bearing loads was very similar, but the variations during the crack growth are negligible at both the cracked hole and at the adjacent one (less than 3%).
The main drivers for both the small variations in the through thickness distribution of the bearing loads and for the difference in patterns of the two sets of resultsare the lateral restraints of the pates. A similar conclusion is presented in reference [7] where it is stated that "it was found that in cases of splice plates with un-symmetric reinforcement relative to the fastener holes, the secondary bending was significantly different on the two sides of the same fastener hole" and that the "approach must consider the distribution of contact forces between the plates and between the fastener heads and the plates (due to fastener tilt)" 7 Simple analogies could be made with the weighting factor applied to an RBE3 when transferring a nodal concentrated force as a uniformly distributed load on the edge of a plate.

Effect of Crack Length on the Far Field and Bypass Strip Loads and on the Running Loads
Because for the single shear tension joint presented, the effects of secondary bending are significant, the results are split in two sub-groups: 1) The effects of the cracking configurations and of the crack lengths on the running loads in the membrane 8 , 2) The effects of secondary bending due to the load path eccentricity.
The results have been processed and they are shown in terms of internal loads for the convenience offered in visualization, coordinates transformations and in interpretations by the vectorial quantities.
1) The left side of Figure 14 shows proximately the same distance ahead or after the holes edges. Table 5 and Table   6 show that the strip approximation is adequate for the average running loads in baseline configuration because when the values are divided by the strip width, the result is very close to the uniformly distributed externally applied load (63 lbf/in).
When the crack grows, the running loads at each fastener are different from    Figure 15 and Figure 16. The curves have been plotted in real force and distance units for clarity and to overcome the "change sign" alerts occurring when building envelope summaries of percent variations from baseline configuration.
All curves are relatively smooth and can be easily represented as formulas (polynomial, splines etc.). Conversion to running loads or average stresses involves only divisions by specific constants dependent on hole location relative to the crack and/or the reference value preferred (strip wide or gross average far-field). For non-dimensional or parametric form, a few more runs with different plate gauges and fastener sizes are required.
For by-pass load the curves are not shown because for the first row the bypass curves are almost identical to the far-field curves in Figure 15 and Figure 16, and the bypass values for the second row are almost zero due to the vicinity of plate's end. It is noted that after the cracked hole, there is significant variation in the strip-wide bypass membrane force.
2) The distribution of moments due to secondary bending at the configurations and crack lengths exemplified are shown in Figure 17. The strip-wide approach detailed before provides the averages of the far-field and bypass bending moment ahead and after each fastener hole. All the fastener level variation diagrams are quite flat in all configurations and therefore they are not included in the results shown.
With the two sets of average internal loads within the band, the components in Oy direction for the far-field and bypass stresses can be easily determined at any plate depth using the classical structural analysis formulas.
The variations functions for shear flows follow similar patterns, but the values are small and therefore, the curves are not shown in this presentation.

Effect of Crack Length on the Far Field and Bypass Stresses at Outer Surfaces of the Plates
Because a large proportion of testing is based on the strains collected by strain gauges mounted on the outer surfaces of the parts, the results this section present the approach and the results for the far-field and bypass stresses on the surfaces. The added complexities in processing tensorial quantities were minimized but not quite removed by using only planes parallel with the reference coordinate system (CID 0) for the classical engineering formulations or by using the 2D min-max. invariant formulations (and principal angles) of the stresses.
The post processing was done only for the outer surfaces of the plate 1 (0.063" thick) and are shown for the same configurations exemplified in the running loads formulations, i.e. the averaging used was also the strip wide approach to reduce the variability in the stress fields, and suggests a pragmatic approach for quantification and for calibrations during laboratory testing.   For the direct stresses σ x , and σ y (strictly in CID0) "average" values were constructed by using virtual equivalent forces developed in a strip of unitary thickness that were divided by the cross sectional area of the strip. For the shear stresses τ xy (in CID 0), the min-max principal stresses 9 and the maximum shear stresses, the "average" value was taken as the statistical median value of the these stresses for the 10 elements in each strip. This is because in the early post processing it was found that sometimes: • the differences in the principal angles of any two elements were too big for, at least a rational, if not entirely rigorous, basis of averaging; • the differences between the maximum and the minimum element values were disproportionately large in comparison with the values for each element.
An example of baseline values and outcomes of preliminary checks for the baseline stresses at the inner surface (Z1) in model C are shown in Table 7 and at the outer surface (Z2) in Table 8.
For the baseline configuration, the results of "averagings" are quite good for σ y (marked as Fy in the tables). A quick averaging between the far-field values stresses on inner and outer surface at each hole in the first row gives around 1000 psi, same as the tension stress of the external load, for both configurations.
It was found that the there are no practical differences in these values for the two types of plate lateral restraints.
The effects of secondary bending are though substantial as illustrated also in Figure 18. At the outer surfaces the σ y components of the secondary bending stresses, are 2.5 times higher than the membrane tension stresses, resulting in compressive stresses for most of the upper surface. The secondary bending effects are however relatively easy to estimate from the strain gauges data on the free face once the secondary bending is taken into consideration. To remove the need of lubrication between the faying surfaces, friction data needs to be included in model D.
Results of a very similar nature were found for the other crack configurations.
However, the variations of the surface axial stresses as functions of crack length weregenerally different for each configuration and fastener location both at the inner and outer surfaces. The secondary bending makes the influence of the cracks to be less noticeable than on the curves for fastener loads or the membrane running forces. The diagrams in Figure 19 to Figure 23 illustrate the above results for the far field stress. 9 Calculated inside MSC Patranv19 end extracted as results for post-processing. Small corrections were done, where necessary, for the min principal values, as in this version, the outputs for the 2D min principal stress were truncated to zero when the max principal was large. Similar for max principal values, when the min principal was large. The principal angles were all correct when checked using the stress values in CID 0.  Table 7. Baseline "average" far-field stresses on a 0.455" wide strip ahead of the fasteners for the inner surface of model C Below max-min differences between the elements in the strip.   The ranges of variation for max-min principal stress values are relatively similar to the stress components in Oy.
All the curves of variation can be converted to a non dimensional parametric format relatively easy once the secondary bending is incorporated distinctly and unambiguously.

Discussion
Besides cold working, the parameter with the most influence on the fatigue strength is the fastener fit. A close tolerance fit that makes the fastener interfere with the deformations at the hole, reduces the local stress variations at the hole when subjected varying external loads. When the clearance increases, the inter- The results for the fasteners load transfer function show that regardless of cracking configuration, at the cracked hole the load decreases by up to 20% of its value in non-cracked conditions and that, after the ligament fracture, the fastener continues to transfer that load until new cracks develop and grow. As the load transfer decreases at the cracked hole, the increases at the adjacent holes are up to 10% -12%. It is especially noticeable that after breakage of the free edge ligament, the load at that fastener returns almost to its baseline value (the "hook effect").
These results also showed small but realistic differences in the fastener load The through thickness distribution of the bearing loads/bearing stresses at the fastener holes was found to have maximum values at the faying surface, and this fact is of practical importance for inaccessible single shear joints and for multi layered joints with high eccentricities in the load path, where the crack can grow undetected for long durations. These results are only preliminary until the more in depth work on the effects of clamp-up forces and friction is completed. In [2] it is stated that testing has shown clamping forces to have a pronounced effect on the fatigue life of joints 10 . The clamping forces appear to be also a main factor in the eye-brow cracking around the periphery of fastener heads that was observed and reported for in-service aircraft.
The results for secondary bending have been found to be potentially very significant for the fatigue strength of the joint with secondary bending ratios σ bend /σ tens = 2.5 at a distance of two hole diameters ahead of the holes 11 . As the main scope of these analytical investigations was to establish load transfer functions for the fasteners after cracking occurs, the externally applied stress was un- 10 The use of ovalized "self locking" nuts with standard csk. bolts installed in clearance fit holes produced less than 20% of the normal life of a two row single shear joint, when torqued in a degreased condition compared with the specified lubrication. The cause of the short fatigue life was insufficient clamping.
In another test program (also spectrum testing) studying fatigue crack growth from 1 mm corner flaws made at the edge of fastener holes, the fatigue life increased by a factor of 4 when torqueing up the clearance fit bolts to the specifications. 11 To correlate with the standardized location for strain gauge measurements for secondary bending and for optimum determination of the bypass forces by the integration method (AGARD working group -"Fatigue Rated Fastener Systems").  These programs account for the crack length through a geometry function largely known as "the beta function". Beta functions incorporate the stress variation in the net section at the crack, using as reference a clearly specified and constant far-field stress. In the case of growing cracks in a joint, the reference stress ahead of the hole (the far field stress) is changing as the fastener load is changing, and consequently the bypass load is changing. Some of the parameters that affect the fatigue strength have also a direct influence on the changes in the 12 Tate & Rosenfeld used this approach in 1946 to calculate the bypass force. With adequate knowledge of the baseline stress environment(s) and of the fasteners load transfer, the use of parametrized curves developed analytically (and calibrated by a small number of laboratory tests), may be used to adjust the spectrum stresses at every several blocks during the crack growth and residual strength analysis. The benefit will be in increased inspection intervals for the critical locations, without any reduction in the confidence levels for the results.
In the design stage, they can be used for predictions on the appearance of multi-site fatigue damage at those details.

Summary
The aim of this work was to develop a method that enables a relatively quick and The method proposed has the following characteristics: • It is accounting for the pronounced 3-dimensional distribution of stresses in single shear joints with multiple rows of fasteners; • It is applicable to cycles with bulk stresses in elastic range only; • It accounts for installation clearance and for bolt tilting effects; • It incorporates and quantifies the effects of secondary bending in a manner that is easy to interrogate and to verify by standard tests for fatigue rated joints; • It provides quantified precision information on the averaged loads and stresses of the same nature; • The computations and post processing are highly automated for analysis of different materials, plate thicknesses and fasteners configurations; • All results are verifiable at all stages and the data reduction is shown in user friendly formats for interpretation and for engineering judgement; • When used in conjunctions with calibration data from laboratory testing, it provides reliable data for design and sizing in new programs or for use in the specification of inspection intervals for structural integrity on existing platforms.
In addition to incorporating the effects of friction for easier calibration of the