# Contrôle de courbes de réponses isolées par optimisation structurelle

Adrien Mélot, Enora Denimal Goy, Ludovic Renson

► To cite this version:

Adrien Mélot, Enora Denimal Goy, Ludovic Renson. Contrôle de courbes de réponses isolées par optimisation structurelle. CSMA 2024 - 16ème Colloque National en Calcul de Structures, CNRS, CSMA, ENS Paris-Saclay, CentraleSupélec, May 2024, Giens, France. http://hal-04610912

HAL Id: hal-04610912

https://hal.science/hal-04610912

Submitted on 19 Jun 2024

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 Attribution 4.0 International License
# Contrôle de courbes de réponses isolées par optimisation structurelle

Adrien Mélot, Enora Denimal Goy, Ludovic Renson

► To cite this version:

Adrien Mélot, Enora Denimal Goy, Ludovic Renson. Contrôle de courbes de réponses isolées par optimisation structurelle. 16ème Colloque National en Calcul de Structures (CSMA 2024), CNRS; CSMA; ENS Paris-Saclay; CentraleSupélec, May 2024, Hyères, France. hal-04610912

HAL Id: hal-04610912

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

Submitted on 19 Jun 2024

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 Attribution 4.0 International License
# Contrôle de courbes de réponses isolées par optimisation structurelle

A. Mélot<fnref n="1" />, E. Denimal Goy<fnref n="2" />, L. Renson<fnref n="3" />

**Résumé** — Nous introduisons un cadre d'analyse numérique permettant de contrôler des courbes de réponse isolées (isolas), c'est-à-dire formant des courbes fermées et non connectées à la branche principale de solutions. La méthodologie repose sur des analyses de suivi de bifurcation afin de suivre l'évolution de bifurcations nœud-col dans un espace de codimension-2. La théorie des singularités est mise en œuvre afin de distinguer les points de formation et de fusion d'isolas des bifurcations de codimension-2. Un problème d'optimisation est défini pour avancer ou retarder la formation ou la fusion d'isolas.

**Mots clés** — Vibrations non-linéaires, Suivi de bifurcations, Singularités, Isolas, Equilibrage harmonique

## 1 Introduction

Les courbes de réponses isolées, ou isolas, sont des branches de solutions de systèmes dynamiques non-linéaires qui forment des courbes fermées dans l'espace des paramètres. Bien qu'elles ne soient pas connectées à la branche principale de solutions, c'est-à-dire la branche liée aux solutions (quasi-)linéaires aux basses énergies, elles peuvent fusionner avec cette dernière lorsqu'un paramètre, par exemple l'amplitude de forçage, est modifié (Fig. 1). Il est également possible d'observer la fusion de plusieurs isolas [7].

<figure><img src="image_2.png" /><figcaption>FIGURE 1 — Exemple d'isola d'une poutre encastrée libre. Les courbes en bleu sont des courbes de réponse forcée pour différentes valeurs d'amplitude de forçage. La courbe rouge est la courbe de bifurcations nœud-col.</figcaption></figure>

De nombreuses études attestent de la présence d'isolas dans la réponse des systèmes non-linéaires, les premières datant des années 1950 [1]. Les isolas sont particulièrement dangereux car ils sont susceptibles d'induire des régions multi-stables dans lesquelles des réponses de grande amplitudes coexistent avec des réponses de faible amplitudes. Il peuvent également conduire à des augmentations significatives des décalages en amplitude et fréquence lorsque la fusion entre un isola et la branche principale survient au niveau d'une résonance non-linéaire. Certains travaux récents ont montré qu'il était également possible de les exploiter. Par exemple, les isolas induits par des brisures de symétrie ont été utilisés à des fins de détection de masses infinitésimales [6]. Bien que leur calcul soit d'un intérêt certain, les isolas ne sont

<footnote n="1">¹ Univ. Gustave Eiffel, Inria, COSYS-SII, I4S Team, Campus Beaulieu, 35042 Rennes, France, adrien.melot@inria.fr</footnote>

<footnote n="2">² Inria, CMAP, Ecole Polytechnique, IPP, Palaiseau, France, enora.denimal@inria.fr</footnote>

<footnote n="3">³ Dynamics Group, Department of Mechanical Engineering, Imperial College London, London, United Kingdom, l.renson@imperial.ac.uk</footnote>
# Contrôle de courbes de réponses isolées par optimisation structurelle

A. Mélot<fnref n="1" />, E. Denimal Goy<fnref n="2" />, L. Renson<fnref n="3" />

**Résumé** — Nous introduisons un cadre d'analyse numérique permettant de contrôler des courbes de réponse isolées (isolas), c'est-à-dire formant des courbes fermées et non connectées à la branche principale de solutions. La méthodologie repose sur des analyses de suivi de bifurcation afin de suivre l'évolution de bifurcations nœud-col dans un espace de codimension-2. La théorie des singularités est mise en œuvre afin de distinguer les points de formation et de fusion d'isolas des bifurcations de codimension-2. Un problème d'optimisation est défini pour avancer ou retarder la formation ou la fusion d'isolas.

**Mots clés** — Vibrations non-linéaires, Suivi de bifurcations, Singularités, Isolas, Equilibrage harmonique

## 1 Introduction

Les courbes de réponses isolées, ou isolas, sont des branches de solutions de systèmes dynamiques non-linéaires qui forment des courbes fermées dans l'espace des paramètres. Bien qu'elles ne soient pas connectées à la branche principale de solutions, c'est-à-dire la branche liée aux solutions (quasi-)linéaires aux basses énergies, elles peuvent fusionner avec cette dernière lorsqu'un paramètre, par exemple l'amplitude de forçage, est modifié (Fig. 1). Il est également possible d'observer la fusion de plusieurs isolas [7].

<figure><img src="image_2.png" /><figcaption>FIGURE 1 — Exemple d'isola d'une poutre encastrée libre. Les courbes en bleu sont des courbes de réponse forcée pour différentes valeurs d'amplitude de forçage. La courbe rouge est la courbe de bifurcations nœud-col.</figcaption></figure>

