Numerical simulation of plasticity in double lap metallic bolted joints

Abstract This article aims to compare three plasticity models, namely Bai-Wierzbicki (BW), von Mises and Drucker Prager, to predict the failure modes in double-lap metallic bolted joints under shear loading. The double-lap bolted joint finite element model has been validated against literature data, which demonstrated the application of BW ductile fracture model to bolt joint failure prediction. Taguchi Design is used to investigate the effects of geometrical parameters, preload and stress triaxiality on the failure modes. Force–displacement response amongst the three models differs up to 20% when failing in shear and a different prediction of the failure mode (shear out vs. net section) under the same loading and geometry conditions is observed.


Introduction
Metallic bolted connections are widely adopted in the automotive industries since they allow disassembly and reassembly of structures.Although they are easy to inspect and repair, bolted joints add weight and are prone to corrosion.Double lap joints are preferred to single lap joints as they avoid secondary bending which promotes plasticity and failure.The failure modes induced in bolted joints are conventionally classified as: net tension, shear-out, bearing, cleavage and pullthrough failure.There are three additional failure modes related to the bolts which are head failure, shank failure and collar failure.Various factors such as joint geometry, preload, materials and loading conditions dictate which failure mode is triggered.Due to stress concentrations, material and geometrical non-linearities the formulation of reliable and accurate closed-form solutions to perform their structural assessment is challenging.
There have been many investigations on the design methods for bolted joints under static loads (Da Silva 2008).These studies have focused on notch effects (Wang, Hung, and Chang 1996), preload (Oskouei and Chakherlou 2009), fastener types (Oskouei, Keikhosravy, and Soutis 2010) and clearance (McCarthy and McCarthy 2005).It is well-known that drilling a hole in a flange causes high stress concentrations near the hole and reduces the strength of the bolted connection.Drilling can also cause a rough surface finish, which may become the site for nucleation and propagation of fatigue cracks (Chakherlou, Alvandi-Tabrizi, and Kiani 2011).Studies have shown that preload (Croccolo, De Agostinis, and Olmi 2013), cold expansion (Chakherlou and Vogwell 2003) or interference fit (Raju et al. 2016) can help to mitigate the notch effect and increase the joint's fracture strength.Keikhosravy et al. (2012) investigated the effect of the geometry on the stress and strain distribution in AA2024-T3 single-lap bolted joints.The failure mode changed from net-tension to bearing when the distance of the hole from the plate edge increased and changed from shear-out to bearing when the width of the plate increased.Kim and Yura (1999) performed an experimental study on the influence of distance between the bolt and the edge on the bearing strength of a double lap bolted joint.The results confirmed that the bearing strength is proportional to the material UTS and the joint can bear higher deformation when the edge distance is increased.Kim et al. (2007) proposed four different FEM techniques to investigate bolted connections: solid bolt, coupled bolt, spider bolt and no-bolt.The solid bolt model with surface-to-surface contact was the most accurate compared to the experimental tests.Ju et al. (2004) studied the structural behavior of the butt-type steel bolted joints.It was concluded that plate deformation and bolt bending effects could be neglected to calculate the ultimate load for the bolt shear failure of the connection.Chakherlou et al. (2009) developed a numerical model of AA7075-T6 bolted plates to estimate stresses and strains under preload and external loads.The results of this study showed that higher preloads decrease the magnitude of the stresses around the hole when subjected to longitudinal tensile loads.Chakherlou et al. (2008) studied with a combined numerical and experimental approach the beneficial effects of the residual compressive stresses induced by a radial preload on the hole, which reduces the magnitude of the stresses generated by tensile loads and improves the fatigue life of the bolted joint.
The common approach adopted to implement finite element model of single and double lap joints uses the classical von Mises or Tresca theory to describe the plastic response of ductile metals.Although these theories may be successful in applications where the plasticity experienced before failure is localized, it has been demonstrated that fracture initiation is often preceded by large plastic deformation with considerable stress and strain gradients around the point of fracture (Chen et al. 2011).
The plastic response and ductile fracture in metals are greatly affected by the stress-state.French and Weinrich (1973) experimentally demonstrated that hydrostatic pressure greatly affected the strain to failure under tension.Hancock and Mackenzie (1976) investigated the connection between the ductility and the stress triaxiality for different steels.They found that the ductility decreased with the stress triaxiality.Richmond and Spitzig (1980) conducted an extensive experimental study on Aluminum alloys, which demonstrated the effect of hydrostatic pressure on yielding.
Extensive theoretical studies on ductile fracture models have been carried out since the 1960s, when McClintock (1968) and Rice and Tracey (1969) discussed the effect of stress triaxiality on void growth rate.Gurson (1977) proposed a porous plasticity model that describes the concept of a shrinking yield surface, showing the vital role of hydrostatic pressure in plastic yield and void growth.This model has been modified by others to include additional dependencies such as void nucleation (Needleman and Tvergaard 1984) (Pardoen and Hutchinson 2000), enhanced strain hardening models (Holte et al. 2019) or plastic anisotropy (Grange, Besson, and Andrieu 2000).
Johnson and Cook developed a ductile fracture criterion (Johnson and Cook 1985) using the results of dynamic tests at different strain rates and temperatures.The experimental results demonstrated that stress triaxiality on the ductility of their specimens was more significant than strain rate and temperature.
Wierzbicki's group (Xue 2007;Wierzbicki et al. 2005) observed that triaxiality was not the only factor that affected ductile failure.As a result of testing several specimens under different loading conditions, it was observed that the torsion/shear test, which has lower triaxiality than the tensile test, had lower fracture strain.A non-linear damage rule which characterized the damage accumulation with respect to the plastic strain was proposed.Bao (2003) and Bai and Wierzbicki (2008) conducted extensive investigations on the strain at failure, covering stress triaxiality ranging from compression to multi-axial tension.Tensile, shear and compression tests on 2024-T351 aluminum alloy were performed and it was demonstrated that the amount of plastic strain accumulated before fracture varies over a large range of stress triaxiality (-1/3 to 1.1).Gao and Kim (2006) showed that the Lode angle should be considered to distinguish the stress state with the same triaxiality ratio.Xue and Wierzbicki (2008) postulated a symmetric ductile fracture criterion that considers stress triaxiality and deviatoric state parameter related to the Lode angle.The accuracy of the criterion had been compared against 7 other fracture models and demonstrated to be the most accurate (Wierzbicki et al. 2005).Bai and Wierzbicki (2008) developed a new plasticity model with correction for the hydrostatic pressure and Lode parameter to address the possible asymmetry of the fracture locus.The same authors also transformed the classical Mohr-Coulomb fracture criterion into the space of stress triaxiality, Lode angle and equivalent plastic strain, and proposed a modified version (MMC) (Bai and Wierzbicki 2010).This criterion was later extended by Ebnoether and Mohr (2013) in order to take anisotropic fracture, softening and crack propagation into consideration.
There are other methods in modeling and simulating crack initiation and propagation for mechanical joints, for example the cohesive zone model (CZM) and extended FEM (XFEM) (Mubashar, Ashcroft, and Crocombe 2014;Stein et al. 2017;Dimitri et al. 2021), both of which are based on crack-tip fracture mechanics.This is different from double lap bolted joints presented in the paper, where there are no predeterminate crack paths.
In this work, the failure modes observed in double lap bolted joints are investigated comparing three different plasticity models.It is shown that under different combinations of load and geometry the three models can predict a different failure and strength.These factors are discussed to provide guidelines on the approach to be used in the identification of the most appropriate plasticity model when modeling this type of joint.

