Finite Elements in Analysis and Design 110 (2016) 32–42

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design

journal homepage: www.elsevier.com/locate/finel

Full Length Article

A reduced fastener model using Multi-Connected Rigid Surfaces

for the prediction of both local stress field and load distribution
between fasteners
Ramzi Askri a,n, Christophe Bois a, Hervé Wargnier a, Julie Lecomte a,b
a
Univ. Bordeaux, I2M, UMR 5295, F-33400 Talence, France
b
ASTF, 8 Avenue du Val d'Or, 33700 Mérignac, France

art ic l e i nf o a b s t r a c t

Article history: This paper describes the development of a reduced model of a fastener using Multi-Connected Rigid
Received 27 March 2015 Surfaces (MCRS). The stiffness of the connectors is determined, based on a physical approach, considering
Received in revised form different deformation modes of the bolt. The reduced model is constructed and identified from a
2 October 2015
numerical simulation of a single lap reference joint under tensile load, with the adherent parts and bolts
Accepted 22 November 2015
represented by 3-D solid elements. A single simulation with a given clearance, axial preload and friction
Available online 11 December 2015
coefficient is used to identify equivalent stiffnesses. The reduced model is then compared with the 3-D
Keywords: solid elements model in a two-fastener configuration for different values of clearance, preload and
Bolted joint friction coefficient. The comparison covers overall response in terms of stiffness and load distribution
Composite material
between fasteners, local response in terms of stress fields and calculation times. Results show that the
Finite elements
reduced model proposed here is able to reduce calculation times while still providing a good estimate of
Rigid surface
Connector the mechanical quantities needed for the study and dimensioning of multi-fastener joints.
& 2015 Elsevier B.V. All rights reserved.

1. Introduction preload due to its low resistance in the out-of-plane direction. In

addition, the materials in contact tend to be damaged as a result of
Because of the costs and delays incurred in experimental stu- bearing stresses, which also leads to a change in stiffness. Gen-
dies, numerical simulation is an essential tool when designing erally, with a single lap joint, the rotation of the parts caused by
assembled structures. The design of fastened joints is based on the secondary bending moment is directed by the bending stiff-
making a good estimate of the distribution of loads between fas- ness of the fastener heads. This is an important phenomenon since
teners. This load distribution depends on the stiffness of the fas- it gives rise to a non-uniform distribution of contact pressure
tener, the stiffness of the parts, the way in which external forces between the pin and the two holes.
are introduced, and also clearances and geometrical defects. Note The quantity and complexity of the phenomena involved lead
that by stiffness of the fastener, we refer to the ratio of the load us naturally to model the fastened joints using a refined mesh of
borne by the fastener to the relative displacement, one against the 3-D solid elements. This type of model was developed in order to
other, of the parts adjacent to the fastener. Depending on whether study the effect of the friction coefficient [1], clearance [2], loca-
the load is transmitted by adherence or by fastener pin-hole tion error [3] or preload [4]. As well as predicting load distribution
contact, this stiffness will be very different. In the former case
between the fasteners, these models also provide accurate results
(load transmitted by adherence), the role of the axial preload is
for the stress fields in the high gradient regions around the hole.
fundamental; in the latter case (load transmitted by contact), the
They are also known for their ability to take material non-
presence of clearance results in contact occur gradually and hence
linearities into account, such as damage [2,5,6]. Several modes of
stiffness evolves. For a composite structure joint, the two modes of
degradation can be introduced, such as delamination, matrix
load transfer combine since the material is unable to bear a high
cracking or fibre failure. Despite these performances, models using
3-D solid elements are still very costly in terms of calculation time
n
Corresponding author. Tel.: þ 33 5 56 84 79 85; fax: þ 33 5 56 84 58 43. and they require specifically adapted calculation strategies, nota-
E-mail addresses: [email protected] (R. Askri),
bly parallel computing by domain decomposition to process
[email protected] (C. Bois),
[email protected] (H. Wargnier), assemblies with several dozen fasteners or to carry out parametric
[email protected] (J. Lecomte). studies on several hundred configurations [7,8].

http://dx.doi.org/10.1016/j.finel.2015.11.004
0168-874X/& 2015 Elsevier B.V. All rights reserved.
 R. Askri et al. / Finite Elements in Analysis and Design 110 (2016) 32–42 33

Thus in order to deal with assemblies with a very large number