De nombreuses études attestent de la présence d'isolas dans la réponse des systèmes non-linéaires, les premières datant des années 1950 [1]. Les isolas sont particulièrement dangereux car ils sont susceptibles d'induire des régions multi-stables dans lesquelles des réponses de grande amplitudes coexistent avec des réponses de faible amplitudes. Il peuvent également conduire à des augmentations significatives des décalages en amplitude et fréquence lorsque la fusion entre un isola et la branche principale survient au niveau d'une résonance non-linéaire. Certains travaux récents ont montré qu'il était également possible de les exploiter. Par exemple, les isolas induits par des brisures de symétrie ont été utilisés à des fins de détection de masses infinitésimales [6]. Bien que leur calcul soit d'un intérêt certain, les isolas ne sont

<footnote n="1">¹ Univ. Gustave Eiffel, Inria, COSYS-SII, I4S Team, Campus Beaulieu, 35042 Rennes, France, adrien.melot@inria.fr</footnote>

<footnote n="2">² Inria, CMAP, Ecole Polytechnique, IPP, Palaiseau, France, enora.denimal@inria.fr</footnote>

<footnote n="3">³ Dynamics Group, Department of Mechanical Engineering, Imperial College London, London, United Kingdom, l.renson@imperial.ac.uk</footnote>
souvent pas détectés car les analyses de continuations employées classiquement sont initialisées à basse énergie sur la branche principale de solutions.

Un certain nombre d'auteurs ont proposé des méthodes analytiques afin de prédire les isolas [8, 12]. D'un point de vue numérique, la méthode la plus populaire est le suivi de bifurcations, qui consiste à suivre l'évolution de la position d'un point de bifurcation dans un espace de codimension-2. Les isolas étant des courbes fermées, ils possèdent au moins deux bifurcations de type nœud-col. Il est donc possible de suivre une bifurcation nœud-col en fonction de la fréquence d'excitation et d'un paramètre additionnel, nommé paramètre de suivi dans ce qui suit, et de détecter un isola lorsque la courbe de bifurcations possède des solutions multiples pour certaines valeurs du paramètre de suivi [4, 2] (Fig. 1). Cependant, bien que les concepts de base des outils numériques soient relativement bien établis, la détection et la classification des points de retournement des courbes de suivi de bifurcations afin de détecter les points de formation et de fusion d'isolas nécessite toujours une interprétation a posteriori incompatible avec une procédure d'optimisation.

Nous proposons ici une méthodologie numérique permettant d'optimiser la conception de systèmes mécaniques afin de contrôler la position des isolas dans l'espace des paramètres. Nous introduisons un problème d'optimisation reposant sur une série d'analyses de continuation en codimension-2 afin de calculer les courbes de bifurcations nœud-col. Les points de retournement de ces courbes sont détectés et les points de formation et de fusion d'isolas, correspondant à des singularités de codimension-1, sont identifiés de l'ensemble des bifurcations de codimension-2 grâce à la théorie des singularités.

## 2 Formulation du problème d'optimisation

Dans ce qui suit, nous considérons des problèmes de la forme générale suivante :

$$ \begin{array}{ll} \underset{\mathbf{x}}{\text{maximiser}} & \mathcal{F} \\ \text{soumise à} & b_i^l \le x_i \le b_i^u \quad \forall i \in [1, p], \end{array} \qquad (1) $$

où $\mathcal{F}$ est la fonction objectif considérée, $\mathbf{x} \in \mathbb{R}^p$ est le vecteur des paramètres d'optimisation, c'est-à-dire des paramètres de conception. Nous considérons ici uniquement des contraintes simples destinées à borner les valeurs des paramètres d'optimisation avec les vecteurs $\mathbf{b}^l$ and $\mathbf{b}^u$ afin de garantir que les valeurs sont réalistes, par exemple avec une masse strictement positive. Néanmoins, sous réserve du choix d'un algorithme d'optimisation approprié, la mise en œuvre de contraintes plus complexes est relativement aisée.

### 2.1 Définition de la fonction objectif

Plusieurs objectifs peuvent être considérés concernant la position d'isolas dans l'espace des paramètres. Nous considérons ici la position au sens de la valeur d'un paramètre de fonctionnement $\mu$ dont la variation peut conduire à la formation ou à la fusion d'un isola. Un exemple classique d'un tel paramètre est l'amplitude du forçage appliqué au système. Notons $\mu_{\text{birth}}$ (resp. $\mu_{\text{merger}}$) les valeurs de paramètres auxquelles un isola est créé (resp. fusionne) avec une autre branche de solution. Ces valeurs sont calculées par la procédure de suivi de bifurcations présentée dans la section suivante.

**Problème 1 : retarder la formation (resp. fusion) d'un isola.** L'objectif le plus simple consiste à retarder l'apparition d'un point de formation d'isola afin d'augmenter la taille de la gamme du paramètre de suivi de bifurcation sans multi-stabilité. Si nous supposons que la formation et la fusion d'isola a lieu pour des valeurs croissantes du paramètre $\mu$, la fonction objectif s'écrit :

$$ \underset{\mathbf{x}}{\text{maximiser}} \quad \mathcal{F} = \mu_{\text{birth}} \quad (\mathcal{F} = \mu_{\text{merger}}) \qquad (2) $$

**Problème 2 : déplacer la formation d'isola à une valeur cible du paramètre de suivi.** Une autre possibilité consiste à déplacer le point de formation d'isola vers valeur $\mu_{\text{target}}$ du paramètre de suivi de bifurcation. Dans ce cas, la fonction objectif s'écrit :

$$ \underset{\mathbf{x}}{\text{minimiser}} \quad \mathcal{F} = \left| \frac{\mu_{\text{birth}} - \mu_{\text{target}}}{\mu_{\text{target}}} \right| \qquad (3) $$
souvent pas détectés car les analyses de continuations employées classiquement sont initialisées à basse énergie sur la branche principale de solutions.

Un certain nombre d'auteurs ont proposé des méthodes analytiques afin de prédire les isolas [8, 12]. D'un point de vue numérique, la méthode la plus populaire est le suivi de bifurcations, qui consiste à suivre l'évolution de la position d'un point de bifurcation dans un espace de codimension-2. Les isolas étant des courbes fermées, ils possèdent au moins deux bifurcations de type nœud-col. Il est donc possible de suivre une bifurcation nœud-col en fonction de la fréquence d'excitation et d'un paramètre additionnel, nommé paramètre de suivi dans ce qui suit, et de détecter un isola lorsque la courbe de bifurcations possède des solutions multiples pour certaines valeurs du paramètre de suivi [4, 2] (Fig. 1). Cependant, bien que les concepts de base des outils numériques soient relativement bien établis, la détection et la classification des points de retournement des courbes de suivi de bifurcations afin de détecter les points de formation et de fusion d'isolas nécessite toujours une interprétation a posteriori incompatible avec une procédure d'optimisation.

