A nonlinear finite element connector for the simulation of bolted assemblies

Fine scale computations of bolted assemblies are generally too costly and hardly tractable within an optimization process. Thus, finite elements (FE) connectors or user-elements are commonly used in FE commercial codes by engineers as substitutes for bolts. In this paper, a non-linear FE connector with its identification methodology is proposed to model the behaviour of a single-bolt joint. The connector model is based on design parameters (bolt prestress, friction coefficient, bolt characteristics...). The axial behaviour of the connector reflects the preload effect and the axial bolt stiffness. The tangential connector behaviour accounts for frictional phenomena that occur in the bolt's vicinity due to preload thanks to an elasto-plastic analogy for friction. Tangential and normal behaviours identification is performed on a generic elementary single bolt joint. The connector has been implemented in ABAQUS through a user-element subroutine. Comparisons of the quasistatic responses between full fine scale 3D computations and 3D simulations with connectors on various bolted assemblies are proposed. Results are in good agreement and a significant gain in terms of CPU time is obtained.

Keywords: Bolted joint; Nonlinear user-element; Connector; Frictional contact; Simplified model

Bolted joints are widely used in many fields of mechanical engineering, such as aerospace or automotive industries. In spite of well-established design rules and assembly procedures, important sources of uncertainty remain due to a high sensitivity of the joint mechanical behaviour to dispersions (tightening procedure, frictional effects, geometrical tolerances, bolt-hole clearance...). Additionally, these structures are subjected to complex loads and are including numerous non-linear effects. The difficulties in taking these effects into account often lead designers to use excessive safety coefficients. For those reasons, the general design method of bolted connections leads to oversized structures, while the weight of aeronautical structures always remains a crucial issue. Thus, the need for enhanced design methods where the emphasis is on mass and cost optimizations must be balanced with predictive simulations that account for all these complex phenomena.

Some of those complicated phenomena in the vicinity of the bolt are studied in the literature, such as the tightening and loosening processes due to shear loading [[1]]. The modelling of the axial behaviour of a bolted assembly subjected to axial load is examined in [[2]] in order to determine an equivalent stiffness depending on various geometrical parameters, material properties and friction coefficients. Moreover, influences of other parameters such as the bolt thread [[3]] or the magnitude and the position of the external load [[4]] on those coefficients were also studied.

Furthermore, the need for a fine description of the geometrical details often requires a very fine FE mesh size and a high number of solver iterations due to non-linear phenomena. A three-dimensional FE method for modelling friction in composite bolted joints is presented in [[5]]. An elaborated numerical analysis of the structural behaviour of single-lap protruding composite joints under tensile loading has been conducted in [[6]]. In this article, a three-dimensional progressive damage FEM model is used, which involves a combination of Hashin's failure Criteria, ply-discount material degradation rules and penalty method. Many methods exist to take into account other modes of degradation in composite bolted joints, such as delamination, matrix cracking or failure in [[7]–[9]]. Assemblies subjected to complex fatigue loads are laid out in [[10]]. Other major phenomena such as secondary bending may have a large impact on fatigue control of assemblies. This phenomenon occurs in single overlapping joints, when the neutral bending line is discontinuous and has eccentricities. Numerical studies of the secondary bending phenomenon on FE models including rivets have been carried out in [[12]]. A 3D FE method to study secondary bending due to the eccentric load path in composite-to-aluminium lap joints has also been developed in [[13]]. Unfortunately, 3D models with very fine FE mesh size are computationally expensive. It may require dedicated computational strategies, such as domain decomposition with parallel computing [[14]], especially in the case of parametric studies with a significant number of configurations [[15]].

Consequently, simplified models are required during the optimization or design phases. In order to model large assemblies with a considerable number of fasteners, alternative computational strategies have been developed. An important part of the simplified models remains one-dimensional models, where bolts and assemblies parts are represented by spring models. They can model load distribution in multi-bolt joints and are able to take into account loss of stiffness due to bearing damage, clearance and location errors [[16]], or adhesion and slip state with plastic deformation-shearing [[17]]. Despite a significant gain in term of computational time in comparison with the initial 3D model, these models are often limited to elementary geometries, and do not take into account non-linear effects in the vicinity of the fastener.

Several approaches studied in literature use equivalent elements to consider the mechanical state of the assembly and reduce the complexity of fasteners in FE models. Common engineering practices make use of two-nodes connectors or user-elements available in commercial FE codes in order to model the behaviour of bolted joints that connect 3D or plate models of the assembled parts. The connector behaviour is generally identified from fine scale 3D computations on an elementary bolted joint and/or from experimental campaigns, which are often costly and where non-linear phenomena remains difficult to identify. A linear connector element which links the six degrees of freedom in each connecting node of element joints is presented in [[18]], where the connector kinematics is based on degrees of freedom decoupling assumption. A static case is described in [[19]], where fastened elements are represented with two-dimensional elastic bodies, and the interaction between the joining element and the elastic body is modelled by normal displacements prescribed on distributed springs. Friction is taken into account with the addition of external tractions in an iterative process. Each fixation is finally represented by a non-linear relation between the transferred force and the relative displacements. The impact of bolt-hole clearance for a quasi-static problem on the load distribution in composite multi-bolt joints is studied in [[20]]. Stiffnesses of the equivalent elements are then analytically determined, assuming the fact that frictional effects are negligible. However, identification is often limited in literature to the elastic part of the behaviour and for a given set of loadings and design parameters. For specific situations, analytical models [[21]–[25]] are also available, but they mainly address linear static or dynamic problems.

More elaborated connector models have been recently developed. Ekh and Schön display structural elements in [[26]], considering the mechanisms governing the load distributing process in the vicinity of the bolt such as out-of-plane distribution, bolt-hole clearance, friction between assembled parts, and fastener flexibility. Comparisons with results obtained experimentally and with detailed FE analyses are in good agreement. The Global Bolted Joint Model (GBJM) presented by Gray et al. in [[27]] puts forward a different modelling strategy for bolted composite joints. Shell elements are employed to form the composite laminates, and the bolt is substituted by a combination of beam elements coupled to rigid contact surfaces, taking into account non-linear behaviours such as bolt-hole clearance, bolt torque, friction between laminates. Despite the robustness and the accuracy of those methods, identification of equivalent element stiffnesses remains difficult, and the physical sense of the identified parameters stays far from corresponding 3D model values.

Graph: Fig. 1 Description of the elementary single lap bolted joint test

An innovative reduced model of a fastener using Multi-Connected Rigid Surfaces (MCRS) has been recently introduced by Askri et al. in [[28]]. The analysis of deformation modes of an assembly in a single lap configuration led the authors to model the fastener with 4 rigid surfaces coupled to 4 nodes with 2 translational and 2 rotational stiffnesses. The equivalent stiffnesses of the elements are based on a single 3D simulation of a bolted joint with a given bolt-hole clearance, bolt preload and friction coefficient. Comparisons between 3D and reduced models for loads leading to contact between the hole and the bolt are in excellent agreement for global and local behaviours of the assembly. However, the suggested reduced model remains in a linear range, and requires a rigidification of surfaces under the screw head and nut.

