Gilles Dowek
email: gilles.dowek@polytechnique.edu

The Undecidability of Unification Modulo σ Alone

Definition 1. A second order context is a context in which all types are of order at most 2.

In λ-calculus, a unification problem a = b is a second order unification problem if its context is second order, the common type of a and b is atomic and all the free variables of a and b have a type of order at most 2.

Proposition 1. There is no algorithm that decides if a second order unification problem has a solution or not.

Proof. See [START_REF] Goldfarb | The Undecidability of the Second-Order Unification Problem[END_REF]. Definition 2. A λσ-term is said to be simple if every of its subterms of the form X[s] is such that s =↑ n for some n. A substitution is said to be simple if for every variable X, θX is simple.

.c p . ↑ n ]) ↓ λσ = (d[c 1 ...c p . ↑ n ]) ↓ σ • if d = (i e 1 .
.. e q ) then the e i 's have an atomic type and we apply the induction hypothesis, Remark 2. We have reduced the second order unification problem to the σ-unification problem. Reducing the full higher order unification problem seems to be much more difficult. Thus σ alone seems to be a formalism that has links with second order languages. In particular in sorts A 1 ...A n ⊢ B, the A i 's and B are types and not sorts and thus the "arrow" ⊢ cannot be nested.

• if d = X[↑ r ]
More precisely σ-calculus seems to be a precise formulation of "second order languages without λ's" as defined for instance in [START_REF] Goldfarb | The Undecidability of the Second-Order Unification Problem[END_REF]. We conjecture the decidability of σ-matching.

Proposition 2 .Remark 1 .

 21 A problem a = b has a solution, if and only if a F = b F has a solution in λ-calculus, if and only if a F = b F has a solution in the image of F , if and only if a F = b F has a simple solution. Proof. See [1]. Proposition 3. Let σ be the substitution mapping every variable X of sort (Γ, A 1 → ... → A n → B) occurring in a = b to the term λ...λY where Y is a new variable of sort (A 1 ...A n .Γ, B). Let ã be the normal form of σa F and b be the normal form of σb F . The problem a = b has a solution if and only if ã = b has a simple solution. Proof. The problem a = b has a solution if and only if a F = b F has a simple solution. If a F = b F has a simple solution mapping X i to λ...λc i then the substitution mapping Y i to c i is a simple solution to ã = b. Conversely if ã = b has a simple solution mapping X i to c i , the substitution mapping X i to λ...λc i is a simple solution to problem a F = b F . The context of the equation ã = b is second order. The contexts of all the variables occurring in this problem are second order and their types is atomic. If X[c 1 ...c p . ↑ n ] occurs in this problem then c i 's have a first order type. Proposition 4. Let a be a normal term well typed of atomic type in a second order context, such that all the variables occurring in a have a second order context and an atomic type and if X[c 1 ...c p . ↑ n ] occurs in this problem then c i 's have a first order type. Let θ be a simple substitution. Then (θa) ↓ λσ = (θa) ↓ σ . Proof. By induction over the structure of a. -If a = (i b 1 ... b n ) then the types of b i 's are atomic. The result follows by induction. -If a = X[c 1 ...c p . ↑ n ], then let d = θX we prove by induction on the structure of d that (d[c 1 ..

  then the result is trivial. Proposition 5. The problem a = b has a solution if and only if ã = b has a solution for σ alone. If a = b has a solution, then ã = b has a simple solution, in λσ and this it has a solution for σ alone. Converselly if ã = b has a solution for σ alone, then it has a solution for λσ and thus a = b has a solution. Corollary 1. σ-unification is undecidable.


Richard Verwaerde

Pierre-Alain Guidault

Pierre-Alain Boucard
email: pierre-alain.boucard@ens-paris-saclay.fr

A non-linear finite element connector model with friction and plasticity for the simulation of bolted assemblies

Keywords: bolted joint, nonlinear user-element, plasticity, frictional contact, connector

Accurate simulation of bolted joints is not always consistent with industrial requirements, since numerous nonlinearities in the vicinity of the bolt can lead to overly expensive calculations. For this reason, commercial finite element (FE) codes preferentially use substitutes for bolts, such as simplified models or connectors. In this paper, a nonlinear FE connector with its identification methodology is proposed to model the behavior of a bolted assembly. The connector model is based on practical design parameters, such as bolt preload, friction coefficients, or plastic material parameters. The connector is based on a separation of the different phenomena governing the macroscopic behavior of a bolted assembly. The axial behavior of the connector reflects the preload effect and the axial stiffness of the assembly, while the tangential behavior of the connector takes into account the friction phenomena that occur in the vicinity of the bolt. At the same time, a plastic behavior law with isotropic hardening for the bolt is proposed. The identification of the connector parameters is performed on a generic elementary one-bolt assembly. The connector is implemented in ABAQUS through a user-element subroutine. Comparisons of quasi-static responses between large-scale full 3D calculations and simulations with connectors on different bolted assemblies are provided. The results obtained are very similar, while observing a significant reduction in CPU time.

Introduction

Over the past few years, the development of new mechanical models and numerical tools has led to improvements in the design of assembled structures. The design of bolted joints is based on making a good estimate of the loads distribution between fasteners. This load distribution depends on the stiffness of the fastener and of the parts, the way in which external forces are introduced and of all physical parameters.

5

The dimensioning of these elements therefore requires rigorous modeling and description. In this context, it is necessary to take into account the preload installed in the assembly, the friction between the parts, the plasticity, etc.

In this work, preload and friction phenomena as well as plasticity within the screw are considered. To take these elements into account, one can imagine using some of the possibilities offered in commercial softwares. [START_REF] Alkatan | Equivalent axial stiffness of various components in bolted joints subjected to axial loading[END_REF] For example, using an elastoplastic beam element would make it possible to represent the screw behavior.

Commercial codes usually enables one to impose bolt preload, generally as a prescribed displacement jump, at an arbitrary bolt cross-section or interface. However, the displacement value corresponding to a desired preload is unknown to engineers if not through an iterative process. Finally, to take friction into account, dedicated connector elements can be added. Modeling many connections in this way is at the cost of tedious, and often cumbersome according to the code, model manipulations. This does not enable one to treat the problem in an integrated way. Moreover, a method of identification for the associated connection model parameters has to be proposed.

The goal of the present work is to fill this gap by proposing a coherent and fully integrated implementation of a nonlinear finite element (FE) connector and its identification methodology. It thus consists in proposing a simplified model that can behaves similarly to a finely discretized bolt and accounts for the nonlinearities that occur within a bolted assembly as previously described. The non-linearities of the initial numerical model are then taken into account within an unique connector.