Nous proposons ici une méthodologie numérique permettant d'optimiser la conception de systèmes mécaniques afin de contrôler la position des isolas dans l'espace des paramètres. Nous introduisons un problème d'optimisation reposant sur une série d'analyses de continuation en codimension-2 afin de calculer les courbes de bifurcations nœud-col. Les points de retournement de ces courbes sont détectés et les points de formation et de fusion d'isolas, correspondant à des singularités de codimension-1, sont identifiés de l'ensemble des bifurcations de codimension-2 grâce à la théorie des singularités.

## 2 Formulation du problème d'optimisation

Dans ce qui suit, nous considérons des problèmes de la forme générale suivante :

$$ \begin{array}{ll} \underset{\mathbf{x}}{\text{maximiser}} & \mathcal{F} \\ \text{soumise à} & b_i^l \le x_i \le b_i^u \quad \forall i \in [1, p], \end{array} \qquad (1) $$

où $\mathcal{F}$ est la fonction objectif considérée, $\mathbf{x} \in \mathbb{R}^p$ est le vecteur des paramètres d'optimisation, c'est-à-dire des paramètres de conception. Nous considérons ici uniquement des contraintes simples destinées à borner les valeurs des paramètres d'optimisation avec les vecteurs $\mathbf{b}^l$ and $\mathbf{b}^u$ afin de garantir que les valeurs sont réalistes, par exemple avec une masse strictement positive. Néanmoins, sous réserve du choix d'un algorithme d'optimisation approprié, la mise en œuvre de contraintes plus complexes est relativement aisée.

### 2.1 Définition de la fonction objectif

Plusieurs objectifs peuvent être considérés concernant la position d'isolas dans l'espace des paramètres. Nous considérons ici la position au sens de la valeur d'un paramètre de fonctionnement $\mu$ dont la variation peut conduire à la formation ou à la fusion d'un isola. Un exemple classique d'un tel paramètre est l'amplitude du forçage appliqué au système. Notons $\mu_{\text{birth}}$ (resp. $\mu_{\text{merger}}$) les valeurs de paramètres auxquelles un isola est créé (resp. fusionne) avec une autre branche de solution. Ces valeurs sont calculées par la procédure de suivi de bifurcations présentée dans la section suivante.

**Problème 1 : retarder la formation (resp. fusion) d'un isola.** L'objectif le plus simple consiste à retarder l'apparition d'un point de formation d'isola afin d'augmenter la taille de la gamme du paramètre de suivi de bifurcation sans multi-stabilité. Si nous supposons que la formation et la fusion d'isola a lieu pour des valeurs croissantes du paramètre $\mu$, la fonction objectif s'écrit :

$$ \underset{\mathbf{x}}{\text{maximiser}} \quad \mathcal{F} = \mu_{\text{birth}} \quad (\mathcal{F} = \mu_{\text{merger}}) \qquad (2) $$

**Problème 2 : déplacer la formation d'isola à une valeur cible du paramètre de suivi.** Une autre possibilité consiste à déplacer le point de formation d'isola vers valeur $\mu_{\text{target}}$ du paramètre de suivi de bifurcation. Dans ce cas, la fonction objectif s'écrit :

$$ \underset{\mathbf{x}}{\text{minimiser}} \quad \mathcal{F} = \left| \frac{\mu_{\text{birth}} - \mu_{\text{target}}}{\mu_{\text{target}}} \right| \qquad (3) $$
Notons que d'autres fonctions objectifs sont envisageables. Par exemple, il est possible de considérer plusieurs objectifs en retardant la formation d'isola tout en imposant son amplitude et/ou sa fréquence. De tels cas ne sont pas considérés ici mais les informations nécessaires à l'évaluation de telles fonctions objectifs peuvent être facilement obtenues avec les sorties de l'analyse de suivi de bifurcation.

# 3 Méthodes de résolution

## 3.1 Continuation des bifurcations nœud-col

Dans le cas général, les équations du mouvement d'un modèle à n degrés de liberté (ddl) s'écrivent sous la forme matricielle suivante :

$$ M\ddot{q}(t) + C\dot{q}(t) + Kq(t) + f_{nl}(q(t)) = f_{ex}(t) \quad (4) $$

où M, C, K sont respectivement les matrices de masse, d'amortissement et de raideur avec q le vecteur déplacement de chaque degré de liberté. $f_{nl}$ et $f_{ex}$ sont les vecteurs d'efforts non-linéaires et d'efforts extérieurs périodiques.

Sans perte de généralités, nous utilisons ici la méthode de l'équilibrage harmonique (HBM) [14, 10] afin de discrétiser les orbites périodiques du problème dynamique. La méthode consiste à exprimer les solutions périodiques q de fréquence fondamentale $\Omega$ sous la forme de séries de Fourier tronquées à un ordre H suffisamment élevé pour garantir une précision suffisante des résultats :

$$ q = \Re \left( \sum_{k=0}^{\infty} \tilde{q}_k e^{ik\Omega t} \right) \approx \Re \left( \sum_{k=0}^{H} \tilde{q}_k e^{ik\Omega t} \right) \quad (5) $$

où le vecteur $\tilde{q}$ contient les coefficients complexes de la transformée de Fourier. Les déplacements, vitesses et accélérations s'expriment alors :

$$ q = \Re([T \otimes I_n] \tilde{q}), \quad (6) $$

$$ \dot{q} = \Re([T\Omega \nabla \otimes I_n] \tilde{q}), \quad (7) $$

$$ \ddot{q} = \Re([T\Omega^2 \nabla^2 \otimes I_n] \tilde{q}), \quad (8) $$

avec $\otimes$ le produit de Kronecker, $I_n$ la matrice identité de taille $n \times n$, T le vecteur des fonctions de base exponentielles et $\nabla$ l'opérateur de dérivation fréquentielle complexe. Les efforts non-linéaires et le forçage périodique extérieur peuvent s'écrire de façon similaire :