In this work, a non-linear FE connector with its identification methodology is proposed to model the behaviour of a bolt joint. Model parameters are based on design parameters (bolt preload, friction coefficient, bolt-hole clearance, bolt dimensions...) [[29]]. Axial connector behaviour models the preload effect and the axial bolt stiffness [[2]]. Tangential connector behaviour accounts for frictional phenomena that occur in the bolt's vicinity due to preload thanks to an elasto-plastic analogy for friction [[30]]. Contributions of the different elements, such as the bolt or the interface between the attached plates, are separated. This allows the introduction of different non-linear behaviours for each element, and the separation of the different mechanisms. Identification of tangential behaviour is performed on a generic elementary single 3D bolted joint [[31]]. Model parameters have a strong mechanical sense. The connector has been developed in ABAQUS through a Fortran user-element subroutine [[33]]. Comparisons between fine scale 3D computations and simulations with connectors on various bolted assemblies are proposed.

The article is structured as follows. A phenomenological study of a bolted connection is presented in Sect. 2. A 3D reference model presenting the different design parameters is proposed, as well as the modelling principles of the proposed connector model. The connector behaviour model is detailed in Sect. 3. The linearization of contact conditions and the numerical integration is detailed. An identification strategy is suggested in Sect. 4, in order to obtain a set of parameters for a given bolt preload, friction coefficient and geometrical dimensions of the reference model. The method is validated in Sect. 5 on a more complex case with a bolted assembly of flanged housings.

This section focuses on the phenomenological description of an elementary bolted assembly. More precisely, the framework of the proposed study is defined, as well as the quantities of interest used in the subsequent work.

Graph: Fig. 2 Phenomenological behaviour of a bolted assembly subjected to a shear cycle

The studied elementary bolted joint is presented in Fig. 1a, as well as boundary conditions. Two plates are assembled using a bolt, and frictional contact conditions between the plates and in the vicinity of the screw head and nut are set up. Friction coefficients are different between the assembled plates (0.1) and in the others frictional contact areas (0.15). All these elements are meshed with quadratic hexahedral elements. The 3D reference model studied is composed of 180, 000 degrees of freedom.

Graph: Fig. 3 Bolt substitution by a two-nodes connector model and kinematic couplings

The assembly is first preloaded and subsequently subjected to a shear loading-unloading cycle. The average tangential displacement jump of the nodes of the assembled elements in the vicinity of the screw head and nut, denoted by ${\mathit{g}}_T$ is recorded during the loading cycle, as it is shown in Fig. 1b. More precisely, ${\mathit{g}}_T$ corresponds to the jump of mean displacements extracted on the bearing surfaces under screw head and nut on the assembled plates as depicted in Fig. 1b. Figure 2a shows the evolution of the displacement jump ${\mathit{g}}_T$ as a function of the applied load. Several stages are identified, including:

• An elastic behaviour and micro-slip phenomena between point A and B;
• A macroscopic sliding between the assembled elements, followed by a stiffened behaviour due to contact between the screw body and the hole between point B and C;
• An elastic unloading between point C and D, with a residual displacement jump when the external loading becomes zero, due to frictional phenomena that occur during the whole cycle.

Figure 2 shows a cross-section view of the deformed assembly at points A, B, C and D illustrating these phenomena. The suggested cutting plan is shown in Fig. 2a. In this work, only the micro-slip part is investigated (i.e. the AB segment in Fig. 2a) since industrial bolted joints are a priori sized so that point B is not reached. The macroscopic sliding stage (i.e. the BC segment in Fig. 2a) would require a quasi-dynamic or even dynamic framework and is therefore not within the scope of our study.

Characteristic elements of a 2-nodes FE connector as a substitute for the bolt are shown in Fig. 3. It is divided into coupling zones, connected to master nodes via kinematic couplings, modelling the interactions between the screw head and the nut with the plates. A connector taking into account all frictional phenomena in the vicinity of the assembly finally connects the master nodes.

The connector normal behaviour (in the bolt axis) takes bolt preload into account. The tangential behaviour in the normal plane of the connector axis reflects the bending stiffness of the bolt as well as the friction phenomena at the interface of the assembled plates. Identification of the connector parameters is presented in Sect. 4. All of them reflect a strong mechanical sense and rely on phenomena involved in the vicinity of the assembly.

The modelling choices and main assumptions are the following:

• Bolt parts will remain in their linear elastic range, under small-perturbations assumption;
• Friction phenomena are only taken into account in the connector in order to limit the number of frictional contact conditions thanks to the Rötscher's cone, which gives an approximation of the contact area and pressure distribution, allowing the effects of friction to be located in the vicinity of the bolted area; consequently, frictionless contact conditions will be used elsewhere;
• Kinematic couplings are defined thanks to classical functionalities available in legacy FE softwares: in an average sense with a Distributing Coupling or in a more rigid way with a Multi-Points Constraints in ABAQUS [[34]];
• Connector behaviour is considered equivalent to a Timoshenko beam behaviour in parallel with a one-dimensional friction model as illustrated in Sect. 2.3.

In the following, bold lowercase is used for vectors. Moreover, a quantity a evaluated at time increment n will be denoted by $a_{(n)}$ . Similarly, a quantity $a_i$ will be denoted by $a_{i,(n)}$ . Connector initial and current configurations are presented in Fig. 4.

Graph: Fig. 4 Connector model and kinematics

The location (or position) of a connector node in the current configuration, ${\mathit{x}}_1$ (resp. ${\mathit{x}}_2$ ) is determined by adding the displacement ${\mathit{u}}_1$ (resp. ${\mathit{u}}_2$ ) to the reference location ${\mathit{x}}_{1,(0)}$ (resp. ${\mathit{x}}_{2,(0)}$ ):

Graph

The rotations of the master nodes are similarly defined, so that in the current configuration, ${\mathit{\varphi }}_1$ (resp. ${\mathit{\varphi }}_2$ ) is determined by adding the rotation increment ${\mathit{\theta }}_1$ (resp. ${\mathit{\theta }}_2$ ) to the reference configuration ${\mathit{\varphi }}_{1,(0)}$ (resp. ${\mathit{\varphi }}_{2,(0)}$ ):

2 $$\begin{array}{c}\hfill {\mathit{\varphi }}_1={\mathit{\varphi }}_{1,(0)}+{\mathit{\theta }}_1\text{and}{\mathit{\varphi }}_2={\mathit{\varphi }}_{2,(0)}+{\mathit{\theta }}_2\end{array}$$

Graph

The normal (or axial) direction $\mathit{n}$ is defined under the small-perturbations assumption by:

3 $$\begin{array}{c}\hfill \mathit{n}\approx {\mathit{n}}_{(0)}={\displaystyle \frac{{\mathit{x}}_{2,(0)}-{\mathit{x}}_{1,(0)}}{\Vert {\mathit{x}}_{2,(0)}-{\mathit{x}}_{1,(0)}\Vert }}\end{array}$$