of fasteners, or to carry out extensive parametric studies, analytical
or semi-analytical models have been developed. The 1-D models
proposed in the literature can predict load distribution between
fasteners in a wide variety of configurations: bolted [3,9–12] or
bolted/bonded [13], with clearance and location errors [3], taking
into account loss of stiffness due to bearing damage [3,9,13,14].
However, these models are limited to one-dimensional plane
geometries and loadings.
Several approaches have been developed in recent years to
simplify the way fasteners are represented in finite element
models. The first consists of representing each fastener by a special
connection element, available in most finite element analysis
software packages, which links a single node from each adherent
Fig. 1. Geometry of the joint.
part, potentially incorporating various non-linear behaviours [15–
17]. In this case, the holes and the contacts are not represented.
This is a very interesting approach in terms of calculation time, but matrix T700GC/M21. The stacking sequence for the composite
it requires a specific identification procedure to include the role of adherent is [90/0/  45/ þ45/0]s and ply-thickness is 0.25 mm. The
clearance, friction and material non-linearities overall in the second adherent part assembled is made of aluminium alloy 2024.
behaviour of the connector. Moreover, the stress field around the The two steel bolts, diameter 6.35 mm, are fixed with a radial
fasteners is not represented, and thus macroscopic failure criteria clearance of 10 mm and axial preload of 3500 N generating an
must be used for the assembled parts [18,19]. average contact pressure under the head of about 75 MPa. All
Gray et al. [20] have modelled the adherent part with shell ele- material properties are listed in Table 1.
ments and the bolts with two rigid cylindrical surfaces connected by a
beam exhibiting elastic behaviour. The heads were replaced by cou- 2.2. Analysis method of deformation modes of bolt
pling the movement of the beam ends with the adherent parts. The
model is validated for a single lap joint with one bolt and three bolts, The model created with Abaqus code uses 3-D solid elements
including clearance and preload. The problem with these approaches with reduced integration (C3D8R). The threaded junction between
remains the identification of the equivalent stiffnesses associated with the nut and the screw is not modelled. The bolt is therefore
the connection elements. A simple estimate of shear stiffness can be represented as a single part. The flexibility resulting from the
made using various analytical models to be found in the literature threaded junction between the nut and bolt is therefore neglected
[21,22], but the values obtained can vary by a factor of 1 to 5. in the reference model and hence in identifying stiffnesses in the
In this paper, a reduced model of a fastener using Multi- MCRS model. According to the literature, this flexibility can be
Connected Rigid Surfaces (MCRS) is proposed. The number of estimated analytically and thus incorporated afterwards into the
rigid surfaces and the stiffnesses introduced between them were MCRS model [23,24]. Preload is modelled by initial penetration of
selected by analysing deformation modes of the fastener in a 15 mm of each head into the adherent parts giving an intended
single-lap configuration obtained from a reference model made up axial preload of 3500 N. A contact algorithm using the penalty
of 3-D solid elements. From this decomposition, 4 rigid surfaces method with a friction coefficient of 0.1 is also introduced to
and 4 stiffnesses were defined. The procedure for identifying these model contact between the bolts and the adherent parts. Defor-
equivalent stiffnesses is also based on the 3-D reference model in a mation modes in the bolt are analysed with a tensile load of about
fixed configuration (clearance, friction coefficient, preload). The 15 kN, which corresponds to the initiation of bearing damage
efficiency of the model was assessed on the basis of global criteria obtained experimentally.
(at the scale of the joint) and local criteria (at the scale of the A section S of a bolt is defined as a set of nodes having initially the
materials). The validity domain of the identification was studied by same coordinate Z. To represent the kinematics of the bolt, for each
varying the clearance, friction coefficient and preload. section S, initially written S0 with centre M0, we define a least-squares
plane P from actual position of the nodes of section S. In the event of
!
uniaxial loading along X , plane P remains perpendicular to plane (XZ).
2. Analysis of deformation modes of the fastener !!!
An orthonormal coordinate system ðM; t ; y ; n Þ can then be
assigned to each deformed section, as shown in Fig. 2. M is defined as
In this section, an analysis of deformation modes of the fastener is the intersection of plane P and the curve representing the deflected
proposed and in particular the relative movement between the func- shape of the bolt axis. For the centre section O (O being the inter-
tional surfaces of the fastener. Studies referred to in the bibliography section of the axis of the bolt in its deformed state and the overlap
[1,3,21,23] show that deformations in the bolt can be decomposed into plane of the joint) the associated coordinate system is written
4 main modes: tension and bending along the fastener pin axis, !!!
ðO; x ; y ; z Þ.
shearing of the pin and bearing–compression of the contact semi- In order to obtain a displacement field that can be exploited to
cylinders. A finite element model where the fastener is represented by analyse the different deformation modes in the bolt, a rigid body field
3-D solid elements is proposed to extract quantitatively these different displacement was subtracted from the calculated displacement field
deformation modes. This analysis is used as a basis for building and !!! !!!
so that coordinate system ðO; x ; y ; z Þ and ðO0 ; X ; Y ; Z Þ coincide. A
identifying the reduced model.
column matrix of small displacements [D] is then defined to represent
!
the displacements of translations ux and uz, and also rotation θy ¼ ð x
2.1. Description of the reference joint !
; t Þ for each section, as shown in Fig. 2:
The reference configuration used to study deformation modes 2 3
ux
in the fastener is a single lap joint made of aluminium/composite 6u 7
½D ¼ 4 z 5: ð1Þ
material with two bolts, as shown in Fig. 1. The composite material
is a carbon fibre laminate with unidirectional ply and a thermoset θy
34 R. Askri et al. / Finite Elements in Analysis and Design 110 (2016) 32–42