The model proposed in [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF] already accounts for friction and preload. In particular, Section 3.3 of [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF] details the process of setting the preload: the preload is introduced by giving a force value, and the algorithm proposed computes at the first iteration of the first computational increment the displacement jump to be imposed, as well as the internal stress due to the preload. As en extension of this model, plasticity of the screw is here introduced. The adopted modeling principle consists in separating the phenomena and mechanisms in order to facilitate the identification of the model parameters. This results in parameters with a strong mechanical meaning, which turns out to be very advantageous to conduct a robust bolted connection optimization process. The whole is designed to be integrated in a user element subroutine which is easily adaptable to any FE codes: a first implementation has been proposed in SAMCEF [START_REF] Guidault | Modelling and identification of a nonlinear finite element connector for the simulation of bolted assemblies[END_REF], the enhanced version presented here has been implemented in ABAQUS. Without appropriate tools to take all these non-linearities into account, simplified generic patterns are usually used that lead to structures that are often oversized, although mass saving is crucial for aeronautical structures.

Historically, a first way to apprehend the mechanical behavior of a bolted assembly naturally lies in the realization of experimental campaigns. In the literature, there are numerous experimental set-ups which aim at characterizing the behavior of bolted joints. For example, scanning electron microscopy techniques are used [START_REF] Kounoudji | Role of third body on bolted joints' self-loosening[END_REF] in order to examine the surfaces of the nut and bolt to analyze the initial contact conditions.

A phenomenon that is much studied in the literature and in industry is the phenomenon of self-loosening of bolted joints. Several approaches are thus detailed in [START_REF] Bhattacharya | An investigation on the anti-loosening characteristics of threaded fasteners under vibratory conditions[END_REF][START_REF] Ksentini | Experimental study of bolted joint selfloosening under transverse load[END_REF] in order to provide elements contributing to the understanding of the phenomenon under transverse loading. When the connection is subjected to vibrations, the anti-loosening capacity of many fastening screws is tested on a test bench. Contributions of [START_REF] Delcher | Nouvelle approche pour la caractérisation du dévissage et desserrage des assemblages vissés[END_REF] try to determine a general methodology for carrying out loosening tests. More generally, the study of the identification of friction properties in bolted connections remains crucial to understand the different interactions between each element. The work proposed in [START_REF] De Crevoisier | Identification of in situ frictional properties of bolted assemblies with digital image correlation[END_REF] thus aims at determining the coefficient of friction in which parameters such as surface condition or bolt preload vary.

In many cases, especially in an industrial context, it is interesting to be able to quickly estimate the performance of assemblies in the presizing phases. For this, it is necessary to have analytical models for the representation of assemblies that quickly provide estimates of the mechanical stresses that pass through the assembly. For the calculation of normal equivalent stiffnesses, two distinct approaches can be identified. The first is based on the work of [START_REF] Rötscher | Die maschinenelemente[END_REF] which focuses essentially on the determination of the value of the opening angle of the Rötscher cone as a function of the geometry of the parts and the type of threaded element. The second approach is the one of [START_REF] Rasmusen | Bolted Joints A two body contact problem with friction[END_REF], which is based on the search for an equivalent section of a cylindrical part of the same length as the assembly, of the same Young's modulus as the real parts, subjected to uniform compression, from the energy of elastic deformation. In [START_REF] Alkatan | Equivalent axial stiffness of various components in bolted joints subjected to axial loading[END_REF], an analytical expression based on this model is proposed.

For the formulation of an analytical formula for the tangential stiffness of an elementary bolted joint, different values are found in [START_REF] Tate | Preliminary Investigation of the loads caried by individual bolts in Bolted Joints[END_REF][START_REF] Swift | Development of the Fail-safe Design Features of the DC-10[END_REF][START_REF] Huth | Influence of Fastener Flexibility on the Prediction of Load Transfer and Fatigue Life for Multiple-Row Joints[END_REF][START_REF] Niu | Aircraft structural design[END_REF][START_REF] Cope | Stress intensity determination in lap joints with mechanical fasteners[END_REF]. By applying these formulas on the same assembly, the difference between the different tangential stiffnesses can reach 40%. Even if part of this gap can be explained by the fact that they are drawn from different experimental assemblies and/or based on different modeling hypotheses, this gap is not satisfactory and shows the great diversity of analytical models in the literature.

However, the majority of analytical models in the literature are one-dimensional, and often valid for only one type of solicitation, which makes their generic use delicate and restricts their scope of application.

The improvement of design methods is then oriented towards precise and reliable numerical simulations, which try to take into account all the complex phenomena highlighted experimentally. The work of [START_REF] Gerendt | On the progressive failure simulation and experimental validation of fiber metal laminate bolted joints[END_REF], for example, proposes a framework, based on finite elements, to predict the progressive static failure behavior of bolted assemblies of fiber-metal laminates. The influence of the number of bolted connections, their geometrical arrangement and the clearance between the screw body and the bore for composite materials are research topics that are still important in the academic and industrial communities. Multi-bolted composite connections have thus been studied in [START_REF] Mccarthy | Finite element analysis of effects of clearance on single shear composite bolted joints[END_REF][START_REF] Mccarthy | Three-dimensional finite element analysis of single-bolt, single-lap composite bolted joints: Part ii--effects of bolt-hole clearance[END_REF]. The development of three-dimensional finite element models then demonstrated that increasing play leads to an increase in bolt rotation, a decrease in the contact area in the vicinity of the bolted connections, and a decrease in joint stiffness. Complementarily, a numerical model of the stress distribution in the vicinity of the bore, including a complete discretization of the fastener model, was provided by [START_REF] Egan | Stress analysis of single-bolt, single-lap, countersunk composite joints with variable bolt-hole clearance[END_REF][START_REF] Egan | Modelling bearing failure in countersunk composite joints under quasi-static loading using 3d explicit finite element analysis[END_REF].

Improving prediction capability through 3D numerical simulations of multi-bolt assemblies generates prohibitive calculation costs. There are nevertheless specific calculation strategies, as in [START_REF] Kapidžić | Finite element modeling of mechanically fastened composite-aluminum joints in aircraft structures[END_REF] or in [START_REF] Vadean | Bolted joints for very large bearings-numerical model development[END_REF]. The work of [START_REF] Roulet | An efficient computational strategy for composite laminates assemblies including variability[END_REF] presents an efficient mixed domain decomposition method based on the LATIN method to study the influence of material parameters for 3D assemblies of composite parts with contact and friction. This facilitates the implementation of parametric studies with a large number of configurations [START_REF] Roulet | A parallel strategy for the multiparametric analysis of structures with large contact and friction surfaces[END_REF].

However, a numerical model describing the fastener at a fine scale cannot be used on an assembly with a large number of fasteners. Taking into account the non-linear phenomena inherent to bolted connections generates far too much calculation time in an industrial context. Thus, the use of simplified models that can account for these phenomena, but that do not require the full discretization of the connections is increasingly popular with the community. Thanks to the savings in terms of computing time, the use of connectors has become commonplace in many areas of numerical simulation. Indeed, connectors allow to report in a robust and efficient way on the behavior of links.

