Mesoscale modelling of extended bearing failure in tension-absorber joints

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

2 applied to "extended bearing failure". In this case, the aim is to follow the bearing response, far beyond peak load, with pin displacements up to 30 mm.
Mechanical joints are key components in fuselage structures and their behaviour has been studied extensively [7][8][9][10][11][12][13][14]. Recently, [15], DLR and Airbus investigated the potential use of joints as energy-absorbing devices in future narrow-body composite aircraft fuselages, to assist with meeting regulatory crashworthiness requirements. Due to the reduced available space compared to wide-body aircraft, not all the energy in a "foreseeable survivable impact event", such as a 30 ft/s (9.14 m/s) vertical drop onto a hard surface, can be absorbed by sub-floor crush beams, so additional mechanisms are needed. The DLR/Airbus "tension-absorber" concept, illustrated in Figure 1(a), involves the modification of joints in areas such as the cargo and passenger crossbeams, which are loaded in tension as the fuselage deforms during impact. The modified joints would behave like normal joints under in-service loads, but in a crash would absorb considerable energy. The key design requirement is to prevent bolt pull-through or fracture, so that "extended bearing failure" occurs, resulting in the absorption of energy through crushing of the material in front of the bolt.  [15,16], (b) pin joint "extended bearing failure" tests for parameter investigation: rig schematic, image at end of test, and typical force-displacement curve [17].
To examine the effects of individual material and geometric parameters, Airbus and DLR have studied a simplified version of the problem, namely a pin being pulled through a composite plate [18][19][20][21]. Recently [17], the current authors used this configuration to study the effects of pin diameter, laminate thickness and stacking sequence on bearing strength and energy absorption at quasi-static loading rates. The test rig used, and global response obtained, are illustrated in Figure 1(b). From a modelling standpoint, this recent study, [17], has the advantage that it employed IM7/8552 carbon fibre/epoxy, a material used in the third world-wide composites failure exercise [22], which has been extensively characterised in the literature [23][24][25][26]. Computed Tomography (CT) scans of interrupted tests were also performed, providing a detailed mapping of the internal damage and failure. Thus the results in [17] provide a useful dataset for testing numerical models on a very challenging problem.
In the current work an advanced 3D FE approach is developed for modelling the extended bearing failure mode illustrated in Figure 1(b). The aim is to predict the ultimate bearing strength (UBS), the mean crushing stress (MCS) after bearing failure, and the mass-specific energy absorption (SEA). The model needs to follow the joint response far beyond the peak load and be capable of predicting the change in response when material and geometric parameters are varied. In order to do so, it must represent the physical internal damage processes as faithfully as possible. It also has to remain numerically stable, at extreme levels of material damage, while adjusting contact conditions between the pin and the laminate, and within the laminate itself, as the simulation progresses. In line with several recent studies on composite joint modelling, [6,33,34], these stability requirements lead to the use of an explicit solver to avoid the convergence issues that plague implicit methods [40], even though the loading rate is quasi-static.
The approach taken is to extend the physically-based damage model of Egan et al. [10,36,44], which addressed bearing failure in normal (not tension-absorbing) countersunk bolted joints. That model included Puck"s criteria for matrix damage [45], a nonlinear law for in-plane shear, a maximum stress criterion for fibre failure, a crack band model to mitigate mesh sensitivity [46], and frictional contact between the fastener and the laminate. In the present work, the model is extended through the use of in-situ ply strengths, stress-based fibre failure criteria that incorporate fibre kinking effects in an efficient manner [47], a cohesive zone model (CZM) for simulating delamination, frictional contact between adjacent plies after they delaminate, and an efficient search algorithm for matrix damage identification. The model predictions are compared in detail to the test results in [17], at both the global and mesoscopic levels, and the fidelity and capability of the model to forecast the effects of variations in joint parameters are assessed.