Table 1
Material properties used in reference model.

Properties of unidirectional ply for T700GC/M21

E11 [GPa] E22 [GPa] E33 [GPa] G12 [GPa] G13 [GPa] G23 [GPa] ν12 ν13 ν23

110 7.6 7.6 4.75 4.75 2.65 0.33 0.33 0.43

Properties of aluminium alloy 2024

E [GPa] ν

70 0.3

Properties of steel

E [GPa] ν

210 0.3

Initial state Deformed state

z n
S
S0

M0
θy
M x

t
Z + uz
Z P
Z
z

ux
X
O0 O
x

Fig. 2. Kinematics of bolt sections.

Fig. 3. (a) Coordinate system in the deformed state and (b) definition of sections of the bolt head–pin junction.

Specific points on the bolt are defined in plane (XZ) as shown in Since the rigid-body field displacement has been subtracted, dis-
Fig. 3a. To take the deformation of the heads into account and to placements and rotations in the bolt are weak. Thus components ux, uz
represent properly the movement of the functional surfaces that make and θy are representative of the deformation modes in the bolt. Dis-
up the surfaces under the head, the sections centred on A and B are placement uz represents strains along the bolt axis, the rotation of
divided into “Pin” section (exponent P) and “Head” section (exponent sections θy is the result of bending strains and du
dz
x
 θy represents the
H) as shown in Fig. 3b. shear slip of the sections.
 R. Askri et al. / Finite Elements in Analysis and Design 110 (2016) 32–42 35

2.3. Tensile and bending deformation modes

The graph in Fig. 4a showing change in θy as a function of the

position on z-axis shows a slight rotation of all sections from the
initial state (θymax  0:121, or 0.0021 rad). However, there is a
steep slope close to the ends of the pin and a significant
deviation between sections SH P
A;B and SA;B is observed. The graph
in Fig. 4b showing change in uz as a function of the position on
z-axis also shows a wide deviation in axial displacements
between sections SH P
A;B and SA;B . These results are due mainly to
the deformation of the bolt heads from pressure applied by the
adherent parts.

Fig. 5. Deflected shapes of segments (AB), (CD) and (EF) extracted from the bolt pin.

2.4. Shear and bearing–compression deformation modes

Results show that angle θy remains small along the pin, com-
pared with dux =dz. Displacement ux and the deflected shape of the
axis (AB) are therefore mainly the result of shear strain. Fig. 5
shows the deflected shape (AB) of the axis of the pin, also the
deflected shape (CD) and (EF) of its two contact lines in the loading
plane (XZ) which accumulate both shear and compression strains
in the pin in areas of contact with the adherent parts.
The proportion of displacements caused by shear loading and
bearing–compression depends for the most part on the geometry of
the bolt. The shorter the bolt in relation to its diameter, the greater the
compression strains in relation to shear strains. Note also that during
transverse loading, the adherents are subjected to shear deformation
due to tangential frictional forces in the overlap surface. If adherent
material stiffnesses are lower than the fastener stiffness and bolt-hole
clearance is low, contact can be established between the pin and
adherents near the overlap plane before sliding. That point and the
effects of bending in the adherent parts tend to concentrate the con-
tact pressure near the overlap plane as shown in Fig. 5. Moreover, on
comparing the deflected shape of (CD) and (EF) it can be seen that the
orientation of the composite plies affects the regularity of the dis-
placement field because of the variation in stiffness of the material
adjacent to the contact.
Gray and McCarthy [20] modelled fasteners with two rigid surfaces
connected by an elastic beam in order to account for shear and
bending deformation. This model therefore overlooks the deformation
induced by bearing–compression.

3. Construction of the MCRS model

3.1. Choice of kinematics for the MCRS model

The purpose behind reducing the model is to do away with the