$$ f_{nl} = \Re \left( \sum_{k=0}^{H} \tilde{f}_{nl,k} e^{ik\Omega t} \right), \quad (9) $$

$$ f_{ex} = \Re \left( \sum_{k=0}^{H} \tilde{f}_{ex,k} e^{ik\Omega t} \right). \quad (10) $$

où $\tilde{f}_{nl}$ et $\tilde{f}_{ex}$ contiennent les coefficients de Fourier des efforts non-linéaires et du forçage extérieur, respectivement.

En substituant ces expressions dans l'équation (4) et en réalisant une projection de Galerkin, nous obtenons un système de $n(H+1)$ équations algébriques non-linéaires :

$$ R(\tilde{q}, \Omega) = Z(\Omega)\tilde{q} + \tilde{f}_{nl}(\tilde{q}) - \tilde{f}_{ex} = 0, \quad (11) $$

avec $Z(\Omega)$ la matrice de rigidité dynamique dupliquée sur chaque harmonique considéré :

$$ Z(\Omega) = \Omega^2 \nabla^2 \otimes M + \Omega \nabla \otimes C + I_{H+1} \otimes K. \quad (12) $$

Afin de réaliser le suivi de points de bifurcation et de dépasser les points de retournements une méthode de continuation est mise en œuvre. Pour se faire, deux paramètres sont traités comme des inconnues. Il s'agit ici de la fréquence $\Omega$ et du paramètre $\mu$. Deux équations additionnelles sont ajoutées
Notons que d'autres fonctions objectifs sont envisageables. Par exemple, il est possible de considérer plusieurs objectifs en retardant la formation d'isola tout en imposant son amplitude et/ou sa fréquence. De tels cas ne sont pas considérés ici mais les informations nécessaires à l'évaluation de telles fonctions objectifs peuvent être facilement obtenues avec les sorties de l'analyse de suivi de bifurcation.

# 3 Méthodes de résolution

## 3.1 Continuation des bifurcations nœud-col

Dans le cas général, les équations du mouvement d'un modèle à n degrés de liberté (ddl) s'écrivent sous la forme matricielle suivante :

$$ M\ddot{q}(t) + C\dot{q}(t) + Kq(t) + f_{nl}(q(t)) = f_{ex}(t) \quad (4) $$

où M, C, K sont respectivement les matrices de masse, d'amortissement et de raideur avec q le vecteur déplacement de chaque degré de liberté. $f_{nl}$ et $f_{ex}$ sont les vecteurs d'efforts non-linéaires et d'efforts extérieurs périodiques.

Sans perte de généralités, nous utilisons ici la méthode de l'équilibrage harmonique (HBM) [14, 10] afin de discrétiser les orbites périodiques du problème dynamique. La méthode consiste à exprimer les solutions périodiques q de fréquence fondamentale $\Omega$ sous la forme de séries de Fourier tronquées à un ordre H suffisamment élevé pour garantir une précision suffisante des résultats :

$$ q = \Re \left( \sum_{k=0}^{\infty} \tilde{q}_k e^{ik\Omega t} \right) \approx \Re \left( \sum_{k=0}^{H} \tilde{q}_k e^{ik\Omega t} \right) \quad (5) $$

où le vecteur $\tilde{q}$ contient les coefficients complexes de la transformée de Fourier. Les déplacements, vitesses et accélérations s'expriment alors :

$$ q = \Re([T \otimes I_n] \tilde{q}), \quad (6) $$

$$ \dot{q} = \Re([T \Omega \nabla \otimes I_n] \tilde{q}), \quad (7) $$

$$ \ddot{q} = \Re([T \Omega^2 \nabla^2 \otimes I_n] \tilde{q}), \quad (8) $$

avec $\otimes$ le produit de Kronecker, $I_n$ la matrice identité de taille $n \times n$, $T$ le vecteur des fonctions de base exponentielles et $\nabla$ l'opérateur de dérivation fréquentielle complexe. Les efforts non-linéaires et le forçage périodique extérieur peuvent s'écrire de façon similaire :

$$ f_{nl} = \Re \left( \sum_{k=0}^{H} \tilde{f}_{nl,k} e^{ik\Omega t} \right), \quad (9) $$

$$ f_{ex} = \Re \left( \sum_{k=0}^{H} \tilde{f}_{ex,k} e^{ik\Omega t} \right). \quad (10) $$

où $\tilde{f}_{nl}$ et $\tilde{f}_{ex}$ contiennent les coefficients de Fourier des efforts non-linéaires et du forçage extérieur, respectivement.

En substituant ces expressions dans l'équation (4) et en réalisant une projection de Galerkin, nous obtenons un système de $n(H+1)$ équations algébriques non-linéaires :

$$ \mathbf{R}(\tilde{q}, \Omega) = \mathbf{Z}(\Omega)\tilde{q} + \tilde{f}_{nl}(\tilde{q}) - \tilde{f}_{ex} = \mathbf{0}, \quad (11) $$

avec $\mathbf{Z}(\Omega)$ la matrice de rigidité dynamique dupliquée sur chaque harmonique considéré :

$$ \mathbf{Z}(\Omega) = \Omega^2 \nabla^2 \otimes \mathbf{M} + \Omega \nabla \otimes \mathbf{C} + \mathbf{I}_{H+1} \otimes \mathbf{K}. \quad (12) $$

Afin de réaliser le suivi de points de bifurcation et de dépasser les points de retournements une méthode de continuation est mise en œuvre. Pour se faire, deux paramètres sont traités comme des inconnues. Il s'agit ici de la fréquence $\Omega$ et du paramètre $\mu$. Deux équations additionnelles sont ajoutées
au résidu (11) afin de fermer le système. Ces équations sont l'équation de longueur d'arc $P(\tilde{q}, \Omega, \mu)$ qui contraint la solution recherchée à se situer à l'intersection entre la branche de solution et une hypersphère de rayon $\Delta s$, $s$ étant l'abscisse curviligne de la branche de solution :

$$P(\tilde{q}, \Omega, s) = (\Delta \tilde{q})^T (\Delta \tilde{q}) + \Delta \Omega^2 + \Delta \mu^2 - \Delta s^2 = 0. \quad (13)$$