Graph

The normal (resp. tangential) displacement jumps $g_N$ (resp. ${\mathit{g}}_T$ ) are defined using the normal contact direction in the current configuration $\mathit{n}$ and the global displacement jump $\mathit{g}$ , defined as follows:

4 $$\begin{array}{c}\hfill \mathit{g}={\mathit{u}}_2-{\mathit{u}}_1=g_N\mathit{n}+{\mathit{g}}_T\end{array}$$

Graph

by introducing under the small perturbation assumption the projection operator $\Pi (·)=({\mathit{I}}_3-\mathit{n}\otimes \mathit{n})(·)$ on the local tangent plane, defined as the normal plane to $\mathit{n}$ by:

5 $$\begin{array}{c}\hfill \left\{\begin{array}{ccc}\hfill g_N& =& \left({\mathit{x}}_2,-,{\mathit{x}}_1\right)·\mathit{n}-g_{N,(0)}\hfill \\ \hfill & \approx & \left({\mathit{u}}_2,-,{\mathit{u}}_1\right)·{\mathit{n}}_{(0)}\hfill \\ \hfill {\mathit{g}}_T& =& \Pi \left({\mathit{x}}_2,-,{\mathit{x}}_1\right)=({\mathit{I}}_3-\mathit{n}\otimes \mathit{n})\left({\mathit{x}}_2,-,{\mathit{x}}_1\right)\hfill \\ \hfill & \approx & ({\mathit{I}}_3-{\mathit{n}}_{(0)}\otimes {\mathit{n}}_{(0)})\left({\mathit{u}}_2,-,{\mathit{u}}_1\right)\hfill \end{array}\right)\end{array}$$

Graph

where ${\mathit{I}}_3$ corresponds to the dimension 3 identity matrix.

The connector behaviour model and the resulting assumptions are presented in this section. In particular, it is proposed to study independently the influence of friction phenomena between the assembled elements on one side, and the bolt itself on the other side. The identification strategy proposed in Sect. 4 results from this suggested behaviour model.

In the following, two types of behaviours are distinguished: the normal behaviour in the axial direction of the connector, and the tangential behaviour in the plane normal to the connector axis.

Normal behaviour takes into account the influence of preload, as well as tensile and compressive stresses in bolt axis direction. It is thus illustrated by a normal stiffness denoted by $c_N$ . It is also characterized by a internal normal force denoted by $t_N$ and the resultant of normal contact loads, denoted by $p_N$ . Links between $c_N$ , $t_N$ and $p_N$ are detailed in Sect. 3.3.

Following the phenomenological study of an elementary bolted connection carried out in Sect. 2, the tangential behaviour is driven by several independent contributions:

• The tangential behaviour of the bolt, accounted by a beam model with a bending stiffness $c_{\mathit{bolt}}$ (see Fig. 4);
• The shear behaviour of the assembled elements, assimilated to a spring stiffness $c_{\mathit{pl}}$ ;
• Interactions due to contact rubbing at the interface of the assembled elements, assimilated to a spring stiffness $c_{\mathit{int}}$ , which may be seen as a regularization term for the friction law from a numerical point of view.

The identification of the numerical value of the parameters $c_{\mathit{bolt}}$ , $c_{\mathit{pl}}$ and $c_{\mathit{int}}$ is proposed in Sect. 4.

The phenomenological approach displays two different states when the assembly is subjected to shear loading. When preload remains sufficiently significant compared to external load, an adherent (or elastic) state occurs, where there is no sliding between the assembled elements. A spring stiffness $c_{\mathit{pl}}+c_{\mathit{int}}$ , denoted by $c_T$ in the following, influencing only the translation displacements of connector nodes in tangential direction, as well as a beam with a stiffness denoted by $c_{\mathit{bolt}}$ influencing the translation and rotation displacements, are acting in parallel, as shown in Fig. 4.

When loads resulting from friction phenomena at the interface are no longer sufficiently substantial to counterbalance external load, a sliding (or plastic) state occurs, where only the beam contributes to the overall behaviour of the assembly.

A threshold function $f_s$ monitoring the transition from the adherent state to the sliding state therefore depends on the resultant of the contact forces $p_N$ , a friction coefficient $\mu$ and the internal tangential force denoted by ${\mathit{t}}_T$ , as shown in Sect. 3.2.

Considering that the contact between the plates has to be friction-free in order to reduce the computational time, the connector has to take into account the friction phenomena and the dissipative effects they cause. The implementation of friction phenomena at the interface of bolted elements is detailed in this section. Coulomb's law and other equations of friction may be formulated within the framework of elasto-plasticity [[30]].

The key idea of the elasto-plastic approach is a split of the tangential slip ${\mathit{g}}_T$ into an elastic (or adhesive) part ${\mathit{g}}_T^e$ and a plastic (or sliding) part ${\mathit{g}}_T^s$ :

6 $$\begin{array}{c}\hfill {\mathit{g}}_T={\mathit{g}}_T^e+{\mathit{g}}_T^s\end{array}$$

Graph

The simplest possible model for the elastic part of the tangential contact is an isotropic linear elastic relation, which yields:

7 $$\begin{array}{c}\hfill {\mathit{t}}_T=c_T{\mathit{g}}_T^e\end{array}$$

Graph

where $c_T$ is defined as in Sect. 3.1.

The plastic part of the tangential slip ${\mathit{g}}_T^s$ is driven by an evolution equation which may be formulated using standard concepts from elasto-plasticity theory [[35]].

For the classical Coulomb's law, the following threshold function may be used:

8 $$\begin{array}{c}\hfill f_s({\mathit{t}}_T)=\Vert {\mathit{t}}_T\Vert -\mu p_N⩽0\end{array}$$

Graph

which is the plastic slip criterion for a given (positive) contact pressure $p_N$ .

The constitutive evolution equation for the plastic part may be stated in a form of a slip rule as follows:

9 $$\begin{array}{c}\hfill {\dot{\mathit{g}}}_T^s=\lambda {\displaystyle \frac{\partial f_s({\mathit{t}}_T)}{\partial {\mathit{t}}_T}}=\lambda {\mathit{n}}_T\text{with}{\mathit{n}}_T={\displaystyle \frac{{\mathit{t}}_T}{\Vert {\mathit{t}}_T\Vert }}\end{array}$$

Graph

In addition, Kuhn-Tucker's conditions must be satisfied, i.e.

10 $$\begin{array}{c}\hfill \lambda \ge 0,f_s({\mathit{t}}_T)\le 0,\lambda f_s({\mathit{t}}_T)=0.\end{array}$$

Graph

Parameter $\lambda$ represents the magnitude of the plastic slip at each increment. It is obtained using a radial return mapping algorithm detailed in "Appendix A".

Normal behaviour is determined before tangential behaviour, since it results in the determination of the connector internal normal force $t_N$ as well as the estimation of the contact pressure $p_N$ required for the evaluation of the threshold function $f_s$ . The connector must also account for the influence of the preload applied on the assembly.