3-D elements of the bolt, which represent a large proportion of the
number of degrees of freedom in the joint. The idea is to use a
limited number of rigid surfaces to represent the functional sur-
faces of the bolt and connect them elastically. For the deformation
mode generating a rotation located at the pin–head junctions, it is
natural to propose a division along the planes passing through A
!
and B accompanied by transverse rotational stiffness along x and
Fig. 4. Displacement of sections of the bolt: (a) rotation around y and !
y . The two rigid surfaces modelling the bolt heads are written RSH
1
(b) translation along z.
and RSH 2 respectively. For the other deformation modes (axial on
36 R. Askri et al. / Finite Elements in Analysis and Design 110 (2016) 32–42

Fig. 6. MCRS model of bolt: (a) kinematic behaviour and (b) stiffnesses associated with connectors.

Fig. 7. (a) Equivalent axial forces and (b) equivalent moments.

! !
z , shear between the normal planes z and compression at the according to fastener axis are not explicitly represented, but only
contact points) as the strains are not localised, the division is not a the global relative displacement of two contact surfaces. The
! geometry and the relative displacements permitted in the differ-
natural one. Concerning the axial deformation on z , the aim of
! ent parts of the MCRS are shown in Fig. 6a.
the model is to predict relative displacement along z of the sur-
faces under the heads as a function of the normal load borne by Stiffnesses and constraints are introduced into Abaqus by 2-
!
the bolt along z . A single division associated with stiffness in node connector element with 3-D behaviour (CONN3D2). They
translation is therefore proposed. The decision to place the divi- enable the master nodes of each rigid surface to be linked and
sion in the centre section O was an arbitrary one. Concerning
! introduce stiffnesses according to certain degrees of freedom. As
shearing between the normal planes z and compression at the
contact points, the aim of the model is to predict the relative the bolt is axisymmetric, the same stiffnesses are used in direc-
! !
transversal displacement of the section portions in contact as a tions x and y . Abaqus proposes a variety of usual links in its
! !
function of the transverse forces along x and y . Normally, a connector element library [25]. The connector C0 is an assembly of
!
model of this kind should be discretised along z . However, as the
two connectors ALIGN and CARTESIAN. The ALIGN connector
ultimate aim is to use the reduced model of the bolt with an
adherent part model using continuum shell elements, it seems ensures that the directions of the local coordinate system for the
judicious to divide the pin into 2 cylindrical rigid surfaces RSP1 and two rigid surfaces remain parallel. The CARTESIAN connector
RSP2 in order to maintain some coherence between the kinematics defines 3 translational stiffnesses (K n , K tx , K t y ). Connector C1 (and
of the bolt model and that of the adherents. It is worth noting that C2 respectively) of the UJOINT type allows for two possible rota-
since in the MCRS model, shear deformations, bearing–compres-
tions with rotational stiffnesses K 1x and K 1y (K 2x and K 2y respec-
sion and bending deformations are gathered and concentrated in a
single section (pin section which is in the overlap plane), the tively) and thus blocks the other degrees of freedom. These dif-
sliding of the pin section due to shear and associated rotations ferent stiffnesses are presented in Fig. 6b.
 R. Askri et al. / Finite Elements in Analysis and Design 110 (2016) 32–42 37

Fig. 8. (a) Evolution of the percentage of strain energy for each stiffness introduced in the MCRS model with the reference configuration and (b) change in bolt strain energy
in the MCRS model at 15 kN as a function of transverse stiffness K tx .

Table 2
Connectors and associated stiffnesses.

Connector C0 C1 C2
Type ALIGN þCARTESIAN UJOINT UJOINT

Stiffness K t x [kN/mm] K ty [kN/mm] Kn [kN/mm] K 1x [kN mm/rad] K 1y [kN mm/rad] K 2x [kN mm/rad] K 2y [kN mm/rad]

368 368 243 4239 4239 4669 4669

!
3.2. Identification of translational stiffness along z and transverse same process is applied to determine Θ2y which represents the
! !
rotational stiffness along x and y equivalent angle of rotation between the rigid surface of the head
RSH P
2 and the surface of the pin RS2 . We therefore obtain
The purpose of the connectors introduced into the MCRS model   
M  M 
 1
is to transmit forces and moments and to position the rigid sur- K 1y ¼  y  ¼  y 
1
ð4Þ
faces in the joint. The division that was decided on and described  Θ1y  θy1
in Section 3.1 results in the kinematic behaviour of sections SH
A and  
SH
B being used to determine stiffness K n and two rotational stiff-
M  M 
 2y   2y 