Experimental set-up
A brief description of the experimental set-up is provided here, with full details available in [17]. The material used was HexPly ® IM7/8552 (EU version: 134 gsm) unidirectional, continuous-fibre carbon fibre/epoxy composite prepreg (nominal ply thickness = 0.125 mm). Specimens with the geometry shown in Figure 2 were extracted from panels manufactured in an autoclave using water jet cutting. As outlined in Table   1, fifteen configurations were tested. All layups were quasi-isotropic, but the stacking sequence varied, as did the pin diameter and laminate thickness. The 2 mm and 3 mm thick specimens are labelled "dispersed" or Each test consisted of pulling a steel pin through the laminate at quasi-static loading speed (1.67x10 -4 m/s), using the specialised rig illustrated in Figure 1(b). The reaction force was obtained by the load cell of the servo-4 hydraulic test machine and the pin displacement was measured by a digital image correlation method. Four repeats of each configuration were tested, with further interrupted tests undertaken for CT analysis. Figure 2 Composite specimen dimensions. Thicknesses tested were 1 mm, 2 mm and 3 mm. A typical response, illustrated in Figure 1(b), exhibited a peak force at bearing failure, followed by a load drop for 5-10 mm, and a transition to a relatively constant crushing force which persisted until the pin exited the end of the laminate. To allow comparisons on a material level, the following performance parameters are defined. The ultimate bearing strength (UBS), ult  is defined in accordance with ASTM standard D 5961/D 5961M [48]: where D is the pin diameter and t is laminate thickness. The mean crushing stress (MCS) is: where F mean is calculated between 5 mm and 30 mm pin displacement. Finally the mass-specific energy absorption (SEA) is defined as: where is the mass of the material involved in energy absorption, ρ is material density and s m is taken here as 30 mm. In [17], s m was taken to be 40 mm, but the shorter distance is used here for comparison with the models, due to considerations of computation time. The difference between the SEA calculated over 30 mm or 40 mm is very small and running the simulations beyond 30 mm provides little additional value. The factor of 1.2 in equation (3) is based on an estimation in [19] that, for materials with brittle fibres, the width of destroyed material is 20% larger than the pin diameter.

3D elastic behaviour and nonlinear shear law
A unidirectional (UD) ply damage model has been implemented in an Abaqus/Explicit VUMAT subroutine.
The damage variable 12 d controls the reduction of the original shear modulus ( 0 12 G ) due to progressive matrix damage occurring under shear loading, Figure 3. Following the method of Donadon et al. [50], it varies with 12  according to equation (8), where a is the slope of 0 12 12 GG versus 12  (the "gradual stiffness reduction curve"), which is determined experimentally. 12 12 da   (8) Inelastic strain is determined as 12 12 12 12 in e ed        and in-plane shear stress is updated as: The parameter 12 f  defines the shape of the nonlinear response. The IM7/8552 carbon fibre/epoxy used here shows a response similar to that shown in Figure 3, and the curve is fitted using two cubic polynomials, one applicable up to a shear strain of 1,max P  (see Figure 3), the other applicable above 1,max P  , i.e.: 12 32 12 1,max 1 12 2 12 3 12 32 12 1,max 1 12 2 12 where i c and i d (i = 1, 2, 3) are the coefficients of the fitted curve and 4 d is the shear stress at 1,max A feature of the model, previously developed by Egan et al. [44], is a purely symmetric shear stress-strain law on load reversal. This aspect is discussed in detail in [44]. , 1 12 f  and , 2 12 f  . 12  is maximum shear strain over time, e  , ed  and in  are elastic, elastic-damage and inelastic parts, respectively.