Graph: Fig. 5 One-dimensional model of normal behaviour

In order to understand how the bolt preload is taken into account, let first consider the simplified one-dimensional example depicted in Fig. 5 where $c_N$ represents the connector normal stiffness and $c_{\mathit{elem}}$ corresponds to the equivalent normal stiffness of the connector environment. The behaviour of the connector environment is assumed to remain in its elastic linear domain in order to alleviate the notations in the following. The appplication to non-linear behaviours is straightforward. In its initial state and before the application of an external load, the assembly is pre-loaded with a preload $P_C$ . The difficulty is the lack of knowledge of the normal stiffness of the assembled elements $c_{\mathit{elem}}$ .

For this purpose, without any external load, the first computation increment is used to estimate this stiffness. A difficulty is that the preload must be imposed through the internal normal force $t_N$ .

The internal normal force in the connector is defined by:

11 $$\begin{array}{c}\hfill t_N=P_C+c_N{g}_N\end{array}$$

Graph

with notations defined in Sect. 2.3. Similarly, the internal force in the assembled elements is equal to $c_{\mathit{elem}}{g}_N$ , assuming a linear behaviour as discussed previously.

The internal force of the problem depicted in Fig. 5 simply reads:

12 $$\begin{array}{c}\hfill f_{N,int}(g_N)=t_N(g_N)+c_{\mathit{elem}}{g}_N.\end{array}$$

Graph

The correction of the displacement jump increment $\delta g_{N,(n)}$ is such that:

13 $$\begin{array}{c}\hfill g_{N,(n)}=g_{N,(n-1)}+\delta g_{N,(n)}\end{array}$$

Graph

and is obtained by solving the following tangent problem (14):

14 $$\begin{array}{c}\hfill \left(c_N,+,c_{\mathit{elem}}\right)\delta g_{N,(n)}=-f_{N,int}(g_{N,(n-1)})\end{array}$$

Graph

In the case of linear elasticity and by considering that $g_{N,(0)}$ is initialized to 0, it converges in one iteration and leads to:

15 $$\begin{array}{c}\hfill \left\{\begin{array}{c}{g}_{N,(1)}=g_{N,(0)}+\delta g_{N,(1)}=-{\displaystyle \frac{P_C}{c_N+c_{\mathit{elem}}}}\hfill \\ t_N(g_{N,(1)})=P_C+c_N{g}_{N,(1)}=P_C{\displaystyle \frac{c_{\mathit{elem}}}{c_N+c_{\mathit{elem}}}}\hfill \end{array}\right)\end{array}$$

Graph

It can therefore be seen that, after this initial step, the normal stress in the bolt $t_N$ is not equal to the preload $P_C$ , but depends on the stiffnesses $c_N$ and $c_{\mathit{elem}}$ . It is nevertheless possible to modify the preload introduced in (11) by multiplying it by the factor ${\displaystyle \frac{c_N+c_{\mathit{elem}}}{c_{\mathit{elem}}}}$ in order to obtain the correct internal force $t_N$ .

Thus, the preload is corrected accordingly in subsequent increments by a proportionality rule, replacing $P_C$ by $P_C{\displaystyle \frac{c_N+c_{\mathit{elem}}}{c_{\mathit{elem}}}}={\displaystyle \frac{P_C^2}{t_{N,(1)}}}$ , where $t_{N,(1)}=t_N(g_{N,(1)})$ , in the definition of the internal force $t_N$ given by (11).

Consequently, for all increments n other than the initial increment $(n>1)$ , $t_N$ is therefore estimated by:

16 $$\begin{array}{c}\hfill t_{N,(n)}={\displaystyle \frac{P_C^2}{t_{N,(1)}}}+c_N{g}_{N,(n)}.\end{array}$$

Graph

With this new definition of $t_N$ , (12) becomes:

17 $$\begin{array}{c}\hfill f_{N,\mathit{int}}=t_N+c_{\mathit{elem}}{g}_N=P_C{\displaystyle \frac{c_N+c_{\mathit{elem}}}{c_{\mathit{elem}}}}+(c_N+c_{\mathit{elem}})g_N\end{array}$$

Graph

In the case of linear elasticity and without external loads, this leads to:

18 $$\begin{array}{c}\hfill \left\{\begin{array}{ccc}\hfill g_{N,(2)}& =& -{\displaystyle \frac{P_C}{c_{\mathit{elem}}}}\hfill \\ \hfill t_{N,(2)}& =& P_C\hfill \end{array}\right)\end{array}$$

Graph

It converges in one iteration to the desired internal normal force $P_C$ .

The resultant of the normal contact force $p_N$ could be determined at each calculation step from the integrated contact pressure on the contact elements located at the interface between the assembled plates in the bolt's vicinity, and, subsequently, returned as an input parameter to the connector routine. For the sake of simplicity, only an estimation of $p_N$ is determined.

The Eq. (15) gives an estimation of the value of $c_{\mathit{elem}}$ by

19 $$\begin{array}{c}\hfill c_{\mathit{elem}}={\displaystyle \frac{t_{N,(1)}}{P_C-t_{N,(1)}}}{c}_N\end{array}$$

Graph

At the initial step, the resultant of the normal contact force $p_N$ is simply defined as follows:

20 $$\begin{array}{c}\hfill p_{N,(1)}=P_C.\end{array}$$

Graph

For the subsequent load increments, one estimates $p_N$ by the following expressions:

21 $$\begin{array}{c}\hfill \begin{array}{ccc}\hfill p_{N,(n+1)}& =& {〈-,c_{\mathit{elem}},g_{N,(n+1)}〉}_{+}\hfill \\ \hfill & =& {〈-,{\displaystyle \frac{t_{N,(1)}}{P_C-t_{N,(1)}}},c_N,g_{N,(n+1)}〉}_{+}\hfill \end{array}\end{array}$$

Graph

where the positive part ${〈·〉}_{+}$ is introduced to account for the fact that the resultant of normal contact force is zeroing if assembled parts separate from each other.

The update of the tangential stress ${\mathit{t}}_{T,(n+1)}$ is carried out by a radial return mapping algorithm based on an integration (backward Euler) of the evolution equation for the plastic slip. The algorithm is given in Algorithm 1 and detailed in "Appendix A". The elasto-plastic tangent modulus computed by the routine is detailed in "Appendix B".

The projection of the trial stresses ${\mathit{t}}_{T,(n+1)}^{\mathit{tr}}$ which do not fulfil the slip condition is graphically depicted in Fig. 6 for a two-dimensional case. As can be seen, the quantity of ${\mathit{t}}_{T,(n+1)}^{\mathit{tr}}$ which exceeds the tangential load allowed according to Coulomb's law $\mu p_{N,(n+1)}$ is used to adjust loads leading to ${\mathit{t}}_{T,(n+1)}$ . Furthermore, the increase of the slip (inelastic) part of the relative tangential motion ${\mathit{g}}_{T,(n+1)}^s$ is shown graphically.