For the case of bolted joints, a Global Bolted Joint Model (GBJM) implemented in the ABAQUS software is proposed in [START_REF] Gray | A global bolted joint model for finite element analysis of load distributions in multi-bolt composite joints[END_REF]. This model consists of two linear beam elements, representing the bolt, and two rigid surfaces, representing the contact surfaces, and allows to take into account the preload of the assembly. In addition, the assembled elements are discretized by shell-type elements. The contact between the assembled components, and between the bore and the bolt, is treated by penalization. Friction is only taken into account for the contact between the plates. An improvement of this model is proposed in [START_REF] Gray | A highly efficient user-defined finite element for load distribution analysis of large-scale bolted composite structures[END_REF]. This time, the model presents only one beam element, and is able to represent the non-linear behavior of bolted composite connections up to the failure of the connection. A third node must be added to the simplified model for a double overlapping configuration.

The study proposed in [START_REF] Askri | A reduced fastener model using multi-connected rigid surfaces for the prediction of both local stress field and load distribution between fasteners[END_REF] highlights that the identification of the equivalent stiffnesses of the GBJM was delicate. Moreover, the model proposed by [START_REF] Gray | A global bolted joint model for finite element analysis of load distributions in multi-bolt composite joints[END_REF] with two rigid surfaces connected by an elastic element accounts for shear and bending deformation, but does not take into account the deformation induced by the bolt tension. An original approach to overcome this problem is that of [START_REF] Askri | A reduced fastener model using multi-connected rigid surfaces for the prediction of both local stress field and load distribution between fasteners[END_REF] which introduces an MCRS model (for Multi Connected Rigid Surfaces) to represent the functional surfaces of the assembly and connects them elastically.

A problem common to all the models described is that they do not take into account all the non-linear phenomena occurring in the vicinity of a bolted connection such as frictional contact, bolt plasticity, or large displacements. The recent developments [START_REF] Guidault | Modelling and identification of a nonlinear finite element connector for the simulation of bolted assemblies[END_REF][START_REF] Guidault | Identification of a finite element connector for the simulation of bolted joints[END_REF][START_REF] Guidault | Modélisation et identification d'un connecteur élément fini non-linéaire pour la simulation d'assemblages boulonnés[END_REF] resulting from the CARAB project (Robust Advanced Design of Bolted Joints) propose a new model of connector allowing to take into account the bolt preload and the friction phenomena between the assembled elements.

The model proposed in [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF] proposes a non-linear connector element with 2 nodes, themselves linked to the rest of the model by kinematic couplings. The identification of the parameters is carried out on a 3D numerical model, including the phenomena at the interfaces of the assembled elements. A formulation of friction laws through an elastoplastic analogy is enhanced. The reduction in calculation time is significant on industrial scale structures.

In this work, an improvement of a non-linear FE connector modeling the behavior of a bolted assembly is detailed. The main novelty of this work compared with [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF] is the introduction of the plastic behavior of the screw combined with the friction phenomena in the vicinity of the assembly. The chosen strategy makes it possible to separate the different phenomena and mechanims in order to identify and illustrate their effects independently.

The connector parameters are based on design parameters : bolt preload, friction coefficient, bolt dimensions and material parameters. The contributions of the different elements, such as the bolt behavior or the frictional phenomena at the interface between the bolted elements, are separated. This makes it possible to introduce a non-linear behavior for the bolt, identified as a Timoshenko beam with an elastoplastic behavior, where the threshold function is based on beam quantities under traction/compression, bending and torsion solicitations. The behavior of axial connectors models the effect of bolt preload and axial stiffness [START_REF] Alkatan | Equivalent axial stiffness of various components in bolted joints subjected to axial loading[END_REF].

Then, friction phenomena that occur in the vicinity of the bolt due to preload are taken into account by an elasto-plastic analogy [START_REF] Wriggers | Constitutive Equations for Contact Interfaces[END_REF]. The identification of the behavior is performed on a generic simple elementary 3D single-bolt assembly [START_REF] Guidault | Modelling and identification of a nonlinear finite element connector for the simulation of bolted assemblies[END_REF][START_REF] Guidault | Identification of a finite element connector for the simulation of bolted joints[END_REF]. The connector is developed in ABAQUS through a Fortran user element subroutine [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF][START_REF] Verwaerde | A nonlinear connector element with physical properties for modelling bolted connections[END_REF]. Comparisons between fine-scale 3D computations and simulations with connectors on various bolted assemblies are provided.

The article is structured as follows. A phenomenological approach based on the simulation of a generic bolted assembly is first proposed in Section 2. The effects of friction and plasticity phenomena are then clearly distinguished. A 3D reference model with the different design parameters is presented, as well as the principles of connector modeling. The connector behavior model is detailed in Section 3. The constitutive plastic law is presented and an efficient local integration of this law and its implementation in the global resolution scheme is detailed. Finally, the method is validated in Section 4 on a complex multi-bolted assembly. The response of the model with connectors is compared to a so-called reference model, meshed with 3D elements.

Study of an elementary bolted assembly

The study of an elementary bolted connection is presented below. The phenomena due to the elastoplastic behavior of the bolt are of particular interest.

Phenomenological approach on a generic model

The view cut of the studied elementary bolted joint is presented in Figure 1. Two plates are assembled using a bolt, consisting of a screw and a nut. The bolt has a plastic behavior. The friction coefficients (0.1 between the assembled plates and 0.15 elsewhere) are chosen such that the proposed study is limited to micro-slip phenomena under the screw head, under the nut and between the assembled plates. Indeed, the situation where macroscopic sliding occur between the screw head and the plate (or the nut and the plate) would require a quasi-dynamic or dynamic study, and therefore does not fall within the scope of the work presented in this paper. Screw head and nut remain in adhesion with the plates. The mesh consists of hexahedral quadratic elements, and is composed of 33, 000 degrees of freedom.

First, a preload is applied to the assembly to establish contact conditions. Tensile loading in y-direction is then applied in order to stress the bolt in shear. The loading as well as the boundary conditions considered are shown in Figure 1.

The quantity of interest studied is the displacement tangential jump g T , of the nodes of the assembled elements located near the screw head and the nut. More specifically, the tangential displacement jump g T is defined by the difference between the average displacement of the nodes under the screw head, and the average displacement of the nodes under the nut, extracted as shown in Figure 1.

In order to be able to observe the influence of the plasticity of the bolt independently from the friction phenomena, different simulations with different plasticity parameters are presented. Is is thus proposed to study two bolt behavior models: an elastic model and an elastoplastic model with isotropic strain-hardening.

Depending on the value of the parameters (and more specifically the yield strength), different scenarios are thus presented. The parameters summarized in Table 1 Depending on the geometrical parameters, materials and coefficients of friction involved in the assembly, it is therefore shown that different non-linear phenomena can occur independently of each other. More precisely, two phenomena are identifed. The transition from an adherent state to a so-called micro-slip state governed Extraction of the average displacement of the nodes of the plate by friction phenomena and represented by a slope break in Figure 2a and Figure 2b plays a major role on the global behavior of the assembly [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF]. Simulaneously, the assembly behavior is influenced by the plastic behavior oh the screw. In this study, we are then interested in taking into account the elastoplastic behavior of the bolt. Friction [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF] and plasticity must then be treated in parallel using separate behavioral models.