Stress state representation and different plasticity/fracture models
A vector in a 3-dimensional space, where the principal stresses are the main axes (Bai 2007), can represent the stress state.
Figures 1-18 show a particular stress state of a point "P:" The hydrostatic axis (z) is the axis in which the three principal stresses are equal r 1 ¼ r 2 ¼ r 3 .The plane containing P and perpendicular to the hydrostatic axis is called Deviatoric Plane.
Figure 1.Three types of coordinate system in the space of principal stresses (Wilson 2002).
P can also be represented geometrically in terms of the three invariants of the stress tensor.Assuming that the three principal stresses have the order as r 1 !r 2 !r 3 : Hydrostatic stress or the first stress invariant is the distance from the origin to the deviatoric plane containing the point P ( OO 0 j j ). (1) Second deviatoric stress invariant characterizes the radius of the cylinder containing the point P ( O 0 P j j ). (3) Third deviatoric stress invariant characterizes the position of the point on the circle (or p-plane) and defines the Lode angle (h).
(5)  The stress state can also be determined by two dimensionless functions of the stress triaxiality (g) and Lode angle parameter ( h).   h The size of the examined yield surfaces is given by the equivalent stress (or J 2 ), the Lode parameter gives the shape in the deviatoric plane whereas the stress triaxiality characterizes the change in direction of the hydrostatic axis.von Mises criterion The von Mises plasticity criterion assumes that yielding begins when the second deviatoric stress invariant reaches a critical value that is a characteristic of each material.Since it is independent from the first and third stress invariant, yielding is independent from the hydrostatic stress and the flow stress is independent from the Lode angle.The von Mises stress r, is defined as a function of the second invariant which gives the radius of the yield surface (Fig. 2) (Wilson 2002).
The FEA results with the implementation of the von-Mises criterion reported and discussed in this work use the material constants reported in Tables 1-12.Table 1.Material properties used in the simulation for the bolt and nut (AISI 1045) (Bai 2007).