avec $\Delta \tilde{q}$, $\Delta \Omega$ et $\Delta \mu$ les composantes du vecteur tangent à la branche de solutions au dernier point calculé $\Delta \tilde{Q} = (\Delta \tilde{q}, \Delta \Omega, \Delta \mu)^T$. La seconde équation sert à caractériser le type de point de bifurcation qui doit être suivi. Nous adoptons ici une approche reposant sur les systèmes augmentés minimaux [3, 4, 11] et considérons une équation scalaire définie, dans le cas des bifurcations nœud-col, par :

$$\begin{bmatrix} \partial_{\tilde{q}} \mathbf{R} & \mathbf{b} \\ \mathbf{d}^{\dagger} & 0 \end{bmatrix} \begin{pmatrix} \mathbf{w} \\ g \end{pmatrix} = \begin{pmatrix} \mathbf{0} \\ 1 \end{pmatrix} \quad (14)$$

où $\partial_{\tilde{q}} \mathbf{R}$ correspond à la dérivée du résidu HBM par rapport aux coefficients de Fourier des déplacements et où les vecteurs $b$ et $d$ sont choisi de façon à rendre le système inversible. Le système complet formé du résidu HBM et des deux équations additionnelles peut ensuite être résolu par une méthode de type Newton afin de calculer la courbe de bifurcations.

## 3.2 Caractérisation des singularités de codimension-1

Conformément aux pratiques habituelles en analyse de bifurcations, nous définissons des fonctions tests permettant de détecter des bifurcations spécifiques. Les points de retournement de la courbe de bifurcation sont facilement identifiables en surveillant le signe de la composante correspondant au paramètre de suivi $\mu$ du vecteur tangent :

$$\varphi_F = \Delta\mu \quad (15)$$

Cette fonction test ne permet cependant pas de discriminer les bifurcations de codimension-2 de type fromce (cusp) et les points de formation et de fusion d'isolas. Il est donc nécessaire de définir d'autres indicateurs afin de les identifier. Les bifurcations de type fromce peuvent être identifiées avec la fonction test suivante :

$$\varphi_{CP} = |\Delta\Omega| + |\Delta\mu| \quad (16)$$

Grâce à la théorie des singularités [5], il est possible de montrer que la courbure de la surface de solutions est localement positive aux points de formation d'isolas et négative aux points de fusion avec une autre branche. La surface de solutions est ainsi elliptique au voisinage d'un point de formation d'isola et hyperbolique au voisinage d'un point de fusion. Dans ces travaux, nous suivons le concept introduit dans [9] afin de déterminer le signe de la courbure de la surface de solutions via la fonction :

$$\varphi_{IS} = - \frac{\psi_0^T \partial_{\tilde{q}\tilde{q}} \mathbf{R} \phi_0 \phi_0 (\partial_s \Delta\Omega \Delta\mu - \partial_s \Delta\mu \Delta\Omega)}{\psi_0^T \partial_\mu \mathbf{R} \Delta\Omega - \psi_0^T \partial_\Omega \mathbf{R} \Delta\mu} \quad (17)$$

où $\partial_{\tilde{q}\tilde{q}} \mathbf{R}$ est le tenseur Hessien de $\mathbf{R}$ approximé par différences finies. $\partial_\mu \mathbf{R}$ and $\partial_\Omega \mathbf{R}$ sont directement obtenables dans la Jacobienne du résidu étendu calculée pendant le suivi de bifurcations. $\partial_s \Delta\Omega$ et $\partial_s \Delta\mu$ sont les dérivées par rapport à l'abscisse curviligne des composantes relatives à $\Omega$ et $\mu$ du vecteur tangent qui sont également calculées par différences finies. Notons que l'équation (17) suppose que le point n'est pas une bifurcation de type fromce (cf. Equation (16)) et donc que ces dernières soient identifiées avant les points de formation et de fusion d'isolas. Le point de retournement de la courbe de bifurcation est un point de formation d'isola si l'équation (17) est strictement positive et un point de fusion d'isola si strictement négative.

## 3.3 Processus général

A chaque itération, un suivi de bifurcation est effectué afin de détecter les singularités de codimension-1. Comme pour toute procédure de continuation, il est nécessaire de fournir un itéré initial au solveur. Dans ce travail, nous effectuons au préalable nous réalisons une analyse de réponse forcée à $\mu$ fixé dans un intervalle fréquentiel le plus petit possible afin de détecter une bifurcation nœud-col appartenant à
au résidu (11) afin de fermer le système. Ces équations sont l'équation de longueur d'arc $P(\tilde{q}, \Omega, \mu)$ qui contraint la solution recherchée à se situer à l'intersection entre la branche de solution et une hypersphère de rayon $\Delta s$, $s$ étant l'abscisse curviligne de la branche de solution :

$$P(\tilde{q}, \Omega, s) = (\Delta \tilde{q})^T (\Delta \tilde{q}) + \Delta \Omega^2 + \Delta \mu^2 - \Delta s^2 = 0. \quad (13)$$

avec $\Delta \tilde{q}$, $\Delta \Omega$ et $\Delta \mu$ les composantes du vecteur tangent à la branche de solutions au dernier point calculé $\Delta \tilde{Q} = (\Delta \tilde{q}, \Delta \Omega, \Delta \mu)^T$. La seconde équation sert à caractériser le type de point de bifurcation qui doit être suivi. Nous adoptons ici une approche reposant sur les systèmes augmentés minimaux [3, 4, 11] et considérons une équation scalaire définie, dans le cas des bifurcations nœud-col, par :

$$\begin{bmatrix} \partial_{\tilde{q}} \mathbf{R} & \mathbf{b} \\ \mathbf{d}^{\dagger} & 0 \end{bmatrix} \begin{pmatrix} \mathbf{w} \\ g \end{pmatrix} = \begin{pmatrix} \mathbf{0} \\ 1 \end{pmatrix} \quad (14)$$

où $\partial_{\tilde{q}} \mathbf{R}$ correspond à la dérivée du résidu HBM par rapport aux coefficients de Fourier des déplacements et où les vecteurs $b$ et $d$ sont choisi de façon à rendre le système inversible. Le système complet formé du résidu HBM et des deux équations additionnelles peut ensuite être résolu par une méthode de type Newton afin de calculer la courbe de bifurcations.

## 3.2 Caractérisation des singularités de codimension-1