Graph: Fig. 6 Projection of trial stress onto the yield/slip surface

The connector behaviour choice has been detailed in Sect. 3. A reminder of the parameters to be identified and a procedure for identifying the connector parameters are presented in this section.

Graph: Fig. 7 Identification of tangential behaviour

The parameters that must be identified must be consistent with the different loading modes of the assembly:

• Quantities related to the connector normal behaviour: a normal stiffness $c_N$ for the tensile/compressive behaviour and a stiffness $c_{\mathit{tor}}$ for torsional loadings in the direction of bolt axis;
• A bending stiffness $c_{\mathit{bolt}}$ related to the Timoshenko beam behaviour;
• A tangential stiffness $c_T=c_{\mathit{int}}+c_{\mathit{pla}}$ associated with friction phenomena at the interface between the assembled plates (regularization paramter $c_{\mathit{int}}$ ), and with the equivalent shear stiffness $c_{\mathit{pla}}$ of plates;
• An equivalent friction coefficient $\mu$ driving the threshold function associated with micro-slip phenomena.

Some of these stiffnesses may be determined analytically [[2], [21], [37]]. The strategy presented in this paper proposes to identify all these parameters from a 3D reference FE model of a single bolt assembly, assumed to be representative of those used in a given bolted structure.

An example of numerical values resulting from the identification strategy is presented in Sect. 4.2.

Assembly loading modes have been separated into three elementary loads which are tension/compression, torsion and shear. The proposed identification strategy is therefore based on three distinct simulations, which are isolating the contribution of each element to the assembly global behaviour with different boundary conditions and loads.

The connector model remains in the linear domain for tensile/compressive and torsional loads. The identification of $c_N$ and $c_{\mathit{tor}}$ is therefore carried out by respectively measuring the average displacement and average rotation of the plates nodes located under the screw head and under the nut as it is shown in Fig. 1b, respectively when a tensile load and torsion load is imposed.

The remaining parameters are determined using a simulation of a single lap joint under shear loading. Boundary conditions applied on an elementary preloaded bolted assembly are shown in Fig. 7a and b. Figure 7a shows the boundary conditions in the axial bolt direction: the four lateral faces of the two assembled plates are locked in bolt axis direction, preventing an overall rotation of the assembly. Figure 7b represents the boundary conditions in the tangential direction: the lower plate is clamped on one of its lateral surfaces, and a tangential displacement is imposed on the opposite lateral surface of the upper plate. This configuration makes it possible to study the behaviour of the assembly in pure shear, preventing rotational movements of the assembly, but leaving free the local rotation of the bolt due to bending phenomena near the screw head and nut.

For this study, the assembled plates are in aluminium $(E=70\text{GPa},\nu =0.3)$ and the bolt is in steel $(E=210\text{GPa},\nu =0.3)$ .

The response to several loading-unloading cycles of increasing amplitude is presented in Fig. 7c, where average displacement and average rotation measured as presented in Fig. 1b are compared to the resultant force in $\mathbf{y}$ direction which works with prescribed displacement $U_d$ , denoted by F.

Two different states are displayed, corresponding respectively to the adherent and to the micro-slip states.

Only the beam (shown by $c_{\mathit{bolt}}$ ) plays a role during the micro-slip state, while both the bolt and the frictional interface influence the behaviour of the adherent state.

In accordance with the beam behaviour reminded in "Appendix C", tangential load of a Timoshenko beam during the micro-sliding state is expressed as a function of the displacement jump ${\mathit{g}}_T$ , the length of the beam L and the sum of the rotations as

22 $$\begin{array}{c}\hfill {\mathit{t}}_T=c_{\mathit{bolt}}\left({\mathit{g}}_T,-,{\displaystyle \frac{L}{2}},({\theta }_1+{\theta }_2)\right).\end{array}$$

Graph

The $c_{\mathit{bolt}}$ coefficient is then obtained using the Fig. 7d by determining the slope of the behaviour during the micro-slip phase.

Parameter $c_T$ is then identified using Fig. 7c. The average of the discharges slopes, corresponding to the elastic adherent state is removed from the influence of the bolt behaviour previously identified, to determine the behaviour of the rubbing interface.

Finally, the friction coefficient $\mu$ driving the threshold function is determined by dividing the load value of the frictional interface at the transition from the adherent state to the micro-slip state, by the resultant of normal actions $p_N$ at this time increment.

The connector parameters have a strong mechanical meaning, and their numerical values are close to the physical phenomena they represent. They depend on several model parameters, such as geometric dimensions of the bolt, material parameters of the bolted elements, or friction coefficient of the 3D reference model. The model used for identification is parameterized, and identification is automated. Any bolt configuration is easily identified with the suggested strategy.

The identification strategy has been applied to an elementary assembly. The results are resumed in Table 1. A comparison between the connector model and the reference model is presented in Fig. 8. Results are very similar. The slight differences for high loads may be explained by the increase of micro-slip phenomena. However, the results are satisfactory for a large loading range. Bolts used in numerical cases of Sect. 5 are identical, allowing comparison between 3D reference model and model with connectors.

Identified values for connector parameters

Graph: Fig. 8 Comparison between reference model and identified connector model

Graph: Fig. 9 Description of the 3D reference model and the connector model

Graph: Fig. 10 Tangential gap for each bolt/connector versus the magnitude of the resulting load F

The identified friction coefficient is very close to the friction coefficient between the assembled elements use in the 3D reference model, which is equal to 0.1 in our case.

The geometry is described in Fig. 9. The bolts are identical to those used in Sect. 4. Thus, the identified connector parameters used are those of Table 1.

During the first computation step, bolts are subjected to a preload equivalent to $10,000N$ . The assembly is then subjected to a shear loading by prescribing a displacement $U_d$ in y direction at the right-hand end of the upper plate in Fig. 9. Left end of the bottom plate is clamped. The assembled plates are in aluminium $(E=70\text{GPa},\nu =0.3)$ and the bolts are in steel $(E=210\text{GPa},\nu =0.3)$ . For the 3D reference model, the friction coefficient between the plates, under screw heads and under nuts is set to 0.1. As a reminder, all contacts are frictionless contacts for the simulation with connectors.

Elements used for the simulation are hexahedral elements (C3D8 in ABAQUS) (Fig. 10). The simulation is divided into 200 time increments.

The deformed shape of both models at the final time step are presented in Fig. 11. There is a very good agreement between the 3D reference model and the connector model.

Graph: Fig. 11 Comparison of the tangential displacement field between the 3D reference model and the connector model (scale factor = 100)

Evolutions of tangential gap at each bold/connector as defined in Fig. 1 are shown in Fig. 10. There is also a good agreement between the two models.

The value of the cumulated slip at each connector is shown in Fig. 12. The local sliding of connectors 1 and 4 is observed before connectors 2 and 3, and explains slope changes which can be seen in the connector model response. The gains in terms of computation time as well as the number of solver iterations are presented in Table 2. Computations were performed with 16 Gb of RAM memory on a node of two Intel Xeon 6148 2.40 GHz processors. The speed-up observed for this example is equal to 8.42. This is a satisfactory result, which is mainly explained by a significant decrease of model number of degrees of freedom (gain of a factor of 4.15), but also by a reduction of the Newton solver iterations number (gain of a factor of 1.98), thanks to the elimination of friction conditions between the assembled elements.