u x = u y = u z = 0 F = F y x y z

Connector model implementation assumptions

The strategy for modeling the connector from the previous FE model is presented The modeling choices and key assumptions are the following: 180

Bolt substitution Kinematic couplings Connector Connector nodes

• Under small perturbations assumption, it is assumed that the two assembled plates remain in a linear elastic behavior range.

• The preload introduced in the first calculation step is taken into account in the normal behavior of the connector. An elastoplastic behavior with isotropic strain-hardening is chosen for the bolt. The plasticity criterion, detailed in Section 3, considers in particular the tensile/compression, bending, and torsional solicitations.

• In the so-called connector modeling (Figure 3), the friction phenomena are treated only in the user element. Frictionless contact conditions are thus used at the interface of the assembled elements instead of the frictional contact. This assumption is based on the study of the contact area size and the pressure distribution in this area, which shows that the effects of friction are localized near the contact area, and, consequently, can be encapsulated in the connector modeling. This amounts to consider that micro-slip phenomena only occur at the interface of the assembled elements, so that there is always adhesion between the screw head/nut and the plates.

• The interactions between the screw head and the nut with the assembled elements are modeled thanks to classical functionalities available in existing FE software: in an average sense with a Distribution

Coupling or in a more rigid way with a Multi-Point Constraints or a Kinematic Coupling in ABAQUS [START_REF]Abaqus-Analysis-User's-Guide, Abaqus Analysis User's Guide[END_REF]. In this study, it is chosen to use rigid connections, considering that the bending under screw and nut is taken into account in the preliminary identification phase.

• The behavior of the bolt is assimilated to the behavior of an elastoplastic Timoshenko beam with linear isotropic strain-hardening. The parameters governing the plastic evolution will thus be determined directly from the material parameters of the bolt. The detail of the integration of the elastoplastic constitutive law is specified in Section 3.

Phenomena separation

In the following, bold lowercase is used for vectors. A quantity a evaluated at time increment k will be denoted by a (k) . Quantities related to the normal behavior (i.e. in the bolt axis) such as the normal displacement jump g N will be noted with an N subscript, and tangential quantities in the plane normal to the bolt axis such as the displacement jump g T will be noted with a T subscript.

The coordinates of the connector nodes x 1 (resp. x 2 ) are then calculated by adding to their initial position

x 1,(0) (resp. x 2,(0) ) the displacements u 1 (resp. u 2 ). The rotations of the master nodes are similarly defined, so that in the current configuration, φ 1 (resp. φ 2 ) is determined by adding the rotation increment θ 1 (resp.

θ 2 ) to the reference configuration φ 1,(0) (resp. φ 2,(0) ).

Moreover, under the hypothesis of small perturbations, the normal direction of the bolt is considered as fixed, and is then calculated from the initial state, as shown in Figure 4. The normal direction n is then defined under the small-perturbations assumption by:

n ≈ n (0) = x 2,(0) -x 1,(0) x 2,(0) -x 1,(0) (1)
Normal and tangential displacement jumps are then defined with respect to the definition of n by

g = u 2 -u 1 = g n n + g T (2) with g T • n = 0
The treatment of the problem of the frictional interface in the plane normal to n, as well as the identification of the characteristic parameters, are detailed in [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF]. The main principles are recalled in Section 3.5. The main contribution of the present article consists in improving the model proposed in [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF] by intro-215 ducting an elastoplastic Timoshenko beam model to account for plasticity phenomena within the bolt, as highlighted in blue in Figure 4. The phenomenological approach presented has indeed highlighted the fact that, according to the assumptions of Section 2.2, plasticity and friction phenomena could be treated separately.

1 2 1 2 (x 2,(0) , 2,(0) ) (x 2 , 2 ) (x 1 , 1 ) n (0) = x 2,(0) x 1,(0) kx 2,(0) x 1,(0) k n = n (0) u 1 u 2 t c T f s F rictional interface (x 1,(0) , 1,(0) ) Initial configuration Deformed configuration Elastoplastic Timoshenko beam

Connector behavior model 220

The elastoplastic behavior law and its local integration are presented in this section. In particular, a plasticity criterion is presented, based on the beam quantities of tension/compression, bending and torsion.

Elastoplastic constitutive relation with isotropic strain-hardening

The beam behavior law links the beam generalized stresses vector t bolt to the vectors of generalized strain and its plastic part p by:

t bolt = H.( -p ) (3)
In the rest of the study, it is assumed that the generalized stresses required to describe the plasticity phenomenon are the normal stress N , the two bending moments M f y and M f z and the torsional moment M x . This hypothesis corresponds to a pure bending phenomenon, which leads to the nullity of the last two terms of p in (4), which means that the plastic shear deformation is thus neglected. The terms of the equation ( 3) are then written as follows in (4), as a classic Timoshenko-type beam element with the obvious notations where x, y and z correspond to the local frame of the beam and where x is collinear to the neutral fiber of the beam.

t bolt =               N M f y M f z M x T y T z               H =               ES 0 0 0 0 0 0 EI y 0 0 0 0 0 0 EI z 0 0 0 0 0 0 GC 0 0 0 0 0 0 k y SG 0 0 0 0 0 0 k z SG               =                 x = ∂u ∂x χ y = ∂θ y ∂x χ z = ∂θ z ∂x χ x = ∂θ x ∂x y = ∂u y ∂x -θ z z = ∂u z ∂x + θ y                 p =               p x χ p y χ p z χ p x 0 0               (4)
The terms contained in the matrix H depict the elastic beam behavior: E represents the Young's modulus, S its cross-section, I y and I z the bending moments, G the shear modulus, C the torsion constant, k y and k z the shear correction coefficients.

Threshold function and isotropic strain-hardening

The threshold function is based on a limit tensile force N p , limit bending moments M py and M pz and a limit torsional moment M px depending only on the section geometry and material used. This modeling choice is explained by the fact that the normal stresses are much higher than the shear stresses following the introduction of a preload, and it was therefore chosen for a sake of simplicity to not include the shear forces T y and T z in the plasticity criterion.

For a beam characterized by a circular cross-section area S, a radius R and made of material with yield strength σ Y , the characteristic values are expressed using the beam quantities as in [START_REF] Ksentini | Experimental study of bolted joint selfloosening under transverse load[END_REF] :

N p = Sσ y M ey = I y σ y R M py = ϕM ey M ez = I z σ y R M pz = ϕM ez M px = Cσ y (5)
The term ϕ is called plastic shape factor and depends only on the shape of the beam section [START_REF] Bažant | Stability of structures: Elastic, inelastic, fracture and damage theories[END_REF]. In the case of a circular section, it is demonstrates in Appendix A that it is equal to 16 3π .

Isotropic strain-hardening is characterized by the Q(p) function dependent on cumulated plasticity p. A linear isotropic strain-hardening is chosen, such that

Q(p) = S(σ y + Hp) (6)
The term H in ( 6) is defined as the strain-hardening modulus. It is constant in the case of linear strainhardening.

