Toward the modelling of riveted assemblies by super-elements in fast dynamics

– The ﬁnite element analysis of the behaviour of airframes subjected to crash or impact loadings requires the use of suitable ﬁnite elements, in particular for the modelling of riveted assemblies. In order to predict the structure survivability, it is indeed necessary to focus on these areas because stress concentrations, and consequently crack initiations, which can lead to catastrophic loss of the airplane, are likely to occur. Because of the local nature of the phenomenon, the disproportion between the aircraft and the assembly scale, and the large number of fasteners in a complete structure (more than 100 000), super-elements for the fasteners and for the perforated sheets have been developed in order to suitably model assemblies in structural calculations. However, these two types of ﬁnite elements can not be currently connected together. The paper presented here focuses on how to link these ﬁnite elements.


Introduction
The numerical study of airframes subjected to extreme loads is tightly related to the modelling of the riveted assemblies behaviour because phenomena such as stress localisation can be observed in their neighbourhood. In order to predict the crack initiation and propagation, and therefore, evaluate the structure survivability, it is essential to take into account the riveted assemblies and include their mechanical behaviour in the structural analysis.
Normally, a fine mesh close to the assembly areas would be necessary to describe accurately stress concentrations, but the large difference between the aircraft and a Corresponding author: claire.hennuyer@onera.fr the assembly scale, together with the substantial number of rivets in a complete airplane (more than 100 000), will imply excessive computing costs. Consequently, for several years, research in this domain has focused on the development of macro-elements, i.e., elements with a smaller number of degrees of freedom compared to that of the fine modelling, but with an accurate mechanical behaviour implemented.
At first, macro-elements (1D beam, spring) dedicated to the modelling of fasteners were formulated. The most advanced formulations are able to describe the non-linear behaviour of the rivet until rupture [1,2]. The assembly is thus represented with such an element connected to two conventional shell elements using kinematic conditions (Fig. 1). However, the experimental and computational results comparison has shown that the structural stiffness was overestimated by the numerical model and the plastic strain within the FE shell elements connected to rivets was not efficiently enough predicted so as to initiate the failure of the sheet metal plates [3,4]. The assembly modelling suffers in fact from a lack of representativeness. The structural embrittlement caused by holes in the sheets is not introduced in the conventional shell elements to which rivets macro-elements are connected. The forces are moreover not correctly transferred from the rivet macro-element to the shell elements. The interaction between the fastener and the perforated sheets is improperly taken into account, and, several failure modes (e.g., bearing, pull-through) can not be handled.