K 2y ¼ ¼ ð5Þ
nesses K 1y and K 2y . The displacement vectors associated with  Θ2y   θy2 
these two surfaces are defined as follows:
2 3 2 3 In principle, the above definitions of stiffnesses Kn, K 1y and K 2y
ux1 ux2
6 uz 7 6 uz 7 can be computed for each level of applied load. Numerical results
½D1  ¼ 4 1 5 and ½D2  ¼ 4 2 5 ð2Þ
show that during the sliding phase and clearance recovery, contact
θy1 θ y2
pressure under the bolt head is distributed uniformly. Once pin-
For stiffness K n , as displacement uz cumulates for the most part hole contact is established, contact pressure under the head on the
in the bolt heads (see Fig. 4b), the representative force F n , which opposite side of pin-hole contact tends to increase because of the
connector C0 must transmit, is defined as the mean of the axial rotation of the heads, thus generating a greater strain on the
forces, written F z1 and F z2 transmitted by sections SPA and SPB of the heads. As a result, the equivalent tension stiffness Kn decreases
3-D model, as shown in Fig. 7a. Equivalent axial displacement U n is with the load applied (Kn ¼368 kN/mm for a load of 1 kN and
the difference in the axial translational displacements of the Kn ¼ 243 kN/mm for a load of 15 kN). K 1y and K 2y are less impacted.
heads. We therefore obtain As when analysing the deformation modes of the fastener, the load
  selected to identify stiffnesses is equal to 15 kN, which corre-
F n 1F z1 þ F z2  sponds to the initiation of bearing damage obtained expe-
Kn ¼ ¼  ð3Þ
U n 2 uz2  uz1  rimentally.
For the transverse rotational stiffnesses, the equivalent angle of
rotation Θ1y between the rigid surface of the head RSH ! !
1 and that of 3.3. Identification of translational stiffness along x and y
the pin RSP1 is equal to θy1 , the rotational component of the dis-
placement vector ½D1 . The equivalent moment M 1y is the external The method we propose here to determine shear/bearing–com-
moment applied by the adherents to the heads and obtained pression stiffness is based on an energy approach, given the difficulty
directly from the 3-D model on section SH A as shown in Fig. 7b. The in estimating displacements and the equivalent forces applied to the
38 R. Askri et al. / Finite Elements in Analysis and Design 110 (2016) 32–42

functional surfaces modelling the pin in the same way as the trans- terms of bolt rotation and deviation of the adherent parts at the
! ends of the overlap region.
verse rotational and translational along z stiffnesses.
As only stiffness K tx ¼ K t y has not yet been identified, we search Fig. 10 shows the force–displacement graphs for the 3-D and
for the value of this stiffness that will equal the strain energy of MCRS models in the different configurations tested. Results show a
the bolt calculated with the 3-D model and the sum of the strain good level of prediction for the joint stiffness after the sliding
energies for all the connectors in the MCRS model. As shown in phase and clearance recovery, whatever the value of the clear-
Fig. 8a, the strain energy associated with transverse stiffness K t x ances. Sliding forces are underestimated by about 30% with the
represents about 85% of the total strain energy of the bolt at a load MCRS model, for all configurations. Indeed, as mentioned in
level of 15 kN (joint displacement equal to 0.7 mm). Deducing the Section 3.2, stiffness Kn depends on the load selected for the
value of K t x from the intersection of the graph shown in Fig. 8b identification of the MCRS model. The deviation on sliding forces
representing the change in strain energy estimated with the MCRS matches the relative gap between the equivalent tension stiffness
model and the level of strain energy calculated with the 3-D model Kn obtained with the 3-D finite element model for a load close to
at a load of 15 kN is therefore a reliable approach. Stiffnesses the sliding force (Kn ¼368 kN/mm for a load of 1 kN), and the load
identified for each connector are listed in Table 2. taken to identify the MCRS model (Kn ¼ 243 kN/mm for a load of
15 kN). Obtaining a better prediction of both sliding forces and
final joint stiffness will require defining stiffness Kn which
depends on the load applied to the considered fastener. Using
4. Results and discussion
stiffness identified with a fairly high load (15 kN) explains why the
sliding forces are under-estimated during the clearance compen-
The purpose of this validation is to show that the MCRS model can
sation phase. We should remember, nevertheless, that for the
be used both for predicting load distribution in a multi-fastener joint
composite parts, load transfer by friction remains low since the
and also for predicting the stress state in assembled parts and thus
material cannot support a large preload as it has low resistance in
apply adequate resistance criteria.
the out-of-plane direction.
Validating the MCRS model is based on comparing global and
A graph of changes in loads transmitted by each bolt in con-
local criteria with the 3-D model. To assess the domain of validity figuration 2, where clearance differs for the two bolts, is shown in
of the MCRS model, the comparison covers different configura- Fig. 11. The good agreement between the two models shows that
tions compliant with Fig. 1 where clearance, friction coefficient the MCRS model estimates overall transverse stiffness of the bolt
and preload have been modified. The different configurations correctly and that it can be used in cases of multi-fastener joints.
tested are listed in Table 3. It should be remembered that only case
no. 1 was used for the identification. 4.2. Local behaviour