Conformément aux pratiques habituelles en analyse de bifurcations, nous définissons des fonctions tests permettant de détecter des bifurcations spécifiques. Les points de retournement de la courbe de bifurcation sont facilement identifiables en surveillant le signe de la composante correspondant au paramètre de suivi $\mu$ du vecteur tangent :

$$\varphi_F = \Delta\mu \quad (15)$$

Cette fonction test ne permet cependant pas de discriminer les bifurcations de codimension-2 de type fromce (cusp) et les points de formation et de fusion d'isolas. Il est donc nécessaire de définir d'autres indicateurs afin de les identifier. Les bifurcations de type fromce peuvent être identifiées avec la fonction test suivante :

$$\varphi_{CP} = |\Delta\Omega| + |\Delta\mu| \quad (16)$$

Grâce à la théorie des singularités [5], il est possible de montrer que la courbure de la surface de solutions est localement positive aux points de formation d'isolas et négative aux points de fusion avec une autre branche. La surface de solutions est ainsi elliptique au voisinage d'un point de formation d'isola et hyperbolique au voisinage d'un point de fusion. Dans ces travaux, nous suivons le concept introduit dans [9] afin de déterminer le signe de la courbure de la surface de solutions via la fonction :

$$\varphi_{IS} = - \frac{\psi_0^T \partial_{\tilde{q}\tilde{q}} \mathbf{R} \phi_0 \phi_0 (\partial_s \Delta\Omega \Delta\mu - \partial_s \Delta\mu \Delta\Omega)}{\psi_0^T \partial_\mu \mathbf{R} \Delta\Omega - \psi_0^T \partial_\Omega \mathbf{R} \Delta\mu} \quad (17)$$

où $\partial_{\tilde{q}\tilde{q}} \mathbf{R}$ est le tenseur Hessien de $\mathbf{R}$ approximé par différences finies. $\partial_\mu \mathbf{R}$ and $\partial_\Omega \mathbf{R}$ sont directement obtenables dans la Jacobienne du résidu étendu calculée pendant le suivi de bifurcations. $\partial_s \Delta\Omega$ et $\partial_s \Delta\mu$ sont les dérivées par rapport à l'abscisse curviligne des composantes relatives à $\Omega$ et $\mu$ du vecteur tangent qui sont également calculées par différences finies. Notons que l'équation (17) suppose que le point n'est pas une bifurcation de type fromce (cf. Equation (16)) et donc que ces dernières soient identifiées avant les points de formation et de fusion d'isolas. Le point de retournement de la courbe de bifurcation est un point de formation d'isola si l'équation (17) est strictement positive et un point de fusion d'isola si strictement négative.

## 3.3 Processus général

A chaque itération, un suivi de bifurcation est effectué afin de détecter les singularités de codimension-1. Comme pour toute procédure de continuation, il est nécessaire de fournir un itéré initial au solveur. Dans ce travail, nous effectuons au préalable nous réalisons une analyse de réponse forcée à $\mu$ fixé dans un intervalle fréquentiel le plus petit possible afin de détecter une bifurcation nœud-col appartenant à
l'isola qui servira de point de départ à la procédure de suivi de bifurcation. Pour des raisons pratiques, nous sélectionnons $\mu$ suffisamment grand pour que l'isola ait fusionné avec la branche principale. Aussi, nous adaptons dynamiquement les bornes de l'intervalle de suivi de bifurcation afin de garantir la fusion de l'isola avec la branche principale et maintenir des temps de calcul raisonnables.

## 4 Exemple de résultat

Nous illustrons notre méthodologie sur le cas d'une poutre encastrée libre de longueur $L$ présentant une non-linéarité cubique à son extrémité libre (voir Fig. 4). Les caractéristiques mécaniques et géométriques initiales sont tirées de [13] et sont résumées dans la Table 1. Le système est discrétisé avec 59 éléments poutre d'Euler-Bernoulli possédant une translation transverse et une rotation dans le plan par nœud. Un forçage harmonique $F$ est appliqué à une distance $d = 0.21$ m de l'encastrement.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Schéma du système étudié.</figcaption></figure>

<table><caption>TABLE 1 – Propriétés géométriques et mécaniques de référence [13].</caption><thead><tr><th>Longueur (m)</th><th>Largeur (m)</th><th>Hauteur (m)</th></tr></thead><tbody><tr><td>0.7</td><td>0.014</td><td>0.014</td></tr><tr><td>Module d'Young (N/m²)</td><td>Masse volumique (kg/m³)</td><td>Coeff. non-linéaire (N/m³)</td></tr><tr><td>2.05e11</td><td>7800</td><td>6e9</td></tr></tbody></table>

Le vecteur des variables d'optimisation s'écrit : $\mathbf{x} = (L, \mathbf{w}_e^T, \mathbf{h}_e^T, k_{nl})^T$ avec $\mathbf{w}_e$ et $\mathbf{h}_e$ les vecteurs contenant les largeurs et hauteurs de chaque élément. Au total le problème d'optimisation possède 120 paramètres. Afin d'accélérer les calculs, un modèle réduit de la poutre est créé via la méthode de Craig-Bampton en retenant les deux ddl liés à la non-linéarité et au forçage extérieur ainsi que cinq modes à interface fixe. Par ailleurs la matrice d'amortissement est également mise à jour afin de maintenir des taux d'amortissement modaux égaux à 1% pour le premier mode et 5% pour tous les autres.

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Projections dans les plans amplitude-forçage (a) et forçage-fréquence (b) de la courbe de bifurcation nœud-col de la translation de l'extréminté libre de la poutre initiale (---) et optimisée (—). Les bifurcations de type froncé et les points de formation et de fusion d'isola sont représentés par des losanges verts, (♦), oranges (◆) et rouges (◆), respectivement.</figcaption></figure>
l'isola qui servira de point de départ à la procédure de suivi de bifurcation. Pour des raisons pratiques, nous sélectionnons $\mu$ suffisamment grand pour que l'isola ait fusionné avec la branche principale. Aussi, nous adaptons dynamiquement les bornes de l'intervalle de suivi de bifurcation afin de garantir la fusion de l'isola avec la branche principale et maintenir des temps de calcul raisonnables.

## 4 Exemple de résultat

