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

Claire Hennuyer, Nicolas Leconte, Bertrand Langrand, Eric Markiewicz

<figure><img src="image_1.png" /></figure>

## To cite this version:

Claire Hennuyer, Nicolas Leconte, Bertrand Langrand, Eric Markiewicz. Toward the modelling of riveted assemblies by super-elements in fast dynamics. Mechanics & Industry, 2014, 15 (2), pp.133-137. 10.1051/meca/2014012. hal-03627075

HAL Id: hal-03627075

https://uphf.hal.science/hal-03627075

Submitted on 15 Apr 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Toward the modelling of riveted assemblies by super-elements in fast dynamics

CLAIRE HENNUYER<fnref n="1" />,<fnref n="2" />,ª, NICOLAS LECONTE<fnref n="2" />, BERTRAND LANGRAND<fnref n="1" /> AND ÉRIC MARKIEWICZ<fnref n="2" />

Received 27 June 2013, Accepted 17 February 2014

**Abstract** – The finite element analysis of the behaviour of airframes subjected to crash or impact loadings requires the use of suitable finite 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 finite elements can not be currently connected together. The paper presented here focuses on how to link these finite elements.

**Key words:** Crashworthiness / riveted assembly / interaction rivet-perforation / perforated super-element

**Résumé** – L'analyse par éléments finis du comportement des structures aéronautiques soumises au crash ou à l'impact requiert l'utilisation d'éléments finis adaptés, en particulier pour la représentation des assemblages rivetés. Afin de prédire la survivabilité de la structure, il est en effet nécessaire de s'intéresser à ces zones d'assemblages puisqu'elles sont propices à des concentrations de contraintes, et donc à l'amorçage de fissures pouvant mener à la ruine de l'appareil. Du fait du caractère local de ces phénomènes, de la disproportion entre l'échelle de l'avion et celle de l'assemblage, et, du nombre important de rivets dans une structure complète (plus de 100 000), des super-éléments pour les fixations et les tôles perforées ont été développés pour une représentation des assemblages adaptée aux calculs de structures. Cependant, ces deux types d'éléments finis ne peuvent être connectés à l'heure actuelle. Les travaux présentés se focalisent donc sur leur association.

**Mots clés** : Résistance à l'impact / assemblage riveté / interaction rivet-perforation / super-élément perforé

## 1 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

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

<footnote n="1">¹ ONERA, The French Aerospace Lab, 59045 Lille, France</footnote>

<footnote n="2">² LAMIH UMR CNRS 8201, University of Valenciennes and Hainaut-Cambrésis, 59313 Valenciennes, France</footnote>

<footnote>ª Corresponding author: claire.hennuyer@onera.fr</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Current modelling of assemblies.</figcaption></figure>

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.

<figure><img src="image_4.png" /><figcaption>Fig. 2. An eight-node perforated super-element.</figcaption></figure>

## 2 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-Muskelishvili (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 $\phi$ and $\psi$. 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 $(\alpha_j, \beta_j)$ parameters. For example, the displacement in the x-direction, for a circular hole of radius $a$, is given by equation (1).

$$u_x^h(R, \theta) = \frac{1}{2\mu} \sum_{j=-N}^{M} \alpha_j [(kR^j + R^{-j}) \cos(j\theta) - j(R^j - R^{-j}) \cos(j - 2)\theta] \quad (1)$$

with $(R, \theta)$ the polar coordinates (the radial coordinate is dimensionless $R = r/a$), $\mu$ the shear modulus of the material and $k$ the Muskelishvili constant ($k = (3 - \nu)/(1 + \nu)$ in plane stress, $k = 3 - 4\nu$ in plane strain, where $\nu$ is the Poisson's ratio).

The fields $u^h$ and $\sigma^h$ can be written more concisely:

$$u^h = Nc \quad \text{and} \quad \sigma^h = Qc \quad (2)$$

where $N$ and $Q$ are matrices of interpolation functions associated with displacements and stresses, respectively, and, $c = (\alpha_j, \beta_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).

$$n_u \ge n_\sigma \ge n_q - r \quad (3)$$

where $n_u$ and $n_\sigma$ 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 = (\alpha_j, \beta_j)$, so $n_u = n_\sigma$. 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-11], it is recommended to take $n_\sigma$ approximately 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-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.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Modelling of the assembly by super-elements.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. A perforated super-element featuring nodes on its perforation.</figcaption></figure>

## 3 Modelling of the rivet/perforation interaction

Loading the perforation of the super-element is only possible by superimposing 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.

### 3.1 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 super-element. 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<fnref n="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).

### 3.2 Results

Here, the influence of the nodes disposition and of the order of the interpolation functions on the accuracy of the K-M's solution is studied. The following five configurations are thus analysed: {$N = M = 4$; 8 pts rect}, {$N = M = 4$; 16 pts rect}, {$N = M = 4$; 8 pts rect + 8 pts perfo}, {$N = M = 8$; 16 pts rect}, {$N = M = 8$; 8 pts rect + 8 pts perfo}.

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 = 74000$ MPa and $ν = 0.3$, for the Young modulus and the Poisson ratio respectively. The value of the shear modulus $μ$ and Muskelishvili constant $k$, required for the computation of the K-M's solution, are $μ = 28461.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

<footnote n="1">¹ The methodology took on this study is similar to the one set up during the development of the original perforated super-element (8 nodes on external boundary) [5].</footnote>
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Comparison between reference and K-M's solutions for a pure shear load case, with the configuration {$N = M = 8$; 8 pts rect + 8 pts perfo}.</figcaption></figure>

found in any load cases<fnref n="2" /> and for each configuration. A little loss of accuracy is still noted when only the number of collection points is increased, or when data are collected on the perforation, but with minor significance. Concerning the increase of the order, it has an influence on the accuracy of the solution in every instance, but with more or less importance according to the analysed configuration. Two particular configurations give unexpected results: {$N = M = 8$; 16 pts rect} for a tension load case (example A) and {$N = M = 8$; 8 pts rect + 8 pts perfo} for a shear load case (example B). The stresses/displacements distribution is globally consistent with the load case, and especially close to the hole. Nevertheless, important differences are noticed in the vicinity of some points located on the boundary of the studied area [-1, 1] x [-1, 1], 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 $\sigma_{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 $\sigma_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 $(\alpha_j, \beta_j)$ parameters which appear to be non significant are responsible for these local differences. In example B, the terms $\beta_{-5} \approx 10^{-4}$ and $\beta_7 \approx 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<fnref n="3" /> but only on numerical solutions.

## 3.3 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 $(\alpha_j, \beta_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 demonstrated 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 $(\alpha_j, \beta_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] x [-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.

## 4 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

<footnote n="2">² In reference [5], it has been showed that each load case can be represented by the K-M formalism where only few parameters take part. For example, the biaxial tension case is completely defined by the parameter $\alpha_1$.</footnote>

<footnote n="3">³ An analytical solution has been generated for each of the four load cases considered, starting from K-M's solution and significant corresponding parameters.</footnote>
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.

*Acknowledgements.* The authors are grateful to the Région Nord-Pas-de-Calais and Onera for their financial support of the DRSC Laboratory (Collisions and Rapid Structural Dynamics) with the participation of Onera/DADS/CRD, UVHC/LAMIH and CNRS.

## References

[1] A. Combescure, F. Delcroix, L. Caplain, S. Espanol, P. Eliot, A finite element to simulate the failure of weld points on impact, Int. J. Impact Eng. 28 (2003) 783-802

[2] B. Langrand, E. Deletombe, E. Markiewicz, P. Drazétic, Riveted joint modelling for numerical analysis of airframe crashworthiness, Finite Elem. Anal. Des. 38 (2001) 21-44

[3] B. Langrand, L. Patronelli, E. Deletombe, E. Markiewicz, P. Drazétic, An alternative numerical approach for full scale characterisation for riveted joint design, Aerosp. Sci. Technol. 6 (2002) 343-354

[4] B. Langrand, A.S. Bayart, Y. Chauveau, E. Deletombe, Assessment of multi-physics FE methods for bird impact modelling – Application to a metallic riveted airframe, Int. J. Crashworthiness 7 (2002) 415-428

[5] N. Leconte, B. Langrand, E. Markiewicz, On some features of a Hybrid-Trefftz displacement element containing a hole, Finite Elem. Anal. Des. 46 (2010) 819-828

[6] N.I. Muskelishvili, Some basic problems of the mathematical theory of elasticity, Noordhoff, Groningen, Holland, 1953

[7] R. Piltner, Special finite elements with holes and internal cracks, Int. J. Numer. Methods Eng. 21 (1985) 1471-1485

[8] R. Piltner, Comment on: "A hybrid-Trefftz element containing an elliptic hole", by Manicka Dhanasekar, Jianjun Han and Qinghua Qin, Finite Elem. Anal. Des. 42 (2006) 1314-1323, Finite Elem. Anal. Des. 43 (2007) 975.

[9] R. Piltner, Some remarks on finite elements with an elliptic hole, Finite Elem. Anal. Des. 44 (2008) 767-772

[10] P. Tong, New displacement hybrid finite element models for solid continua, Int. J. Numer. Methods Eng. 2 (1970) 73-83

[11] T.H.H. Pian, D. Chen, On the suppression of zero energy deformation modes, Int. J. Numer. Methods Eng. 19 (1983) 1741-1752

[12] G.C. Sih, Stress distribution around an ovaloid hole under arbitrary concentrated forces, Appl. Sci. Res. 12 (1964) 378-390

[13] G. Kirsch, Die theorie der elastizität und die bedürfnisse der festigkeitslehre, Z. VDI 42 (1898) 797-807

[14] C. Hennuyer, B. Langrand, Modélisation par superéléments finis de rupture de tôles perforées et rivetées soumises au choc et à l'impact : Bilan de l'année de thèse, Onera Report no. RA 2/19554 DADS, 2012

[15] M. Dhanasekar, J. Han, Q.H. Qin, A hybrid-Trefftz element containing an elliptic hole, Finite Elem. Anal. Des. 42 (2006) 1314-1323

[16] F. Cesari, F.M. Furgiuele, D. Martini, Analisi di strutture piane in materiale eterogeneo con l'elemento finito ibrido, in: XXX AIAS, Alghero, Italy, 2001, pp. 1093-1101

[17] C.H. Yang, A.K. Soh, Modeling of voids/cracks and their interactions, Theor. Appl. Fract. Mech. 38 (2002) 81-101.