The plasticity criterion P defining the elasticity surface depends on the generalized stresses bolt t bolt and internal variables characterizing plasticity and is inspired by [START_REF] Flejoux | Plasticite dans les poutres[END_REF]. It is then written:

P (t bolt , p) = F (t bolt , p) -Q(p) (7) = N p N N p 2 + M f y M py 2 + M f z M pz 2 + M x M px 2 -Q(p) (8)
The chosen plasticity model is an associated plasticity model: P = 0 defines both the yield surface and the plastic flow. Hill's principle leads to the law of normality, which implies that the vector of generalized plastic strain velocities ˙ p is normal to the yield surface. Moreover, it was assumed that the plasticization was not progressive: the section is either fully elastic or fully plastic.

The constitutive evolution equations can then be written as follows:

˙ p = ṗ ∂P ∂t bolt (9)
In addition, the complementarity conditions (Kuhn-Tucker's conditions) must be satisfied, i.e. ṗ ≥ 0 ; P (t bolt , p) ≤ 0 ; ṗ.P (t bolt , p) = 0 [START_REF] Alkatan | Equivalent axial stiffness of various components in bolted joints subjected to axial loading[END_REF] The p parameter represents the magnitude of the plastic strain increment. It is obtained using a radial return mapping algorithm detailed in Appendix B.

Numerical integration of the plastic law

This section presents the computation of generalized stresses from a generalized strain increment, using a classic backward Euler implicit scheme. It lies in the determination of the plastic multiplier ṗ in [START_REF] Rasmusen | Bolted Joints A two body contact problem with friction[END_REF]. 245

Definition of a plastic flow direction

The diagonal matrix A is introduced such that:

A =                 1 0 0 0 0 0 0 N p M py 0 0 0 0 0 0 N p M pz 0 0 0 0 0 0 N p M px 0 0 0 0 0 0 0 0 0 0 0 0 0 0                 (11)
and A.t bolt is denoted by t bolt , such that:

P (t bolt , p) = t bolt -Q(p) (12)
One can show that:

∂P ∂t bolt = ∂ t bolt ∂t bolt = ∂ A.t bolt ∂t bolt = A T t bolt t bolt (13)
Thus, by finally noting n plas = t bolt t bolt the plastic flow direction, one gets:

∂P ∂t bolt = A T n plas (14)

Trial state and resolution with Newton secant method

The integration on the element is carried out using 3 Gauss points. The generalized constraints t i,(n+1)

are thus evaluated at each Gauss point i (with i ∈ [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF][START_REF] Kounoudji | Role of third body on bolted joints' self-loosening[END_REF]) from the nodal displacement vector of the element.

It is proposed to carry out an elastic prediction and to write the associated trial state, assuming that the considered strain increment remains elastic, i.e.

t bolt,tr i,(n+1) = H( i,(n+1) -p i,(n) ) (15)
The update of generalized constraints vector t bolt i,(n+1) is then achieved from t bolt,tr i,(n+1) . The value of ∆p i for each Gauss point is then computed respectively thanks to a search of zeros by secant Newton of the scalar functions f (∆p i ) detailed in Appendix B. The resolution of this system concludes the local integration at a Gauss point by obtaining the values of ∆p i for the considered generalized strain increment.

Tangent matrix

The algorithm must then return to the global solver a tangent matrix and a residue associated with the bolt element. The following section details the calculations.

The tangent matrix K T of the element for a generalized strain increment is obtained using the evaluation of the elastoplastic tangent operator H ep i at each Gauss point. It is shown in Appendix C that for a pure elastic strain increment at a Gauss point (without any plasticization), the elastoplastic tangent operator is then written

H ep i = H (16)
On the other hand, when plasticization occurs at a Gauss point, i.e. for P (t bolt,tr i,(n+1) , p i,(n) ) > 0, the elastoplastic tangent operator depends on the plastic flow direction by :

H ep i = H.     I d - ∂P ∂t bolt ⊗ ∂P ∂t bolt ∂P ∂t bolt H ∂P ∂t bolt - ∂P ∂p H     (17)
The tangent matrix of the element is then computed by integration over the length of the beam using 3 Gauss points. In a similar way, the nodal forces are computed by integration on the Gauss points. Details of these calculations are recalled in Appendix C .

The local integration algorithm for the plasticity law for a time increment, and the computation of the tangent matrix and the residue for the global solver can be summarized as in Algorithm 1.

Coupling with friction phenomena

In addition to the plastic behavior of the bolt, the friction phenomena in the vicinity of the bolt are taken into account with a radial return mapping algorithm. It is indeed possible to formulate Coulomb's laws (or other behavioral equations for friction) within an algorithm similar to the one used to solve an elastoplastic problem.

The key idea of the elastoplastic approach detailed in [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF] is a split of the tangential slip g T into an elastic (or adhesive) part g e T and a plastic (or sliding) part g s T . A threshold function f s monitoring the transition from the adherent state to the sliding state therefore depends on the resultant of the normal contact forces p N , a friction coefficient µ and the internal tangential force denoted by t T . The update of the tangential stress 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. More precisely, the quantity of t tr T,(n+1) which exceeds the tangential load allowed according to Coulomb's law µ p N,(n+1) is used to adjust loads leading to t T,(n+1) and leads to the increase of the slip (inelastic) part of the relative tangential motion g s T,(n+1) . The global tangent matrix to be returned by the routine is the addition of the tangent matrix K bolt T given in Algorithm 1 for the bolt, and the one presented in [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF] for the frictional interface, while the residual, similar to internal forces, is the sum of the two residuals resulting from these algorithms.

) = t bolt,tr i,(n+1) -Q(p i ) if P i ≤ 0 then p i,(n+1) = p i,(n) p i,(n+1) = p i,(n) H ep i = H else Computation of ∆p i with Appendix B Computation of n plas = t bolt t bolt p i,(n+1) = p i,(n) + ∆p i n plas p i,(n+1) = p i,(n) + ∆p i H ep i = H.     I d - ∂P ∂t bolt ⊗ ∂P ∂t bolt ∂P ∂t bolt H ∂P ∂t bolt - ∂P ∂p H     end if
Compute the tangent matrix for the global solver with Appendix D

K bolt T = L 0 B T HB dx = L 2 3 i=1 ω i B T 1 + ξ i 2 L .H ep i .B 1 + ξ i 2 L
Compute residual for the global solver with Appendix D

F bolt int = L 2 3 i=1 ω i B T L 1 + ξ i 2 H ep i ( i -p i ) 13
The overall algorithm of the connector is thus summarized in Algorithm 2. The identification of the 280 connector parameters introduced in Algorithm 2 is also detailed in [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF].

Algorithm 2 Connector algorithm

Require: p i,(n) , p i,(n) , g s T,(n)

Compute bolt state

Use Algorithm 1

Get K bolt T , F bolt int
Compute frictional interface state (see [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF] for details)

Compute trial state

t tr T,(n+1) = c T g T,(n+1) -g s T, (n+1)
K f r T = c T   P T -P T     I 3 -n T ⊗ n T 0 3 0 3 0 3   P T -P T end if F f r int =         -t T,(n+1) 0 3 t T,(n+1) 0 3        
Compute connector quantities One can observe a micro-sliding phase of bolts 1 and 4 preceding the one of bolts 2 and 3 for the same 300 preload value. The same asynchronous evolution of the sliding of each associated connector is similarly observed. This evolution has been studied analytically in the work of [START_REF] Guillot | Calcul des assemblages vissés assemblages de pièces planes de faibles épaisseurs. partie 1, Techniques de l'ingénieur Assemblages et fixations mécaniques base documentaire[END_REF] and was numerically highlighted in [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF].

K T = K bolt T + K f r T F int = F bolt int + F f r
The importance of taking into account the elasto-plastic behavior of bolts is also illustrated on Figure 6.

The value of the yield strength (σ Y = 800 MPa) was chosen so that the bolts do not plasticize during 305 the preload phase. On the other hand, they plasticize when the prescribed shear loading reaches 1800 N.

Comparisons between models show excellent agreement for both elastic and elastoplastic simulations.

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 5.13 for the elastic case, and equal to 5.08 for the elastoplastic case. This is a satisfactory result, which is mainly explained by a significant decrease of the number of degrees of freedom (gain of 4.2), but also by a reduction of the number of solver iterations thanks to the reduced number of contact interactions between the assembled plate (frictionless contact) as well as the fast integration of frictional and plasticity phenomena at the connector/user-element level compared with the 3D reference model. The test-case presented in this section is a multi-bolted assembly of a manifold. The assembly contains 55 bolts to secure the 4 tubes constituting the manifold. The friction coefficient at each contact surface is chosen equal to 0.1 for the 3D reference simulation, and all contacts are frictionless for the simulation using the connector model. To appreciate the influence of the elastoplastic behavior, results for two different connector models are proposed: a connector model whose beam is purely elastic, and a second model taking into account plasticity phenomena. As a reminder, the assembled elements have a linear elastic behavior.

Model

In each simulation, the elements used are C3D20 elements, which are hexahedral elements of order 2.

Meshes are conformal at the contact surfaces. A 50 kN force F is applied on the top tube of the manifold, and the 4 lower tubes are clamped as shown in Figure 7.

Result of the simulations

Comparisons between the reference model and connector models are made for the proposed load case.

The results are proposed at two different instants: after the preload step, and then, after the applied shear Figure 10a shows the evolution of the axial displacement jump g N after the preload step, and after the load application. In order to highlight the contribution of the connector elasto-plastic behavior, the load with elastoplastic behavior, a so-called elastic model using connectors that only take into account friction phenomena but not those related to plasticity, and a model including connectors capable of reflecting the elasto-plastic behavior of the bolt and the frictional interface contribution. The differences obtained between the two connector models underline the contribution of the implementation of the elasto-plastic behavior in the connector. The normal and the norm of the tangential displacement jumps of the elastic connector model are very small, which explains why green curves are very similar in Figure 10. Thus, it can be observed that for all bolts that have plasticized, the elastic connector model is not able to reflect the load distribution in the structure.

X Y Z F = F y X Y Z Clamping
The local behavior difference of each bolt/connector can be explained by the complex geometry of the model under consideration. For example, it is interesting to note that bolts/connectors 7 and 33 of Figure 8 are the most distanced from the centre of the assembly, and are the most solicited for the two models. In the same way, the bolts/connectors 20 and 47 are more central and are less stressed in both configurations. The influence of the plasticization, and thus a loss of stiffness, of a few bolts on the overall structure is observed.

The figure on the right-hand side of Figure 10a shows, for example, that bolts 2 to 8 have not plasticized because the elastic and plastic curves coincide, whereas bolt 32, for example, has plasticized to a large extent because of its location as shown in Figure 8. The results given by the two models are in good agreement, because load transfers among bolts are thus adequately represented. On the other hand, the purely elastic model is not able to adequately represent the phenomena for bolts where the plastic behavior occurs, and significant differences from the reference model are then observed.

Comparison of the performances

Simulations were performed with 128 GB RAM on 32 cores of a node of two 2.40 GHz Intel Xeon 6148 processors.

Table 3 presents the relative differences on the previous quantities of interest for the applied load. It can be seen in Figure 9 that these maximum values are reached near the clamping for the Von Mises stress, and near the applied loading for the displacement magnitude. The speed-up factors can change significantly, depending on the reference mesh. For example, if the bolts were meshed in a much more refined manner, the gain in terms of the number of degrees of freedom would have been greater, thus increasing the speed-up factors. It is also interesting to note that the connector model does not require a mesh as fine as the reference model to converge, thanks to the substitution of frictional contact conditions by frictionless contact. An additional acceleration factor could then be obtained by meshing the bolted structure more coarsely. However, it was chosen not to modify the initial mesh of the reference model for the presented results.

Quantity

The global Newton solver requires fewer convergence iterations for the connector model than for the reference model. This is due to the fact that some non-linearities (frictional contact, plasticity) are encapsulated and integrated locally at the connector level, and not in the global model.

For the same reason, it can also be observed that the elasto-plastic connector does not require much more CPU time or iterations than the elastic connector. The reason is that the integration of the plastic behavior by searching for zeros by secant Newton is performed locally in the user-element subroutine, and requires only a negligible amount of CPU time.

Profits obtained in terms of CPU time are very satisfactory. The small differences the reference and connector models highlighted, coupled with the profits in computing time attest to the interest and efficiency of the proposed connector.

• σ xx (y) = σ Y for -R ≤ y ≤ 0 • σ xx (y) = -σ Y for 0 ≤ y ≤ R By definition, M f z = S σ xx (y)y dS. (21)
Thus by symmetry,

M pz = 2 π θ=0 R r=0 σ Y r sin(θ) r dθ dr = 4 3 σ Y R 3 (22)
It finally comes

M pz M ez = ϕ = 16 3π
which depends only on the geometrical characteristics of the section.

Thus, with the rule of summation of index i:

n T Dn = d i (Q(p + ∆p) + ∆p d i ) -2 ( t bolt i ) 2 (33) t bolt,tr .n = (Q(p + ∆p) + ∆p d i ) -1 ( t bolt i ) 2 (34)
The scalar equation [START_REF] Verwaerde | A nonlinear connector element with physical properties for modelling bolted connections[END_REF] to be solved in ∆p may then be written:

∆p d i (Q(p + ∆p) + ∆p d i ) -2 ( t bolt,tr i ) 2 = (Q(p + ∆p) + ∆p d i ) -1 ( t bolt,tr i ) 2 -Q(p + ∆p) (35)
By factoring the terms in ( t bolt,tr i

) 2 , one gets:

d i ∆p (Q(p + ∆p) + ∆p d i ) 2 - 1 Q(p + ∆p) + ∆p d i ( t bolt,tr i ) 2 = -Q(p + ∆p). (36)
Thus,

- Q(p + ∆p) (Q(p + ∆p) + ∆p d i ) 2 ( t bolt,tr i ) 2 = -Q(p + ∆p) (37)
One finally obtains the expression of f (∆p) that we have to zero:

f (∆p) = 1 - ( t bolt,tr i ) 2 (Q(p + ∆p) + ∆p d i ) 2 (38)
The function f (∆p) of unknown ∆p has thus as explicit expression

f (∆p) = 0 = 1 - N tr 2 (Q(p + ∆p) + ES∆p) 2 - N 2 p M 2 py M tr f y 2 Q(p + ∆p) + EI y N 2 p M 2 py ∆p 2 - N 2 p M 2 pz M tr f z 2 Q(p + ∆p) + EI z N 2 p M 2 pz ∆p 2 - N 2 p M 2 px M tr x 2 Q(p + ∆p) + GC N 2 p M 2 px ∆p 2 (39)
The resolution scheme is finally based on a zero search of the function f (∆p). This equation is solved by a secant Newton method, which requires two initial points to be given. ∆p is positive, and the maximum value The computation of the tangent operator consists of linking the quantity ṫbolt to the quantity ˙ . In the elastic case, the relationship is simply:

ṫbolt = H ˙ ( 41
)
This relationship is no longer exact in the case where the strain increment is plastic.

In addition, the coherence conditions make it clear that if P (t bolt , p) = 0, then Ṗ (t bolt , p) = 0, which leads The elastoplastic tangent modulus is therefore finally written as follows:

H ep = H. I d - A T nn T A n T AHA T n + SH H (51)
The tangent matrix of the element is then computed by integration over the length of the beam using 3

Gauss points by:

K T = L 0 B T HB dx = L 2 3 i=1 ω i B 1 + ξ i 2 L T .H ep i .B 1 + ξ i 2 L ( 52
)
where ω i and ξ i respectively represent the Gauss weights and Gauss points, and where B symbolizes the matrix linking the nodal displacements U to the generalized strain i at each Gauss point, detailed in Appendix D such that:

i = B i U (53)
In a similar way to get the tangent matrix of the element, the nodal forces are computed by integration on the Gauss points by:

F int = L 2 3 i=1 ω i B T L 1 + ξ i 2 H ep i ( i -p i ) (54)
The functions used to approximate the displacement field within the element are order-one polynomials for traction/compression and torsion, and order-three polynomials for bending.

The test functions to be chosen must satisfy the equilibrium equations without right-hand side, i.e. for example in the case of bending in the plane (xOz):

Ψ + ω = 0 EI y ω -k z SG(Ψ + ω) = 0
where Ψ represents the test functions for the displacements, and ω those for the rotations. The test functions 440 will therefore explicitly depend on the geometrical and material characteristics of the beam. The value of the constants of each shape function is then determined using the boundary conditions leading to the expression of the following shape functions.

By denoting 12EI y k z SGL 2 by Φ y , let us choose the following pairs of test functions (N i (x), N i+4 (x)) corresponding respectively to the degrees of freedom u 1,z , θ 1,y , u 2,z , θ 2,y of the element.

445

Similarly, the pairs of test functions chosen for the bending phenomenon in the plane (xOy) are identical, except for the sign: (N 3 , N 7 ), (-N 4 , N 8 ), (N 5 , -N 9 ) and (-N 6 , N 10 ). For the displacements:

N 3 (x) = 1 1 + Φ y 2 x L 3 -3 x L 2 -Φ y x L + (1 + Φ y ) N 4 (x) = L 1 + Φ y - x L 3 + 4 + Φ y 2 x L 2 - 2 + Φ y 2 x L N 5 (x) = 1 1 + Φ y -2 x L 3 + 3 x L 2 + Φ y x L N 6 (x) = L 1 + Φ y - x L 3 + 2 -Φ y 2 x L 2 + Φ y 2 x L
For the rotations:

N 7 (x) = 6 L(1 + Φ y ) x L 1 - x L N 8 (x) = 1 1 + Φ y 3 x L 2 -(4 + Φ y ) x L + (1 + Φ y ) N 9 (x) = -6 L(1 + Φ y ) x L 1 - x L N 10 (x) = 1 1 + Φ y 3 x L 2 + (-2 + Φ y ) x L
One deduces the expression of the displacement u of a point located on the neutral fiber at the abscissa Finally, by writing the generalized strain vector as:

=                 x = ∂u x ∂x χ y = ∂θ y ∂x χ z = ∂θ z ∂x χ x = ∂θ x ∂x ∂u y ∂x -θ z ∂u z ∂x + θ y                 , (57)
and thanks to (56), the matrix B is then obtained by writing:

=                 ∂ ∂x 0 0 0 0 0 0 0 0 0 ∂ ∂x 0 0 0 0 0 0 ∂ ∂x 0 0 0 ∂ ∂x 0 0 0 ∂ ∂x 0 0 0 -1 0 0 ∂ ∂x 0 1 0                 u = BU (58)

  lead to different states: a plasticization of the bolt during the adhesion phase of the bolted assembly elements called Scenario 1, and plasticization of the bolt during the micro-slip phase called Scenario 2. The evolution of the applied load F as a function of the tangential displacement gap norm ||g T || for the different scenarios are then presented in Figure 2.

Figure 1 :

 1 Figure 1: View cut of the elementary single lap bolted joint

  170

Figure 3 .Figure 2 :

 32 Figure 2: Evolution of the applied load F as a function of the tangential displacement gap norm ||g T || for different values of bolt plasticity parameters

Figure 3 :

 3 Figure 3: Bolt substitution by a two-node connector model and kinematic couplings

Figure 4 :

 4 Figure 4: Connector behavior model with plasticity in the bolt where fs and c T represent the threshold function and the interface stiffness of the frictional behavior introduced in[1]

  if f tr s,(n+1) ≤ 0 then g s T,(n+1) = g s,trT,(n+1)t T,(n+1) = t tr T,(n+1) n+1)µp N,(n+1) c T g s T,(n+1) = g s T,(n) + λn T,(n+1) t T,(n+1) = c T g T,(n+1)g s T,(n+1)