Nous illustrons notre méthodologie sur le cas d'une poutre encastrée libre de longueur $L$ présentant une non-linéarité cubique à son extrémité libre (voir Fig. 4). Les caractéristiques mécaniques et géométriques initiales sont tirées de [13] et sont résumées dans la Table 1. Le système est discrétisé avec 59 éléments poutre d'Euler-Bernoulli possédant une translation transverse et une rotation dans le plan par nœud. Un forçage harmonique $F$ est appliqué à une distance $d = 0.21$ m de l'encastrement.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Schéma du système étudié.</figcaption></figure>

<table><caption>TABLE 1 – Propriétés géométriques et mécaniques de référence [13].</caption><thead><tr><th>Longueur (m)</th><th>Largeur (m)</th><th>Hauteur (m)</th></tr></thead><tbody><tr><td>0.7</td><td>0.014</td><td>0.014</td></tr><tr><td>Module d'Young (N/m²)</td><td>Masse volumique (kg/m³)</td><td>Coeff. non-linéaire (N/m³)</td></tr><tr><td>2.05e11</td><td>7800</td><td>6e9</td></tr></tbody></table>

Le vecteur des variables d'optimisation s'écrit : $\mathbf{x} = (L, \mathbf{w}_e^T, \mathbf{h}_e^T, k_{nl})^T$ avec $\mathbf{w}_e$ et $\mathbf{h}_e$ les vecteurs contenant les largeurs et hauteurs de chaque élément. Au total le problème d'optimisation possède 120 paramètres. Afin d'accélérer les calculs, un modèle réduit de la poutre est créé via la méthode de Craig-Bampton en retenant les deux ddl liés à la non-linéarité et au forçage extérieur ainsi que cinq modes à interface fixe. Par ailleurs la matrice d'amortissement est également mise à jour afin de maintenir des taux d'amortissement modaux égaux à 1% pour le premier mode et 5% pour tous les autres.

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Projections dans les plans amplitude-forçage (a) et forçage-fréquence (b) de la courbe de bifurcation nœud-col de la translation de l'extréminté libre de la poutre initiale (---) et optimisée (—). Les bifurcations de type froncé et les points de formation et de fusion d'isola sont représentés par des losanges verts, (♦), oranges (◆) et rouges (◆), respectivement.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Géométrie de la poutre optimisée. Vues de côté (a) et de dessus (b). La géométrie initiale est tracée en rouge.</figcaption></figure>

Une optimisation est effectuée afin de retarder la formation de l'isola, induit par une résonance interne 3 :1 entre le premier et second mode de la poutre. L'itéré initial pour l'optimisation sont les paramètres reportés dans la Table 1. Le critère d'arrêt est fixé à 100 itérations. Les résultats sont présentés en Figure 3 qui montre les projections des courbes de bifurcations nœud-col dans les plans amplitude-forçage et forçage-fréquence. Nous pouvons voir que le point de formation d'isola, symbolisé par le losange rouge (♦) est décalé vers des amplitudes de forçage plus importantes avec un gain de plus de 200 N. On remarque aussi un fort décalage en fréquence, de 15 Hz environ.

La figure 4 montre la géométrie de la poutre optimisée avec une vue de dessus et une vue de côté. Nous pouvons voir que la longueur de la poutre est inchangée. Par ailleurs, elle présente une augmentation de 25% de sa largeur et de sa hauteur au niveau de l'encastrement et de son milieu, ce qui est cohérent avec le fait que l'isola est associé à une interaction modale entre les premier et second modes de flexion. Ceci se traduit par une augmentation en masse de 28%. A titre indicatif, pour ce système, le temps de calcul global est d'environ 2.5 h sur un PC portable (i7-12700H @ 2.30 GHz, 16 Go RAM) opérant sous Linux avec Julia 1.9.3.

La méthodologie numérique proposée ici permet bien de décaler de façon importante le point de naissance de l'isola pour une structure mécanique en modifiant sa forme.

# 5 Conclusion

Les courbes de réponses isolées sont un phénomène non-linéaire important pouvant avoir d'importantes conséquences sur les performances et la tenue en service des systèmes mécaniques. Nous avons proposé une méthodologie générale permettant de contrôler leur position dans un espace de codimension-2. Cette dernière met en œuvre la méthode de l'équilibrage harmonique et des systèmes augmentés minimaux afin de suivre l'évolution des points de bifurcation nœud-col. Des indicateurs ont également été définis grâce à la théorie des singularités afin de détecter et classifier les bifurcations de type fronce ainsi que les points de formation et de fusion d'isolas.

Bien que nous ayons choisi la méthode de l'équilibrage harmonique pour ses bonnes performances en calcul de structures, la méthodologie présentée est tout à fait générale et peut être facilement adaptée
<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Géométrie de la poutre optimisée. Vues de côté (a) et de dessus (b). La géométrie initiale est tracée en rouge.</figcaption></figure>

Une optimisation est effectuée afin de retarder la formation de l'isola, induit par une résonance interne 3 :1 entre le premier et second mode de la poutre. L'itéré initial pour l'optimisation sont les paramètres reportés dans la Table 1. Le critère d'arrêt est fixé à 100 itérations. Les résultats sont présentés en Figure 3 qui montre les projections des courbes de bifurcations nœud-col dans les plans amplitude-forçage et forçage-fréquence. Nous pouvons voir que le point de formation d'isola, symbolisé par le losange rouge (♦) est décalé vers des amplitudes de forçage plus importantes avec un gain de plus de 200 N. On remarque aussi un fort décalage en fréquence, de 15 Hz environ.

La figure 4 montre la géométrie de la poutre optimisée avec une vue de dessus et une vue de côté. Nous pouvons voir que la longueur de la poutre est inchangée. Par ailleurs, elle présente une augmentation de 25% de sa largeur et de sa hauteur au niveau de l'encastrement et de son milieu, ce qui est cohérent avec le fait que l'isola est associé à une interaction modale entre les premier et second modes de flexion. Ceci se traduit par une augmentation en masse de 28%. A titre indicatif, pour ce système, le temps de calcul global est d'environ 2.5 h sur un PC portable (i7-12700H @ 2.30 GHz, 16 Go RAM) opérant sous Linux avec Julia 1.9.3.

La méthodologie numérique proposée ici permet bien de décaler de façon importante le point de naissance de l'isola pour une structure mécanique en modifiant sa forme.

# 5 Conclusion