A hybrid-Trefftz perforated super-element
In order to enhance the modelling of riveted assemblies, Onera and LAMIH performed research jointly to develop a perforated shell super-element that is able to localise mechanical fields around the hole during fast dynamic loads. At this stage, an eight-node perforated super-element (Fig. 2), able to describe accurately stress concentrations due to the hole in plane elasticity, has been developed [5].
This finite element (FE) is based on a hybrid-Trefftz displacement variational principle. It permits to ensure the compatibility with neighbouring conventional FE, and to reduce the interior domain integrations to boundary integrations. The interpolation functions associated with this super-element come from the Kolosov-Muskhelishvili (K-M) formalism [6] and do account for the influence of the hole inside the element. In this formulation, an analytic expression is obtained for displacements and stresses in accordance with two complex potentials, denoted by φ and ψ. Their expressions depend on the considered problem (i.e., boundary conditions applied to the hole), and these two functions are generally expanded in a Laurent series. Here, considering a perforated plate whose hole is a load-free boundary, the displacements and stresses (known as homogeneous) are expressed in the form of a Laurent series, and depend on (α j , β j ) parameters. For example, the displacement in the x-direction, for a circular hole of radius a, is given by equation (1). The fields u h and σ h can be written more concisely: where N and Q are matrices of interpolation functions associated with displacements and stresses, respectively, and, c = (α j , β j ) represents the vector of the generalized degrees of freedom (dofs) of the special element. The upper and lower bounds of series appearing in equations (1, 2) are fixed in accordance with the number of nodes of the element thanks to the relation (3).
where n u and n σ represent the sufficient number of parameters to define displacements and stresses, n q is the number of physical dofs, and r is the number of rigid body modes. The displacement and stress fields considered in equation (2) are both expressed with the vector c = (α j , β j ), so n u = n σ . The problem considered being bidimensional, n q is twice the number of nodes, and r = 3. With this relation, the minimum number of required parameters is set relatively to the number of nodes of the formulated element. In the literature [7,[9][10][11], it is recommended to take n σ approximatively equal to n q − r. This minimum number does not always guarantee a stiffness matrix with full rank. Full rank may be achieved by suitably increasing the number of parameters. However, too many parameters would make the element become too stiff [9][10][11]. The optimal value of parameters for a given type of element should be found by numerical experimentation. For now, following these theoretical requirements, N and M need to be set to 4 for the formulation of an eight-node perforated element (Fig. 2).
This super-element is implemented into a FE code developed by Onera, and gives very accurate results with low computing costs (few degrees of freedom) in comparison with a fine FE mesh [5]. However, this special FE is formulated in such a way that the perforation is a free analytic boundary. Thus, it is obvious that it is not suitable for the modelling of riveted assemblies, since the perforated element is intended to be linked with a rivet macro-element (Fig. 3), as mentioned in Section 1. Moreover, it does not allow the hole boundary to be subjected to any FE-like load.

Modelling of the rivet/perforation interaction
Loading the perforation of the super-element is only possible by superimposing to the calculated perforated homogeneous solution (1, 2) a particular analytical solution corresponding to a particular perforation loading kind. However, the loads that can be considered are quite restricted [6,7,12]: uniform pressure applied on a part or on the whole hole boundary, constant shear load on the perforation, or concentrated forces applied along the perforation. They are not sufficient to model the generality and the complexity of the fastener/perforated plate system interaction. In order to establish a constant (data) exchange between the fastener and the perforation, it is proposed to develop a finite element internal perforation permitting the interaction between the rivet and the hole boundary of the super-element presented in Section 2.

A perforated super-element featuring nodes on its perforation
A suitable solution to do so is to "materialize" the perforation of the element by nodes formulated on this boundary, and consequently, enable to stress these perforation nodes with the load provided by the rivet superelement. Therefore, it is proposed to study the formulation of a perforated super-element featuring 8 nodes on its external boundary (as in Sect. 2) plus 8 additional nodes placed on its perforation (Fig. 4). Since the order of the Laurent series is linked to the number of element nodes, and, that the proposed FE features 16 nodes, N and M need to be fixed to 8.
A necessary but not sufficient step for the theoretical development of this new FE, is to make sure that the interpolation functions are able to build the displacement and stress fields inside the element when additional nodes are placed on the perforation. The numerical study proposed here 1 consists in collecting information (e.g., displacements) on a reference solution (i.e., an analytical or a numerical solution (FEA)) at given points placed like the 16 element nodes. Then, the K-M's solution and the reference solution are considered equal at these 16 particular points, and consequently, a linear system whose unknowns are (α j , β j ) parameters is formed. All fields are then rebuilt by introducing the identified parameters within the analytical expressions (1, 2).

Results
Here, the influence of the nodes disposition and of the order of the interpolation functions on the ac- In the first part of this study, a perforated plate featuring a centred circular hole of radius a = 0.2 subjected to a far-field uniaxial tensile load σ ∞ = 100 MPa is considered. The material properties are E = 74 000 MPa and ν = 0.3, for the Young modulus and the Poisson ratio respectively. The value of the shear modulus μ and Muskhelishvili constant k, required for the computation of the K-M's solution, are μ = 28 461.5 MPa and k = 2.0769 (plane stress). A reference analytical solution of this problem is provided by Kirsch [13]. By computing the K-M's solution from data collected on Kirsch's solution, it appeared that the value of significant parameters was the same for all studied configurations (α −1 = 10, α 1 = 5). A little decrease of the accuracy is observed when only the order N has been modified (i.e., fixed to 8 instead of 4). It seems to be counter intuitive, but, in fact, it can be suggested that increasing the order means increasing the number of parameters, which can make new terms become significant, even if they are close to zero. A little bias is consequently introduced for each contribution added to the solution (1), which generates differences between the reference solution and K-M's solution. The accuracy of K-M's solution also decreases in a more significant way when displacement values are also collected at points placed on the perforation (and not only on the external boundary), but the results remain accurate enough. More details about the value of these differences can be found in [14].
In the second part of this study, the geometrical and material properties are kept and several load cases are now considered (uniaxial tension, biaxial tension, pure shear, and, simple shear). FE models featuring refined meshes in particular close to the hole, have been generated in order to build a reference numerical solution. Globally, the results obtained through K-M's solution are rather correct. The value of significant parameters is always , and in particular for the computation of stresses (Fig. 5). In example B, when a pure shear load is applied to the plate, a difference of more than 35 MPa is measured for σ xy close to the four corners of the domain: the reference value is 91.4 MPa, and the value obtained is 126.5 MPa (Figs. 5a-5c). More important errors are observed for example A, for instance the calculated value of σ x at point (1,1) is equal to 176.5 MPa instead of 331 MPa, the reference value, the difference is important but very localised. The origin of these errors has been identified. Firstly, some (α j , β j ) parameters which appear to be non significant are responsible for these local differences. In example B, the terms β −5 ≈ 10 −4 and β 7 ≈ 10 −6 are the causes of the observed differences; besides, the reference solution is recovered when these parameters are set to zero (Fig. 5d). Secondly, it has been verified that the accuracy of input data (displacement values at 8 or 16 points) is important. Indeed, these differences are never observed when collecting data on analytical solutions 3 but only on numerical solutions.

Comparison with the literature
The analysis presented above shows clearly problems of accuracy of the K-M's solution which are linked to the increase of the number of (α j , β j ) parameters. It appeared that this subject had been discussed by Dhanasekar et al. [15], and Piltner [8,9]. Piltner was the first to formulate finite elements featuring a hole [7]. Dhanasekar et al. announced that Piltner's solution was unstable because the associated interpolation functions would be sensitive when the order of the series grows. They proposed an enhanced version (i.e., the variational principle is the same, the improvement suggested is related to the interpolation functions expression), and made comparisons. Piltner disproved the results obtained by Dhanasekar et al. and indicated that mistakes had been made in the implementation of its element.
In any case, one can notice that each proposed solution is evaluated and compared at only one or two points near the hole (e.g., calculation of hoop stress). Piltner assures that there is no loss of accuracy with the increase with the number of (α j , β j ) parameters [9]. In fact, he presented more accurate results with a 24-node element than with a 8-node element (cf. Tabs. 2-5 in [9]), but only at two points.
In our study, we endeavoured to analyse each obtained solution in the whole domain [−1, 1]×[−1, 1], and not only at some points, in such a way to anticipate the reconstruction of fields inside the future FE from nodal values. Close to the perforation, all configurations give good results, as in Piltner's articles, but problems related to the order of the Laurent series appeared close to the domain boundary.

Conclusions
A new formulation of a perforated super-element has been proposed to link it with a rivet macro-element, in order to model a riveted assembly suitably in structural calculations. The original perforated element has been developed in such a way that the boundary that represents the hole is analytical, preventing any connection with a rivet super-element. Consequently, the authors propose to "materialize" the hole of the element by nodes. A necessary numerical study, before the theoretical development of this new FE, has been presented in order to evaluate if the associated interpolation functions are able to describe the displacement and stress fields inside the element when additional nodes are placed on perforation.
The results of this study show that the mechanical fields are generally correctly rebuilt, and particularly close to the perforation. Some cases highlight significant discrepancies, but these are very localised, and their origin is now understood (additional series terms and numerical solution). Moreover, it should be noticed that the presented work allows to elucidate a disagreement appearing in the literature. Indeed, it is now clear that this disagreement comes from the locations on which reconstructed data are compared to the reference solution: the rebuilt solution is still accurate when increasing the order in the vicinity of the hole, but it becomes less accurate, due to spurious terms of the series becoming significant, close to the boundary of the domain. It is shown that spurious terms become significant only when the reference solution, on which data are collected, is of numerical nature. Note that the contribution of each mode of the series will be deepened to explain their mechanical meaning. Moreover, a systematic parametric analysis will be included to select meaningful modes.
The works in progress concern the formulation of the new super-element, i.e. the development of the variational principle and the stiffness matrix of the element: additional nodes are placed on the perforation (taking into consideration the corresponding dofs), and the displacement compatibility is only ensured along the external boundary. Because of the discordance between these two requirements, the study of stationary conditions and the analysis of finite elements from the literature (totally compatible FE [16], or, on the contrary, FE without compatibility [17]) lead us to focus on the development of a partially compatible perforated FE. Once the new super-element with internal nodes and FE boundary implemented and fully validated for elastic material behaviour, its extension to non-linear material behaviours will be studied and addressed.