int 4 .Figure 5 :

 45 Figure 5: Description of the 3D reference model and the connector model for the single lap four-bolt joint under shear loading

Figure 6 :

 6 Figure 6: Evolution of the tangential displacement gap for each bolt/connector

Figure 7

 7 shows the loading and the boundary conditions for the test case. The dimensions of the bolts used in this test case are identical to those used in Section 4.1. The bolts numbering for the presentation of the results is shown Figure 8. Each bolt is first preloaded with a value of 10 kN before the application of the external load. The manifold is made of an elastic material with Young's modulus E = 70 GPa and Poisson's ratio ν = 0.3, and the bolts are made of an elasto-plastic material (E = 210 GPa, ν = 0.3), and a yield strength σ y equal to 800 MPa.

Figure 7 :Figure 8 :

 78 Figure 7: Description of the manifold geometry and applied boundary conditions

FigureFigure 9 :

 9 Figure9aintroduces the comparison of the Von Mises stress fields at the end of the loading step. The reference model is presented by masking the bolts, in order to be able to appreciate the behavior of the nodes under the screw head and under the nut. Both models provide very similar results. A similar comparison for the displacement fields is presented Figure9b. The very small differences between the two models validate 340

( b )Figure 10 :

 b10 Figure 10: Normal and norm of the tangential displacement jumps for the 3D reference model and the connector model