The test case presented in this section represents a simplified model of a bolted turbo-jet turbine casing. The assembly contains 120 bolts. The bolts are identical to the one identified in Sect. 4, and no additional identification is required. Using the problem cyclic symmetry, only one fifteenth of the problem is modelled, corresponding to a sector of 8 bolts differently preloaded.

Graph: Fig. 12 Cumulated slip at each connector

Comparison of CPU time for the 4-bolts assembly

Graph: Fig. 13 Boundary conditions for turbine casings test-case

For the sake of simplicity, the same material properties as those given in Sect. 5.1 are chosen, i.e. aluminium for plates $(E=70\text{GPa},\nu =0.3)$ and steel for bolts $(E=210\text{GPa},\nu =0.3)$ . Consequently, the same values of connector parameters are retained as given in Table 1. The friction coefficient in each frictional contact area is still equal to 0.1 for the 3D reference simulation, and all contacts are still frictionless for the connector model simulation. The elements used are still C3D8 elements. Tensile and torsional loading-unloading are prescribed with 200 time increments. Simulations are performed with the same computer configuration than the one of example presented in Sect. 5.1.

The lower housing is clamped on its lower surface, and two different loads are applied, as it is represented in Fig. 13: a tensile loading, where the upper housing is subjected to loading along the axis of symmetry, and a torsional loading, where the upper housing is subjected to a torque along the same axis of symmetry. Different preload values within the range of $10\text{kN}\pm 30\%$ are considered for the 8 bolts, as given in Table 3.

It is worth noting at the fact that even though different preload values are here used, connector parameter identification does not have to be performed for each preload value. Consequently, the same values of connector parameters identified in Sect. 4 for a preload of $10\text{kN}$ have be used regardless of the preload used subsequently.

Preloads applied to the bolts

Comparisons between the 3D reference model and the connector model are made for the two different load cases. Figure 14 shows the evolution of the axial displacement jump for Bolt and Connector 1 as a function of the resultant force in $\mathbf{z}$ direction which works with prescribed displacement $u_d$ . Only the result for Bolt 1 is presented since the responses for the other bolts are similar. Figure 15 shows the evolution of tangential gap $∥{\mathit{g}}_T∥$ of each bolt as a function of the magnitude of the torque resulting from the prescribed rotation ${\theta }_d$ . The results given by the two models are in good agreement.

Graph: Fig. 14 Normal gap for Bolt/Connector 1 compared to the magnitude of axial load resulting from prescribed displacement ud

Graph: Fig. 15 Tangential gap for each bolt/connector versus magnitude of the prescribed torque

Performances in terms of computation time are presented in Table 4. An interesting speed-up of about $2-2.5$ for both loading cases is obtained.

The phenomenological analysis of an elementary bolted assembly and the non-linear phenomena governing its behaviour has been investigated in order to propose a non-linear FE connector model, as well as a dedicated identification methodology.

Non-linear phenomena such as frictional contact at the interface of the bolted elements are taken into account only in the proposed connector model, allowing a significant computational time reduction, while ensuring a satisfactory quality of the results.

CPU time comparison for the turbine casing assembly

The comparison between a model of a flanged housing portion with connectors and an entire 3D reference model validates the suggested approach. Encouraging computation time gains, ranging from 2 to 8 were observed, while maintaining a satisfactory quality of results.

Future work will focus on the industrial need to reduce calculation times even further, in order to simulate the addition of plasticity in bolts. A sensitivity analysis of the connector parameters to several input factors (discretization, preload, plate thickness...) will be also investigated. The case of bolted assemblies undergoing large displacements or rotations will be also addressed in future works.

This work was performed using HPC resources from the "Mésocentre" computing center of CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Île-de-France (http://mesocentre.centralesupelec.fr/).

• Increment of tangential gap within the time step $\Delta t_{n+1}$
• As specified in (5), the expression of the tangential displacement jump ${\mathit{g}}_T$ is given by
• 23 $$\begin{array}{c}\hfill \Delta g_{T,(n+1)}=({\mathit{I}}_3-{\mathit{n}}_{(0)}\otimes {\mathit{n}}_{(0)})\left(\Delta ,u_{2,(n+1)},-,\Delta ,u_{1,(n+1)}\right)\end{array}$$

Graph

• in the case of small pertubations.
• Computation of the elastic trial state from (7) and evaluation of the slip criterion (8) at time $t_{n+1}$ . The trial state consists in assuming that the considered increment is purely elastic, and then correcting it according to the value of the obtained threshold function. Since the total slip ${\mathit{g}}_{T,(n+1)}={\mathit{g}}_{T,(n)}+\Delta {\mathit{g}}_{T,(n+1)}$ is decomposed into an elastic and a plastic part, one has:
• 24 $$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\mathit{t}}_{T,(n+1)}^{\mathit{tr}}& =& c_T\left({\mathit{g}}_{T,(n+1)},-,{\mathit{g}}_{T,(n)}^s\right)\hfill \\ \hfill & =& {\mathit{t}}_{T,(n)}+c_T\Delta {\mathit{g}}_{T,(n+1)}\hfill \end{array}\end{array}$$

Graph

• Here the vector ${\mathit{t}}_{T,(n)}=c_T\left({\mathit{g}}_{T,(n)},-,{\mathit{g}}_{T,(n)}^s\right)$ is the tangential force at time $t_n$ . A value of the slip criterion which fulfills $f_{s,(n+1)}^{\mathit{tr}}⩽0$ indicates stick. For $f_{s,(n+1)}^{\mathit{tr}}>0$ sliding occurs in the tangential direction, and a return mapping of the trial tractions to the slip surface has to be performed. For Coulomb's model, one simply has:
• 25 $$\begin{array}{c}\hfill f_{s,(n+1)}^{\mathit{tr}}=\Vert {\mathit{t}}_{T,(n+1)}^{\mathit{tr}}\Vert -\mu p_{N,(n+1)}\end{array}$$

Graph

• where $p_{N,(n+1)}=c_N|g_{N,(n+1)}|$ (since $g_{N,(n+1)}<0$ when contact occurs) and where $c_N$ is the normal connector stiffness.
• Return mapping procedure is derived from the time integration algorithm. In case the implicit Euler scheme is applied to approximate the evolution equation (9), one obtains:
• 26 $$\begin{array}{c}\hfill {\mathit{g}}_{T,(n+1)}^s={\mathit{g}}_{T,(n)}^s+\lambda {\mathit{n}}_{T,(n+1)}\end{array}$$

Graph

• with

• 27 $$\begin{array}{c}\hfill {\mathit{n}}_{T,(n+1)}={\displaystyle \frac{{\mathit{t}}_{T,(n+1)}}{\Vert {\mathit{t}}_{T,(n+1)}\Vert }}\end{array}$$