Les courbes de réponses isolées sont un phénomène non-linéaire important pouvant avoir d'importantes conséquences sur les performances et la tenue en service des systèmes mécaniques. Nous avons proposé une méthodologie générale permettant de contrôler leur position dans un espace de codimension-2. Cette dernière met en œuvre la méthode de l'équilibrage harmonique et des systèmes augmentés minimaux afin de suivre l'évolution des points de bifurcation nœud-col. Des indicateurs ont également été définis grâce à la théorie des singularités afin de détecter et classifier les bifurcations de type fronce ainsi que les points de formation et de fusion d'isolas.

Bien que nous ayons choisi la méthode de l'équilibrage harmonique pour ses bonnes performances en calcul de structures, la méthodologie présentée est tout à fait générale et peut être facilement adaptée
à d'autres méthodes de résolution comme les méthodes de tir ou de collocation orthogonales.

Références

[1] H.N. Abramson. *Response curves for a system with softening restoring force*, Journal of Applied Mechanics **22**, 434–435, 1955.

[2] R. Alcorta, S. Baguet, B. Prabel, P. Piteau, G. Jacquet-Richardet. *Period doubling bifurcation analysis and isolated sub-harmonic resonances in an oscillator with asymmetric clearances*. Nonlinear Dynamics **98**, 2939–2960, 2019.

[3] E.L. Allgower, H. Schwetlick. *A general view of minimally extended systems for simple bifurcation points*, ZAMM J. Appl. Math. Mech. Zeitschrift für Angewandte Mathematik und Mechanik **77**(2), 83–97, 1997.

[4] T. Detroux, L. Renson, L. Masset, G. Kerschen. *The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems*. Computer Methods in Applied Mechanics and Engineering **296**, 18–38, 2015.

[5] M. Golubitsky, D. G. Schaeffer. *Singularities and Groups in Bifurcation Theory*, Springer New York, NY, 1984.

[6] C. Grenat, S. Baguet, C.H. Lamarque, R. Dufour. *Mass sensing by symmetry breaking and localization of motion in an array of electrostatically coupled nonlinear MEMS resonators*. International Journal of Non-Linear Mechanics **140**, 103903, 2022.

[7] G. Habib, G.I. Cirillo, G. Kerschen. *Isolated resonances and nonlinear damping*, Nonlinear Dynamics **93**, 979–994, 2018.

[8] T. Hill, S. Neild, A. Cammarano. *An analytical approach for detecting isolated periodic solution branches in weakly nonlinear structures*, Journal of Sound and Vibration **379**, 150–165, 2016.

[9] A. Jepson, A. Spence. *Folds in solutions of two parameter systems and their calculation. part I*. SIAM Journal on Numerical Analysis **22**, 347–368, 1985.

[10] M. Krack, J. Gross. *Harmonic Balance for Nonlinear Vibration Problems*, Springer, 2019.

[11] A. Mélot, E. Rigaud, J. Perret-Liaudet. *Bifurcation tracking of geared systems with parameter-dependent internal excitation*, Nonlinear Dynamics **107**, 413–431, 2023.

[12] S. Ponsioen, T. Pedergnana, G. Haller. *Analytic prediction of isolated forced response curves from spectral submanifolds*, Nonlinear Dynamics **98**, 2755–2773, 2019.

[13] F. Thouverez. *Presentation of the ECL benchmark*, Mechanical Systems and Signal Processing **17**, 195–202, 2003.

[14] M. Urabe. *Galerkin's procedure for nonlinear periodic systems*, Archive for Rational Mechanics and Analysis **20**, 120–152, 1965.
à d'autres méthodes de résolution comme les méthodes de tir ou de collocation orthogonales.

Références

[1] H.N. Abramson. *Response curves for a system with softening restoring force*, Journal of Applied Mechanics **22**, 434–435, 1955.

[2] R. Alcorta, S. Baguet, B. Prabel, P. Piteau, G. Jacquet-Richardet. *Period doubling bifurcation analysis and isolated sub-harmonic resonances in an oscillator with asymmetric clearances*. Nonlinear Dynamics **98**, 2939–2960, 2019.

[3] E.L. Allgower, H. Schwetlick. *A general view of minimally extended systems for simple bifurcation points*, ZAMM J. Appl. Math. Mech. Zeitschrift für Angewandte Mathematik und Mechanik **77**(2), 83–97, 1997.

[4] T. Detroux, L. Renson, L. Masset, G. Kerschen. *The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems*. Computer Methods in Applied Mechanics and Engineering **296**, 18–38, 2015.

[5] M. Golubitsky, D. G. Schaeffer. *Singularities and Groups in Bifurcation Theory*, Springer New York, NY, 1984.

[6] C. Grenat, S. Baguet, C.H. Lamarque, R. Dufour. *Mass sensing by symmetry breaking and localization of motion in an array of electrostatically coupled nonlinear MEMS resonators*. International Journal of Non-Linear Mechanics **140**, 103903, 2022.

[7] G. Habib, G.I. Cirillo, G. Kerschen. *Isolated resonances and nonlinear damping*, Nonlinear Dynamics **93**, 979–994, 2018.

[8] T. Hill, S. Neild, A. Cammarano. *An analytical approach for detecting isolated periodic solution branches in weakly nonlinear structures*, Journal of Sound and Vibration **379**, 150–165, 2016.

[9] A. Jepson, A. Spence. *Folds in solutions of two parameter systems and their calculation. part I*. SIAM Journal on Numerical Analysis **22**, 347–368, 1985.

[10] M. Krack, J. Gross. *Harmonic Balance for Nonlinear Vibration Problems*, Springer, 2019.

[11] A. Mélot, E. Rigaud, J. Perret-Liaudet. *Bifurcation tracking of geared systems with parameter-dependent internal excitation*, Nonlinear Dynamics **107**, 413–431, 2023.

[12] S. Ponsioen, T. Pedergnana, G. Haller. *Analytic prediction of isolated forced response curves from spectral submanifolds*, Nonlinear Dynamics **98**, 2755–2773, 2019.

[13] F. Thouverez. *Presentation of the ECL benchmark*, Mechanical Systems and Signal Processing **17**, 195–202, 2003.

[14] M. Urabe. *Galerkin's procedure for nonlinear periodic systems*, Archive for Rational Mechanics and Analysis **20**, 120–152, 1965.