4.1. Global behaviour The strength of the assembled parts is based on the distribution
and amplitude of the stresses around the hole. Only the results for the
As shown in Fig. 9, at the scale of the overall structure, com- composite adherent part are shown here. For composite materials, the
paring the deflected shape of the joint gives good agreement in model should correctly predict stress in the direction of the fibres

Table 3
Configurations of joints tested.

Case Penetration gen- Friction Clearance in Clearance in

erating preload coefficient bolt 1 [mm] bolt 2 [mm]
[mm]

1 15 0.1 10 10
2 15 0.1 10 125
3 15 0.1 80 80
4 15 0.3 80 80
5 5 0.1 80 80

Bolt 1 Bolt 2
3-D Model

MCRS Model

Fig. 9. Comparison overall of deflected shape of 3-D and MCRS models (configuration 1) at a load of 15 kN: (a) overall view and (b) close-up on bolt 2.
 R. Askri et al. / Finite Elements in Analysis and Design 110 (2016) 32–42 39

Fig. 10. Force–displacement graphs for the different configurations tested.

40 R. Askri et al. / Finite Elements in Analysis and Design 110 (2016) 32–42

local response of the MCRS model is thus promising for dealing with
bearing failure.
The maximum tensile stress in the direction of the fibres
(positive values of σ11) is very well estimated by the MCRS model,
and thus this model can be used to predict net-section failure.

4.3. Calculation time

First, it is worth noting that the objective with the proposed

reduced model is to study joints with a large number of fasteners
and to perform parametric studies by running tens or hundreds of
calculations. For a given case defined by adherent materials and
thicknesses, only one calculation with a two-bolt 3-D model has to
be performed in order to identify equivalent connector stiffnesses.
The calculation time using the reference 3-D model will be neg-
ligible compared to total calculation time.
Although the aim of this study was ultimately to use the MCRS
bolt model with adherents modelled by shell elements, it is
nevertheless interesting to evaluate the calculation time saved by
simply replacing the 3-D elements of the bolt by the MCRS model.
Numerical simulations of the different models were carried out
with an Abaqus/Standard implicit scheme using 4 cores of an Intel
Xeon 16-cores 2 GHz processor and 64 GB RAM computer. With the
MCRS model, the saving in total calculation time compared with the
Fig. 11. Load distribution between the two bolts when clearances differ
(configuration 2).
3-D bolt model varies from 10% to 68% depending on configuration, as
shown in Table 4. Reducing the number of degrees of freedom by
under tension (failure initiation in net-section) and under compres- using rigid analytical surfaces means that not only can the size of the
sion (bearing damage initiation). Note that many studies [26–28] have stiffness matrix of the FE model be reduced but also the number of
shown that the failure of fibres under compression also depends on degrees of freedom involved in the contact algorithm. However, this
shear stresses and transverse stresses. Therefore, in order to assess the last point may have a negative effect since by making the bolt rigid,
efficiency of the MCRS model at the scale of the material, it is vital to resolving the contact law may become more difficult. As a result, the
number of iterations per increment of load may increase. Thus
compare the field of the 6 components of the stress tensor at the edge
although the average time saving per iteration is about 40%, the time
of the hole in the ply coordinate system: direction 1 is the direction
saving overall is markedly less.
of the fibres, direction 2 is transverse direction in the stratification
plane and direction 3 is transverse direction perpendicular to the stra-
tification plane.
The comparison is made between the 3-D and MCRS models 5. Conclusion
using configuration 1 with a load of 15 kN. Results for the other
This article describes a model of a fastener based on rigid surfaces
configurations are similar. The images in Fig. 12a–g show the edge
connected elastically. The construction of the model is based on a
of the hole in the composite adherent part in contact with bolt 2.
reference finite element model made up entirely of solid 3-D ele-
Results for bolt 1 are virtually identical. Note that the resultants for
ments. The method can be used to define the number, the position
the forces transmitted by the pin of bolt 2 are very similar for the
and the type of stiffness (translational or rotational) to introduce
two models: 6864 N for the 3-D model and 6934 N for the MCRS
between these rigid surfaces based on analysis of the different
(see results in Section 4.1). To facilitate analysis of the results, the
deformation modes of fastener. Using this approach, it will be possible
contact pressure field is also drawn (Fig. 12a).
to make choices linked with physical phenomena and to understand
The pressure fields are fairly similar in appearance, however, the
the limitations associated with these choices. As identifying the stiff-
contact pressure is located closer to the overlap plane (Fig. 12a) for the
nesses associated with each deformation mode is partly dissociated
MCRS model. This accounts for the deviations of about 25% observed
and based only on stress and strain fields in bolt, it therefore results in
for the maximum values of stresses linked directly with the pressure a unique set of parameters, which are not much dependent on the
field (negative values of σ11 and σ22, positive and negative values of behaviour of the assembled parts, contacts between parts and fastener
σ12). The same deviations were noticed for the aluminium plate. The and level of preload. This method therefore provides very satisfactory
main explanation is that the flexibility of the pin in the 3-D model is results not only for global stiffness and the distribution of loads
more distributed than in the MCRS model, which generates a more between fasteners, but also for the local response assessed from
spread contact pressure field whether the material used is layered contact pressure fields and stress fields.
composite or isotropic. However, whatever the FE model used, it is With this reduced model of a fastener combined with a con-
complicated to estimate correctly the stress field at the edge of the tinuum shell elements model for the assembled parts, the number
hole because of the high stress gradient in this area which generates of degrees of freedom can be reduced by about 85%, a saving in
high mesh sensitivity. Practically, in order to prevent mesh sensitivity calculation time of about 80% and it is possible to envisage accu-
and also introduce physical material characteristic length, failure cri- rate simulations of complex structures with hundreds of fasteners.
teria, especially for composite materials, are based on point stress or In addition, the fact that the true shape of the parts is explicitly
average stress approaches, with the result that we do not only use the represented with the holes means that failure criteria or non-linear
stress field at the edge of the hole [29,30]. Results show that the gap laws resulting from plasticity or damage at the material scale can be
between the MCRS model and the 3-D FE model tends to decrease applied directly, thus avoiding a stage of local re-analysis or the use of
quickly as soon as one moves away from hole edge. Even so, the use of macroscopic criteria in a narrow validity domain.
 R. Askri et al. / Finite Elements in Analysis and Design 110 (2016) 32–42 41