Graph

• With the standard arguments regarding the projection schemes [[36]], one obtains:
• 28 $$\begin{array}{c}\hfill \left\{\begin{array}{ccc}\hfill {\mathit{t}}_{T,(n+1)}& =& {\mathit{t}}_{T,(n+1)}^{\mathit{tr}}-\lambda c_T{\mathit{n}}_{T,(n+1)}\hfill \\ \hfill {\mathit{n}}_{T,(n+1)}& =& {\mathit{n}}_{T,(n+1)}^{\mathit{tr}}\hfill \end{array}\right)\end{array}$$

Graph

• with

• 29 $$\begin{array}{c}\hfill {\mathit{n}}_{T,(n+1)}^{\mathit{tr}}={\displaystyle \frac{{\mathit{t}}_{T,(n+1)}^{\mathit{tr}}}{\Vert {\mathit{t}}_{T,(n+1)}^{\mathit{tr}}\Vert }}\end{array}$$

Graph

• The multiplication of (28) by ${\mathit{n}}_{T,(n+1)}$ yields the condition from which $\lambda$ may be computed:
• 30 $$\begin{array}{c}\hfill \kappa (\lambda )=\Vert {\mathit{t}}_{T,(n+1)}^{\mathit{tr}}\Vert -\mu p_{N,(n+1)}-c_T\lambda =0\end{array}$$

Graph

• since $\Vert {\mathit{t}}_{T,(n+1)}\Vert -\mu p_{N,(n+1)}=0$ on the yield surface. Thus, one gets:
• 31 $$\begin{array}{c}\hfill \lambda ={\displaystyle \frac{1}{c_T}}\left(\Vert ,{\mathit{t}}_{T,(n+1)}^{\mathit{tr}},\Vert -\mu ,p_{N,(n+1)}\right)\end{array}$$

Graph

• Once $\lambda$ is known, one obtains the explicit results for Coulomb's model:
• 32 $$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\mathit{t}}_{T,(n+1)}& =& \mu p_{N,(n+1)}{\mathit{n}}_{T,(n+1)}^{\mathit{tr}}\hfill \\ \hfill {\mathit{g}}_{T,(n+1)}^s& =& {\mathit{g}}_{T,(n)}^s+\lambda {\mathit{n}}_{T,(n+1)}^{\mathit{tr}}\hfill \end{array}\end{array}$$

Graph

• This update completes the local integration algorithm for the frictional interface law.

Let us consider the yield surface described in (25):

33 $$\begin{array}{c}\hfill f_s({\mathit{t}}_T,p_N)=\Vert {\mathit{t}}_T\Vert -\mu p_N\end{array}$$

Graph

Thus:

34 $$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\dot{f}}_s({\mathit{t}}_T,p_N)& =& {\displaystyle \frac{\partial f_s}{\partial {\mathit{t}}_T}}{\dot{\mathit{t}}}_T+{\displaystyle \frac{\partial f_s}{\partial p_N}}{\dot{p}}_N\hfill \\ \hfill & =& {\displaystyle \frac{\partial f_s}{\partial {\mathit{t}}_T}}{c}_T\left({\dot{\mathit{g}}}_T,-,{\dot{\mathit{g}}}_T^s\right)+{\displaystyle \frac{\partial f_s}{\partial p_N}}{\dot{p}}_N\hfill \\ \hfill & =& {\displaystyle \frac{\partial f_s}{\partial {\mathit{t}}_T}}{c}_T\left({\dot{\mathit{g}}}_T,-,\lambda ,{\displaystyle \frac{\partial f_s({\mathit{t}}_T)}{\partial {\mathit{t}}_T}}\right)+{\displaystyle \frac{\partial f_s}{\partial p_N}}{\dot{p}}_N\hfill \end{array}\end{array}$$

Graph

The consistency condition ( $\lambda {\dot{f}}_s({\mathit{t}}_T,p_N)=0$ if $f_s({\mathit{t}}_T)=0$ ) enables one to determine $\lambda$ :

35 $$\begin{array}{c}\hfill \lambda ={\displaystyle \frac{c_T{\displaystyle \frac{\partial f_s}{\partial {\mathit{t}}_T}}{\dot{\mathit{g}}}_T+{\displaystyle \frac{\partial f_s}{\partial p_N}}{\dot{p}}_N}{c_T{\displaystyle \frac{\partial f_s}{\partial {\mathit{t}}_T}}{\displaystyle \frac{\partial f_s}{\partial {\mathit{t}}_T}}}}\end{array}$$

Graph

Since ${\dot{\mathit{g}}}_T^s=\dot{\gamma }{\displaystyle \frac{\partial f_s({\mathit{t}}_T)}{\partial {\mathit{t}}_T}}=\dot{\gamma }{\mathit{n}}_T=\lambda {\mathit{n}}_T$ , one deduces that:

36 $$\begin{array}{c}\hfill \lambda ={\displaystyle \frac{c_T{\mathit{n}}_T^T{\dot{\mathit{g}}}_T+{\displaystyle \frac{\partial f_s}{\partial p_N}}{\dot{p}}_N}{c_T{\mathit{n}}_T^T{\mathit{n}}_T}}={\displaystyle \frac{c_T{\mathit{n}}_T^T{\dot{\mathit{g}}}_T+{\displaystyle \frac{\partial f_s}{\partial p_N}}{\dot{p}}_N}{c_T}}\end{array}$$

Graph

For the Coulomb's friction model, one obtains:

37 $$\begin{array}{c}\hfill \lambda ={\displaystyle \frac{c_T{\mathit{n}}_T^T{\dot{\mathit{g}}}_T-\mu {\dot{p}}_N}{c_T}}\end{array}$$

Graph

This provides the expression for ${\dot{\mathit{t}}}_T$ as a function of ${\dot{\mathit{g}}}_T$ and ${\dot{g}}_N$ :

38 $$\begin{array}{c}\hfill \begin{array}{ccc}\hfill {\dot{\mathit{t}}}_T& =& c_T\left({\dot{\mathit{g}}}_T,-,{\dot{\mathit{g}}}_T^s\right)\hfill \\ \hfill & =& c_T{\dot{\mathit{g}}}_T-c_T\lambda {\mathit{n}}_T\hfill \\ \hfill & =& c_T{\dot{\mathit{g}}}_T-{\displaystyle \frac{c_T{\mathit{n}}_T\left(c_T,,{\mathit{n}}_T^T,,{\dot{\mathit{g}}}_T\right)-\mu c_T{\dot{p}}_N{\mathit{n}}_T}{c_T}}\hfill \\ \hfill & =& c_T\left({\mathbf{I}}_{\mathbf{3}},-,{\mathit{n}}_T,\otimes ,{\mathit{n}}_T\right){\dot{\mathit{g}}}_T+\mu c_N\mathrm{sign}\left(g_N\right){\dot{g}}_N{\mathit{n}}_T\hfill \end{array}\end{array}$$

Graph