Drucker-Prager
The classical Drucker-Prager yield criterion, which was introduced to model soil mechanics, is a modification of the von Mises criterion that introduces the hydrostatic pressure (f(I 1 )) into the von Mises equation (Wilson 2002).The yield surface is conical with the radius that changes throughout the hydrostatic axis.The effect of the Lode angle is not taken in consideration and the mathematical representation is shown in Eq. ( 9) where J 2 is the second invariant of the Cauchy stress tensor, I 1 is the first invariant of the Cauchy stress tensor and A and B are material constants to be found experimentally.ffiffiffi ffi In this study, the Drucker-Prager criterion, implemented as an user-defined material model (VUMAT on Abaqus), was calibrated with the experimental results from Bai and Wierzbicki (2008) and Bai (2007).The constant values are reported in Table 3. C g g 0 C s h C c h m are four material parameters to be calibrated, e y is the equivalent plastic strain, E is the Young's modulus and is the Poisson's ratio.

Bai-Wierzbicki
Bai and Wierzbicki (2008) developed a model (Eq.10) of metal plasticity and fracture that consider both dimensionless hydrostatic pressure g and Lode angle parameter h.This model assumes homogeneous, isotropic, elastic-plastic with isotropic hardening material and plastic incompressibility (with deviatoric associated flow rule).The mathematical description of the yield criterion using the normalized Lode angle h is: where C ax h are defined as The first part of the equation represents the stress-strain curve at a reference loading condition of stress triaxiality g 8 , and it corresponds to the material strain hardening function.The second part represents the effect of the hydrostatic pressure on the yield surface and g o is the value of the stress triaxiality from the tests performed to calibrate the material strain hardening function.The stress triaxiality in the original model has been further modified (Ghazali et al. 2020) is the yield stress at the reference loading of g 8 : The third part of the yield criterion c À c mþ1 mþ1 is introduced into the deviatoric function to make the yield surface smooth and differentiable with respect to Lode angle h around c ¼ 1 characterizes the dependence on the Lode angle.The values of C t h , C s h and C c h are relative, and at least one of them is equal to unity.This depends on which type of reference test is used to calibrate the strain hardening function rð p Þ: The constants used to model the double lap bolted joint were calibrated using experimental results from Bai and Wierzbicki (2008) and from Bai (2007).They are summarized in Table 4.