Fibre failure
In line with previous researchers [51][52][53][54][55][56], for tensile fibre failure, a simple maximum stress criterion is used: 11 11 For 0 : 1 where t X is the tensile fibre-direction strength. For compressive fibre failure, some researchers use a fibrekinking model [6,[57][58][59], on the basis that longitudinal compressive failure is caused by the formation of fibre kink bands following degradation of the supporting matrix. The model of Pinho et.al [59] for example, successfully predicted the increase in compressive strength due to hydrostatic pressure shown experimentally in [60]. Fibre-kinking models are however, computationally expensive, and Egan et al. [44] justified using a simple maximum stress criterion instead ( compression test would be lower than the in-situ compressive strength in the torqued-up countersunk joints in [44], so predictions would be conservative. Unfortunately in the current pin-loaded specimens, no torque exists, so this argument cannot be made. However, Raimondo et.al [47] derived a fibre compression failure criterion from polynomial fitting of experimental failure envelopes. They proposed that longitudinal and shear stresses contribute to shear fracture of fibres, followed by fibre rotation and in-plane matrix shearing at the crack tip, which in turn promotes kink band development [61]. The failure criteria for quasi-static loading was given as: , 1 12 f   S are the quasi-static compressive fibre, in-plane shear and out-of-plane shear strengths respectively. This fitted experimental data well for 1   . Given the high level of detail in the other parts of our model, and the computational expense of including an explicit fibre kinking model, the approach of using equation (12) was followed instead ( is selected here to be 2.5). Although this is a compromise it is still an improvement over the criterion used in [44].

Matrix failure
To predict matrix failure, the stress tensor is rotated about the 1-direction by variable angle  , 0   , as shown in Figure 4 where T Y is the transverse tensile strength, 12 S is the longitudinal shear strength, and 23 S is the transverse shear strength. The transverse friction coefficient t  is defined using Mohr-Coulomb theory as: where 0  is the fracture surface orientation for pure transverse compressive failure which typically has a value of about 53 for unidirectional polymer matrix composites. The longitudinal friction coefficient is given by: The angle at which failure is detected,  , determines the orientation of the fracture surface. This search process can be computationally expensive, so to efficiently find the fracture surface a "Golden Search Algorithm", as outlined in [62], was implemented. (a)

Crack band model and definition of characteristic lengths
To mitigate mesh sensitivity, a crack-band model [46] where G is the fracture energy and o  is the stress at failure onset. Damage variable growth is described by [59]: (17) which causes nonlinear softening of stress components, as illustrated in Figure 4 where x = max(0, ) x and stress and strain components are obtained from the stress tensor rotated by the fracture plane angle,  . As in [59], the strain tensor prior to rotation contains elastic in-plane shear strain. This ensures that only elastic internal energy contributes to energy absorption associated with localised fracture. The strain measure used to grow the matrix damage variable ( t t m   ) after the damage onset is given in equations (20) and (21), depending on whether failure is tensile or compressive.
where Ic G and IIc G are the transverse fracture energies, and the stress components are the tractions at the onset of matrix failure.
The characteristic element length is provided by ABAQUS to VUMAT routines as the variable "charLength". However, this value is only appropriate for models with cubic elements in which fracture is perpendicular to an element side. In simulations of bearing failure, near-hole elements are inevitably of noncubic shape and non-standard orientation, while the ply orientations and Mohr-Coulomb material behaviour promote cracks which are angled with respect to element sides. Hence instead of using ""charLength"" an approach developed by Egan et al. [44] is used. A python script pre-processes the mesh to compute element geometries from nodal coordinates. These are used to compute characteristic lengths, f L and m L , which accurately account for crack orientation, element orientation and element shape. Full details are given in [44].

Damage enforcement and element deletion
When fibre failure is detected, all stress tensor components are softened according to equation (27), on the assumption that fibre rupture damages the supporting matrix. The resulting Cauchy stress tensor, tt Excessive element distortion can cause explicit simulations to abort, so element deletion criteria are defined in equations (28) and (29). After extensive trial simulations, 1  and 2  were set to be 0.8.

Interlaminar damage
The uncoupled, mixed-mode cohesive zone model (CZM), already implemented in ABAQUS, is used to model interlaminar failure [63]. The CZM follows a bi-linear traction-separation law with linear softening for opening (mode I) and shearing (mode II and III) modes, as shown in Figure 5(a) and (b). The initial linear behaviour is defined by: where the Macauley bracket  is used to ensure that a purely compressive stress state does not initiate damage,  [65] is used to propagate the delamination: c  c  c  c  II  III  I  II  I  cc  II  I c II G and c III G are the critical fracture energies in the normal and shear directions, while  is a cohesive property parameter.

Input parameters for the VUMAT
UD carbon fibre/epoxy IM7/8552 has been used by many researchers and has been fully characterised at various loading rates [23][24][25][26][66][67][68][69]. The elastic and damage properties at quasi-static loading rates are summarised in Table 2 to Table 7. In [6] it was proposed that in-situ ply strengths should be used in the simulation of pin-bearing problems, and comparisons with experiments validated this approach. As shown in Table 5, transverse and shear strengths were made dependent on the location and thickness of a ply e.g. thin outer ply (to), thin embedded ply (te) or thick ply (thick).     Table 5 IM7/8552 in-situ shear and transverse tension strengths [70,71]. Single plies are "thin", blocked plies (two or more plies with same orientation) are "thick". Strengths depend on whether the ply is a thin outer ply (to), a thin embedded ply (te), or a thick ply (thick).

Model details
The pin-bearing tests were modelled in ABAQUS/Explicit. Meshing was performed as shown in Figure 6, using a python script for consistency among the various pin diameters and laminate thicknesses. Three-dimensional, reduced integration, 8-noded solid elements (C3D8R) were used. In displacement-based non-linear FE models, the use of full integration elements over-estimates the stiffness matrix and the elements may suffer from shear locking. To mitigate this effect, the use of reduced integration elements (C3D8R) is recommended. Furthermore, using these elements lowers the computational cost [63]. One downside of using reduced integration elements is that a fine mesh must be used near structural boundaries to capture the stress concentrations at such locations [63,76]. However, this disadvantage was more than compensated for by the advantages listed above for the current problem. These elements are also prone to hourglassing due to the presence of several spurious zeroenergy modes, so to ameliorate these effects "pure stiffness" section control was applied [63]. For interface elements, an 8-noded three-dimensional cohesive elements (COH3D8) were used, which had a thickness of six microns. This element type is typically used to capture delamination behaviour in composite structures [6,77].
With these elements, it is assumed that the cohesive layers are only subjected to direct through-thickness ( 33  ) strains and transverse shear strains ( 23 13 ,  ). The other two direct strains, 11 22 ,  and shear strain 12  are assumed to be zero. These assumptions are appropriate in situations where a relatively thin and compliant cohesive layer bonds two relatively rigid parts [63].
Each ply was modelled with a single element through the thickness. To accurately predict bearing strength, a variable element size was used over the first 2.5 mm of crushing length (from 0.3125 mm × 0.2 mm to 0.3125 mm × 0.5 mm), see Figure 6. The rest of the crushing length had a constant element size (0.5 mm × 0.5 mm), as can be seen in Figure 6. The reason for selecting such a fine mesh near the hole edge is to capture the stress concentration at the pin-hole contact effectively using the chosen reduced integration elements. In the material model, the crack band model is implemented to ameliorate mesh sensitivity, as discussed in section 3.4. Elastic, isotropic material properties were used for the hardened steel pins (E = 210 GPa, ν = 0.25), while the laminate was modelled with elastic properties in the grip region, where no damage was expected, and via the VUMAT described in the previous sections, in regions predicted to damage, see Figure 7(a). Non-linear geometrical effects were included due to the finite strains developed during crushing of the elements.

13
Contact between the laminate and pin was defined using the "general contact" algorithm [63] with the choice of friction co-efficient discussed in section 4.1. Continuous laminate crushing results in deletion of elements, exposing new surfaces for contact with the pin. To ensure continuity of the simulation, the contact surface needs to be updated by creating interior surfaces. ABAQUS/CAE does not currently support the creation of interior surfaces, so the input files had to be modified manually. Firstly an element set was defined containing elements lying in the crushing zone, see Figure 6, then interior surfaces were created using the *SURFACE command. Additional contact pairs were then defined in the ABAQUS input file using these surfaces. The leftmost (as per Figure 7(b)) 45 mm of the laminate was given a fixed boundary condition (imitating the clamps of the fixed end of the test machine), while a velocity boundary condition was applied to the centreline of the pin (see Figure 7(b)). To aid in achieving an efficient quasi-static solution, and to provide an accurate and repeatable measure of the peak force and bearing strength, a sigmoid function was used to increase the velocity gradually from zero to the applied velocity. After this, the velocity was kept constant. 14 cohesive elements and one element per ply in the thickness direction. It was decided not to use mass scaling based on past experience. To examine the feasibility of using an increased velocity, a series of simulations were performed and the kinetic to internal energy ratio was examined until the peak load, as this is important for quasi-static analysis with an explicit solver [78]. Three loading rates were trialled on a BK_D4_T2 configuration (see Table 1), 0.1 m/s, 0.5 m/s and 1 m/s. Note that no strain-rate dependency was included in the material model in the VUMAT, so any change in behaviour is not due to material properties. The numerical reaction force was calculated from the fixed end of the laminate while the displacement was extracted from the centreline of the pin. Reaction force is plotted against pin displacement for the different loading rates in Figure   8(a), and the ratio of kinetic to internal energy in each model is given in Figure 8(b). As a rule of thumb, the kinetic energy should not exceed 5-10% of the internal energy for a quasi-static solution [63]. In the 1 m/s model, the ratio of kinetic-to-internal energy exceeds 15% early in the analysis, indicating a quasi-static solution is not achieved. The load-deflection response for the other loading rates converge reasonably well onto one response, Figure 8(a), and the energy ratio is well within bounds. Based on this trial study, 0.5 m/s was used in all the quasi-static analyses, as it gave good agreement with experimental data without showing an detrimental increase in kinetic-to-internal energy ratio, and it was also computationally efficient. Simulations were run on a single node of a high-performance computing system, where each node has 2 x 20-core 2.4 GHz Intel Xeon Gold 6148 (Skylake) processors.

Model calibration
In this section, the model is calibrated in three iterations. In iteration 1, the contact conditions and friction coefficients are chosen. In iteration 2, a deeper analysis of the fracture energies used for the cohesive elements between plies is undertaken. In iteration 3, a simplification is made in the interests of saving computational time.
The experimental study, [17], showed that delamination plays an important role in the joint response.
Delamination is modelled here with the in-built ABAQUS bi-linear traction-separation law, as explained in Section 3.6. Now, as the bearing load increases, cohesive elements between the plies start to delete. Once that happens, if contact is not defined between the plies, they are free to pass through each other, which is nonphysical. Thus, we defined contact between plies via interior contact between ply elements. Having done this, friction coefficients between the pin and the laminate, and between the plies themselves, have to be chosen.
Iteration 1 addresses these issues, with configuration BK_D4_T2 used as a sample case.  (1). The following observations can be made: (i) Figure 9(a) shows that UBS increases with increasing friction coefficient between the pin and laminate, (ii) Figure 9 (iii) Figure 9(c) shows the case where interply contact has been defined, and (iv) Figure 9(d) shows the longitudinal section when interply friction is defined. Penetration of the ply elements is avoided, and a reasonable damage morphology is obtained. The second iteration addresses the under-prediction of MCS. It was decided to look more closely at the selection of fracture toughness values (G Ic and G IIc ) for the cohesive elements. Up to this point, all cohesive elements were assigned the material properties of a 0/0 interface, see Table 7. However, it is known from the literature for carbon fibre/epoxy systems [80,81], that mode II fracture toughness for +θ/-θ and 0/θ interfaces, where θ is an arbitrary ply orientation other than 0°, are higher than for a 0/0 interface. The effect of interface ply orientations on mode I fracture toughness is relatively small. Table 8 summarises G Ic and G IIc values found in the literature for IM7/8552, for the various interfaces in our models. As can be seen the G IIc values for the non-0/0 interfaces are higher than for the 0/0 interface. Using the 0/0 value for all interfaces in the model (as done so far) is likely to lead to early deletion of cohesive elements, when they should remain in the model, thus reducing the bending stiffness of the laminate and the resistance to pin movement. Thus, the FE model was modified to include the appropriate toughness values for each interface. Figure 10 shows the dramatic improvement this made in the prediction of the MCS. Meanwhile the already good prediction of UBS is virtually unaffected.  In the experiments, [17], it was found that most delaminations occurred between plies of different orientation, with very few existing within blocks of similarlyoriented plies. So to save computational cost, it is reasonable to place cohesive elements only between plies of different orientation, which reduces the number of delamination sites in 2 mm thick laminates from 15 to seven, and in 3 mm thick laminates from 23 to seven. This approach has also been used recently in [83], and it reduced the CPU time for our models from 98 hours to 43 hours for the 3 mm thick specimens. Figure 11  Experimental average UBS Experimental average MCS 0/0 interface properties everywhere Different properties at each interface load) were always between blocks of similarly-oriented plies, whereas later in the crushing process, some delaminations appeared within such ply blocks. Overall, the difference in MCS values due to use a reduced number of cohesive layers is relatively small, so this approach was used for the remainder of the paper.
In summary, the model calibration consisted of choosing appropriate values for friction coefficients, correct G IIc values at each interface, and a simplification to reduce computation time for blocked laminates. No other calibration or tuning was applied in the remainder of this study.

Prediction of global response, including the effects of varying joint parameters
In this section, the model is tested for its ability to predict the global response, through comparison with experimental results in [17]. In Figure 12(a), the movement of the pin, after 30 mm displacement, is shown for the model and the experiment for the BK_D4_T2 configuration. Because the tension rods in the experiment (shaded gold in Figure 1(b)) are only pinned at one end, they are free to rotate. Consequently the pin is free to follow the path of least resistance, and it sometimes veers away from straight-line motion, as shown in the example in Figure 12(a). The pin is also not restricted from sideways motion in the model, and it can be seen that the pin also veers away from straight-line motion. The exact direction of travel is unpredictable, as it depends very sensitively on the sequence of damage events in the laminate. One thing the model does not fully capture is the peeling off of surface 45° plies which can be seen in the experiment in Figure 12(a) and is discussed in [17]. However, this peeling process has a relatively minor effect on energy absorption. In Figure   12(b), the damage from a 4 mm pin is compared to that of a 12 mm pin. It is noticeable that the composite material on either side, and in front of, the pin, splays out-of-plane to a greater extent for the 4 mm pin than for the 12 mm pin. This will be discussed further below. To assess the capability of the model to predict the effects of varying pin diameter, laminate thickness and stacking sequence, all test configurations in Table 1 were simulated. Table 9 shows the predicted values of UBS, The predictions for MCS and SEA are not as good, which is not surprising, given the extremely complex damage and failure events that occur during crushing over such a long distance. Nonetheless, eight and 10 configurations are predicted within 10% of the experimental mean, for MCS and SEA respectively, with the remainder being within 20%. The experimental standard deviation for MCS and SEA is seen to be up to 14.8%.
The experimental and numerical rankings for SEA again differ by no more than two for any configuration. The good prediction of UBS and SEA for blocked laminates suggests the successful representation of the in-situ effect for thick plies (Table 5). Figure 13 shows a selection of experimental and numerical stress-displacement curves. The experimental curves are just one sample from among four repeat tests, so one should not expect exact agreement with the simulation. Figure 13(a) shows that the model correctly predicts that increasing pin diameter leads to decreased MCS/SEA. Similarly, Figure 13(b) shows that the model correctly predicts that increased thickness leads to increased MCS/SEA, and Figure 13(c) shows that, for 12 mm diameter pins, the model forecasts that a dispersed stacking sequence results in a higher MCS/SEA than a blocked stacking sequence, in accordance with the experiments. The reasons why pin diameter, laminate thickness and stacking sequence have these effects on the performance parameters are discussed in detail in [17].
One of the major findings in [17] was the strong correlation between UBS, MCS and SEA and the ratio of pin diameter to laminate thickness, D/t. In Figure 14, the predictions for the relationships between UBS, MCS and SEA and D/t ratio are plotted together with the experimental results, for dispersed stacking sequences (on the left) and blocked stacking sequences (on the right). Results for 1 mm thick laminates are included in the dispersed stacking sequence plots. The predictions of UBS, MCS and SEA are seen to be excellent for dispersed stacking sequences, while the prediction of UBS versus D/t is very good for blocked stacking sequences. The trend line equation for MCS (and consequently SEA) versus D/t is somewhat off for blocked laminates, but the overall direction of the trend is more or less correct. Possibly, dispensing with the simplification in iteration 3 above, i.e. placing cohesive elements at all ply interfaces for blocked laminates, might give better results.
Overall though, the model is shown to be capable of predicting the global response and performance parameters, UBS, MCS and SEA over a wide variation of geometric and material parameters. Table 9 UBS, MCS and SEA mean experimental values versus model predictions, with percentage that the model value is above or below the experimental mean. Also shown are the relative standard deviation, RTSD, (standard deviation as a percentage of mean) in the experiments, and the rankings (1 -15) of the configurations, found experimentally and numerically.

Validation of mesoscopic response
In [17], tests were performed in which the pin displacement was halted at 0.75 mm (about 0.3 -0.4 mm beyond peak load). The specimens were then scanned using three-dimensional computed tomography (3D CT), details of which can be found in [17]. In Figure 15 As noted in [17], in the experiments, Figure 15(a), (c), (e) and (g), the outer plies (red labels) delaminate and splay outwards, while the centre plies (black labels) stay more or less aligned with the load and undergo crushing. Blocks of apparently undelaminated central plies at a distance of one laminate thickness (2 mm) from the hole edge are indicated. The size of these blocks and their degree of alignment with the loading direction were found to be good indicators of SEA in [17]. The alignment of the central blocks with the loading direction depends on the lateral support provided by the outer plies. Deeper delaminations tend to weaken this support.
In the simulations, Figure 15(b), (d), (f) and (h), the overall shape of the predicted deformation is correct.
Outer plies delaminate and splay outwards, while the centre plies stay more or less aligned with the load.
Furthermore, for the 4 mm pin, Figure 15   and (f), in both the experiments and the simulations. These differences are due to the difference in width (4 mm versus 8 mm) of the crush zone in front of the pin. This zone tears away from the rest of the laminate, and the outer plies delaminate and buckle outwards. The support against ply buckling is less for the wider strip (larger pin), so the delaminations extend further into the laminate. As a consequence of these deeper delaminations, the support provided by the outer plies to the central plies is less for the larger pin, which allows the central plies to bend more easily out of the way of the oncoming pin. As discussed in [17], this is one of the key reasons why the performance parameters drop off with increasing pin diameter, and the model captures this effect well. Another observation is that the maximum delamination depth is greater for the blocked laminates, Figure   15(e), (f), (g) and (h), than for the dispersed laminates, Figure 15(a), (b), (c) and (d), which is a contributing factor towards the lower performance of the blocked laminates. Again this predicted by the simulations.
In Figure 16, the stress state in the cohesive elements in front of the pin is examined in detail for the DS_D4_T2 case. Figure 16(a) shows the location at which cohesive zone tractions were extracted (first cohesive element in front of the pin in each layer) and the coordinate system used. The normal traction, 33  , relates to opening Mode I when positive, while 13 31 /  relate to sliding Mode II, and 23 32 /  relate to tearing Mode III.

 
(see Table 7) is shown as a vertical dashed line on the shear traction plots.
far exceeds one at interfaces 1 and 15 indicating delamination initiation. Additionally though, it also slightly exceeds one at interfaces 5 and 11 (also 45°/-45° interfaces). Thus at peak load, delamination initiates at interfaces 1 and 15, followed by interfaces 5 and 11. Sliding Mode II is the biggest contributor at both sets of interfaces, followed by tearing Mode III. Opening Mode I plays no part.
After a 10% drop in load (post-peak), the tractions have changed significantly. 33  is still negative for all plies so still plays no part in delamination. 13  is now largest at interfaces 3 and 13, indicating Mode II-driven delamination at these 0°/90° interfaces, while 23  is highest at interfaces 5 and 11. The bottom-right figure shows that most of the outer plies have delaminated (note that in some cases, the delamination initiation criterion value has fallen back to less than it was at peak load, but this does not indicate the absence of delamination, since delamination is non-recoverable). A centre block of five plies (90/0 2 /90/-45) is undelaminated, which correlates quite well with the CT scan in Figure 15(a).
In Figure 16(c), compressive fibre damage is shown at peak load, and at 10% drop in peak load (post-peak).
It can be seen that at peak load, fibre compressive damage is already well under way, particularly in the 0° plies, but also to a lesser extent in the ±45° plies. However, none of the elements are shaded red, indicating that no elements have completely failed due to fibre compressive damage yet. By the time the load drops by 10%, multiple elements are indicating complete fibre damage failure. The first drop in load from its peak value thus appears to be triggered by delamination (not the initiation of fibre damage), which most likely then accelerates the accumulation of fibre damage.  Fibre compressive damage at peak load  Figure 17 shows a similar set of plots for a blocked stacking sequence, BK_D4_T2. Once again, the normal tractions are negative for all plies at both load levels, so delamination is entirely shear-driven. This time, the tractions do not change nearly as dramatically between peak load and post-peak (10% load drop), as they did for the dispersed stacking sequence. Delamination occurs at interfaces 2 and 14 (45°/-45° interfaces) and interfaces 4 and 12 (-45°/90° interfaces). No cohesive elements exist between plies of the same orientation, so no delamination is possible at odd-numbered interfaces. A central block of eight undelaminated plies (90 2 /0 4 /90 2 ) is predicted after 10% load drop (bottom right plot), which correlates exactly with the CT scan in Figure 15(e).
Once again fibre compressive damage has initiated before the peak load, but no elements have fully failed due to fibre damage at peak load. So the initial drop in load is delamination driven.   Figure 16 for BK_D4_T2 case, at peak load and 10% load drop (after peak load), (b) fibre compressive damage at peak load and at 10% load drop In Figure 18(a) and (b) scans of the bearing plane are shown for 3 mm thick laminates, tested with a 12 mm diameter pin until 0.75 mm displacement, for dispersed (DS_D12_T3) and blocked (BK_DS_T3) stacking sequences respectively. The dispersed laminate contains a large number of delaminations, extending up to 6 mm into the laminate, and leaving a central block of about nine undelaminated plies at one laminate thickness from the hole edge. The blocked laminate separates at the interfaces between blocks of three similarly oriented plies, resulting in fewer but deeper delaminations (up to 9 mm), and central block of only six undelaminated plies. The blocked laminate has significantly lower UBS and SEA than the dispersed laminate (see Table 9) due to the deeper delaminations and the smaller block of undelaminated plies. Below the CT scans are plots of fibre compression damage at 0.25 mm, 0.5 mm and 0.75 mm pin displacement, and also a 3D view. Fibre compression damage begins in the 0° plies, but spreads to all plies eventually. The predicted deformation agrees well with the experiments, and in the 3D views, it can be seen that fibre damage in the central plies spreads more widely to either side of the hole for the blocked laminate than for the dispersed laminate.
In Figure 19, a view of the pin-laminate contact region shows the progression of fibre compression damage (on the left) and fibre matrix damage (on the right) for the DS_D8_T2 configuration, with a comparison with a CT scan. It can be observed that at 0.75 mm pin displacement, the FE model shows good agreement with the CT scan in terms of the spread of damage and level of brooming of the outer plies.

Conclusions
A 3D finite element approach, incorporating a mesoscale damage model, has been developed for the prediction of extended bearing failure, which occurs in tension-absorber joints. The model is compared with experimental data covering a wide range of parameter variations, and very good agreement is found at both the global and mesoscopic scale. Using material parameters for IM7/8552 readily available in the literature, the model is capable, without any special tuning, of predicting the effects of variable pin diameter, laminate thickness and stacking sequence on the key performance parameters: ultimate bearing strength, mean crushing force and specific energy absorption. The model was able to rank 15 different joint configurations in terms of performance with a high degree of accuracy. Predictions of bearing strength were mostly within 5%, and all within 12%, of the experimental mean. Predictions of specific energy absorption were mostly within 10%, and all within 20%, of the experimental mean. The meso-scale response also closely followed the experiments, as shown by comparison with CT scans. The general shape of the ply deformations, depth of delaminations, size of the central block of undelaminated plies, and evolution of fibre and matrix damage, matched the experiments well. Overall, the model is found to be genuinely predictive for an extremely challenging problem. Although the model is computationally expensive, the approach will become more practical as computer hardware continues to develop.
In future work, the model will be extended to include strain-rate dependent material properties and tested on dynamic pin-bearing tests which are currently under review for publication. In addition, another test series involving non-quasi-isotropic layups has also been performed, and the model will be tested on those layups. A numerical optimisation study will then be performed to find the highest performing stacking sequences for tension-absorber joints.