Table 3 :

 3 Relative differences on maximum values between the reference and connector modelsIt highlights the fact that the connector model adequately depicts not only the local behavior of each substituted bolt, but also the global behavior of the assembly. The most important stresses are observed in the same areas and are very close for each model. The maximum of the displacement magnitude is also in adequacy for each model. The design of bolted structures with complex geometry and various loads is thus simplified with the proposed connector model. Taking into account the elasto-plastic behavior in the connector is essential when the loads become important and strongly impact the local behavior of the bolted joint.

6 .

 6 of ∆p corresponds to a full plastic increment. The research interval for the (positive) solution of f (∆p) = 0 is thus [0, ∆p 1 ] with ∆p 1 = |F bolt,tr (n+1) ) -Q(p (n) Appendix C : Elastoplastic tangent modulus435

  can finally replace the partial derivatives by their respective expressions, noting that:∂P ∂t bolt = A T n(48)∂P ∂p = -SH for a linear isotropic strain-hardening (49) in order to get:ṫbolt = H. I d -A T nn T A n T AHA T n + SH H ˙(50)

Table 1 :

 1 Bolt material parameters

  Algorithm 1 Integration algorithm for the plasticity law of the bolt for a time increment

	Require: p i,(n) , p i,(n)
	For each Gauss point i
	Compute trial state : t bolt,tr i,(n+1) = H( i,(n+1) -p i,(n) ) Compute trial threshold function : P i (t bolt,tr i,(n+1) , p i,(n)

Table 2 :

 2 Comparison of computation time

		CPU time (s) Number of dof Nb. of Newton iterations
	3D Reference -Elastoplastic	1.21 10 5	1.03 10 6	490
	Connector -Elastoplastic	2.38 10 4	2.44 10 5	396
	3D Reference -Elastic	1.18 10 5	1.03 10 6	471
	Connector -Elastic	2.30 10 4	2.44 10 5	289
	4.2. Manifold test-case
	4.2.1. Presentation

Table 4

 4 presents a comparison of the performance of the different models used.

		Reference model Elasto-plastic connector Elastic connector
	Number of d.o.f.	5.6 M	3.7 M	3.7 M
	Number of Newton solver iterations	256	166	126
	CPU time (s)	8.46 10 5	2.45 10 5	2.35 10 5
	Speed-up	-	3.45	3.6

Table 4 :

 4 CPU time comparison of elasto-plastic reference and connector models simulations for a manifold assembly

Table 5

 5 describes the performances of the connector model in terms of memory space. It must be compare with other connector models, such as the BARC connector resulting from the work of[START_REF] Soulé De Lafont | Conception d'un connecteur élément fini pour la simulation des assemblages boulonnés[END_REF] implemented in SAMCEF[START_REF] Paleczny | Développement et validation d'un connecteur pour les assemblages boulonnés, in: 3e Colloque Supméca Alumni, Les assemblages mécaniques, évolutions récentes et perspectives[END_REF]. The performances of the connector proposed in this work are thus interesting, in particular in terms of data setting and memory required to carry out the simulations.

		Reference model Connector model Gain
	Input file size	169 MB	110 MB	53%
	Output file size	10 GB	9.1 GB	9%
	Memory required	64 GB	43 GB	49%
	Pre-processing time	85 s	50 s	70%

Table 5 :

 5 Performance comparison between the reference and the connector models for the manifold test-case

Acknowledgments

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

Appendix A : Plastic shape factor for a circular beam section

Conclusion

The phenomenological analysis of an elementary bolted assembly and the non-linear phenomena governing its behavior has been studied in order to propose a non-linear FE connector model.

If the consideration of friction phenomena in the connector had already been explained in [START_REF] Verwaerde | A nonlinear finite element connector for the simulation of bolted assemblies[END_REF], plasticity phenomena were not treated. The separation of the different phenomena has been proposed in this article, allowing the integration of the bolt elasto-plastic behavior, in addition to the treatment of the frictional contact at the interfaces. These are treated using a threshold function related to beam quantities. The integration of constitutive laws is introduced with dedicated fast algorithms.

The comparison on bolted structures representative of industrial assemblies between a full 3D reference model and a model using the proposed connector validates the proposed approach. Physical phenomena studied are thus highlighted on models containing several million of degrees of freedom, including complex geometries and loading cases. A significant reduction of the calculation time is observed, while guaranteeing a satisfactory quality of the results, especially in terms of normal and tangential gaps at the bolt level.

The implementation of the connector model in ABAQUS has made it possible to highlight the excellent compromise between calculation time and the accuracy of the obtained results. The relative deviations of the computed displacement and stress fields are indeed less than 2% between models. Future work will focus on the need for industry to further reduce calculation times. The treatment of bolted joints undergoing large displacements or rotations will be addressed in future work. A sensitivity analysis of connector parameters to several input factors (discretization, preload, plate thickness...) is also envisaged. Finally, other non-linear phenomena, such as bolt failure or contact between the screw body and the bore as a result of bolt macro-slip, can be also investigated.

For a beam subject to pure bending by a bending moment M f z , the stress σ xx varies linearly in a section proportionnaly to quadratic bending moment I y according to the relation:

Plasticity initiates when the limit value σ Y is reached. Thus, for a circular beam cross section with radius R, one has:

or

In addition, the section is fully plasticized when :

Appendix B : Radial return mapping algorithm

One proposes to carry out an elastic prediction and to write the associated trial state, i.e.

and therefore:

= t bolt,tr (n+1) -H∆ p

= t bolt,tr (n+1) -H∆p

t bolt (n+1) = t bolt,tr (n+1) -∆p HA T n plas,(n+1) [START_REF] Askri | A reduced fastener model using multi-connected rigid surfaces for the prediction of both local stress field and load distribution between fasteners[END_REF] with ∂P ∂t bolt = A T n plas by definition of n plas . Consequently, one gets:

By multiplying (29) by n plas,(n+1) , one gets:

.n plas,(n+1) -∆p n T plas,(n+1) AHA T n plas,(n+1) [START_REF] Wriggers | Constitutive Equations for Contact Interfaces[END_REF] By introducing the diagonal matrix D = AHA T , and by noting that At bolt (n+1) = Q(p + ∆p), one finally gets:

Let us express the quantities n T Dn and t bolt,tr .n as a function of ∆p in the case of a linear isotropic hardening.

The equation [START_REF] Guidault | Modélisation et identification d'un connecteur élément fini non-linéaire pour la simulation d'assemblages boulonnés[END_REF] gives an expression of n (n+1) (which will be denoted by n in the following to lighten the notations):

Appendix D : Behavior matrix of a Timoshenko beam

The vector of the nodal displacements of the beam element gathers the nodal displacements and rotations of the two master nodes, and is written: