# Simulation dynamique non linéaire de structures assemblées par jonctions boulonnées par une méthode de balance harmonique adaptative

Vincent Jaumouillé, Jean-Jacques Sinou, Benoît Petitjean

► To cite this version:

Vincent Jaumouillé, Jean-Jacques Sinou, Benoît Petitjean. Simulation dynamique non linéaire de structures assemblées par jonctions boulonnées par une méthode de balance harmonique adaptative. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01421854

HAL Id: hal-01421854

https://hal.science/hal-01421854v1

Submitted on 23 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Simulation dynamique non linéaire de structures assemblées par jonctions boulonnées par une méthode de balance harmonique adaptative

V. Jaumouillé¹,², J-J. Sinou¹, B.Petitjean²

¹ Laboratoire de Tribologie et Dynamique des Systèmes
Ecole Centrale de Lyon
36 Avenue Guy de Collongue, F-69134 Ecully Cedex
{vincent.jaumouille, jean-jacques.sinou}@ec-lyon.fr

² EADS Innovation Works
12 Rue Pasteur, BP76, F-92152 Suresnes Cedex
benoit.petitjean@eads.net

**Résumé** — La simulation dynamique d'un assemblage boulonné intègre des phénomènes de frottement complexes que l'on appréhende dans cet article par une méthode non linéaire de balance harmonique adaptative qui ajuste le nombre d'harmoniques retenus à chaque pas de fréquence. Cette méthode associée à une réduction sur les degrés de liberté d'interface utilise comme crête de sélection des harmoniques une variation de l'approximation de l'énergie de déformation en fonction du nombre d'harmoniques. L'application à un système de deux poutres boulonnées montre un nombre d'harmoniques grandissant au passage des résonances et met en évidence le caractère non linéaire de la liaison (amortissement et assouplissement des modes).

**Mots clés** — balance harmonique, jonction boulonnée, frottement, modèle non linéaire, dynamique des structures.

## 1 Introduction

La prédiction du comportement dynamique des structures spatiales assemblées par jonctions boulonnées ou rivetées est aujourd'hui un enjeu important et nécessite de prendre en compte des phénomènes dissipatifs complexes de micro-glissement [3]. L'amortissement de liaisons spatiales industrielles a déjà été étudiée lors de précédents travaux [1, 2] et a nécessité parfois une modélisation fine des liaisons. Pour réaliser une analyse dynamique, des modèles plus simples doivent être utilisés tout en conservant le caractère non linéaire des jonctions. Pour simuler de tels systèmes, la méthode de balance harmonique (HBM) développée dans les travaux de Pierre, Ferri et Dowell [5] permet de calculer des réponses périodiques de systèmes non linéaires. La méthode se base sur une décomposition en série de Fourier tronquée de la solution cherchée. Le nombre d'harmoniques nécessaires pour une précision donnée peut varier fortement sur une plage de fréquence étudiée. Cet article vise donc à développer une balance harmonique adaptative en ordre qui limite le nombre d'harmoniques retenus en calculant une énergie de déformation approchée et en observant sa saturation. Le cas d'application est un assemblage de deux poutres boulonnées utilisant un modèle de liaison intégrant un modèle de LuGre [3] inspiré des travaux de Song [7]. Une condensation sur les degrés de liberté (ddls) non-linéaires est réalisée. Dans un premier temps
la formulation de la méthode et du critère de sélection est détaillée. Enfin le modèle de jonction est
précisé et l'analyse des fonctions de réponse en fréquence (FRF) évalue la validité de la méthode.

2 Formulation de la HBM adaptative

2.1 Formulation générale et condensation

On considère un système mécanique discret à nddl ddls décrit par ses matrices de masse
**M**, de raideur **K** et d'amortissement **D**. Une force externe périodique F<sub>L</sub>(Ω, t) est appliquée sur
le système à la pulsation Ω. Les non-linéarités du système sont considérées comme des forces
extérieures F<sub>NL</sub>(X, ˙X, Ω, t) dépendantes des déplacements X, des vitesses ˙X, de la pulsation Ω et
du temps t. L'équation du mouvement peut ainsi s'écrire :

$$
\mathbf{M}\ddot{X} + \mathbf{D}\dot{X} + \mathbf{K}X = F(X, \dot{X}, \Omega, t) = F_L(\Omega, t) + F_{NL}(X, \dot{X}, \Omega, t) \quad (1)
$$

On suppose la réponse périodique et on la développe en série de Fourier tronquée :

$$
X(t) = B_0 + \sum_{k=1}^{m} \left( A_k \sin\left(\frac{k}{v}\Omega t\right) + B_k \cos\left(\frac{k}{v}\Omega t\right) \right) \quad (2)
$$

avec m le nombre d'harmoniques retenus, Z = [B₀ Aₖ Bₖ ...]<sup>T</sup> le vecteur des coefficients de Fourier de taille (2m + 1)nddl × 1 et v un entier permettant de représenter d'éventuels harmoniques non entiers. Si l'on note b(Z, Ω) la fonction qui associe les coefficients de Fourier de F(X, ˙X, Ω, t) = F<sub>L</sub>(Ω, t) + F<sub>NL</sub>(X, ˙X, Ω, t) à ceux de X(t), le problème revient à annuler une fonction H(Z) de ℝ<sup>(2m+1)×nddl</sup> → ℝ<sup>(2m+1)×nddl</sup> dépendant des coefficients du développement en série de X(t).

$$
H(Z) = \mathbf{A}(\Omega)Z - b(Z, \Omega) \quad \text{avec } \mathbf{A} = \text{diag} \left( \mathbf{K}, \dots, \left[ \begin{array}{cc} \mathbf{K} - \left(\frac{k}{v}\Omega\right)^2 \mathbf{M} & -\frac{k}{v}\Omega \mathbf{D} \\ \frac{k}{v}\Omega \mathbf{D} & \mathbf{K} - \left(\frac{k}{v}\Omega\right)^2 \mathbf{M} \end{array} \right], \dots \right) \quad (3)
$$

Si l'expression des coefficients des efforts d'excitation est généralement une donnée du problème,
il n'en est pas de même pour les efforts non linéaires. Dans le cas où il n'existe pas d'expression
analytique liant les coefficients de X à ceux de F<sub>NL</sub>, la méthode consiste à calculer dans le domaine
temporel les signaux approchés X(t) et ˙X(t) afin de calculer F<sub>NL</sub>(X, ˙X, Ω, t) sur une période. On
redécompose ensuite en série de Fourier le signal pour en déduire b(Z, Ω).

*Condensation* - Il est possible de réaliser une substitution dans les équations de façon à ne conser-
ver que les degrés de liberté dits non-linéaires, sur lesquels s'appliquent les efforts. Dans un pre-
mier temps, on réorganise les ddls par la transformation linéaire X = [P_p P_q] [X^p X^q]^T (avec X^p
et X^q les vecteurs des p et q ddls linéaires et non linéaires). P est une matrice de passage boo-
léennes, P_p et P_q étant ses p et q premières et dernières colonnes. Il s'ensuit une réorganisation
des coefficients de Fourier Z de X et b de F via la relation :

$$
\left[ \begin{array}{cc} {}^{T}\mathbf{N}_{\mathbf{P}_p} \mathbf{A} \mathbf{N}_{\mathbf{P}_p} & {}^{T}\mathbf{N}_{\mathbf{P}_p} \mathbf{A} \mathbf{N}_{\mathbf{P}_q} \\ {}^{T}\mathbf{N}_{\mathbf{P}_q} \mathbf{A} \mathbf{N}_{\mathbf{P}_p} & {}^{T}\mathbf{N}_{\mathbf{P}_q} \mathbf{A} \mathbf{N}_{\mathbf{P}_q} \end{array} \right] \left[ \begin{array}{c} Z_p \\ Z_q \end{array} \right] = \left[ \begin{array}{c} b_p \\ b_q \end{array} \right] \qquad (4)
$$

avec $\mathbf{N}_{\mathbf{P}_p} = \mathrm{diag}(\mathbf{P}_p, \mathbf{P}_p, \dots)_{(2m+1)nddl \times (2m+1)p}$ et $\mathbf{N}_{\mathbf{P}_q} = \mathrm{diag}(\mathbf{P}_q, \mathbf{P}_q, \dots)_{(2m+1)nddl \times (2m+1)q}$. Par substitution, on peut ainsi éliminer la dépendance en $Z_p$ et, si la décomposition est faite de telle sorte qu'aucun effort n'est présent sur les ddls linéaires ($b_p = 0$), on ramène le problème à l'annulation d'une fonction $H_q : \mathbb{R}^{(2m+1) \times q} \to \mathbb{R}^{(2m+1) \times q}$ :

$$
H_q(Z_q) = A_q(\Omega)Z_q - b_q(Z_q, \Omega) \tag{5}
$$

avec $\mathbf{A}_q(\Omega) = {}^T\mathbf{N}_{\mathbf{P}_q}\mathbf{A}\mathbf{N}_{\mathbf{P}_q} - {}^T\mathbf{N}_{\mathbf{P}_q}\mathbf{A}\mathbf{N}_{\mathbf{P}_p}({}^T\mathbf{N}_{\mathbf{P}_p}\mathbf{A}\mathbf{N}_{\mathbf{P}_p})^{-1}{}^T\mathbf{N}_{\mathbf{P}_p}\mathbf{A}\mathbf{N}_{\mathbf{P}_q}$.
## 2.2 Critère développé

Pour une fréquence d'excitation donnée, aucun outil théorique ne permet de déterminer quels sont les harmoniques prédominants dans la réponse. De plus, le nombre d'harmoniques nécessaires peut varier fortement sur l'intervalle de fréquence étudiée. Quelques outils numériques ont déjà été développés, notamment dans les travaux de Laxalde [4], qui se basent sur le degré d'approximation de la force non linéaire. Ce critère ne considère pas le comportement global du système du système et si les efforts non linéaires sont prédominants devant les efforts linéaires. Le critère présenté ici se base sur une saturation de l'énergie de déformation approchée.

*Energie de déformation approchée* - L'expression de l'énergie de déformation est obtenue en utilisant le développement en série de $X(t)$ et en calculant sa valeur moyenne sur une période. L'orthogonalité des fonctions trigonométriques mènent au résultat :

$$ \langle U \rangle = \langle \frac{1}{2}^T X(t) K X(t) \rangle = \frac{1}{2} Z L N_K Z \quad \text{avec} \quad L = \frac{1}{2} \begin{bmatrix} 2I & 0 \\ I & 0 \\ 0 & \ddots \end{bmatrix} \quad (6) $$

avec $I$ la matrice identité de taille $nddl \times nddl$.

*Energie de déformation approchée dans le cas d'une condensation* - Dans le cas d'une condensation, l'expression de l'énergie approchée doit être adaptée pour être calculée à partir des coefficients de Fourier des ddls non linéaires $Z_q$ uniquement, évitant ainsi la reconstitution des coefficients des termes linéaires $Z_p$ à chaque évaluation du critère. Pour ce faire, on utilise la relation liant $Z_p$ à $Z_q$ par substitution dans l'équation (5) :

$$ Z_p = ({}^T N_{P_p} A N_{P_p})^{-1} (b_p - {}^T N_{P_p} A N_{P_p} Z_q) \quad (7) $$

On suppose que les coefficients linéaires et non linéaires sont liés dans un cas statique ($\Omega = 0$ dans l'équation précédente) et que le découpage est tel que $b_p = 0$. En utilisant l'équation (5) dans l'équation (6), on obtient :

$$ \langle U \rangle = \frac{1}{2} Z_q L_q N_{K_{red}} Z_q \quad \text{avec} \quad K_{red} = {}^T P_q K P_q - {}^T P_q K P_p ({}^T P_p K P_p)^{-1} {}^T P_p K P_q \quad (8) $$

*Critère* - Le critère développé est calculé pour une valeur de fréquence donnée et commence par réaliser un calcul avec un harmonique. Ensuite un calcul avec un harmonique supplémentaire (1,3,5 si les harmoniques répondant sont impairs) est réalisé et la différence relative entre les deux valeurs approchées de l'énergie donne le critère. L'algorithme est stoppé lorsque la variation relative est inférieure à la valeur seuil déterminée par l'utilisateur. Enfin, les matrices $L$ et $L_q$ étant des matrices diagonales par blocs et constantes, observer la saturation de $\langle U \rangle$ est équivalent à observer la saturation d'une grandeur $\langle \hat{U} \rangle$ sans le terme $1/2$ et ces deux matrices.

# 3 Simulation et Résultats

## 3.1 Description du modèle

Un assemblage de deux poutres encastrées libres reliées entre elles par une jonction boulonnée est simulé. Les poutres en Al 7075 et de section $5.1 \text{cm} \times 2.5 \text{cm}$ mesurent $34.7 \text{cm}$ et $84.7 \text{cm}$. Les poutres sont modélisées par éléments finis sous Abaqus (formulation de Timoshenko) avec 10 éléments et 25 éléments. L'élément représentant la jonction boulonnée associe une matrice de
masse d'un élément poutre et une matrice de raideur non linéaire considérée comme un effort ex-
térieur. Un amortissement de Rayleigh est utilisé en utilisant les matrices de masse et de raideur
d'une poutre monobloc. L'amortissement modal est ainsi de 0.1% sur les deuxième (0.24 kHz) et
cinquième (1.14 kHz) modes. Par analogie avec les travaux de Song [7], la raideur non linéaire de

<figure><img src="image_2.png" /><figcaption>Figure 1 – Elément ALBE</figcaption></figure>

la jonction utilise un modèle noté Adjusted LuGre Beam Element (ALBE) présenté en Figure 1.
Il se base sur un remplacement des raideurs d'un élément poutre linéaire par une association d'un
modèle de LuGre et d'une raideur résiduelle $k_{a,i}$, $i \in 1,2$, caractéristique des jonctions boulon-
nées. L'élément possède 2 ddls de rotation $R1$ et $R2$ et deux ddls de translation $T1$ et $T2$. $h$ et $L$
sont respectivement la hauteur de section et la longueur de l'élément.

Pour exprimer la relation entre l'effort $F_{NL,ALBE}$ non linéaire et les ddls de l'élément, il faut considérer les allongements $\Delta_1 = \frac{L}{2}(R_1 + R_2) + (T1 + T1)$ et $\Delta_2 = \frac{h}{2}(R_1 - R_2)$ de chaque ressort. L'effort de chaque modèle de LuGre $f_{LuGre,i}$, $i \in 1,2$ s'écrit donc en fonction de l'allongement $\Delta_i$ mais aussi en fonction de la valeur de sa variable interne $\zeta_i$ et de sa dérivée $\dot{\zeta}_i$. Le modèle s'écrit :

$$f_{LuGre,i}(\Delta_i, \dot{\Delta}_i, \zeta_i, \dot{\zeta}_i) = \sigma_{0i}\Delta_i + \sigma_{1i}\dot{\zeta}_i + \alpha_{2i}\dot{\Delta}_i \quad (9)$$

$$\dot{\zeta}_i = \dot{\Delta}_i - \frac{\sigma_{0i}}{\alpha_{0i} + \alpha_{1i}e^{-(\frac{\Delta_i}{v_{0i}})^2}} |\dot{\Delta}_i| \zeta_i \quad (10)$$

L'effort de chaque couple (LuGre,raideur) est donc :

$$f_i(\Delta_i, \dot{\Delta}_i, \zeta_i, \dot{\zeta}_i) = f_{LuGre,i}(\Delta_i, \dot{\Delta}_i, \zeta_i, \dot{\zeta}_i) + k_{a,i}\Delta_i \quad (11)$$

Pour représenter la diminution de raideur observée lors des phases de microslip, un coefficient $\gamma_i \in [0;1]$ permet de lier la raideur $k_{a,i}$, le paramètre $\sigma_{0i}$ et la raideur de l'élément $k_i$ correspondant à un élément linéaire, par les relations $\sigma_{0i} = (1 - \gamma_i)k_i$ et $k_{a,i} = \gamma_i k_i$. Les efforts $F_1, F_2$ et les moments résultant $M_1, M_2$ s'écrivent donc :

$$F_{NL,ALBE} \begin{bmatrix} \Delta_1 \\ \Delta_2 \\ \dot{\Delta}_1 \\ \dot{\Delta}_2 \end{bmatrix}, \begin{bmatrix} \zeta_1 \\ \zeta_2 \\ \dot{\zeta}_1 \\ \dot{\zeta}_2 \end{bmatrix}) = \begin{bmatrix} F1 \\ M1 \\ F1 \\ M2 \end{bmatrix} = \begin{bmatrix} f_1(\Delta_1, \dot{\Delta}_1, \zeta_1, \dot{\zeta}_1) \\ \frac{L}{2}f_1(\Delta_1, \dot{\Delta}_1, \zeta_1, \dot{\zeta}_1) + \frac{h}{2}f_2(\Delta_2, \dot{\Delta}_2, \zeta_2, \dot{\zeta}_2) \\ -f_1(\Delta_1, \dot{\Delta}_1, \zeta_1, \dot{\zeta}_1) \\ \frac{L}{2}f_1(\Delta_1, \dot{\Delta}_1, \zeta_1, \dot{\zeta}_1) - \frac{h}{2}f_2(\Delta_2, \dot{\Delta}_2, \zeta_2, \dot{\zeta}_2) \end{bmatrix} \quad (12)$$

A partir des expressions des allongements $\Delta_i$, on obtient donc l'expression des efforts $F_{NL,ALBE}$ en fonction des ddls de l'élément $T_1, R_1, T_2, R_2$ et des variables internes $\zeta_1, \zeta_2$.

Pour les valeurs numériques, les deux raideurs $k_i$ correspondant à l'élément linéaire sont obtenues via $k_1 = 12EI/L^3 = 1,43.10^6$ N/mm et $k_2 = 4EI/Lh^2 = 8,92.10^5$ N/mm. Les autres paramètres sont déduits par analogie des travaux de [6], à savoir $\gamma_1 = \gamma_2 = 0.1078, \alpha_{01} = \alpha_{02} = 81.9$ N, $\sigma_{11} = \sigma_{12} = \alpha_{11} = \alpha_{12} = \alpha_{21} = \alpha_{22} = 0$. Finalement $\sigma_{01} = 1,27.10^6$ N/mm, $\sigma_{02} = 7,96.10^5$ N/mm, $k_{a,1} = 1,55.10^5$ N/mm, $k_{a,2} = 9.62.10^4$ N/mm.

Le système est ensuite excité sur une plage de fréquence [0 – 2,3] kHz par un effort harmonique $F_L$ d'amplitude $42N$ et de pulsation $\Omega$ appliqué sur la plus grande poutre. En reprenant la
formulation définie dans l'équation (1), on considère les efforts non linéaires $F_{NL}$ comme étant l'opposé de ceux provenant de la jonction $F_{NL,ALBE}$.

A cette équation, on doit ajouter les deux équations (10) liant les variables internes que l'on doit réadapter à la formulation de la balance harmonique. En notant $b_{\zeta_j}(Z, Z_{\zeta_j}, \Omega)$ les coefficients de Fourier du terme non linéaire correspondant au terme de droite de l'équation (10), le problème revient donc à annuler une fonction $HH(Z, Z_{\zeta})$ dépendant des coefficients du développement en série respectivement de $X$ et de $\zeta = [\zeta_1 \ \zeta_2]$. Ces coefficients sont notés respectivement $Z$ et $Z_{\zeta} = [Z_{\zeta_1} \ Z_{\zeta_2}]$. $HH$ est une fonction de $\mathbb{R}^{(2m+1) \times (nddl+2)} \to \mathbb{R}^{(2m+1) \times (nddl+2)}$.

$$HH(Z, Z_{\zeta}) = \begin{cases} H(Z, Z_{\zeta}) = \mathbf{A}(\Omega)Z - b(Z, Z_{\zeta}, \Omega) \\ C(Z, Z_{\zeta}) = \mathbf{N}_{I_{2 \times 2}} \nabla Z_{\zeta} - b_{\zeta}(Z, Z_{\zeta}, \Omega) \end{cases} \quad \text{avec} \quad b_{\zeta} = \begin{bmatrix} b_{\zeta_1} \\ b_{\zeta_2} \end{bmatrix} \quad (13)$$

Dans le cas d'une condensation on définit de façon analogue une fonction $HH_q(Z_q, Z_{\zeta})$ de $\mathbb{R}^{(2m+1) \times (q+2)}$ vers $\mathbb{R}^{(2m+1) \times (q+2)}$ qui associe les fonctions $H_q(Z_q, Z_{\zeta}) = \mathbf{A}_q(\Omega)Z_q - b_q(Z_q, Z_{\zeta}, \Omega)$ et $C(Z_q, Z_{\zeta}) = \mathbf{N}_{I_{2 \times 2}} \nabla Z_{\zeta} - b_{q,\zeta}(Z_q, Z_{\zeta}, \Omega)$. Enfin le calcul de l'énergie approchée de déformation $\langle \hat{U} \rangle$ reste identique à l'équation (6) et (8). Seuls les coefficients de Fourier des ddls physiques $Z$ et $Z_q$ sont utilisés.

<figure><img src="image_3.png" /><figcaption>Figure 2 – Réponses en fréquence.</figcaption></figure>

## 3.2 Résultats

Les différentes réponses en fréquence (FRF) obtenues pour le cas linéaire, les cas non-linéaires à 1 et 5 harmoniques et le cas non-linéaire pour la HBM adaptative sont présentées en Figure 2 pour un des ddls de translation du modèle ALBE. La variation relative de l'énergie de déformation approchée a été fixée à 3%. Tout d'abord l'impact de la non-linéarité est très marqué si l'on observe la réponse pour 1 harmonique et se traduit par deux effets : une diminution de l'amplitude des pics de résonance caractéristique de l'amortissement générée par la liaison et un assouplissement des modes traduisant une baisse de raideur dans la jonction. Certains modes ont ainsi des décalages fréquentiels importants voire des distortions, notamment les deuxième, cinquième et sixième pics. En outre un calcul à l'ordre 1 semble satisfaisant pour approcher le comportement global car les réponses pour différents ordres semblent très proches. En ce qui concerne la réponse pour la HBM adaptative on observe sur le graphe gauche de la Figure 3 qu'au passage des résonances, le nombre d'harmoniques varie de 3 à 15 marquant le caractère non linéaire du système lors des pics. L'évolution de l'énergie de déformation approchée $\langle \hat{U} \rangle$ en fonction du nombre d'harmoniques retenus sur le graphe droit de la Figure 3 confirme également cette analyse. Un zoom sur le pic proche de 1kHz a été effectué pour plus de clarté.
<figure><img src="image_4.png" /><figcaption>Figure 3 – Evolution du nombre d’harmoniques dans le cas adaptatif.</figcaption></figure>

## 4 Conclusion

Nous avons présenté dans cet article une méthode de balance harmonique adaptative appliquée sur un modèle de poutres reliées par une jonction boulonnée. Le critère développé se base sur une variation relative d'une approximation de l'énergie de déformation du système. L'analyse des résultats montre qu'un seul harmonique suffit à approcher correctement le mouvement en dehors des résonances. Enfin la HBM adaptative met en évidence la nécessité de devoir considérer plusieurs harmoniques lors des résonances car les non linéarités du système y jouent un rôle plus important.

## Références

[1] A. Caignot, P. Ladevèze, D. Néron, V. Le Gallo, and L. Gonidou. Virtual testing for the prediction of damping in joints. ICED 2007 International Conference in Engineering Dynamics, 2007.

[2] A. D. Crocombe, R. Wang, G. Richardson, and C. I. Underwood. Estimating the energy dissipated in a bolted spacecraft at resonance. *Computers & Structures*, 84(5-6):340-350, 1 2006.

[3] L. Gaul and R. Nitsche. The role of friction in mechanical joints. *Applied Mechanics Reviews*, 54(2):93-106, 2001.

[4] D. Laxalde. *Etude d'amortisseurs non-linéaires appliqués aux roues aubagées et aux systèmes multi-étages*. PhD thesis, Ecole Centrale de Lyon, 2007.

[5] C. Pierre, A. A. Ferri, and E. H. Dowell. Multi-harmonic analysis of dry friction damped systems using an incremental harmonic-balance method. *Journal of Applied Mechanics-Transactions of the Asme*, 52(4):958-964, 1985.

[6] O. V. Shiryayev, S. M. Page, C. L. Pettit, and J. C. Slater. Parameter estimation and investigation of a bolted joint model. *Journal of Sound and Vibration*, 307(3-5):680-697, 2007.

[7] Y. Song, C. J. Hartwigsen, D. M. McFarland, A. F. Vakakis, and L. A. Bergman. Simulation of dynamics of beam structures with bolted joints using adjusted iwan beam elements. *Journal of Sound and Vibration*, 273(1-2):249-276, 2004.