where ${\mathit{n}}_T={\displaystyle \frac{{\mathit{t}}_T}{\Vert {\mathit{t}}_T\Vert }}$ . Note that $\text{sign}(g_N)=-1$ .

Note that the frictional term (second term in right-hand side) makes of the matrix which links ${\dot{\mathit{t}}}_T$ to the normal gap increment ${\dot{g}}_N$ , non-symmetric. This is because the Coulomb's law of friction can be viewed as a non-associative constitutive model.

Tangential slip increment and gap increment may be written in matrix form. One can express the variation $\delta g_N{\mathit{n}}_T$ and $\delta {\mathit{g}}_T$ in terms of connector nodes displacement vector variation $\delta \mathit{g}$ as expressed in (5) as:

39 $$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \delta g_N{\mathit{n}}_T& =& \left(\left(\delta ,{\mathit{x}}_2,-,\delta ,{\mathit{x}}_1\right),·,\mathit{n}\right){\mathit{n}}_T\hfill \\ \hfill & =& \left(\begin{array}{cc}-{\mathit{n}}_T\otimes \mathit{n}& {\mathit{n}}_T\otimes \mathit{n}\end{array}\right)\left(\begin{array}{c}\delta {\mathit{x}}_1\\ \delta {\mathit{x}}_2\end{array}\right)\hfill \\ \hfill & =& \left(\begin{array}{cccc}-{\mathit{n}}_T\otimes \mathit{n}& {\mathbf{0}}_{3\times 3}& {\mathit{n}}_T\otimes \mathit{n}& {\mathbf{0}}_{3\times 3}\end{array}\right)\delta \mathit{g}\hfill \\ \hfill & =& {\mathit{N}}_T\delta \mathit{g}\hfill \end{array}\end{array}$$

Graph

where ${\mathbf{0}}_{n\times m}$ represents a n by m null matrix. With the definition of tangential gap given in (5), one gets:

40 $$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \delta {\mathit{g}}_T& =& ({\mathbf{I}}_{\mathbf{3}}-\mathit{n}\otimes \mathit{n})\left(\delta ,{\mathit{x}}_2,-,\delta ,{\mathit{x}}_1\right)\hfill \\ & =& \left(\begin{array}{cc}-\left({\mathbf{I}}_{\mathbf{3}},-,\mathit{n},\otimes ,\mathit{n}\right)& {\mathbf{I}}_{\mathbf{3}}-\mathit{n}\otimes \mathit{n}\end{array}\right)\left(\begin{array}{c}\delta {\mathit{x}}_1\\ \delta {\mathit{x}}_2\end{array}\right)\hfill \\ & =& \left(\begin{array}{cccc}-\left({\mathbf{I}}_{\mathbf{3}},-,\mathit{n},\otimes ,\mathit{n}\right)& {\mathbf{0}}_{3\times 3}& {\mathbf{I}}_{\mathbf{3}}-\mathit{n}\otimes \mathit{n}& {\mathbf{0}}_{3\times 3}\end{array}\right)\delta \mathit{g}\hfill \\ & =& \mathit{M}\delta \mathit{g}\hfill \end{array}\end{array}$$

Graph

One deduces from (38) that:

41 $$\begin{array}{c}\hfill \begin{array}{ccc}\hfill \Delta {\mathit{t}}_T& =& c_T\left({\mathbf{I}}_{\mathbf{3}},-,{\mathit{n}}_T,\otimes ,{\mathit{n}}_T\right)\Delta {\mathit{g}}_T+\mu c_N\mathrm{sign}\left(g_N\right)\Delta g_N{\mathit{n}}_T\hfill \\ \hfill & =& \left[c_T,\left({\mathbf{I}}_{\mathbf{3}},-,{\mathit{n}}_T,\otimes ,{\mathit{n}}_T\right),\mathit{M},+,\mu ,,c_N,,\mathrm{sign},\left(g_N\right),{\mathit{N}}_{\mathit{T}}\right]\Delta \mathit{d}\hfill \end{array}\end{array}$$

Graph

and, finally:

42 $$\begin{array}{cc}& \Delta {\mathit{t}}_T·\delta {\mathit{g}}_T=\delta {\mathit{g}}_T^T\Delta {\mathit{t}}_T=\hfill \\ \hfill & \delta {\mathit{d}}^T\left({\mathit{M}}^T,\left[c_T,\left({\mathbf{I}}_{\mathbf{3}},-,{\mathit{n}}_T,\otimes ,{\mathit{n}}_T\right),\mathit{M},+,\mu ,,c_N,,\mathrm{sign},\left(g_N\right),{\mathit{N}}_T\right]\right)\Delta \mathit{d}\hfill \\ \hfill & =\delta {\mathit{d}}^T{\mathit{K}}_T\Delta \mathit{d}\hfill \end{array}$$

Graph

Note that matrix ${\mathit{K}}_T$ is nonsymmetric.

As a reminder, a quantity a denoted by $a_{(n)}$ corresponds to the value of a evaluated at time increment n. However, to alleviate the notations, the normal ${\mathbf{n}}_{(n+1)}$ is noted only by $\mathbf{n}$ and stays equal to ${\mathbf{n}}_{(0)}$ under the small perturbations assumption.

The different steps of the connector model integration are summarized in Algorithm 1. ${\mathbf{P}}_T$ represents the matrix for passing from the global frame to the local frame of the connector.

Graph

For a two-dimensional problem, the stiffness matrix of a Timoshenko beam ${\mathbf{K}}_{\mathit{bolt}}$ is expressed as a function of the normal stiffness $c_N$ , the bending stiffness $c_{\mathit{bolt}}$ , its length L and the shear coefficient $\Phi$ depending on the section geometry as:

$$\begin{array}{c}\hfill {\mathbf{K}}_{\mathit{bolt}}=\left(\begin{array}{cccccc}{c}_N& 0& 0& -c_N& 0& 0\\ 0& c_{\mathit{bolt}}& c_{\mathit{bolt}}\frac{L}{2}& 0& -c_{\mathit{bolt}}& c_{\mathit{bolt}}\frac{L}{2}\\ 0& c_{\mathit{bolt}}\frac{L}{2}& c_{\mathit{bolt}}\frac{(4+\Phi )L^2}{12}& 0& -c_{\mathit{bolt}}\frac{L}{2}& c_{\mathit{bolt}}\frac{(2-\Phi )L^2}{12}\\ -c_N& 0& 0& c_N& 0& 0\\ 0& -c_{\mathit{bolt}}& -c_{\mathit{bolt}}\frac{L}{2}& 0& c_{\mathit{bolt}}& -c_{\mathit{bolt}}\frac{L}{2}\\ 0& c_{\mathit{bolt}}\frac{L}{2}& c_{\mathit{bolt}}\frac{(2-\Phi )L^2}{12}& 0& -c_{\mathit{bolt}}\frac{L}{2}& c_{\mathit{bolt}}\frac{(4+\Phi )L^2}{12}\end{array}\right)\end{array}$$

Graph

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.