Fig. 12. Comparison of pressure fields and stress fields at the edge of the hole in the composite adherent part in contact with bolt 2 at 15 kN: (a) contact pressure; (b) σ11;
(c) σ22; (d) σ33; (e) σ12; (f) σ13; and (g) σ23.
42 R. Askri et al. / Finite Elements in Analysis and Design 110 (2016) 32–42

Table 4 [8] V. Roulet, L. Champaney, P.-A. Boucard, A parallel strategy for the multi-
Comparison of calculation times for the different configurations tested. parametric analysis of structures with large contact and friction surfaces, Adv.
Eng. Softw. 42 (2011) 347–358.
Configuration Model Calculation time Total number Saving in calcula- [9] C.T. McCarthy, P.J. Gray, An analytical model for the prediction of load dis-
CPU [s] of iterations tion time CPU [%] tribution in highly torqued multi-bolt composite joints, Compos. Struct. 93
(2011) 287–298.
[10] J. Andriamampianina, F. Alkatan, P. Stéphan, J. Guillot, Determining load dis-
Case 1 3-D 2502 99
tribution between the different rows of fasteners of a hybrid load transfer
MCRS 2041 133 18
bolted joint assembly, Aerosp. Sci. Technol. 23 (2012) 312–320.
[11] J. Chakhari, A. Daidié, Z. Chaib, J. Guillot, Numerical model for two-bolted
joints subjected to compressive loading, Finite Elem. Anal. Des. 44 (2008)
Case 2 3-D 3121 135 162–173.
MCRS 2786 254 10 [12] J. Ekh, J. Schön, Finite element modeling and optimization of load transfer in
multi-fastener joints using structural elements, Compos. Struct. 82 (2008)
245–256.
Case 3 3-D 5041 211 [13] C. Bois, H. Wargnier, J.-C. Wahl, E. Le Goff, An analytical model for the strength
MCRS 3457 235 31 prediction of hybrid (bolted/bonded) composite joints, Compos. Struct. 97
(2013) 252–260.
[14] J. Ekh, J. Schön, D. Zenkert, Simple and efficient prediction of bearing failure in
Case 4 3-D 4099 108 single shear, composite lap joints, Compos. Struct. 105 (2013) 35–44.
MCRS 3054 120 25 [15] Z. Kapidžić, L. Nilsson, H. Ansell, Finite element modeling of mechanically
fastened composite-aluminum joints in aircraft structures, Compos. Struct.
109 (2014) 198–210.
[16] M. Bérot, J. Malrieu, F. Bay, An innovative strategy to create equivalent ele-
Case 5 3-D 8205 350 ments for modelling assembly points in joined structures, Eng. Comput.
MCRS 2562 166 68 (Swansea, Wales) 31 (2014) 453–466.
[17] P.J. Gray, C.T. McCarthy, A highly efficient user-defined finite element for load
distribution analysis of large-scale bolted composite structures, Compos. Sci.
Furthermore, the physical approach proposed to construct the Technol. 71 (2011) 1517–1527.
[18] L.J. Hart-Smith, Bolted joint analysis for composite structures, in: Joining and
reduced bolt model could be extended to different pin-based joining
Repair of Composites Structures, Kansas City, MO, 2004.
technologies. For a double-lap joint it seems appropriate to use [19] W.D. Nelson, B.L. Bunin, L.J. Hart-Smith, Critical joints in large composite air-
3 cylindrical rigid surfaces to represent a bolt pin. For a countersunk craft structure, 1983.
head bolt, a thorough re-analysis of deformation modes of the fas- [20] P.J. Gray, C.T. McCarthy, A global bolted joint model for finite element analysis
of load distributions in multi-bolt composite joints, Compos. Part B: Eng. 41
tener will be required in order to define relevant rigid surface geo- (2010) 317–325.
metries and connector types. [21] M.B. Tate, S.J. Rosenfeld, Preliminary Investigation of the Loads Carried by
Individual Bolts in Bolted Joints, NACA, 1947.
[22] H. Huth, Influence of Fastener Flexibility on the Prediction of Load Transfer
and Fatigue Life for Multiple-row Joints, ASTM, Charleston, SC, USA (1986), p.
References 221–250 , ASTM Special Technical Publication.
[23] F. Alkatan, P. Stephan, A. Daidie, J. Guillot, Equivalent axial stiffness of various
components in bolted joints subjected to axial loading, Finite Elem. Anal. Des.
[1] C.T. McCarthy, M.A. McCarthy, W.F. Stanley, V.P. Lawlor, Experiences with 43 (2007) 589–598.
modeling friction in composite bolted joints, J. Compos. Mater. 39 (2005) [24] T. Sawa, K. Maruyama, On the deformation of the bolt head and nut in a bolted
1881–1908. joint, Bull. JSME 19 (1976) 203–211.
[2] C.T. McCarthy, M.A. McCarthy, Three-dimensional finite element analysis of [25] Abaqus 6.12 Analysis User's Manual, Simulia, 2012.
single-bolt, single-lap composite bolted joints: Part II—effects of bolt-hole [26] C.W. Weaver, J.G. Williams, Deformation of a carbon-epoxy composite under
clearance, Compos. Struct. 71 (2005) 159–175. hydrostatic pressure, J. Mater. Sci. 10 (1975) 1323–1333.
[3] J. Lecomte, C. Bois, H. Wargnier, J.-C. Wahl, An analytical model for the pre- [27] T.V. Parry, A.S. Wronski, Kinking and compressive failure in uniaxially aligned
diction of load distribution in multi-bolt composite joints including hole- carbon fibre composite tested under superposed hydrostatic pressure, J. Mater.
location errors, Compos. Struct. 117 (2014) 354–361. Sci. 17 (1982) 893–900.
[4] V. Caccese, K.A. Berube, M. Fernandez, J. Daniel Melo, J.P. Kabche, Influence of [28] T.J. Vogler, S.-Y. Hsu, S. Kyriakides, Composite failure under combined com-
stress relaxation on clamp-up force in hybrid composite-to-metal bolted pression and shear, Int. J. Solids Struct. 37 (2000) 1765–1791.
joints, Compos. Struct. 89 (2009) 285–293. [29] J.M. Whitney, R.J. Nuismer, Stress fracture criteria for laminated composites
[5] B. Egan, M.A. McCarthy, R.M. Frizzell, P.J. Gray, C.T. McCarthy, Modelling containing stress concentrations, J. Compos. Mater. 8 (1974) 253–265.
bearing failure in countersunk composite joints under quasi-static loading [30] C. Hochard, N. Lahellec, C. Bordreuil, A ply scale non-local fibre rupture cri-
using 3D explicit finite element analysis, Compos. Struct. 108 (2014) 963–977. terion for CFRP woven ply laminated structures, Compos. Struct. 80 (2007)
[6] L. Adam, C. Bouvet, B. Castanié, A. Daidié, E. Bonhomme, Discrete ply model of 321–326.
circular pull-through test of fasteners in laminates, Compos. Struct. 94 (2012)
3082–3091.
[7] V. Roulet, P.-A. Boucard, L. Champaney, An efficient computational strategy for
composite laminates assemblies including variability, Int. J. Solids Struct. 50
(2013) 2749–2757.