Fracture locus
The damage evolution variable D can be calculated by means of a ductile fracture model (coupled model) or predict the failure without any link to ductile behavior (uncoupled model) (Baltic et al. 2020).There are several uncoupled models reported in literature including Rice and Tracey (1969) based on micro void growth and Cockcroft and Latham (Cockcroft 1968) on damage mechanism controlled by maximum (positive) principal stress.A two-dimensional fracture locus with the equivalent plastic strain to fracture and stress triaxiality both controlling the material ductility has also been proposed (see for instance Johnson and Cook (1985).Softening mechanisms based on voids growth have been formulated to relate the damage evolution to the material behavior (coupled model).The damage evolution parameter in coupled models is used to modify the effective stress tensor and model the decrease in strength caused by the damage.Xue and Wierzbicki (2008) developed a phenomenological model based on a power law damage rule to predict crack initiation and propagation based on a weakening factor.
Extensive studies and experimental evidence have been used to develop an uncoupled fracture criteria which uses a three-dimensional fracture locus which relates the strain at fracture with stress triaxiality and lode angle parameter.An incremental damage accumulation rule has also been developed to define the onset of fracture and also capture the nonlinear strain paths.Bai and Wierzbicki (2008) developed a linear relationship to relate the damage evolution factor (D), the fracture locus function f and the equivalent plastic strain e p D e p ð Þ ¼ The initial loss of stiffness and the onset of the damage occurs when D ¼ 1, and between 1 < D < D c the yield surface is modified to consider material softening according to the equation It is assumed that materials (or an element) total fracture occurs when the damage indicator reaches a critical value D c : This theory does not directly model the void growth, nucleation and linkage behaviors, but it assumes that, when a critical value of the equivalent plastic strain is reached, fracture occurs.e f is a general fracture locus function of the stress state described by the stress triaxiality g and the normalized Lode angle h: Six parameters D i describe the asymmetric fracture locus: The damage parameters calibrated by Bai and Wierzbicki (2008) and by Bai (2007) and used in this work for the fracture loci of AA2024-T351 and AISI 1045 steel are reported in Tables 5  and 6.

Validation of the modeling approach
The validation of the double lap bolted joint modeling approach adopted in this work was performed comparing the experimental tests from Abu-Sinna et al. ( 2018) against the results obtained with the FE model described in this work.The experimental study investigated the effect of the clamping force and geometrical parameters on the stiffness and strength of AISI 1006 double lap bolted joints.The experimental results were also used to tune the friction coefficients, bolt to flange clearances and the contact algorithm parameters.
The double lap bolted joint geometry is shown in Fig. 4, where the model is composed of three flanges with identical shape and dimensions.The force-displacement characteristic curves (Abu-Sinna et al. 2018) show slippage, which indicates that there is clearance even though the experimental test was done with a nominal value of the clearance equal to zero.A calibrated clearance of 0.1 mm between the shank of the bolt and the hole of the flanges was modeled to match the experimental results.
The bolt and nut were simplified by having a circular head instead of the hexagonal shape with the same contact area with the flanges.The threads were not modeled to reduce computational time and the complexity associated with the contact modeling strategy.
The mechanical properties used for the AISI 1006 and the M8 bolt are shown in Table 9, Fig. 6 and Table 10, respectively.The preload was applied as a thermal load defining the bolt shank with non-zero orthotropic thermal expansion properties in the axial direction (a y ¼ 10 À4 C À1 ) and null expansion coefficients for the transverse directions (a x ¼ a z ¼ 0 C À1 ).This approach even if is not representing the real condition of the joint if it was subjected to a thermal load is a numerical approximation to induce a given strain corresponding to the axial load and apply the required torque.
Taking advantage on the symmetry of the joint, only half of the geometry is modeled.The simulation is implemented in two steps with the preload of the bolt followed by the shear load applied at one end of the joint as a linear velocity (Fig. 6).
The preload torque on the bolt was set to 25Nm, which corresponds to an axial force on the shank of F preload ¼ 6976 N (Abu- Sinna et al. 2018).The reduction in length of the bolt due to the induced cooling results in a compressive stress in the flanges and tensile stress in the bolt via the internal contact defined between the parts.
The initial temperature value is estimated considering the stiffness equation (Eq.17) and is iteratively tuned via numerical simulations to achieve the expected preload.The assembly consists of three main parts (bolt, flanges and nut) and Discretized with circa 173000 nodes and circa 150500 3 D solid elements (C3D8R).This element formulation is suitable for complex nonlinear analysis that involves contact, plasticity, bending and large deformations.The mesh around the hole is refined to capture the near-hole stress gradients and improve result quality without severely increasing computational time.
The contact between all components is defined as General Contact and modified with interaction properties allowing for self-contact, finite-sliding formulation and edge-to-edge contact, which enforces contact during erosion that could not be detected.Individual property assignments are specified in the general contact interaction to simulate the friction between the surfaces in contact.For the flanges, a frictional tangential behavior and a hard-normal interaction behavior allowing subsequent separation are defined.Finally, an element set "erode" is created to define the contact domain of the elements that could be eroded during crack propagation.The General Contact algorithm can span several unattached bodies; thus, it is not limited to the contact of a single body with itself.For that reason, the element-set-erode consists of 5 parts: 3 flanges, nut and bolt are defined.The nodes of the nut and the lower shank of the bolt are coupled with constraint equations.
Figure 7(a) shows the contour plot of the principal stress under the thermal load applied on the shank and Fig. 7(b) shows the force-temperature characteristic curves.The stress distribution across the whole interface between the bolt and nut does not show any local stress increase that is one of the issues associated to the use of constraint equations.
An iterative numerical procedure using the force-displacement curve from the experimental results is used to find the friction coefficient between the flanges.The optimum value is l ¼ 0.31.The AISC (Manual 2005) suggests that the mean slip coefficient for steel is l ¼ 0.35, which agrees with the optimum value used in this study.
In Fig. 8 the force displacement characteristic curve for the double lap bolted joint from Abu-Sinna et al. ( 2018) is compared against the one predicted with the numerical model developed in this work.The static friction forces generated by the preload in the flanges allows the shear force to be transferred between the flanges up until the shear force limit (circa 4 kN) is reached.Once the frictional capacity is exceeded, the flanges begin to slip until the bolt is engaged and the stiffness suddenly decreases.After this point, the double lap bolted joint stops behaving as a slipresistant joint and acts as a bearing-type joint.Referring to Fig. 8, prior to Pt.1, the curve is almost linear as the load is transferred by bearing forces against the bearing faces of the hole where the flanges begin to yield as the load further increases.The next stages in the damage evolution of the joint are characterized by a gradual decrease in the joint stiffness (Pt. 1 À 2) that continues with bearing deformation (Pt. 2 À 4) and finally the middle flange suffers shear failure (Pt.5).
In order to assess the accuracy of the modeling strategy implemented in this work, the five stages of the damage evolution of the joint (Fig. 8) together with the force-displacement curve were compared against the experimental results.

Design of experiments
Using the modeling approach validated with experimental data from the literature, two parametrized studies are performed.To investigate the influence of geometrical parameters such as flange thickness, edge distance and width of the flanges as well as preload and flange grip, a Taguchi's L9 DOE has been implemented.It adopts orthogonal array and considers equal weight for each parameter.This approach allows investigation on the influence of each parameter with the minimum number of runs even if only the main effects can be observed without the possibility of studying the interactions.
The first set of parameters selected to change the stiffness of the joint and produce different effects on the dynamic characteristics of the joint.Contact pressure is modified by changing the preload of the bolt and the grip length.The stiffness of the joint is modified by changing the geometry of the flanges.These alterations in the model induced non-linearities in the bolted structure causing different failure modes.The hole diameter (U Flanges ) together with width (W ¼ 30 mm), flange length (L ¼ 60 mm) and bolt diameter (MD ¼ e ¼ 8 mm) are kept constant, whereas the amount of preload, grip length and flange thickness ratio are varied.Three different values are carefully selected for each factor.The schematic representation of the double lap bolted joint investigated is shown in Fig. 9.
Two different materials are considered in each test: AA2024-T351 for the flanges and AISI 1045 steel for the bolt and nut.Table 11 and Fig. 9 show a summary of the different DOE combinations analyzed and the parameters considered.D represents the bolt diameter in mm, T1 identifies the top and bottom flanges whilst T2 is the middle flange of the double lap bolted joint.
In order to perform a comprehensive understanding of the structural behavior of double lap bolted joints, additional cases are investigated.Four particular geometries selected from the DOE are further modified.The preload, edge distance and width of the flanges are varied, whilst the flange thickness is kept constant to 2T1 ¼ T2 and the grip length to 8 mm.The four combinations are repeated for Bai-Wierzbicki and von Mises plasticity models only as it will be observed in the DOE that they are those with the greatest difference in terms of force-displacement response and failure mode.

Results and discussion
The nine tests described in Table 11 have been repeated for each plasticity model.The representative force displacement curve with the most important steps observed are shown in Fig. 10.The first portion of the curve is linear (Stage 1) up to the shear force limit where the slip occurs (Stage 2) and the bolt engages with the internal faces of the hole.From this point the stiffness decreases up to the maximum force (Stage 3) where the flanges experience extensive plastic deformation, and the onset of the damage is observed (Stage 4).Due to the propagation of the damage the force decreases up to a point where the crack length reaches the critical value, and it becomes unstable (Stage 5) with the catastrophic failure of the joint.
There is a difference in the crack path followed by the damage across the different models depending on the plasticity models, the fracture locus and the geometrical parameters.
The analysis of the force displacement response confirmed that the preload has a major influence on the slip resistance (Fig. 11a) whilst the variation in grip length results in a non-linear change in the stiffness of the joint (Fig. 11b).
In Table 10(a) summary of the angle between the loading direction and the plane of the crack (a) and the maximum force is reported.Considering Bai-Wierzbicki as datum the maximum difference expressed in percentage is observed for Test 1 and Test 7.For Test 1 the maximum force estimated with von-Mises and Drucker-Prager is 10% and 15% lower respectively (Fig. 12a).For Test 7 the difference increases up to 22% for Drucker-Prager (Fig. 12b).For the other configurations, the difference observed is below 5%.
The American AISC (Manual 2005), the British Standard BS5950 (BSI 2000), and the EuroCode3 (de Normalizaci on 1996), have all identified minimum values for the distance between the edge of the plate and the axis of the bolt in order to control bearing failure.Fig. 13 shows the failure path usually assumed by these standards (red dotted line) since it is the one characterized by the lower strength.
The angle (a) at which actual failure occurs depends ON SEVERAL factors including material properties, preload and thickness of the main plate, amongst others.The nominal resistance of the two shear planes is related to the edge distance, the thickness of the plate and the ultimate shear strength of the main plate (0.6 times the UTS) by the equation: Using this equation, the nominal resistance is calculated for each configuration investigated in DOE and reported in Table 11 in the last two columns for the side flanges together and the middle flange in the last column.
Comparing the nominal resistance with the maximum load predicted with the numerical models, von Mises plasticity model is the closest ranging from 7% to 14% except for TEST 7 and TEST 8 with 74% and 45% respectively.The other two models Drucker Prager and Bai-Wierzbicki have a higher difference ranging from 3% to 73%.The reason for the higher difference with the latter two models is that the equation assumes a shear out failure, characterized by a lower load.The models with von Mises plasticity predict a failure close to the equation, whilst the other two plasticity models are predicting a crack path with a greater angle that results in greater load.
When the difference in the prediction of the angle a is considered, von-Mises and Drucker-Prager differ as high as 90%.The type of failure predicted with von Mises in all cases but Test 7 is shear-out whilst for the other two plasticity models the angle ranges between 9 and 55 for Drucker Prager and 5 and 74 for Bai-Wierzbicki.
In Fig. 14, the comparison of the crack path observed with the same geometry (Test 1) and different plasticity models.The contour plot represents the stress contour map.Shear out failure is predicted with von Mises whilst Bai-Wierzbicki is predicting a greater angle that is very close to net-section failure.Drucker Prager model is predicting a shear out failure with an angle greater than von-Mises.The location of the onset of the failure is the same for the three models.
The response surfaces in Fig. 15 show the effect of the thickness of the flange and the preload on the failure mode represented by the angle a. Considering the results obtained with Bai-Wierzbicki (Fig. 15a), increasing the joint thickness and the preload the angle a first decreases and after increases.This can be explained with the combined effect of the stress triaxiality induced by the increase in thickness and the increase in the compressive stress produced by the preload.For small value of the preload, increasing the thickness the transition from plane stress to plane strain promotes the brittle failure of the flange.For high value of the preload, the compressive stress effect in the region surrounding the hole is dominant and produces an increase in the stress triaxiality producing a ductile failure regardless of the thickness of the flanges.This effect is confirmed by the observation that the failure is pure shear (a ¼ 0 ) when the Lode angle and the hydrostatic pressure is not accounted in the identification of the strain at failure (von Mises) in all geometrical configuration investigated (Fig. 15c).The characteristics of the failure mode observed with Drucker Prager (Fig. 15b) are in between the other two plasticity models.For low preload the failure is shear-out regardless the value of the flange thickness like von-Mises, whilst increasing the preload the failure moves toward net-tension like observed with Bai-Wierzbicki.
In Fig. 16, the force-displacement curve obtained from the three additional cases are reported.The results have been obtained using von Mises and Bai-Wierzbicki plasticity models as previously explained.
The results from this set of simulations are summarized in Table 12 where maximum force is reported for each configuration It was observed that if the width is much larger than the distance of the hole from the edge the failure is tear-out in both plasticity models (Fig. 17 shows a contour plot of the plastic strain).The angle followed by the crack path is close to 45 in case the plasticity model adopted is Bai and Wierzbicki with a greater value of the maximum force.This result is consistent with the experimental evidence reported in literature (Bouchaïr, Averseng, and Abidelah 2008).
The results from Test 3 refers to the condition where the width is comparable to the distance from the edge and in this case bearing failure is observed consistently with the experimental evidence provided in (Bouchaïr, Averseng, and Abidelah 2008).The main difference produced by the plasticity models is that the middle flange would fail when Bai-Wierzbicki is used, whereas top and bottom flanges would fail when von Mises plasticity model is adopted (Fig. 18 shows the plastic strain contour).
When failing in net section (Fig. 18a) force--displacement curve obtained with both plasticity models are very close due to the negligible effects of the hydrostatic pressure and Lode Angle.When failing in shear (Test 1) von Mises plasticity model would predict a tear-out failure whereas Bai-Wierzbicki would predict a higher angle a.

Conclusions
Double lap bolted joints are usually subjected to complex stress states and stress concentrations.For this reason, material properties and contact must be studied and validated against experimental data to achieve reliable and robust modelling strategies.
It has been observed that Bai-Wierzbicki plasticity and ductile fracture model is capable of numerically predict different types of failure modes for double-lap bolted joints a higher a compared to the other two models here considered.Drucker-Prager and von Mises models produce a failure mode closer to an idealised pure shear plane whilst von Mises and Bai Wierzbicki plasticity models have close force-displacement response when net section failure is observed.When a higher preload is introduced a beneficial compressive stress field around the edge of the hole and a reduced shear stresses due to the increasing triaxiality of the stress state is observed resulting in an improved force displacement response.The effect of the width and the distance between the hole and the edge induces different triaxiality in the flanges for different plasticity models resulting in failure in different locations and path.

Figure 7 .
Figure 7. (a) Pressure cone and (b) contact force vs temperature.

Figure 8 .
Figure 8.Comparison of the experimental results vs the ABAQUS simulation.Experimental pictures from (Abu-Sinna et al. 2018).

Figure 10 .
Figure 10.Phases of a double lap bolted joint under shear loading.

Figure 17 .
Figure 17.Failure mode observed with (a) Bai-Wierzbicki and (b) von Mises corresponding to width larger than edge distance.

Figure 18 .
Figure 18.Failure predicted with the use of (a) Bai-Wierzbicki, and (b) von Mises when the width is comparable to the edge distance.

Table 6 .
Damage parameters used for 1045 steel.

Table 10 .
Results from DOE.

Table 11 .
Analytical bearing strength of the different DOE tests performed.

Table 12 .
Comparison of maximum force values for three different tests.