<figure><img src="image_94.png" /><figcaption>FIGURE C.12. Same numerical simulations as in (a) Figs. C.10(f) and (b) C.10(g) depicting the trajectory of the AMPD in the (b, c, a)-space (blue) superimposed on the critical manifold $M_{01}$ (attracting : red solid line; saddle-type : red dashed line) and $M_{02}$ (attracting : green solid line; saddle-type : green dashed line). Unstable fixed points of the AMPD on $M_{01}$ are also represented by black •.</figcaption></figure>

methodologies could be associated to the present results to analyze the behavior of a mechanical system with one unstable mode and coupled to a BNES.

Finally, because of its general nature, the present study could be also used in the future to investigate passive attenuation of free and forced vibrations. It could also be generalized to the case of a multistable NES or to a primary mechanical system with more than one unstable mode.

## C.7 APPENDIX A. APPROXIMATE ANALYTICAL SOLUTION OF THE VAN DER POL OSCILLATOR

The equation of motion of the VdP oscillator alone is

$$ \ddot{x} + \epsilon h(x, \dot{x}) + x = 0. \qquad (C.58) $$

with $h(x, \dot{x}) = -\sigma\dot{x} + \lambda\dot{x}^2$.

An approximate analytical solution of Eq. (C.58) is found using the Krylov-Bogoliubov method of averaging (see e.g. [143]). For that, an amplitude phase representation of the motion is first introduced as

$$ x = r(t) \cos (t + \varphi(t)) \qquad (C.59a) $$

$$ \dot{x} = -r(t) \sin (t + \varphi(t)) \qquad (C.59b) $$

The desired form of Eq. (C.59) requires that

$$ \dot{r} \cos \gamma - r \dot{\varphi} \sin \gamma = 0, \qquad (C.60) $$

with $\gamma = t + \varphi$ which yields

$$ \dot{\varphi} = \frac{\dot{r} \cos \gamma}{r \sin \gamma} \quad \text{and} \quad \dot{r} = r \dot{\varphi} \frac{\sin \gamma}{\cos \gamma}. \qquad (C.61) $$

The time derivative of Eq. (C.59) is

$$ \ddot{r} = -r \cos \gamma - \dot{r} \sin \gamma - r \dot{\varphi} \cos \gamma. \qquad (C.62) $$

The substitution of (C.59) and (C.62) into (C.58) and the use of (C.61) leads to

$$ \dot{r} = h (r \cos \gamma, -r \sin \gamma) \sin \gamma \qquad (C.63a) $$
<figure><img src="image_95.png" /><figcaption>FIGURE C.13. Numerical solutions of the fast subsystem (C.31) for (a) $a = 0.07$, (b) 0.14, (c) 0.25 and (d) 0.4. Parameters are the same as in Fig. C.10(f) and the following initial conditions are used: $b(0) = -0.2$, $c(0) = 0.3$ and $\delta(0) = \pi$.</figcaption></figure>

$$ \dot{\varphi} = h(r \cos \gamma, -r \sin \gamma) \frac{\cos \gamma}{r}. \qquad (C.63b) $$

For $0 < \epsilon \ll 1$, $\dot{r}$ and $\dot{\gamma}$ are small. Therefore, $r$ and $\varphi$ vary much more slowly with $t$ than $\gamma = t + \varphi$. Consequently, the Krylov-Bogoliubov approximation is used, i.e., Eq. (C.63) is averaged over the period $2\pi$ considering $r$, $\varphi$, $\dot{r}$ and $\dot{\varphi}$ are constant. That yields the following system

$$ \dot{r} = \epsilon \frac{1}{8} r (4\sigma - \lambda r^2) \qquad (C.64a) $$

$$ \dot{\varphi} = 0 \qquad (C.64b) $$

that approximates Eq. (C.63).

In this approximation, the phase $\varphi$ is constant (and equal to the initial phase denoted as $\varphi_0$) and uncoupled from the amplitude $r$. The amplitude equation (C.64a) can be solved analytically as

$$ r(t) = \frac{2\sqrt{\sigma}e^{\frac{\sigma t\epsilon}{2}}}{\sqrt{\frac{4\sigma}{r_0^2} + \lambda(e^{\sigma t\epsilon} - 1)}}, \qquad (C.65) $$

where $r_0 = r(t = 0)$.

## C.8 APPENDIX B. ANALYTICAL EXPRESSIONS OF SOME QUANTITIES THAT APPEAR IN SECTION C.4

$$ c_1^{\text{LF}} = \frac{2}{3} \sqrt{\frac{2\beta + 2 - \sqrt{(\beta + 1)^2 - 3\mu^2}}{\alpha}}, \qquad (C.66a) $$
<figure><img src="image_96.png" /><figcaption>FIGURE C.14. Numerical solutions of the fast subsystem (C.31) for (a) $a = 0.22$, (b) 0.35 and (c) 0.45. Parameters are the same as in Fig. C.10(g) and the following initial conditions are used : $b(0) = -0.2$, $c(0) = 0.3$ and $\delta(0) = \pi$.</figcaption></figure>

$$c_1^{\text{RF}} = \frac{2}{3} \sqrt{\frac{2\beta + 2 + \sqrt{(\beta + 1)^2 - 3\mu^2}}{\alpha}} \quad (C.66b)$$

$$c_2^{\text{LF}} = \frac{2}{3} \sqrt{\frac{4\beta - 2 - \sqrt{(1 - 2\beta)^2 - 3\mu^2}}{5\alpha}}, \quad (C.67a)$$

$$c_2^{\text{RF}} = \frac{2}{3} \sqrt{\frac{4\beta - 2 + \sqrt{(1 - 2\beta)^2 - 3\mu^2}}{5\alpha}} \quad (C.67b)$$

$$c_1^D = \frac{2}{3} \sqrt{\frac{2\beta + 2 - \sqrt{((\beta + 1)^2 - 3\mu^2)^{3/2} - \frac{(\beta+1)^2 - 3\mu^2}{3\sqrt{((\beta+1)^2 - 3\mu^2)^{3/2}}}}{\alpha}}, \quad (C.68a)$$

$$c_1^U = \frac{2}{3} \sqrt{\frac{2\beta + 2 + \sqrt{((\beta + 1)^2 - 3\mu^2)^{3/2} + \frac{(\beta+1)^2 - 3\mu^2}{3\sqrt{((\beta+1)^2 - 3\mu^2)^{3/2}}}}{\alpha}}, \quad (C.68b)$$

$$c_2^D = \frac{2}{3} \sqrt{\frac{4\beta - 2 - \sqrt{((1 - 2\beta)^2 - 3\mu^2)^{3/2} - \frac{(1 - 2\beta)^2 - 3\mu^2}{3\sqrt{((1 - 2\beta)^2 - 3\mu^2)^{3/2}}}}{5\alpha}}, \quad (C.69a)$$
$$c_2^U = \frac{2}{3} \sqrt{\frac{4\beta - 2 + \sqrt{((1-2\beta)^2 - 3\mu^2)^{3/2}} + \frac{((1-2\beta)^2 - 3\mu^2)^{3/2})^{2/3}}{5\alpha}} \quad (C.69b)$$

$$F_1 = \frac{c \left( \sigma (3\alpha c^2 - 4(\beta + 1))^2 + 16\mu(\mu\sigma - 1) \right)}{8\sqrt{(3\alpha c^2 - 4(\beta + 1))^2 + 16\mu^2}} \quad (C.70a)$$

$$F_2 = \frac{c \left( \sigma (-8\beta + 15\alpha c^2 + 4)^2 + 16\mu(\mu\sigma - 1) \right)}{8\sqrt{(-8\beta + 15\alpha c^2 + 4)^2 + 16\mu^2}} \quad (C.70b)$$

$$c_{1,1}^* = \sqrt{\frac{1 + \beta - \sqrt{\mu(\frac{1}{\sigma} - \mu)}}{3\alpha}} \quad (C.71a)$$

$$c_{1,2}^* = \sqrt{\frac{1 + \beta + \sqrt{\mu(\frac{1}{\sigma} - \mu)}}{3\alpha}} \quad (C.71b)$$

$$c_{2,1}^* = 2\sqrt{\frac{2\beta - 1 - \sqrt{\mu(\frac{1}{\sigma} - \mu)}}{15\alpha}} \quad (C.72a)$$

$$c_{2,2}^* = 2\sqrt{\frac{2\beta - 1 + \sqrt{\mu(\frac{1}{\sigma} - \mu)}}{15\alpha}} \quad (C.72b)$$

## C.9 APPENDIX C. EXPRESSION OF THE REGULAR NON NORMALLY HYPERBOLIC FIXED POINT OF THE BRANCH $\mathcal{M}_{02}$ OF THE CRITICAL MANIFOLD

Substituting the explicit expressions of the functions $H_{01}$, $H_{02}$ and $H_{03}$ and their derivative into Eq. (C.48) we obtain

$$\mathbf{J}_g(p_2) = \begin{pmatrix} \frac{3\alpha c^2 - 2\beta}{\mu} & \frac{3\sqrt{\alpha c} \sqrt{2\beta - 3\alpha c^2}}{\sqrt{2\mu}} & 0 \\ \frac{3\sqrt{2\sqrt{\alpha c}\mu \sqrt{2\beta - 3\alpha c^2}}}{\mu^2 + 4} & \frac{\mu(-8\beta + 9\alpha c^2 - 4)}{4(\mu^2 + 4)} & -\frac{c(-8\beta + 15\alpha c^2 + 2\mu^2 + 4)}{2(\mu^2 + 4)} \\ \frac{6\sqrt{2\sqrt{\alpha} \sqrt{2\beta - 3\alpha c^2}}}{\mu^2 + 4} & \frac{-8\beta + 9\alpha c^2 + 2\mu^2 + 4}{2c\mu^2 + 8c} & \frac{\mu(-8\beta + 15\alpha c^2 - 4)}{4(\mu^2 + 4)} \end{pmatrix} \quad (C.73)$$

whose third-order characteristic polynomial $p_{\mathbf{J}_g(p_2)}(z)$ has the following form

$$p_{\mathbf{J}_g(p_2)}(z) = a_0 z^3 + a_1 z^2 + a_2 z + a_3. \quad (C.74)$$

For such a third-order system, the Routh-Hurwitz stability criterion (see e.g. [32]) states that the roots of (C.74) have negative real parts if and only if the following inequalities

$$a_0 > 0, \quad a_1 > 0, \quad a_2 > 0, \quad a_3 > 0, \quad d_2 = a_1 a_2 - a_0 a_3 > 0 \quad (C.75)$$

are satisfied.

The coefficient $a_i$ ($i = 0, \dots, 3$) and $d_2$ are plotted in Fig. C.15 as functions of $c$ for the same parameters as used in Fig. C.6. Comparing Figs. C.7 and C.15 we deduce that the value of $c$ corresponding to the regular non normally hyperbolic fixed point of the branch $\mathcal{M}_{02}$ of the critical manifold is one of the solutions of $d_2 = 0$.

The coefficient $d_2$ is a third-order polynomial with respect to $x = c^2$ with the following form

$$d_2 = b_1 x^3 + b_2 x^2 + b_3 x + b_4. \quad (C.76)$$

The roots of $d_2$ are computed by means of the Cardano's method (see e.g. [185]). First, the following parameters are introduced

$$p = -\frac{b_2^2}{3b_1^2} + \frac{c_3}{b_1}, \quad q = \frac{b_2}{27b_1} \left( \frac{2b_2^2}{b_1^2} - \frac{9b_3}{b_1} \right) + \frac{b_4}{b_1}.$$

The discriminant $\Delta$ is defined as $\Delta = -(4p^3 + 27q^2)$. Then:
<figure><img src="image_97.png" /><figcaption>FIGURE C.15. The coefficients $a_i$ ($i = 0, \dots, 3$) and $d_2$ are plotted in Fig. C.15 as functions of $c$ for the same parameters as used in Fig. C.6.</figcaption></figure>

1. If $\Delta < 0$, one root is real and two are complex conjugate.

2. If $\Delta = 0$, all roots are real and at least two are equal.

3. If $\Delta > 0$, all roots are real and unequal.

For the parameters used in this paper one has $\Delta > 0$. In this case, the three real roots $x_k$ ($k = 0, 1, 2$) are given by

$$x_k = 2\sqrt{\frac{-p}{3}} \cos\left(\frac{1}{3}\arccos\left(\frac{3q}{2p}\sqrt{\frac{3}{-p}}\right) + \frac{2k\pi}{3}\right) - \frac{b_2}{3b_1}, \quad k = 0, 1, 2. \qquad (C.77)$$

The wanted root is $x_1$. Then for the parameters used in Fig. C.6 we obtain $c = \sqrt{x_1} = 0.325914$.
Journal of Sound and Vibration 503 (2021) : 116109 [a2]

## SCALING LAW FOR THE SLOW FLOW OF AN UNSTABLE MECHANICAL
SYSTEM COUPLED TO A NONLINEAR ENERGY SINK

Baptiste BERGEOT

### ABSTRACT

In this paper one first shows that the slow flow of a mechanical system with one unstable mode coupled to a Nonlinear Energy Sink (NES) can be reduced, in the neighborhood of a fold point of its critical manifold, to a normal form of the dynamic saddle-node bifurcation. This allows us to then obtain a scaling law for the slow flow dynamics and to improve the accuracy of the theoretical prediction of the mitigation limit of the NES previously obtained as part of a zeroth-order approximation. For that purpose, the governing equations of the coupled system are first simplified using a reduced-order model for the primary structure by keeping only its unstable modal coordinates. The slow flow is then derived by means of the complexification- averaging method and, by the presence of a small perturbation parameter related to the mass ratio between the NES and the primary structure, it appears as a fast-slow system. The center manifold theorem is finally used to obtain the reduced form of the slow flow which is solved analytically leading to the scaling law. The latter reveals a nontrivial dependence with respect to the small perturbation parameter of the slow flow dynamics near the fold point, involving the fractional exponents 1/3 and 2/3. Finally, a new theoretical prediction of the mitigation limit is deduced from the scaling law. In the end, the proposed methodology is exemplified and validated numerically using an aeroelastic aircraft wing model coupled to one NES.

**Keywords:** Passive vibration control, Nonlinear energy sink, Center manifold reduction, Scaling law, Aeroelastic instability

### D.1 INTRODUCTION

The nonlinear energy sinks (NESs) are nowadays well-known devices used for passive shock and vibration mitigation of undesired oscillations caused by either external, parametric or self-excitations of a primary structure. An NES is classically defined as a nonlinear dynamical attachment consisting of a light mass (compared to the total mass of the primary structure), an essentially nonlinear spring (most of the time purely cubic) and a viscous linear damper. The operating of the NES is based on the phenomenon of targeted energy transfer (TET) through which a properly designed strongly nonlinear oscillator can be tuned to any frequency in order to perform an irreversible energy transfer from the primary structure towards itself. In their seminal papers [76, 203] Gendelman, Vakakis and co-workers explain the TET phenomenon by the interaction between two nonlinear modes of vibrations of the system producing a 1:1 resonance capture. Reviews of these concepts can be found in [206] and more recently in a part of [129].

The use of NES to mitigate limit cycle oscillations (LCOs) stemming from dynamic instabilities has been widely studied in the past. The first work reported in this framework concerns the mitigation of LCOs of the Van der Pol oscillator [122]. By means of a perturbation analysis, Gendelman and Bar [75] predicted the response regimes of the same system. This work has been extended to a Van der Pol-Duffing oscillator coupled to one NES by Domany and Gendelman [59]. Numerous papers have been dedicated to the problem of mitigation, by means of one or several NESs, of LCOs due to aeroelastic instabilities. The seminal papers are those of Lee et al. [121, 123, 124] and concerned aeroelastic instabilities in aircraft wing. In these works, the problem has been investigated both numerically and experimentally. The same aircraft wing model coupled to one NES has been studied theoretically using the so-called complexification-averaging method [135] together with a singular perturbation approach [78] or a perturbation algorithm based on a mixed multiple scale-harmonic balance method (MSHBM) [131]. Vaurigaud et al. [207] investigated a problem of passive nonlinear TET between a two degrees of freedom long span bridge model prone to coupled flutter and an NES. The mitigation, using an NES, of vortex-induced vibrations resulting from the nonlinear interaction of a laminar
flow and a rigid circular cylinder has been first studied by Tumkur et al. [201] constructing a two-DOF reduced-order model of the system and validating it by means of a comparison with a full order finite-element model. On the same issue, an improved and experimentally validated reduced model has been proposed by Dai et al. [47]. The passive control of helicopter ground resonance instabilities by a means of NESs has been theoretically and numerically studied by Bergeot et al. considering an NES attached to the helicopter fuselage [a10] and to the blades [a12]. Chatter control in machine tool vibrations has been studied by Gourc et al. [82] considering a vibro-impact NES and by Nankali et al. [142] considering a purely cubic NES. In [a14], mitigation of friction-induced vibrations due to mode coupling instabilities in a friction system has been investigated. The use of several NESs to mitigate LCOs has been analytically investigated by Bergeot et Bellizzi. First, a network of parallel NESs coupled to a Van der Pol oscillator has been analyzed in [a7] and then the prediction of the steady-state regimes of a multi-DOF dynamical system having one unstable mode and coupled to a set of NESs has been performed in [a8]. Finally, in a recent paper, the possibility of mitigating simultaneously two unstable modes of a linear multi-DOF primary system using a single NES has been investigated by Bergeot et al. [a9]. The asymptotic analysis proposed in the latter reference reveals that the NES attachment can produce several and complex mitigated responses which results from the presence of several stable solutions and from the competition between them.

In general, when an NES is attached to a primary structure, the resulting coupled model is analyzed by introducing a small perturbation parameter related to the mass ratio between the NES and the primary structure. Under the assumption of a 1 : 1 resonance capture, the system is first averaged over one period corresponding to a natural frequency of the primary structure. The resulting averaged system, called slow flow, is then analyzed by means of singular perturbation techniques (multiple scales method [144] or geometric singular perturbation theory [107]). The first key point of these analytical treatments is the partition of the slow flow into two time scales (one slow and one fast $^{1}$). In this representation the time evolution of the slow flow is thus described as a succession of slow and fast epochs which are described independently. The second key point is the definition of the critical manifold whose system trajectories approaches during slow epochs. In general, the analysis is carried out within the zeroth-order approximation, i.e. within the limit case in which the perturbation parameter is equal to zero. Consequently, the obtained analytical results depreciate for the largest values of the perturbation parameter by thus limiting their predictive power. Luongo and Zulli [131] used a perturbation technique to express the invariant manifold as a power series with respect to the perturbation parameter in which the first term is the critical manifold $^{2}$. As stated by the Fenichel theory [62], this representation is relevant for normally hyperbolic branch of the critical manifold but fails, at fold points (the junction between slow and fast epochs takes place in the vicinity of these points), to approach the actual trajectory of the system. Indeed, the literature on dynamical systems shows that the trajectory passes through the region near these fold points with a nontrivial scaling behavior with respect to the perturbation parameter that it is not possible to obtain with the classical perturbation methods (see e.g. [40] Sect. 5.4).

By means of a sophisticated analysis, the originality of the present work is to provide an analytical descrip-
tion of the system dynamic behavior which reports on this nontrivial dependence on the small perturbation
parameter. This analytical result is hereafter referred as the scaling law of the slow flow, this notion is cla-
rified below within the main body of the present manuscript. The main objective is to obtain an analytical
description of the system with a high predictive power compared to what has been proposed in the past in
the literature and then to hope in the future of including them in a NES design procedure.

The equations of motion of the system under study are presented in Section. E.2. The full order model equations of a mechanical system coupled to one NES are derived in Section. D.2.1. Assuming a primary structure with only one unstable mode, the modal reduction technique, leading to the reduced-order dimensionless model, is presented in Section. D.2.2. Section. D.3 summarizes and extends some previous analytical results of the literature obtained within the zeroth-order approximation. Section. D.4 expounds the main results of the paper. In Section. D.4.1, the center manifold theorem is used to reduce the dynamics of the slow flow of the system. This reduced form is solved analytically in Section. D.4.2 leading to the scaling law from which a new analytical prediction of the mitigation limit is derived in Section. D.4.3. In Section. D.5, the proposed methodology is exemplified and validated numerically using an aeroelastic aircraft wing model coupled to one NES. Finally, concluding remarks are formulated in Section. F.8.

<footnote n="1">1. Sometimes the terms “slow” and “super-slow” [59, 75] are used to leave the term “fast” for the oscillations on which the average was carried out and therefore to be coherent with term slow flow. In the present work the terms “slow” and “fast” are preferred to be in agreement with the vocabulary usually used in the literature on dynamical systems.</footnote>

<footnote n="2">2. The difference between critical and invariant manifolds will be calrified in the manuscript in Section. D.4.</footnote>
## D.2 EQUATIONS OF THE MODEL

The dynamical system under study in this paper is presented in the present section. The following developments are similar to those presented in [a8].

### D.2.1 The initial model

The primary structure under consideration in this work is described by the following general set of nonlinear differential equations

$$ \tilde{M}\ddot{\tilde{x}} + \tilde{C}\dot{\tilde{x}} + \tilde{K}\tilde{x} + \tilde{g}^{\text{NL}}(\tilde{x}) = 0, \quad (\text{D.1}) $$

where $\tilde{x} = (\tilde{x}_1, \dots, \tilde{x}_N)^T$ with $(\cdot)^T$ the transpose operator, the dot represents time-differentiation and $\tilde{M}, \tilde{C}$ and $\tilde{K}$ are the mass matrix, the damping matrix and the stiffness matrix respectively. The nonlinear vector-valued function $\tilde{g}^{\text{NL}}$ represents the nonlinearity of the primary system. Moreover, cubic nonlinearities are assumed, i.e. each component of $\tilde{g}^{\text{NL}}$ is a linear combination of monomial terms of order 3. In these first equations, the tilde symbol is used before the rescaling which will be carried out below (to obtain variables without tilde after rescaling).

System (D.1) can undergo a dynamic instability of the trivial solution through a Hopf bifurcation, i.e. it has a pair of complex conjugate eigenvalues which cross the complex plane imaginary axis. The loss of stability of the trivial solution comes with the production of a periodic solution and the nonlinear function $\tilde{g}^{\text{NL}}$ allows the existence of Limit Cycle Oscillations (LCOs) on which the system can saturate. For the purpose of mitigating the LCOs, one purely cubic ungrounded NES with mass $\tilde{m}_h$, damping coefficient $\tilde{c}_h$ and cubic stiffness $k_h^{\text{NL}}$ is used. Taking into account the NES displacement $\tilde{h}(t)$, the equations of the coupled system are

$$ \tilde{M}\ddot{\tilde{x}} + \tilde{C}\dot{\tilde{x}} + \tilde{K}\tilde{x} + \tilde{g}^{\text{NL}}(\tilde{x}) + \tilde{B}(\tilde{c}_h(\mathbf{A}\dot{\tilde{x}} - \dot{\tilde{h}}) + \tilde{k}_h^{\text{NL}}(\mathbf{A}\tilde{x} - \tilde{h})^3) = 0 \quad (\text{D.2a}) $$

$$ \tilde{m}_h \ddot{\tilde{h}} - \tilde{c}_h (\mathbf{A} \dot{\tilde{x}} - \dot{\tilde{h}}) - \tilde{k}_h^{\text{NL}} (\mathbf{A} \tilde{x} - \tilde{h})^3 = 0, \quad (\text{D.2b}) $$

where $\mathbf{A} = (a_1, \dots, a_N)$ and $\tilde{\mathbf{B}} = (\tilde{b}_1, \dots, \tilde{b}_N)^T$ are the influence coefficient vectors which depend on the position of the NES.

One assumes now that the mass $\tilde{m}_h$ and the damping coefficient $\tilde{c}_h$ of the NES are small compared to the mass of the primary structure (D.1) and also that the latter is weakly nonlinear. We therefore introduce a small perturbation dimensionless parameter $\epsilon$ (with $0 < \epsilon \ll 1$) and the associated rescaled coefficients $m_h$ and $c_h$ and vector function $\mathbf{g}^{\text{NL}}$ as $\tilde{m}_h = \epsilon m_h$, $\tilde{c}_h = \epsilon c_h$ and $\tilde{\mathbf{g}}^{\text{NL}} = \epsilon \mathbf{g}^{\text{NL}}$. Then, the variables $\tilde{x}$ and $\tilde{h}$ are also rescaled through $\epsilon$ as $\mathbf{x} = \frac{\tilde{x}}{\sqrt{\epsilon}}$ and $\tilde{h} = \frac{\tilde{h}}{\sqrt{\epsilon}}$. Due to the latter rescaling and because cubic nonlinearities are assumed one has $\mathbf{g}^{\text{NL}}(\tilde{x}) = \epsilon^{3/2} \mathbf{g}^{\text{NL}}(\mathbf{x})$. Inserting $m_h, c_h, \mathbf{g}^{\text{NL}}, \mathbf{x}$ and $\tilde{h}$ into (D.2), multiplying Eq. (D.2a) by $\tilde{M}^{-1}$ and Eq. (D.2b) by $1/(\epsilon m_h)$ and finally neglecting terms of order strictly higher than 1 in $\epsilon$ lead to

$$ \ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{K}\mathbf{x} + \epsilon\mathbf{B}\dot{\tilde{h}} = 0 \quad (\text{D.3a}) $$

$$ \ddot{h} - \eta_h (\mathbf{A}\dot{\mathbf{x}} - \dot{\mathbf{h}}) - k_h^{\text{NL}} (\mathbf{A}\mathbf{x} - \mathbf{h})^3 = 0. \quad (\text{D.3b}) $$

with $\mathbf{C} = \tilde{M}^{-1}\tilde{\mathbf{C}}$, $\mathbf{K} = \tilde{M}^{-1}\tilde{\mathbf{K}}$, $\mathbf{B} = m_h\tilde{M}^{-1}\tilde{\mathbf{B}}$, $\eta_h = c_h/m_h$ and $k_h^{\text{NL}} = \tilde{k}_h^{\text{NL}}/m_h$.

Note that because of the previous assumptions and rescaling, the nonlinearities of the primary structure are neglected from now. Therefore, nonlinear modal interactions in the primary structure such as internal resonances are not taken into account in this work.

### D.2.2 Reduction of the dynamics

The system (D.3) is now reduced taking into account that the primary system (D.1) undergoes a single instability of the trivial solution. The objective is to obtain a system with two degrees of freedom (one for the unstable mode of the primary structure and one for the NES) on which one will be able to carry out an analysis which would be impossible on a the original system.

To achieve that, it is first convenient to introduce new coordinates $\mathbf{v} = \mathbf{x} + \epsilon\mathbf{B}\mathbf{h}$ and $\mathbf{w} = \mathbf{h} - \mathbf{A}\mathbf{x}$ giving reciprocally

$$ \mathbf{x} = (\mathbf{I} + \epsilon\mathbf{B}\mathbf{A})^{-1}(\mathbf{v} - \epsilon\mathbf{B}\mathbf{w}) \approx (\mathbf{I} - \epsilon\mathbf{B}\mathbf{A})(\mathbf{v} - \epsilon\mathbf{B}\mathbf{w}) \quad (\text{D.4a}) $$

$$ \mathbf{h} = (1 + \epsilon\mathbf{A}\mathbf{B})^{-1}(\mathbf{w} + \epsilon\mathbf{A}\mathbf{v}) \approx (1 - \epsilon\mathbf{A}\mathbf{B})(\mathbf{w} + \epsilon\mathbf{A}\mathbf{v}) \quad (\text{D.4b}) $$
where **I** is the identity matrix of size *N*.

Using Eq. (D.4), Eq. (D.3) is transformed into the following form

$$
\begin{align}
\ddot{\mathbf{v}} + \mathbf{C}\dot{\mathbf{v}} + \mathbf{K}\mathbf{v} - \epsilon \left[ \mathbf{C}\mathbf{B}(\mathbf{A}\dot{\mathbf{v}} + \dot{\mathbf{w}}) + \mathbf{K}\mathbf{B}(\mathbf{A}\mathbf{v} + \mathbf{w}) \right] &= \mathbf{0} \tag{D.5a} \\
\ddot{\mathbf{w}} + \eta_h \dot{\mathbf{w}} + k_h^{\mathrm{NL}} \mathbf{w}^3 - \mathbf{A}\mathbf{C}\dot{\mathbf{v}} - \mathbf{A}\mathbf{K}\mathbf{v} + \nonumber
\end{align}
$$

$$
\epsilon \left[ \mathbf{A}\mathbf{B}\eta_h \dot{\mathbf{w}} + \mathbf{A}\mathbf{C}\mathbf{B}(\mathbf{A}\dot{\mathbf{v}} + \dot{\mathbf{w}}) + \mathbf{A}\mathbf{K}\mathbf{B}(\mathbf{A}\mathbf{v} + \mathbf{w}) + \mathbf{A}\mathbf{B}k_h^{\mathrm{NL}}\mathbf{w}^3 \right] = 0. \quad (\mathrm{D.5b})
$$

where only terms of order lower or equal to 1 in $\epsilon$ were kept.

To retain only the essential features of the single instability, the reduction of the dynamics is performed on Eq. (D.5a). To achieve this, the latter, is written in state-space form as follows

$$
\dot{\mathbf{y}} = D\mathbf{y} + \epsilon [\mathbf{D}_1\mathbf{y} + \mathbf{D}_1\mathbf{w} + \mathbf{D}_3\dot{\mathbf{w}}] \quad (\text{D.6})
$$

where $\mathbf{y} = (v_1, \dots, v_N, \dot{v}_1, \dots, \dot{v}_N)^T$,

$$
\mathbf{D} = \left( \begin{array}{c|c} \mathbf{0} & \mathbf{I}_N \\ -\mathbf{K} & -\mathbf{C} \end{array} \right), \quad \mathbf{D}_1 = \left( \begin{array}{c|c} \mathbf{0} & \mathbf{0} \\ \mathbf{K}\mathbf{B}\mathbf{A} & \mathbf{C}\mathbf{B}\mathbf{A} \end{array} \right), \quad \mathbf{D}_2 = \left( \begin{array}{c} \mathbf{0} \\ \mathbf{K}\mathbf{B} \end{array} \right) \quad \text{and} \quad \mathbf{D}_3 = \left( \begin{array}{c} \mathbf{0} \\ \mathbf{C}\mathbf{B} \end{array} \right). \qquad (\text{D.7})
$$

Because the matrix **D**, which characterizes the linear dynamics of the primary system, is not symmetric, its
diagonalization must be performed using biorthogonality property of the right eigenvectors **r**<sub>*i*</sub> (*i* = 1, ..., 2*N*)
and the left eigenvectors **l**<sub>*j*</sub> (*j* = 1, ..., 2*N*) of **D**. When the eigenvectors are properly normalized, the bior-
thogonality properties consist in **L**<sup>*T*</sup>**D****R** = **Λ** where **R** = [**r**<sub>1</sub> **r**<sub>1</sub><sup>−</sup> ... **r**<sub>N</sub> **r**<sub>N</sub><sup>−</sup>], **L** = [**l**<sub>1</sub> **l**<sub>1</sub><sup>−</sup> ... **l**<sub>N</sub> **l**<sub>N</sub><sup>−</sup>] and
**Λ** = **diag**(*λ*<sub>1</sub>, **λ**<sub>1</sub><sup>−</sup>, ..., *λ*<sub>N</sub>, **λ**<sub>N</sub><sup>−</sup>) are the right and left eigenvector matrices and the eigenvalue diagonal matrix
respectively.

The so-called biorthogonal coordinates are then introduced. They are constituted of *N* pairs of complex
conjugates, *q*<sub>*n*</sub> and *q*<sub>*n*</sub><sup>−</sup> (*n* = 1, ..., *N*), and defined by the following relations

$$
\mathbf{y} = \mathbf{R}\mathbf{q} \quad \text{and} \quad \mathbf{q} = \mathbf{L}^T\mathbf{y}, \tag{D.8}
$$

where $\mathbf{q} = (q_1, \bar{q}_1, \dots, q_N, \bar{q}_N)^T$.

Substituting Eq. (D.8) into Eq. (D.6), the equations of motion take the form of the following system

$$
\dot{\mathbf{q}} = \mathbf{\Lambda}\mathbf{q} + \epsilon\mathbf{L}^T[\mathbf{D}_1\mathbf{R}\mathbf{q} + \mathbf{D}_2\mathbf{w} + \mathbf{D}_3\dot{\mathbf{w}}], \quad (\text{D.9})
$$

and therefore Eq. (D.5) is equivalent to

$$
\begin{align}
\dot{\mathbf{q}} - \mathbf{\Lambda}\mathbf{q} - \epsilon\mathbf{L}^T[\mathbf{D}_1\mathbf{R}\mathbf{q} + \mathbf{D}_2\mathbf{w} + \mathbf{D}_3\dot{\mathbf{w}}] &= 0 \tag{D.10a} \\
\ddot{\mathbf{w}} + \eta_h \dot{\mathbf{w}} + k_h^{\mathrm{NL}} \mathbf{w}^3 &- \nonumber \\
\mathbf{A}\mathbf{C}(\mathbf{R}^{\mathrm{DL}}\mathbf{q}^\mathrm{U} + \mathbf{R}^{\mathrm{DR}}\mathbf{q}^\mathrm{D}) - \mathbf{A}\mathbf{K}(\mathbf{R}^{\mathrm{UL}}\mathbf{q}^\mathrm{U} + \mathbf{R}^{\mathrm{UR}}\mathbf{q}^\mathrm{D}) + \nonumber \\
\epsilon[\mathbf{A}\mathbf{B}\eta_h \dot{\mathbf{w}} + \mathbf{A}\mathbf{C}\mathbf{B}(\mathbf{A}(\mathbf{R}^{\mathrm{DL}}\mathbf{q}^\mathrm{U} + \mathbf{R}^{\mathrm{DR}}\mathbf{q}^\mathrm{D}) + \dot{\mathbf{w}}) + \nonumber \\
\mathbf{A}\mathbf{K}\mathbf{B}(\mathbf{A}(\mathbf{R}^{\mathrm{UL}}\mathbf{q}^\mathrm{U} + \mathbf{R}^{\mathrm{UR}}\mathbf{q}^\mathrm{D}) + \mathbf{w}) + \mathbf{A}\mathbf{B}k_h^{\mathrm{NL}}\mathbf{w}^3] &= 0, \tag{D.10b}
\end{align}
$$

where the matrix **R** and the vector **q** have been split into a *N* × *N*-block matrix and a *N* × 1-block respectively
as

$$
\mathbf{R} = \begin{bmatrix} \mathbf{R}^{\mathrm{UL}} & \mathbf{R}^{\mathrm{UR}} \\ \mathbf{R}^{\mathrm{DL}} & \mathbf{R}^{\mathrm{DR}} \end{bmatrix} \quad \text{and} \quad \mathbf{q} = ((\mathbf{q}^{\mathrm{U}})^{\mathrm{T}}, (\mathbf{q}^{\mathrm{D}})^{\mathrm{T}})^{\mathrm{T}} \qquad (\mathrm{D}.11)
$$

where “U”, “D”, “L” and “R” in previous superscripts must be understood as the initial letters of up, down,
left and right respectively.

As previously mentioned, one mode of the primary system can become unstable through a Hopf bifurcation. Without loss of generality, the first mode is chosen to be the potential unstable mode. That means that when the chosen physical bifurcation parameter, denoted $\sigma$, crosses the particular parameter value $\sigma_{\text{Hopf}}$, called *Hopf bifurcation point*, a pair of complex conjugate eigenvalues crosses the complex plane imaginary axis (without loss of generality this is first pair $(\lambda_1, \bar{\lambda}_1)$) whereas all the other eigenvalues $(\lambda_n, \bar{\lambda}_n)$ (with $n = 2, \dots, N$) have strictly negative real parts.

Neglecting the coupling term $\epsilon \mathbf{L}^T [\mathbf{D}_1 \mathbf{R} \mathbf{q} + \mathbf{D}_2 \mathbf{w} + \mathbf{D}_3 \dot{\mathbf{w}}]$, Eq. (D.10a) reduces to the diagonalized system $\dot{\mathbf{q}} - \mathbf{\Lambda} \mathbf{q} = \mathbf{0}$. Therefore, in this case, the variables $q_n$ and $\bar{q}_n$ ($n = 2, \dots, N$), related to the stable modes, tend to zero. When the coupling term (which is of the order of magnitude of $\epsilon$) is taken into account, the variables
$q_n$ and $\bar{q}_n$ do not vanish but their contributions are very small and can be neglected [a8, a9, 78]. Consequently, all terms related to $q_n$ and $\bar{q}_n$ (with $n = 2, \dots, N$) are omitted from further consideration and Eq. (D.9) is reduced to

$$ \dot{q}_1 = \lambda_1 q_1 + \epsilon l_1^T (\mathbf{D}_1(\mathbf{r}_1 q_1 + \bar{\mathbf{r}}_1 \bar{q}_1) + \mathbf{D}_2 w + \mathbf{D}_3 \dot{w}). \quad (D.12) $$

Finally, grouping Eq.s (D.10b) and (D.12), splitting $\lambda_1$ into its real and imaginary parts as $\lambda_1 = \lambda_R + j\lambda_I$ (with $j^2 = -1$) and switching the time scale as $t \to \lambda_I t$ yield

$$ \dot{q}_1 - \left( \frac{\lambda_R}{\lambda_I} + j \right) q_1 - \frac{\epsilon}{\lambda_I} l_1^T (\mathbf{D}_1(\mathbf{r}_1 q_1 + \bar{\mathbf{r}}_1 \bar{q}_1) + \mathbf{D}_2 w + \mathbf{D}_3 \dot{w}) = 0 \quad (D.13a) $$

$$ \ddot{w} + \frac{\eta_h}{\lambda_I} \dot{w} + \frac{k_h^{NL}}{\lambda_I^2} w^3 - $$

$$ \frac{1}{\lambda_I^2} (\mathbf{A}\mathbf{C}(\mathbf{r}_1^{\text{DL}} q_1 + \bar{\mathbf{r}}_1^{\text{DL}} \bar{q}_1) + \mathbf{A}\mathbf{K}(\mathbf{r}_1^{\text{UL}} q_1 + \bar{\mathbf{r}}_1^{\text{UL}} \bar{q}_1)) + $$

$$ \epsilon \left[ \mathbf{A}\mathbf{B} \frac{\eta_h}{\lambda_I} \dot{w} + \frac{1}{\lambda_I^2} \mathbf{A}\mathbf{C}\mathbf{B} \left( \mathbf{A}(\mathbf{r}_1^{\text{DL}} q_1 + \bar{\mathbf{r}}_1^{\text{DL}} \bar{q}_1) + \lambda_I \dot{w} \right) + \right. $$

$$ \left. \frac{1}{\lambda_I^2} \mathbf{A}\mathbf{K}\mathbf{B} \left( \mathbf{A}(\mathbf{r}_1^{\text{UL}} q_1 + \bar{\mathbf{r}}_1^{\text{UL}} \bar{q}_1) + w \right) + \mathbf{A}\mathbf{B} \frac{k_h^{NL}}{\lambda_I^2} w^3 \right] = 0. \quad (D.13b) $$

where the vectors $\mathbf{r}_1^{\text{UL}}$ and $\mathbf{r}_1^{\text{DL}}$ correspond to the first columns of the matrices $\mathbf{R}^{\text{DL}}$ and $\mathbf{R}^{\text{UL}}$ respectively (see Eq. (D.11)).

The system of equations (D.13) is the system under study in the following mathematical developments.

## D.3 ASYMPTOTIC ANALYSIS: THE ZEROth-ORDER APPROXIMATION OF THE SLOW FLOW

The present section summarizes and extends previous results of the literature [a7, a8, 78]. In particular, an analytical expression of the mitigation limit (hereafter called the zeroth-order approximation of the mitigation limit) is derived at the end of the section under several assumptions. The first is given just below.

In Eq. (D.13), the quantities $\lambda_R$, $\lambda_I$, the vectors $l_1^T$, $\mathbf{r}_1$, $\mathbf{r}_1^{\text{UL}}$ and $\mathbf{r}_1^{\text{DL}}$ and the matrices $\mathbf{C}$, $\mathbf{K}$ (and those which depend on them) depend in general on the bifurcation parameter under consideration $\sigma$. In this work one assumes that the NES effects are not too far from the Hopf bifurcation point $\sigma_{\text{Hopf}}$ and therefore the latter are evaluated at $\sigma_{\text{Hopf}}$ (using the notation $\mathbf{r}_{\text{Hopf}}$) expect $\lambda_R$ which is by definition equal to zero at $\sigma_{\text{Hopf}}$. $\lambda_R$ is evaluated at the actual value of $\sigma$ and the dimensionless eigenvalue real part $\lambda_R/\omega$ (where $\omega = \lambda_{I,\text{Hopf}}$, the imaginary part of $\lambda_1$ evaluated at $\sigma_{\text{Hopf}}$) is supposed to be small (i.e. of the order of $\mathcal{O}(\epsilon)$). The parameter $\rho$ is then finally introduced as

$$ \rho = \frac{\lambda_R}{\epsilon\omega} \quad (D.14) $$

where $\rho$ is of the order of $\mathcal{O}(1)$.

The rescaled dimensionless eigenvalue real part $\rho$ of the unstable mode is considered as a generalized bifurcation parameter (that results from the physical bifurcation parameter $\sigma$ of the initial mechanical system) from which one will follow the changes in the system dynamic behavior.

### D.3.1 Governing equations of the slow flow

It is known [76, 203] that TET appears in the vicinity of the 1:1 resonance between two nonlinear modes of the coupled structure, the latter resulting from the interaction between a linear mode of the primary structure (here this is the considered unstable mode) and the nonlinear mode of the NES. This phenomenon is called 1:1 resonance capture and occurs at a frequency close to the natural frequency of the linear primary structure. To investigate the solution in the neighborhood of this 1:1 resonance capture the complexification-averaging method [135, 206] is used. The resulting complex averaged dynamics is called slow flow.

The first step consists in introducing the following complex variables

$$ \psi = \dot{w} + jw \quad (D.15) $$

or equivalently

$$ w = \frac{\psi - \bar{\psi}}{2j}, \quad \dot{w} = \frac{\psi + \bar{\psi}}{2} \quad \text{and} \quad \ddot{w} = \dot{\psi} - \frac{j}{2}(\psi + \bar{\psi}) \quad (D.16) $$
and expressing the complex variables $q_1$ and $\psi$ as

$$q_1 = \tilde{\phi}_1 e^{j t} \quad \text{and} \quad \psi = \phi_2 e^{j t}. \qquad (D.17)$$

Substituting (D.17) into (D.13) and averaging over the period $T = 2\pi$ leads to the following system of complex slow modulated amplitudes

$$\dot{\tilde{\phi}}_1 = \epsilon \left[ (\rho + l_{1, \text{Hopf}}^T \mathbf{D}_1, \text{Hopf} \mathbf{r}_1, \text{Hopf}) \tilde{\phi}_1 + \frac{1}{2} (-j l_{1, \text{Hopf}}^T \mathbf{D}_2, \text{Hopf} + l_{1, \text{Hopf}}^T \mathbf{D}_3, \text{Hopf}) \phi_2 \right] \quad (D.18a)$$

$$\begin{aligned} \dot{\phi}_2 = & \frac{\mathbf{A}\mathbf{K}_{\text{Hopf}}\mathbf{r}_{1, \text{Hopf}}^{\text{UL}} + \mathbf{A}\mathbf{C}_{\text{Hopf}}\mathbf{r}_{1, \text{Hopf}}^{\text{DL}}}{\omega^2} \tilde{\phi}_1 - \frac{1}{2} \left( j + \mu - j \frac{3\alpha}{4} |\phi_2|^2 \right) \phi_2 \\ & - \epsilon \left[ \frac{1}{\omega^2} \left( \mathbf{A}\mathbf{K}_{\text{Hopf}}\mathbf{B}\mathbf{Tr}_{1, \text{Hopf}}^{\text{UL}} + \frac{1}{2}\mathbf{A}\mathbf{C}_{\text{Hopf}}\mathbf{B}\mathbf{Tr}_{1, \text{Hopf}}^{\text{DL}} \right) \tilde{\phi}_1 \right. \\ & \left. + \frac{1}{2} \left( \frac{\mathbf{A}\mathbf{C}_{\text{Hopf}}\mathbf{B}}{\omega} - j \frac{\mathbf{A}\mathbf{K}_{\text{Hopf}}\mathbf{B}}{\omega^2} + \mathbf{A}\mathbf{B} \left( \mu - j \frac{3\alpha}{4} |\phi_2|^2 \right) \right) \phi_2 \right] \end{aligned} \qquad (D.18b)$$

with

$$\mu = \frac{\eta_h}{\omega} \quad \text{and} \quad \alpha = \frac{k_h^{NL}}{\omega^2} \qquad (D.19)$$

and recalling that the notation $()_{\text{Hopf}}$ means that the considered quantity is evaluated at the Hopf bifurcation point $\sigma_{\text{Hopf}}$.

Introducing the variable $\phi_1 = 2(\mathbf{A}\mathbf{K}_{\text{Hopf}}\mathbf{r}_{1, \text{Hopf}}^{\text{UL}} + \mathbf{A}\mathbf{C}_{\text{Hopf}}\mathbf{r}_{1, \text{Hopf}}^{\text{DL}})\tilde{\phi}_1/\omega^2$, and ignoring the terms of the order of magnitude of $\epsilon$ in (D.18b)<fnref n="3" />, Eq. (D.18) reduces to the following synthetic form

$$\dot{\phi}_1 + \epsilon [(a - \rho) \phi_1 + b \phi_2] = 0 \qquad (D.20a)$$

$$\dot{\phi}_2 + \frac{1}{2}\phi_1 + \frac{1}{2}(\mu + j)\phi_2 + j\frac{3}{8}\alpha\phi_2|\phi_2|^2 = 0. \qquad (D.20b)$$

with the complex coefficients $a$ and $b$ defined as

$$\begin{aligned} a &= -l_{1, \text{Hopf}}^T \mathbf{D}_1, \text{Hopf} \mathbf{r}_1, \text{Hopf} \\ b &= \frac{1}{\omega^2} (-j l_{1, \text{Hopf}}^T \mathbf{D}_2, \text{Hopf} + l_{1, \text{Hopf}}^T \mathbf{D}_3, \text{Hopf}) (\mathbf{A}\mathbf{K}_{\text{Hopf}}\mathbf{r}_{1, \text{Hopf}}^{\text{UL}} + \mathbf{A}\mathbf{C}\mathbf{r}_{1, \text{Hopf}}^{\text{DL}}). \end{aligned} \qquad (D.21)$$

Then, substituting $\phi_1$ and $\phi_2$ in Eq. (D.18) by their polar coordinates defined as $\phi_1 = re^{j\theta_1}$ and $\phi_2 = se^{j\theta_2}$, new equations of motion for the real amplitudes $r$ and $s$ and phase difference $\Delta = \theta_2 - \theta_1$ can be obtained as

$$\dot{r} = \epsilon f(r, s, \Delta) \qquad (D.22a)$$

$$\dot{s} = \tilde{g}_1(r, s, \Delta, \epsilon) \qquad (D.22b)$$

$$\dot{\Delta} = \tilde{g}_2(r, s, \Delta, \epsilon) \qquad (D.22c)$$

with

$$f(r, s, \Delta) = (\rho - a_R) r - (b_R \cos \Delta - b_I \sin \Delta) s, \qquad (D.23a)$$

$$\tilde{g}_1(r, s, \Delta, \epsilon) = -\frac{1}{2}(\mu s + r \cos \Delta), \qquad (D.23b)$$

$$\tilde{g}_2(r, s, \Delta, \epsilon) = \frac{3}{8}\alpha s^2 - \frac{1}{2} + \frac{r}{2s}\sin\Delta + \epsilon \left[ a_I + b_I \frac{\cos\Delta}{r} + b_R \frac{\sin\Delta}{r} \right] \qquad (D.23c)$$

where the complex parameters $a$ and $b$ were split into real and imaginary parts as $a = a_R + ja_I$ and $b = b_R + jb_I$.

Because of the presence of the small parameter $\epsilon$, the slow flow is governed by different time scales. More precisely, within the framework of the geometric singular perturbation theory [107], Eq. (D.22) appears as a (2, 1)-fast-slow system where $s$ and $\Delta$ are the fast variables and $r$ the slow variable. The latter is related to the modal component $q_1$ whereas the fast variables are associated with the relative displacement $w$ between the NES and the primary structure.

<footnote n="3">3. These terms do not affect qualitatively the dynamic behavior of the slow flow.</footnote>
### D.3.2 Critical manifold and fold points of the slow flow

The real form of the slow flow, i.e. Eq. (D.22), is written with respect to the dimensionless slow time $τ = ϵt$ as follows

$$r' = f(r, s, \Delta) \qquad (D.24a)$$

$$\epsilon s' = \tilde{g}_1(r, s, \Delta, \epsilon) \qquad (D.24b)$$

$$\epsilon \Delta' = \tilde{g}_2(r, s, \Delta, \epsilon) \qquad (D.24c)$$

where $(\cdot)'$ denotes the derivative with respect to the slow time $τ$. Considering $ϵ = 0$ yields the slow subsystem

$$r' = f(r, s, \Delta) \qquad (D.25a)$$

$$0 = \tilde{g}_1(r, s, \Delta, 0) \qquad (D.25b)$$

$$0 = \tilde{g}_2(r, s, \Delta, 0) \qquad (D.25c)$$

which is a differential-algebraic equation, and the fast subsystem

$$\dot{r} = 0 \qquad (D.26a)$$

$$\dot{s} = \tilde{g}_1(r, s, \Delta, 0) \qquad (D.26b)$$

$$\dot{\Delta} = \tilde{g}_2(r, s, \Delta, 0) \qquad (D.26c)$$

The critical manifold of the slow flow is the solution of the algebraic part of Eq. (E.30) and it is expressed as follows

$$\mathcal{M}_0 = \{ (r, s, \Delta) \in \mathbb{R}^{+2} \times [-\pi, \pi] \mid \tilde{g}_1(r, s, \Delta, 0) = 0 \text{ and } \tilde{g}_2(r, s, \Delta, 0) = 0 \}. \qquad (D.27)$$

Combining Eqs. (E.30b) and (E.30c) leads to the following modulus and argument equations

$$r = s \sqrt{\mu^2 + \left(1 - \frac{3\alpha s^2}{4}\right)^2} = H(s) \qquad (D.28a)$$

$$\tan \Delta = \frac{3\alpha s^2 - 4}{4\mu}. \qquad (D.28b)$$

Stability analysis of the critical manifold is now carried out considering the fast subsystem (E.31). It can be shown that a fixed point of (E.31) is stable if $d_s H(s) > 0$ and unstable if $d_s H(s) < 0$. Hence, the subset of $\mathcal{M}_0$ satisfying $d_s H(s) < 0$ defines the attracting zone for the fast dynamics whereas the subset of $\mathcal{M}_0$ satisfying $d_s H(s) > 0$ define the repelling zone.

Exploiting the polynomial properties of $H$, it can be shown that the local extrema of $H$ (i.e $d_s H(s) = 0$) occur at

$$s^{\text{LF}} = \frac{2}{3\sqrt{\alpha}} \sqrt{2 - \sqrt{1 - 3\mu^2}} \quad \text{and} \quad s^{\text{RF}} = \frac{2}{3\sqrt{\alpha}} \sqrt{2 + \sqrt{1 - 3\mu^2}} \qquad (D.29)$$

if the following relation holds

$$\mu < \frac{1}{\sqrt{3}}. \qquad (D.30)$$

In Eq. (E.49), $s^{\text{LF}}$ and $s^{\text{RF}}$ are the abscissa values in the $(s, r)$-plane of a maximum and the minimum of the function $H$ respectively. The superscripts ()$^{\text{LF}}$ and ()$^{\text{RF}}$ refer to left fold point and right fold point respectively. Indeed, in the $(s, r, \Delta)$-space, the two points $(s^{\text{RF}}, r^{\text{RF}}, \Delta^{\text{RF}})$ and $(s^{\text{LF}}, r^{\text{LF}}, \Delta^{\text{LF}})$ (where $r^{\text{LF}}, r^{\text{RF}}, \Delta^{\text{LF}}$ and $\Delta^{\text{RF}}$ can be deduced from $s^{\text{LF}}$ and $s^{\text{RF}}$ using Eq. (E.48)) are generally called fold points.

When the condition (E.50) is not satisfied the function $H$ no longer has local extrema. In the remaining of the paper, one considers situations in which (E.50) holds.

A typical critical manifold (more precisely its real part defined by Eq. (E.48a)) is depicted in Fig. E.9. One can see that the fold points connect attracting parts to repelling part of the critical manifold $\mathcal{M}_0$. The two scalars $s^\text{D}$ and $s^\text{U}$, which are the horizontal projection of the fold points on the critical manifold, are defined by $H(s^{\text{RF}}) = H(s^\text{D})$ and $H(s^{\text{LF}}) = H(s^\text{U})$ giving

$$s^\text{D} = \frac{2\sqrt{2}}{3\sqrt{\alpha}} \sqrt{1 - \sqrt{1 - 3\mu^2}} \quad \text{and} \quad s^\text{U} = \frac{2\sqrt{2}}{3\sqrt{\alpha}} \sqrt{1 + \sqrt{1 - 3\mu^2}}. \qquad (D.31)$$
<figure><img src="image_98.png" /><figcaption>FIGURE D.1. Typical example of the critical manifold in the (s, r)-plane given by Eq. (E.48a) for $\mu = 0.25$ and $\alpha = 5$.</figcaption></figure>

### D.3.3 Fixed points and fold singularities of the slow flow

Because $0 < \epsilon \ll 1$, the fixed points of the slow flow (D.22), which characterize periodic solutions of (D.13), can be approximated by those of the slow subsystem (E.30). From that zeroth-order approximation an analytical expression of these fixed points can be obtained. Note that within this approximation, the fixed points are necessarily located on the critical manifold.

First, substituting Eq. (E.48a) into Eq. (E.30a) results in $(H(s))' = f(H(s), s, \Delta)$ which reduces to

$$d_s H(s)s' = f_s(s), \qquad (D.32)$$

where the explicit expression of the function $f_s$ is not given here.

Then, a (regular) fixed point of Eq. (D.32) is defined as the roots of $f_s(s) = 0$ with $d_s H(s) \neq 0$. Solving $f_s(s) = 0$ gives two fixed points (in addition to the trivial solution) whose analytical expressions are

$$s_1^e = \left( \frac{2}{3\alpha(a_R - \rho)} \left( 2a_R - b_I - 2\rho + \sqrt{b_I^2 - 4\mu(a_R - \rho)(\mu a_R - b_R - \mu\rho)} \right) \right)^{1/2}, \qquad (D.33a)$$

$$s_2^e = \left( \frac{2}{3\alpha(a_R - \rho)} \left( 2a_R - b_I - 2\rho - \sqrt{b_I^2 - 4\mu(a_R - \rho)(\mu a_R - b_R - \mu\rho)} \right) \right)^{1/2}, \qquad (D.33b)$$

with $s_1^e < s_2^e$. As usual, the stability of the fixed points are found looking for the sign of $d_s \left( \frac{f_s(s)}{H'(s)} \right)$.

Both fixed points $s_1^e$ and $s_2^e$ no longer exist when they become complex for

$$\rho^S = \frac{2\mu a_R + \sqrt{b_I^2 + b_R^2} - b_R}{2\mu} \qquad (D.34)$$

where the superscript $^S$ has no particular meaning here, it is just there to give a name to this particular value of $\rho$.

A singular fixed points of Eq. (D.32), also called folded singularites, are defined as the roots of the following nonlinear equations $f_s(s) = 0$ and $d_s H(s) = 0$. Folded singularities are hints of particular solutions of such fast-slow systems, called canard cycles [16]. These solutions are not investigated in this paper.

### D.3.4 The possible responses

Previous analysis provides a qualitative description of the slow flow dynamics in the limit case where $\epsilon = 0$, i.e. within a zeroth-order approximation. Within this approximation, the slow flow evolves at two time scales : a slow time scale in which the slow flow is on the critical manifold and described by the slow subsystem (E.30) and a fast time scale in which the slow flow is outside the critical manifold and described by
the fast subsystem (E.31). The particular S-shape of the critical manifold together with the stability analysis of the fixed points of the slow flow allows us to explain and predict its different responses and consequently those of the initial full order system (D.13). As has been widely discussed in the literature (see e.g. [75, 78]), four scenarios are possible. In previous works by the authors [a7, a8] these responses are classified into two categories. In the first category, the NES acts, resulting in the three responses called *harmless situations*. First, the *Complete suppression* in which the trivial fixed point of the slow flow (D.22) is stable and then reached. In this case, the system is stabilized because of the additional part of the NES including mass and damping. In general, because a light mass and a small damping are considered for the NES, this stabilization effect is negligible, in term of bifurcation parameter range, compared to the mitigation effects described hereinafter. In other words, the Hopf bifurcation point of the primary structure without NES is usually very slightly smaller than the Hopf bifurcation point of the coupled system. Of course, in a context of nonlinear vibrations absorption, the complete suppression is not the desired effect. Then, one can observe *mitigation through periodic response*. In this case, a nontrivial stable fixed point of the slow flow (D.22) is reached, that corresponds to a periodic regime for the initial system (D.13). The last harmless situations correspond to *mitigation through Strongly Modulated Responses (SMRs)*. SMRs correspond to relaxations oscillations of the slow flow (D.22), that corresponds to a quasi-periodic (amplitude and phase modulated) regime for the initial system (D.13). In the second category, called *harmful situation*, the NES is not able to produce small amplitude responses. For the initial nonlinear system one observes in this case a limit cycle with an amplitude close to that of the system without NES. Note that if the nonlinearity of the primary structure is neglected in the governing equations of the slow flow (this the case in the present work), the harmful situation corresponds for it to an unbounded regime in which its trajectory growths to infinity.

In this paper, we are particularly interested in the transition from SMR to no mitigation which is also, in general, the transition from harmless to harmful situations. Therefore one recalls here briefly what relaxation oscillations of the slow flow are. Note that the following description holds in the limit case for which $\epsilon = 0$. In the $(s, r)$-plane, from an initial condition near zero and outside the critical manifold, the trajectory of the slow flow evolves rapidly and horizontally to the left attracting branch of $M_0$ (this fast epoch is described by the fast subsystem (E.31)). Then, the slow flow evolves slowly on this branch (this slow epoch is described by the slow subsystem (E.30)) and, if there is no stable fixed points on it, the trajectory reached the left fold point $(s^{\text{LF}}, r^{\text{LF}})$ at which $M_0$ becomes repelling. From the fold point, the slow flow undergoes a fast horizontal jump to $(s^{\text{U}}, r^{\text{LF}})$ on the right attracting part of $M_0$ and then a second slow epoch to the right fold point $(s^{\text{RF}}, r^{\text{RF}})$. A second horizontal jump occurs from $(s^{\text{RF}}, r^{\text{RF}})$ to $(s^{\text{D}}, r^{\text{RF}})$ and the trajectory returns to the left attracting part of the $M_0$. A third slow epoch occurs on $M_0$ to $(s^{\text{LF}}, r^{\text{LF}})$ and so on.

### D.3.5 Zeroth-order analytical prediction of the mitigation limit

The mitigation limit is defined below.

**DEFINITION D.1.** Considering a set of initial conditions (for the slow flow) as a small perturbation of the trivial solution, the **mitigation limit** is defined as the value of the generalized bifurcation parameter $\rho$ which separates harmless situations from harmful situation.

As discussed in the literature [a7, 75, 78], the mitigation limit (denoted here $\rho_0^*$) is defined, within the zeroth-order approximation, as the value of $\rho$ for which the fixed point $s_2^e$ becomes smaller than $s^U$ defined by Eq. (E.53) or for which $s_2^e$ no longer exists (see Eq. (D.34)). Indeed, when the largest unstable fixed exists (with, in addition, $s_2^e > s^U$), it prevents the phase trajectory to grow to infinity on the right branch of the critical manifold.

From Eqs. (E.53) and (D.33a) the solution of $s_2^e = s^U$ is derived as

$$ \rho_0^U = a_R + \frac{2\sqrt{1-3\mu^2}b_I + b_I - (2\sqrt{1-3\mu^2}+7)\mu^2 b_I}{(\mu^2+1)^2} + \frac{(3\mu^2-4\sqrt{1-3\mu^2}-5)\mu b_R}{(\mu^2+1)^2}. \quad (D.35) $$

Then using (D.34), the following conditional expression for $\rho_0^*$ is obtained:

$$ \rho_0^* = \begin{cases} \rho_0^U \text{ (Eq. (D.35))}, & \text{if } \mu < \mu^* \\ \rho_0^S \text{ (Eq. (D.34))}, & \text{if } \mu > \mu^*, \end{cases} \quad (D.36a) $$

with $\mu^*$ the special value of $\mu$ solution of $\rho_0^U = \rho^S$ which is in general smaller than $1/\sqrt{3}$.

We can notice that within the zeroth-order approximation, the mitigation limit does not depend on the nonlinear parameter $\alpha$ of the NES but only its damping coefficient $\mu$ (together with the parameters of the
primary system). Therefore, for given parameters of the primary structure, it exists an optimal value of $\mu$, denoted $\mu_0^{\text{opt}}$, which maximizes the mitigation limit. Consequently, assuming that $\mu_0^{\text{opt}} < \mu^*$, $\mu_0^{\text{opt}}$ is the solution of $\partial_{\mu}\rho_0^U = 0$. Solving the latter equation yields

$$ \mu_0^{\text{opt}} = -\frac{1}{\mathcal{B}} \left[ b_{\mathcal{R}} b_{\mathcal{I}} + \sqrt{\mathcal{B} (|b|^4 b_{\mathcal{R}}^2)^{1/3} - 3 b_{\mathcal{R}}^2 (b_{\mathcal{I}}^2 + b_{\mathcal{R}}^2)} \right. \\ \left. - \left( -\mathcal{B} (|b|^4 b_{\mathcal{R}}^2)^{1/3} - 2 \mathcal{B} b_{\mathcal{R}}^2 + 2 b_{\mathcal{I}}^2 b_{\mathcal{R}}^2 - \frac{2 b_{\mathcal{I}} b_{\mathcal{R}} (b_{\mathcal{I}}^2 + b_{\mathcal{R}}^2) (8 b_{\mathcal{I}}^2 + 9 b_{\mathcal{R}}^2)}{\sqrt{\mathcal{B} (|b|^4 b_{\mathcal{R}}^2)^{1/3} - 3 b_{\mathcal{R}}^2 (b_{\mathcal{I}}^2 + b_{\mathcal{R}}^2)}} \right)^{1/2} \right] \quad (\text{D.37}) $$

with $\mathcal{B} = 3b_{\mathcal{R}}^2 + 4b_{\mathcal{I}}^2$ and where the subscript 0 is used to highlight that the expression (D.37) is obtained within the zeroth-order approximation.

## D.4 SCALING LAW FOR THE SLOW FLOW NEAR A FOLD POINT

The Fenichel theory [62] guarantees the persistence of the relaxation oscillations scenario, described at the end of Section. D.3.4, for $0 < \epsilon \ll 1$ by stating that Eq. (E.29) has an invariant manifold $M_\epsilon$ in the $O(\epsilon)$-vicinity of $M_0$ and with the same stability properties with respect to the fast variables as $M_0$ (attracting or repelling). In this case, concerning the first jump described previously, the jump point at which the trajectory leaves the attracting branch of $M_\epsilon$ to undergo the jump is no longer the fold point ($s^{\text{LF}}, r^{\text{LF}}$). The same is true for the arrival point at which the trajectory reaches the right attracting part of $M_\epsilon$ after the jump which is no longer ($s^U, r^{\text{LF}}$). The goal of this section is to determine these jump and arrival points and consequently to obtain an analytical prediction of the mitigation limit which, unlike Eq. (D.36), takes into account $\epsilon$.

To achieve that, one must determine beforehand the scaling law of the slow flow dynamics in the neighborhood of the fold point ($s^{\text{LF}}, r^{\text{LF}}$), i.e. a law which describes the $\epsilon$-dependence of the distance between $M_0$ and the actual trajectory of the slow flow for $0 < \epsilon \ll 1$. The jump and arrival points are then deduced from this scaling law.

The Fenichel theorem is true for normally hyperbolic branch of $M_0$ (i.e. outside the fold points) and fails to approach the actual trajectory of the system at fold points. Indeed, the literature on dynamical systems shows that the trajectory still passes through the region near these fold points but with a nontrivial scaling behavior involving the fractional exponents $1/3$ and $2/3$ for the perturbation parameter $\epsilon$ (see e.g. [40] Sect. 5.4).

### D.4.1 Center manifold reduction of the slow flow at the left fold point

The center manifold reduction technique (as presented for example in [27, 88]) is used in this section to reduce the slow flow (E.29), in the vicinity of the left fold point ($r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}}$), to the normal form of the dynamic saddle-node bifurcation which is then solved to deduce the nontrivial scaling law previously mentioned.

First, assuming $\epsilon = 0$ in the right-hand sides of Eqs. (E.29b) and (E.29c), one considers the following simplification of the slow flow

$$ r' = f(r, s, \Delta) \quad (\text{D.38a}) $$

$$ \epsilon s' = \tilde{g}_1(r, s, \Delta, 0) = g_1(r, s, \Delta), \quad (\text{D.38b}) $$

$$ \epsilon \Delta' = \tilde{g}_1(r, s, \Delta, 0) = g_2(r, s, \Delta). \quad (\text{D.38c}) $$

The goal in using the center manifold theorem is to reduce the dynamics with respect to the fast variables $s$ and $\Delta$, i.e. obtain an equivalent system which only depends on one variable. To achieve that, the Jacobian matrix $\mathbf{J}_g$ of (D.38b-D.38c) evaluated at the fold point ($r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}}$) is computed as

$$ \mathbf{J}_g(r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}}) = \begin{pmatrix} -\frac{\mu}{2} & \frac{1}{9} \left(\sqrt{1-3\mu^2}+1\right) \sqrt{\frac{2-\sqrt{1-3\mu^2}}{\alpha}} \\ -\frac{3\left(\sqrt{1-3\mu^2}-1\right)}{4\sqrt{\frac{2-\sqrt{1-3\mu^2}}{\alpha}}} & -\frac{\mu}{2} \end{pmatrix}. \quad (\text{D.39}) $$

which has two real eigenvalues $\lambda_a = 0$ and $\lambda_b = -\mu$.

Then to continue one must obtain polynomial nonlinearities. Indeed, the center manifold reduction requires isolating the linear terms. The functions $g_1$ and $g_2$ are therefore Taylor expanded around the fold point
$(r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}})$, up to order 2 for the variables $s$ and $\Delta$ and up to order 1 for the variable $r$. Remembering that by definition $g_1(r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}}) = 0$ and $g_2(r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}}) = 0$ that yields

$$g_1(r, s, \Delta) \approx (s - s^{\text{LF}}) \partial_s g_1(r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}}) + (\Delta - \Delta^{\text{LF}}) \partial_{\Delta} g_1(r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}}) + \tilde{G}_1(r, s, \Delta) \quad (\text{D.40a})$$

$$g_2(r, s, \Delta) \approx (s - s^{\text{LF}}) \partial_s g_2(r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}}) + (\Delta - \Delta^{\text{LF}}) \partial_{\Delta} g_2(r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}}) + \tilde{G}_2(r, s, \Delta). \quad (\text{D.40b})$$

The polynomial functions $\tilde{G}_1$ and $\tilde{G}_2$ contain nonlinear terms with respect to the $(s - s^{\text{LF}})$ and $(\Delta - \Delta^{\text{LF}})$ (only quadratics terms are taken) and linear terms in $(r - r^{\text{LF}})$.

Denoting $\mathbf{t} = (s, \Delta)^T$ and $\mathbf{t}^{\text{LF}} = (s^{\text{LF}}, \Delta^{\text{LF}})^T$ and using Eq. (D.40), Eqs. (D.38b) and (D.38c) can be written in matrix form as

$$\epsilon \mathbf{t}' = (\mathbf{t} - \mathbf{t}^{\text{LF}}) \mathbf{J}_{\mathbf{g}}(r^{\text{LF}}, \mathbf{t}^{\text{LF}}) + \tilde{\mathbf{G}} (r - r^{\text{LF}}, \mathbf{t} - \mathbf{t}^{\text{LF}}) \quad (\text{D.41})$$

where $\tilde{\mathbf{G}} = (\tilde{G}_1, \tilde{G}_2)^T$. Now, to diagonalize the linear part of Eq. (D.41) the biorthogonal transformation is used again$^4$. For that, the right and left eigenvector matrices $\mathbf{R} = [\mathbf{r}_a \ \mathbf{r}_b]$ and $\mathbf{L} = [\mathbf{l}_a \ \mathbf{l}_b]$ are considered corresponding to the following right and left eigenvalue problems with respect to $\mathbf{J}_{\mathbf{g}}(r^{\text{LF}}, \mathbf{t}^{\text{LF}})$

$$\mathbf{J}_{\mathbf{g}}(r^{\text{LF}}, \mathbf{t}^{\text{LF}})\mathbf{R} = \mathbf{R}\mathbf{\Lambda} \quad \text{and} \quad \mathbf{J}_{\mathbf{g}}(r^{\text{LF}}, \mathbf{t}^{\text{LF}})^T\mathbf{L} = \mathbf{L}\mathbf{\Lambda}. \quad (\text{D.42})$$

where $\mathbf{\Lambda} = \text{diag}(\lambda_a, \lambda_b) = \text{diag}(0, -\mu)$ denotes the diagonal matrix of the eigenvalues. Again the biorthogonal transformation consists of stating

$$\mathbf{t} - \mathbf{t}^{\text{LF}} = \mathbf{R}\mathbf{q} \quad \text{with} \quad \mathbf{q} = \mathbf{L}^T (\mathbf{t} - \mathbf{t}^{\text{LF}}) \quad (\text{D.43})$$

where $\mathbf{q} = (q_a, q_b)^T$. Note that here, on a system undergoing a saddle-node bifurcation, the eigenvalues, the eigenvectors and the modal coordinates are real.

The variable $q_a$ is associated to $\lambda_a = 0$ and the variable $q_b$ is associated to $\lambda_b = -\mu$. Moreover, the matrices $\mathbf{R}$ and $\mathbf{L}$ have been chosen such that $\mathbf{L}^T\mathbf{R} = \mathbf{I}$ and $\mathbf{L}^T\mathbf{J}_{\mathbf{g}}(r^{\text{LF}}, \mathbf{t}^{\text{LF}})\mathbf{R} = \mathbf{\Lambda}$. Eq. (D.41) is then expressed with respect to the variables $q_a$ and $q_b$ as

$$\epsilon \mathbf{q}' = \mathbf{\Lambda} \mathbf{q} + \mathbf{G} (u, \mathbf{q}) \quad (\text{D.44})$$

where $\mathbf{G}(u, \mathbf{q}) = \mathbf{L}^T \tilde{\mathbf{G}} (r - r^{\text{LF}}, \mathbf{R}\mathbf{q})$ and $u = r - r^{\text{LF}}$. Because of the latter change of variable the linear terms in $q_a$ and $q_b$ vanish in $\mathbf{G}$.

The Jacobian matrix $\mathbf{J}_{\mathbf{g}}(r^{\text{LF}}, \mathbf{t}^{\text{LF}})$ has a null eigenvalue and a negative one. Consequently, the center manifold theorem states that, in the neighborhood of the fold point $(r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}})$, studying the following reduced system

$$\epsilon q_a' = G_1(q_a, \ell(q_a), u) \quad (\text{D.45})$$

is equivalent of studying Eq. (D.44). The function $\ell$ cannot be obtained explicitly and it must therefore be approximated. For sake of simplicity one chooses the tangent space approximation (see for example [214], Chap. 18), i.e. $q_b = \ell(q_a) = 0$.

To obtain a normal form of the dynamic saddle-node bifurcation final approximations are needed. First, because, we study the dynamics near the fold point $(r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}})$ we assume that $f(r, s, \Delta) = f(r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}}) = f^{\text{LF}}$. Then, only terms of order 2 in $q_a$ and of order 1 in $u$ are kept in $G_1$ leading to

$$G_1 \approx a_1 q_a^2 + a_2 u \quad (\text{D.46})$$

with

$$a_1 = \frac{3(\sqrt{1 - 3\mu^2} - 1)\mu^2 + \sqrt{1 - 3\mu^2} + 1}{18\mu^2} \quad \text{and} \quad a_2 = \frac{3\alpha\sqrt{\frac{2 - \sqrt{1 - 3\mu^2}}{\alpha}}}{4\sqrt{2}\sqrt{3\mu^2 + \sqrt{1 - 3\mu^2} + 1}} \quad (\text{D.47})$$

An approximated formulation of the slow flow (D.38) as the normal form of the dynamic saddle-node bifurcation is then finally obtained as

$$\hat{\epsilon} q_a' = h(q_a, v) = q_a^2 + v, \quad (\text{D.48a})$$

$$v' = 1 \quad (\text{D.48b})$$

with

$$v = \frac{a_2}{a_1} u, \quad \hat{\epsilon} = \epsilon f^{\text{LF}} \frac{a_2}{a_1^2} \quad (\text{D.49})$$

<footnote n="4">4. To avoid having to introduce new notations, same symbols as in Section. D.2.2 are used in this section to describe the biorthogonal transformation except for the subscripts which are now letters ($a, b, ...$) instead of numbers ($1, 2, ...$).</footnote>
and the following time switching $\tau \to \frac{a_2 f^{\text{LF}}}{a_1} \tau$ and where $h(q_a, v)$ satisfies the bifurcation conditions $h(0, 0) = 0$ and $\partial_{q_a} h(0, 0) = 0$ and the conditions of the saddle-node bifurcation $\partial_{q_a q_a} h(0, 0) \neq 0$ and $\partial_v f(0, 0) \neq 0$.

The bifurcation is said to be *dynamic* because unlike the *static* saddle-node bifurcation the bifurcation parameter *v* is a slowly varying parameter (because of (D.48b)).

The critical manifold $M_0$ of (D.48) consists of an attracting branch $M_{0,a} = \{(x, y) \in \mathbb{R}^- \times \mathbb{R}^- | x = -\sqrt{-y}\}$ and a repelling branch $M_{0,r} = \{(x, y) \in \mathbb{R}^+ \times \mathbb{R}^- | x = \sqrt{-y}\}$. This critical manifold is hyperbolic except in $(0, 0)$. Finally, it may be noted that the center manifold reduction amounts to stating that locally the critical manifold is approximated by a parabolic function.

## D.4.2 Analytical solution for the normal form of the dynamic saddle-node bifurcation

The method to solve analytically Eq. (D.48) is presented for example in [27]. It begins by dividing Eq. (D.48a) by Eq. (D.48b), that leads to

$$ \frac{\dot{\epsilon} dq_a}{dv} = q_a^2 + v. \qquad (D.50) $$

The subsequent scaling $q_a = \dot{\epsilon}^{1/3}z$ and $v = -\dot{\epsilon}^{2/3}s$ is then introduced to obtain the following equation

$$ \frac{dz}{ds} = -z^2 + s. \qquad (D.51) $$

Setting $z(s) = \varphi'(s)/\varphi(s)$ in (D.51) yields the linear second order equation $\varphi''(s) = s\varphi(s)$, whose solution can be expressed in terms of Airy functions as $\varphi(s) = C_1 \text{Ai}(s) + C_2 \text{Bi}(s)$ where $C_1$ and $C_2$ are constants. Assuming that $x(-\infty) = -\sqrt{-y}$, only the contribution of $\text{Ai}(s)$ is kept (this is not proved here). Therefore the solution with respect to the original variable $q_a$, denoted $q_a^*(v)$, is then given by

$$ q_a^*(v) = \dot{\epsilon}^{1/3} \frac{\text{Ai}'(-\dot{\epsilon}^{-2/3}v)}{\text{Ai}(-\dot{\epsilon}^{-2/3}v)}. \qquad (D.52) $$

Note that from Eq. (D.48) to Eq. (D.52) no asymptotic approach has been used.

The first zero of $q_a^*(v)$ is given by the first zero of the Airy function derivative $\text{Ai}'$ and the first singularity by the first zero of $\text{Ai}$ which are tabulated values (see e.g. [1]). Therefore, one has

$$ q_a^*(v) = 0 \quad \text{for} \quad v = 1.01879 \times \dot{\epsilon}^{2/3} \qquad (D.53) $$

and

$$ q_a^*(v) \to \infty \quad \text{for} \quad v = 2.33810 \times \dot{\epsilon}^{2/3}. \qquad (D.54) $$

An illustration of the normal form of the dynamic saddle-node bifurcation is shown in Fig. D.2. The direct numerical integration of Eq. (D.48) with initial condition ($q_a(0) = -1, v(0) - 0.5$) is compared to the analytical solution $q_a^*(y)$ given by (D.52) for $\epsilon = 0.01$. The first zero, the first singularity of $q_a^*(v)$ and the corresponding critical manifold $M_0$ (which consists of an attracting part $M_{0,a}$ and a repelling part $M_{0,r}$) are also shown.

From Eq. (D.43) and (D.49), an approximated analytical solution of the slow flow (D.38) in the neighborhood of the fold point ($r^{\text{LF}}, s^{\text{LF}}, \Delta^{\text{LF}}$) is obtained as

$$ s^*(r) = s^{\text{LF}} + \mathbf{r}_{a,1} \dot{\epsilon}^{1/3} \frac{\text{Ai}'(-\dot{\epsilon}^{-2/3} \frac{a_2}{a_1} (r - r^{\text{LF}}))}{\text{Ai}(-\dot{\epsilon}^{-2/3} \frac{a_2}{a_1} (r - r^{\text{LF}}))} \qquad (D.55) $$

where $a_1, a_2$ and $\dot{\epsilon}$ are given by (D.47) and (D.49) respectively and $\mathbf{r}_{a,1}$ is the first coordinate of the vector $\mathbf{r}_a$.

The previous analysis reveals the nontrivial $\epsilon$-dependence of the slow flow near the fold point, involving the fractional exponents $1/3$ and $2/3$. Consequently, near the fold point, the trajectory no longer tracks the attracting branch of $M_{0,a}$ of the critical manifold at a distance of order $\epsilon$ as it is the case when the trajectory evolves slowly far from the fold point (as predicted by the Fenichel theory in the normally hyperbolic branches of $M_0$).

Two particular values of $r$ can be deduced from Eq. (D.55): the value at which $s^*(r) = s^{\text{LF}}$, denoted $r^0$, and the value for which $s^*(r)$ tends toward infinity, denoted $r^\infty$. These values of $r$ correspond to the first zeros of $\text{Ai}'$ and $\text{Ai}$ respectively (see Eqs. (D.53) and (D.54)). Therefore, one has

$$ s^*(r) = s^{\text{LF}} \quad \text{for} \quad r = r^0 = r^{\text{LF}} + K^0 \dot{\epsilon}^{2/3} \qquad (D.56) $$

$$ \text{where } K^0 = 1.01879 \times \frac{a_1}{a_2} \left( f^{\text{LF}} \frac{a_2}{a_1} \right)^{2/3} \text{ and} $$

$$ s^*(r) \to \infty \quad \text{for} \quad r = r^\infty = r^{\text{LF}} + K^\infty \dot{\epsilon}^{2/3} \qquad (D.57) $$

$$ \text{where } K^\infty = 2.33810 \times \frac{a_1}{a_2} \left( f^{\text{LF}} \frac{a_2}{a_1} \right)^{2/3}. $$
<figure><img src="image_99.png" /><figcaption>FIGURE D.2. Illustration of the normal form of the dynamic saddle-node bifurcation. Result of the numerical integration of Eq. (D.48) with initial condition ($q_a(0) = -1, v(0) = -0.5$) (red dashed line) compared to the analytical scaling law $q_a^*(y)$ given by (D.52) (blue line, the dashed parts are the horizontal asymptotes of $q_a^*(y)$ corresponding to the zeros Airy function) for $\hat{\epsilon} = 0.01$. The first zero and the first singularity of $q_a^*(v)$ (orange and green lines respectively) and the corresponding critical manifold $M_0$ (attracting part $M_{0,a}$ in black and repelling part $M_{0,r}$ in gray) are also shown.</figcaption></figure>

### D.4.3 New theoretical prediction of the mitigation limit

From the results obtained in the previous section, especially Eq. (D.57), it is possible to obtain a more accurate prediction of the mitigation limit than (D.36) which takes into account the value of the small parameter $\epsilon$.

As seen previously (see Section. D.3.5), the mitigation limit is the value of $\rho$ for which the fixed point $s_2^e$ becomes smaller than $s^A$ (the abscissa of the arrival point of the trajectory on the critical manifold) or for which $s_2^e$ no longer exists (see Eq. (D.34)). Within the zeroth-order approximation the arrival point is $s^A = s^U$ defined by (E.53).

It is not easy to derive a new expression of $s^A$ from the scaling law (D.55). However, one can assume that when the trajectory reaches the right attracting par of $M_0$, the limit value $r^\infty$ has already been reached. Therefore, the solution of $s_2^e = s^A$, denoted $\rho_\epsilon^U$, is also the solution of

$$H(s_2^e(\rho_\epsilon^U)) = r^\infty, \qquad (D.58a)$$

$$= r^{\text{LF}} + K^{\infty}(\rho_{\epsilon}^{\text{U}})\epsilon^{2/3}, \qquad (D.58b)$$

$$= H(s^{\text{LF}}) + K^{\infty}(\rho_{\epsilon}^{\text{U}})\epsilon^{2/3}. \qquad (D.58c)$$

Except for very small values of $\epsilon$ for which $\rho_\epsilon^U \approx \rho_0^U$ and the zeroth-order approximation is sufficient, the special value of $\mu$ solution of $\rho_\epsilon^U = \rho^S$ is in general larger than $1/\sqrt{3}$, the upper limit of the study interval (see Eq. (E.50)). Therefore, we consider here that the mitigation limit is directly $\rho_\epsilon^U$.

From Eq. (D.58) a new and more accurate value of the optimal value of the NES damping coefficient can be derived. Indeed, denoting $H(s_2^e(\rho)) = f_1(\rho, \mu)$, $H(s^{\text{LF}}) = f_2(\mu)$ and $K^\infty(\rho)\epsilon^{2/3} = f_3(\rho, \mu)$ one has

$$d\rho_{\mu} = \frac{d_{\mu}f_2 + \partial_{\mu}f_3 - \partial_{\mu}f_1}{\partial_{\mu}f_1 - \partial_{\mu}f_3}. \qquad (D.59)$$

Therefore solving $d_\mu\rho = 0$ together with Eq. (D.58c) gives the coordinates of the maxima of $\rho_\epsilon^U(\mu)$ and then optimal value of the NES damping coefficient, denoted $\mu_\epsilon^{\text{opt}}$.

Introducing a new perturbation parameter $\nu = \epsilon^{2/3}$, Eq. (D.58) takes the following form

$$H(s_2^e(\rho_\epsilon^U)) = H(s^{\text{LF}}) + \nu K^{\infty}(\rho_\epsilon^U). \qquad (D.60)$$
Because $0 < \epsilon \ll 1$ one has also $0 < \nu \ll 1$, the latter equation is a perturbation of $H(s_2^e) = H(s^{\text{LF}})$ which is equivalent to $H(s_2^e) = H(s^U)$ and consequently to $s_2^e = s^U$ whose solution is $\rho_0^U$ given by Eq. (D.35).

Therefore, to find an explicit formulae for the solution of (D.58), a regular perturbation approach is used by expressing $\rho_{\epsilon}^U$ as a power series of $\nu$

$$ \rho_{\epsilon}^{U} = \rho_{\epsilon,0}^{U} + \nu \rho_{\epsilon,1}^{U} + \dots \qquad (D.61) $$

Then, substituting (D.61) into (D.60) and equating the coefficients of identical power of $\nu$ on both sides of the equation lead to the following set of equations (up to the first order)

$$ \nu^0: \qquad H(s_2^e(\rho_{\epsilon,0}^U)) = H(s^{\text{LF}}) \qquad (D.62a) $$

$$ \nu^1: \qquad \rho_1^* d_s H(s_2^e(\rho_{\epsilon,0}^U)) d_\rho s_2^e(\rho_{\epsilon,0}^U) = K^\infty(\rho_{\epsilon,0}^U) \qquad (D.62b) $$

As expected, from (D.62a) one has $\rho_{\epsilon,0}^U = \rho_0^U$ (see Eq. (D.35)) and therefore Eq. (D.62b) yields

$$ \rho_{\epsilon,1}^{U} = \frac{K^{\infty}(\rho_{0}^{U})}{d_{s}H(s_{2}^{e}(\rho_{0}^{U}))d_{\rho}s_{2}^{e}(\rho_{0}^{U})} \qquad (D.63) $$

Finally, to a first approximation, one has the following asymptotic expression for $\rho_{\epsilon}^U$

$$ \rho_{\epsilon}^{U} = \rho_{0}^{U} + \epsilon^{2/3} \frac{K^{\infty}(\rho_{0}^{U})}{d_{s}H(s_{2}^{e}(\rho_{0}^{U}))d_{\rho}s_{2}^{e}(\rho_{0}^{U})} \qquad (D.64) $$

The asymptotic series (D.61) converges if the functions of $\rho$ in Eq. (D.60) (i.e. $s_2^e(\rho)$ and $K^\infty(\rho)$) can be expanded as a convergent Taylor series at $\rho = \rho_0^U$. Unfortunately $s_2^e(\rho)$ has an infinite derivative at $\rho = \rho^S$ (see Eq. (D.34)). Because of that, at $\mu = \mu^*$ where $\rho_0^U = \rho^S$ the term $\rho_{\epsilon,1}^U$ in the series Eq. (D.61) is equal to zero but the following terms grow to infinity making the series divergent. However, except at $\mu = \mu^*$ and for sufficiently small values of $\epsilon$, the expression (D.64) shows again a nontrivial $\epsilon$-dependence stating that $\rho_{\epsilon}^U$ is a $O(\epsilon^{2/3})$-perturbation of $\rho_0^U$.

Equating to zero the derivative, with respect to $\mu$, of the right-hand side of Eq. (D.64) and performing again a perturbative technique with respect to $\nu = \epsilon^{2/3}$, an approximative expression of $\mu_{\epsilon}^{\text{opt}}$ is obtained as

$$ \mu_{\epsilon}^{\text{opt}} = \mu_{0}^{\text{opt}} + \epsilon^{2/3} \frac{d_{\mu} \left( \frac{K^{\infty}(\rho_{0}^{U})}{d_{s}H(s_{2}^{e}(\rho_{0}^{U}))d_{\rho}s_{2}^{e}(\rho_{0}^{U})} \right)_{\mu=\mu_{0}^{\text{opt}}}}{d_{\mu\mu}(\rho_{0}^{U})_{\mu=\mu_{0}^{\text{opt}}}} \qquad (D.65) $$

where $\mu_0^{\text{opt}}$ is given by Eq. (D.37). Obviously, the nontrivial $\epsilon$-dependence is found again as, if $\epsilon$ is sufficiently small, $\mu_{\epsilon}^{\text{opt}}$ is a $O(\epsilon^{2/3})$-perturbation of $\mu_0^{\text{opt}}$. Of course the limitations previously mentioned for Eq. (D.64) are also valuable for Eq. (D.65).

## D.5 APPLICATION TO AN AEROELASTIC AIRCRAFT WING SYSTEM COUPLED TO AN NES

The methodology presented in previous sections is illustrated here on an aeroelastic aircraft wing system coupled to an NES (see Fig. D.3) which has already been used in the past in several works about LCO mitigation by means of one or several NESs [a8, 78, 123, 124]. First, equations of motion the wing-NES system are presented in Section. D.5.1. Then, the comparison between theoretical results and numerical simulations is performed in Section. D.5.2.

### D.5.1 Equations of the aircraft wing model coupled to one NES

The equations of motion of the model presented in Fig. D.3 can be derived for example by means of the Lagrange equations. Assuming a small angle between the wing and the horizontal, the Lagrange equations yield the following equations of motion with respect to the heave $z$ and the angle of attack (pitch) $\tilde{\phi}$ of the wing and to the NES displacement $y$ (see Fig. D.3)

$$ \begin{aligned} m\ddot{z} + S_{\varphi}\ddot{\tilde{\phi}} + K_z z + C_z \dot{z} + qP_a\Gamma(\tilde{\phi} + \dot{z}/U) + K_z^{\text{NL}}\dot{z}^3 + \\ C_y(\dot{z} - d\tilde{\phi} - \dot{y}) + K_y^{\text{NL}}(z - d\tilde{\phi} - y)^3 = 0 \end{aligned} \qquad (D.66a) $$
<figure><img src="image_100.png" /><figcaption>**FIGURE D.3.** Sketch of the two DOFs aircraft wing coupled to one NES. *z* and *z̃* are respectively the heave and the angle of attack (pitch) of the wing and *y* is the displacement of the NES. *A* is the aerodynamic center, *B* the elastic axis, *G* the center of gravity of the aircraft wing. *e* is the location aerodynamic center *A* measured from *B* (positive ahead of *B*). *K*<sub>*z*</sub> and *K*<sub>*φ*</sub> are the linear heave and pitch stiffnesses respectively whereas *K*<sub>*z*</sub><sup>NL</sup> and *K*<sub>*φ*</sub><sup>NL</sup> are the cubic heave and pitch stiffnesses. *C*<sub>*z*</sub> and *C*<sub>*φ*</sub> are the heave and pitch damping coefficients. *U* is the constant and uniform flow speed around the wing and *d* is the offset attachment of the NES to the wing, also measured from *B*.</figcaption></figure>

$$ I_{\varphi} \ddot{\tilde{\varphi}} + S_{\varphi} \ddot{z} + K_{\varphi} \tilde{\varphi} + C_{\varphi} \dot{\tilde{\varphi}} - q e P_a \Gamma (\tilde{\varphi} + \dot{z}/U) + K_{\varphi}^{\text{NL}} \tilde{\varphi}^3 + d C_y (\dot{z} - d \dot{\tilde{\varphi}} - \dot{y}) + d K_y^{\text{NL}} (z - d \tilde{\varphi} - y)^3 = 0 \quad (\text{D.66b}) $$

$$ m_y \ddot{y} + C_y (d\dot{\tilde{\varphi}} + \dot{y} - \dot{z}) + K_y^{\text{NL}} (d\tilde{\varphi} + y - z)^3 = 0 \quad (\text{D.66c}) $$

where again the time derivative is denoted “·”. *b* = *L*/2 is the semichord length. *A* is the aerodynamic center, *B* the elastic axis, *G* the center of gravity of the aircraft wing. *e* is the location aerodynamic center *A* measured from *B* (positive ahead of *B*). The quantities *m* and *I*<sub>φ</sub> are respectively the mass of the wing and its moment of inertia with respect to *B*. *S*<sub>φ</sub> = *mBG* is the mass unbalance in the wing. *K*<sub>*z*</sub> and *K*<sub>*φ*</sub> are the linear heave and pitch stiffnesses respectively whereas *K*<sub>*z*</sub><sup>NL</sup> and *K*<sub>*φ*</sub><sup>NL</sup> are the cubic heave and pitch stiffnesses. *C*<sub>*z*</sub> and *C*<sub>*φ*</sub> are the heave and pitch damping coefficients. *U* is the constant and uniform flow speed around the wing and *q* = 1/2 *ρ*<sub>∞</sub>*U*<sup>2</sup> is the dynamic pressure where *ρ*<sub>∞</sub> is the density of the flow. *P*<sub>a</sub> is the planform of the wing, *Γ* is the lift curve slope and *d* is the offset attachment of the NES to the wing. Finally, *m*<sub>y</sub>, *C*<sub>y</sub> and *K*<sub>y</sub><sup>NL</sup> (*m* = 1, ..., *M*) are the mass, the damping coefficient and the cubic stiffness of the NES respectively.

For convenience, Eq. (D.66) is written in a dimensionless form as follows

$$ \tilde{x}'' + s_{\varphi} \tilde{\varphi}'' + \Omega^2 \tilde{x} + \zeta_x \tilde{x}' + \gamma \Theta (\Theta \tilde{\varphi} + \tilde{x}') + \tilde{\xi}_x \tilde{x}^3 + \tilde{m}_h \tilde{h}'' = 0 \quad (\text{D.67a}) $$

$$ r_{\varphi}^2 \tilde{\varphi}'' + s_{\varphi} \tilde{x}'' + r_{\varphi}^2 \tilde{\varphi} + \zeta_{\varphi} \tilde{\varphi}' - \eta \gamma \Theta (\Theta \varphi + \tilde{x}') + \tilde{\xi}_{\varphi} \tilde{\varphi}^3 + \delta \tilde{m}_h \tilde{h}'' = 0 \quad (\text{D.67b}) $$

$$ \tilde{m}_h \tilde{h}'' + \tilde{\zeta}_h (\delta \tilde{\varphi}' + \tilde{h}' - \tilde{x}') + \tilde{\xi}_h (\delta \tilde{\varphi} + \tilde{h} - \tilde{x})^3 = 0 \quad (\text{D.67c}) $$

where the time *t* has been replaced by the dimensionless time *t*′ = *ω*<sub>*φ*</sub>*t* (with *ω*<sub>*φ*</sub> = √(*K*<sub>*φ*</sub>/*I*<sub>*φ*</sub>) and time derivative *d*<sub>*t*</sub>' is denoted "·". The dimensionless displacements are defined by *x* = *z*/*b*, *ħ* = *y*/*b*. Moreover, *δ* = *d*/*b*, *s*<sub>*φ*</sub> = *S*<sub>*φ*</sub>/(*mb*), *Ω* = *ω*<sub>*z*</sub>/*ω*<sub>*φ*</sub> (with *ω*<sub>*z*</sub> = √(*K*<sub>*z*</sub>/*m*)) and *η* = *e*/*b*. *r*<sub>*φ*</sub> = √(*I*<sub>*φ*</sub>/(*mb*<sup>2</sup>)) is the radius of
gyration of the cross section of the wing. The dimensionless nonlinear stiffnesses and damping coefficients are
$\tilde{\xi}_x = K_z^{\text{NL}}b^2/(m\omega_\varphi^2)$, $\xi_x = C_z/(m\omega_\varphi)$, $\tilde{\xi}_\varphi = K_\varphi^{\text{NL}}/(mb^2\omega_\varphi^2)$, $\zeta_\varphi = C_\varphi/(mb^2\omega_\varphi)$, $\tilde{\zeta}_h = K_y^{\text{NL}}b^2/(m\omega_\varphi^2)$, $\tilde{\zeta}_h = C_y/(m\omega_\varphi)$.
The mass and density ratios are defined by $\tilde{m}_h = m_y/(m\omega_\varphi^2)$ and $\gamma = bP_a\rho_\infty\Gamma/(2m)$ respectively.
Finally, the bifurcation parameter under consideration is the reduced speed of the flow $\Theta = U/(b\omega_\varphi)$.

Eq. (D.67) has the general form given by Eq. (D.2) with $\tilde{\mathbf{x}} = (\tilde{x}, \tilde{\varphi})^T$ and

$$
\begin{equation}
\begin{aligned}
\tilde{\mathbf{M}} &= \begin{pmatrix} 1 & s_{\varphi} \\ s_{\varphi} & r_{\varphi}^2 \end{pmatrix}, & \tilde{\mathbf{C}} &= \begin{pmatrix} \zeta_x + \gamma\Theta & 0 \\ -\eta\gamma\Theta & \zeta_{\varphi} \end{pmatrix}, & \tilde{\mathbf{K}} &= \begin{pmatrix} \Omega^2 & \gamma\Theta^2 \\ 0 & r_{\varphi}^2 - \eta\gamma\Theta^2 \end{pmatrix}, \\
\tilde{\mathbf{A}} &= (1, -\delta), & \tilde{\mathbf{B}} &= (1, -\delta)^T, & \tilde{\mathbf{G}}^{\text{NL}} &= (\tilde{\xi}_x \tilde{x}^3, \tilde{\xi}_\varphi \tilde{\varphi}^3)^T.
\end{aligned}
\tag{D.68}
\end{equation}
$$

The expression of the nonlinear vector function is $\tilde{\mathbf{g}}^{\text{NL}} = (\tilde{\xi}_x \tilde{x}^3, \tilde{\xi}_\varphi \tilde{\varphi}^3)^T$. Moreover, one defines rescaled parameters and variables as $\tilde{\xi}_x = \epsilon \xi_x$, $\tilde{\xi}_\varphi = \epsilon \xi_\varphi$, $\tilde{m}_h = \epsilon m_h$, $\tilde{\zeta}_h = \epsilon \zeta_h$, $x = \frac{\tilde{x}}{\sqrt{\epsilon}}$, $\varphi = \frac{\tilde{\varphi}}{\sqrt{\epsilon}}$ and $h = \frac{\tilde{h}}{\sqrt{\epsilon}}$ to obtain the final form of the equations of motion

$$
x'' + s_{\varphi}\varphi'' + \Omega^2 x + \zeta_x x' + \gamma\Theta(\Theta\varphi + x') + \epsilon m_h h'' = 0 \quad (D.69a)
$$

$$
r_{\varphi}^2 \varphi'' + s_{\varphi} x'' + r_{\varphi}^2 \varphi + \zeta_{\varphi} \varphi' - \eta \gamma \Theta (\Theta \varphi + x') + \epsilon \delta m_h h'' = 0 \quad (D.69b)
$$

$$
h'' + \zeta_h (\delta\varphi' + h' - x') + \xi_h (\delta\varphi + h - x)^3 = 0 \quad (D.69c)
$$

with $\zeta_h = \tilde{\zeta}_h/m_h$ and $\xi_h = \tilde{\xi}_h/m_h$.

The physical bifurcation parameter under consideration is the reduced speed of the air flow $\Theta$ and the
following parameters for the primary system are used

$$
\begin{align}
& \zeta_x = 0.01, \quad \zeta_\varphi = 0.01, \quad s_\varphi = 0.2, \quad r_\varphi = 0.5, \tag{D.70} \\
& \Omega = 0.5, \quad \eta = 0.4, \quad \text{and} \quad \gamma = 0.2. \nonumber
\end{align}
$$

With the parameters (D.70), the aircraft wing model undergoes an Hopf bifurcation at $\Theta_{\text{Hopf}}^{\text{wo}} = 0.933$ at which the real part of the first eigenvalue of Eq. (D.69) vanishes, i.e. Re($\lambda_1$) = 0, and then becomes positive. It is recalled here that the pulsation $\omega$ and the generalized bifurcation parameter $\rho$ are defined from the eigenvalue $\lambda_1$ of the unstable mode by Eq. (D.14).

The NES setup is defined by

$$
m_h = 1, \quad \delta = -0.8, \quad \xi_h = 5 \tag{D.71}
$$

and $\zeta_h$ and $\epsilon$ which are specified below depending on the considered example.

The several response regimes described in Section. D.3.4 are illustrated here by means of numerical si-
mulations of the wing-NES system (D.67) (see Fig. D.4) with, in addition to (D.70) and (D.71), $\epsilon = 0.01$,
$\zeta_h = 0.3$, $\xi_x = 4$ and $\xi_\varphi = 8$. Four values of the reduced speed are used, namely : $\Theta = 0.94$ (see Fig. D.4(a)),
$\Theta = 0.945$ (see Fig. D.4(b)), $\Theta = 0.97$ (see Fig. D.4(c)) and $\Theta = 0.98$ (see Fig. D.4(d)) corresponding respec-
tively to complete suppression, mitigation through a periodic regime, mitigation through a SMR and finally
no mitigation of the harmful LCOs of the primary structure. For comparison purposes, for each situation, a
direct numerical integration of the aircraft wing system without NES is also shown.

D.5.2 Comparison between theoretical results and numerical simulations

The slow flow dynamics of Eq. (D.69), its critical manifold $M_0$ (E.32), the zeroth-order analytical expres-
sions of the mitigation limit (D.36) and of the optimal value of NES damping coefficient (D.37) are derived
following the method described in Sections. D.3.1, D.3.2 and D.3.5 respectively. Note that through Eq. (D.19),
the parameters $\mu$ and $\alpha$ are here defined as $\mu = \zeta_h/\omega$ and $\alpha = \xi_h/\omega^2$. Simultaneously, the scaling law of the
trajectory in the neighborhood of the left fold point of $M_0$, i.e. Eq. (D.55), is computed. Afterwards, solving
Eq. (D.58) and finding the root of Eq. (D.59), $\rho_{\epsilon}^{\text{U}}(\mu)$ and $\mu_{\epsilon}^{\text{opt}}$ are obtained. Finally, the asymptotic analytical
expressions of the latter, given by Eqs. (D.64) and(D.65) respectively, are computed.

To illustrate the proposed analytical procedure, Fig. D.5(a) shows the critical manifold M₀ (E.32), super-
imposed to the direct numerical integration of the slow flow (D.38) (see also Fig. D.5(b) for the times series
of s(τ) and r(τ) focused on the first relaxation cycle) and its scaling law given by Eq. (D.55). The values of
rLF = H(sLF) (see Eq. (E.49)) and r∞ = rLF + K∞ε²/³ (see Eq. (D.57)) are depicted by horizontal lines.
The numerical integration of the slow flow (see red dashed line in Fig. D.5(a)) shows a relaxation oscillations
scenario in (s, r)-plane. Indeed, the trajectory of undergoes a succession of almost horizontal fast parts and
slow parts evolving in the O(ε)-vicinity of the critical manifold M₀ (here, at the scale of the figure and apart
from the fold points, the trajectory and M₀ seem superimposed). The first horizontal fast part is from the
initial condition (s(0) = 0.15, r(0) = 0.04) to the left attracting part of M₀ and the other are fast jumps
<figure><img src="image_101.png" /></figure>

**FIGURE D.4.** Direct numerical integration of the wing-NES system (D.67) (in blue) depicting the several response regimes described in Section. D.3.4 with, in addition to (D.70) and (D.71), $\epsilon = 0.01$, $\zeta_h = 0.3$, $\xi_h = 4$, $\xi_x = 4$ and $\xi_\phi = 8$. Four values of the reduced speed are used, namely : (a) $\Theta = 0.94$, (b) $\Theta = 0.945$, (c) $\Theta = 0.97$ and (d) $\Theta = 0.98$ corresponding respectively to complete suppression, mitigation through a periodic regime, mitigation through a SMR and finally no mitigation of the harmful limit cycle oscillations of the primary structure. For comparison purposes, for each situation, a direct numerical integration of the wing system without NES is also shown (in green).

between the left and right attracting parts of $M_0$ as described in Section. D.3.4 and in the introduction of Section. D.4.

The figure reveals a good agreement between the numerical simulation of the slow flow and its analytical scaling law in the neighborhood of the left fold point. One can see also that choosing $r^\infty$ as the ordinate of the point of intersection between the trajectory of the slow flow and the right attracting branch of $M_0$ is a more accurate approximation than choosing $r^{LF}$ as this is the case within the zeroth-order approximation.

In Fig. D.6, the comparison is now performed in terms of mitigation limit. The theoretical predictions, the one obtained within the zeroth-order approximation (D.36) and those derived from the scaling law (solving numerically Eq. (D.58) and the one given by Eq. (D.64)) are compared with mitigation limits measured on direct numerical integration of the approximated slow flow (D.38) and on the full order system (D.67). In these last two cases, direct numerical integrations are performed for increasing values of $\rho$ (by increasing $\Theta$) and
<figure><img src="image_102.png" /><figcaption>**FIGURE D.5.** Illustration of the proposed analytical procedure. (a) The critical manifold $M_0$ (E.32) (black line) superimposed to the direct numerical integration of the slow flow (D.38) (red line) and its scaling law near the left fold point given by Eq. (D.55) (blue line, the dashed parts are the horizontal asymptotes of $s^*(r)$ due to the zeros of the Airy function). The values of $r^{\text{LF}} = H(s^{\text{LF}})$ (see Eq. (E.49)) and $r^{\infty} = r^{\text{LF}} + K^{\infty}\epsilon^{2/3}$ (see Eq. (D.57)) are depicted by magenta and green horizontal lines respectively. (b) Times series of $s(\tau)$ (blue line) and $r(\tau)$ (green line) obtained from the numerical integration of the slow flow (D.38). The times series are focused on the first relaxation cycle. Parameters (D.70) and (D.71) are used with $\zeta_h = 0.25$, $\epsilon = 0.005$ and $\Theta = 0.95$.</figcaption></figure>

the mitigation limit corresponds to the last value of $\rho$ allowing a harmless situation. In general one observes at this value the transition from a SMR to a “no mitigation regime" as shown for example in Fig. D.4. All the mitigations limits are plotted with respect to the NES damping coefficient $\mu$ (with $\mu \in [0, 1/\sqrt{3}]^5$) and for four values of the perturbation parameter $\epsilon$. One takes first a very little value ($\epsilon = 0.001$) to be sure to respect the assumptions of the asymptotic analysis (i.e. $0 < \epsilon \ll 1$) and thus validate the method under these conditions. Then, three larger values are taken ($\epsilon = 0.005, 0.02$ and $0.1$) to evaluate the robustness of the method when one deviates from the assumption of a small $\epsilon$. To relate these results to the physical bifurcation parameter $\Theta$, the latter is plotted with respect to the generalized bifurcation parameter $\rho$ in Fig. D.7 for the same values of $\epsilon$ as in Fig. D.6.

First, in Fig. D.6, one can see that the zeroth-order approximation gives an inaccurate prediction, even for the smallest value of $\epsilon = 0.001$. Moreover, one observes that, on the contrary, solving Eq. (D.58) provides a very good prediction of the mitigation limits measured on numerical simulations. For the two following values of $\epsilon$ (0.005 and 0.02), these theoretical prediction is closer to the mitigation limits obtained from numerical simulations of the full order system (D.67) than to those measured on numerical integrations of the approximated slow flow (D.38). This is surprising because Eq. (D.58) is derived from Eq. (D.38) but, in general, in these kind of unstable systems coupled to NES, the mitigation limits measured on the full order system are slightly smaller than those measured on the slow flow. Here, when $\epsilon$ increases the prediction depreciates becoming smaller than mitigation limits obtained from numerical simulations of Eq. (D.38) and getting closer to those measured on full order system (D.67). For the largest value $\epsilon = 0.1$, the theoretical mitigation limit obtained solving Eq. (D.58) (in green on the figure) depreciates again but it is still satisfactory in a qualitative point of view. Moreover one can see that the theoretical prediction (in green) no longer has a local maximum and the optimal value of $\mu$ is equal to $1/\sqrt{3}$, its limit value for the existence of the SMR regimes (see Eq. (E.50)).

Finally, the prediction provided by Eq. (D.64) deteriorates enormously for large values of $\epsilon$. This may be due to the non-convergent nature of the series (D.61), as already mentioned above (see comments after Eq. (D.64)).

<footnote n="5">5. Remembering that in this work only situations in which Eq. (E.50) holds are considered.</footnote>
<figure><img src="image_103.png" /><figcaption>**FIGURE D.6.** Theoretical mitigation limits: the one obtained within the zeroth-order approximation (D.36) (blue line) and those derived from the scaling law (solving numerically Eq. (D.58) (green line) and the one given by Eq. (D.64) (orange line)) compared with mitigation limits measured on direct numerical integration of the approximated slow flow (D.38) (gray dots) and on the full order system (D.67) (black dots). Parameters (D.70) and (D.71) are used and ((a) $\epsilon = 0.001$, (b), $\epsilon = 0.005$, $\epsilon = 0.02$ and (d) $\epsilon = 0.1$.</figcaption></figure>

Finally, in Fig. D.8, the theoretical optimal value of the NES damping coefficient obtained from the zeroth-order approximation $\mu_0^{\text{opt}}$ (see Eq. (D.37)) is compared to those derived from the scaling law (D.55): the one obtained solving numerically $d_{\mu}\rho = 0$ (with $d_{\mu}\rho$ given by Eq. (D.59)) together with Eq. (D.58c) and the asymptotic expression given by Eq. (D.65). All this is plotted with respect to the perturbation parameter $\epsilon$. Assuming a good agreement between the theoretical prediction of the mitigation limit obtained solving Eq. (D.58) and the actual value (as shown by Fig. D.6), only theoretical values are presented in Fig. D.8. The figure highlights the relevance of the proposed results especially for large values of the parameter $\epsilon$. We can also see that the asymptotic expression is here able to describe qualitatively the evolution of $\mu^{\text{opt}}$ with respect to the perturbation parameter $\epsilon$ for $\epsilon \lesssim 0.04$. For $\epsilon \gtrsim 0.04$, the prediction obtained from Eqs. (D.59) and (D.58c) tends to $1/\sqrt{3}$ whereas the prediction given by Eq. (D.65) continues to grow.

## D.6 CONCLUSION

In this work a mechanical system with one unstable mode coupled to a Nonlinear Energy Sink (NES) has been studied. The equations of motion of the system have been first simplified using a reduced-order model for the primary structure by keeping only its unstable modal coordinates. Introducing a small perturbation parameter related to the mass ratio between the NES and the primary structure, the slow flow has been then derived by means of the complexification-averaging method. Because of the presence of the small perturbation parameter, the slow flow is governed by two different time scales. More precisely, within the framework of the geometric singular perturbation, in its real form it appears as a (2, 1)-fast-slow system. The critical manifold of this slow flow has been computed. Within the zeroth-order approximation, i.e. the analysis within the limit case in which the perturbation parameter is equal to zero, the solution of the slow flow is approximated by a
<figure><img src="image_104.png" /><figcaption>**FIGURE D.7.** The physical bifurcation parameter $\Theta$ with respect to the generalized bifurcation parameter $\rho$. As in Fig. D.6, parameters (D.70) are used and (a) $\epsilon = 0.001$, (b), $\epsilon = 0.005$, $\epsilon = 0.02$ and (d) $\epsilon = 0.1$.</figcaption></figure>

<figure><img src="image_105.png" /><figcaption>**FIGURE D.8.** The theoretical optimal value of the NES damping coefficient with respect to the small perturbation parameter $\epsilon$. The one obtained from the zeroth-order approximation $\mu_0^{\text{opt}}$ given by Eq. (D.37) (blue line) is compared to those derived from the scaling law (D.55) : the one obtained solving numerically $d_{\mu}\rho = 0$ (with $d_{\mu}\rho$ given by Eq. (D.59)) together with Eq. (D.58c) (green line) and the asymptotic expression given by Eq. (D.65) (orange line). As previously, parameters (D.70) and (D.71) are used.</figcaption></figure>

discontinuous phase trajectory in which the fast and slow epochs are described independently of each other. In particular, the slow epoch are supposed to occur on attracting parts of the critical manifold. Moreover, in this context a theoretical prediction of the mitigation has been obtained. The mitigation limit being defined as the value of the generalized bifurcation parameter, which is obtained rescaling the eigenvalue real part of the primary structure corresponding to the unstable mode, which separates harmless situations (in which the NES acts) from the harmful situation (in which the NES does not act). The results obtained within the zeroth-order approximation (both phase trajectory and mitigation limit) depreciate for the largest values of the perturbation parameter (the latter remaining small in the context of this work).

The novelty of the present work is to use the center manifold theorem to reduce the slow flow, in the neighborhood of a fold point of its critical manifold, to a normal form of the dynamic saddle-node bifurcation. This allowed us to obtain a scaling law for the slow flow in this neighborhood, i.e. a law which describes the
dependance, with respect to the small perturbation parameter, of the distance between the critical manifold
and the actual trajectory of the slow flow. The law reveals that the slow flow scales in a nontrivial way with
respect to the perturbation parameter, involving the fractional exponents 1/3 and 2/3. Finally, a theoretical
prediction of the mitigation limit has been deduced from the scaling law. With certain assumptions, it was
shown that this prediction is a correction with respect to the prediction obtained within the zeroth-order
approximations, involving also a fractional order dependence with respect to the perturbation parameter.

The last section of the paper has been devoted to illustrate and validate the proposed methodology on a an
aeroelastic aircraft wing model coupled to one NES. The accuracy of the present theoretical results compared
to those obtained on numerical simulations of the full order system together with the general nature of the
method suggest that the latter may constitute a tool for the practical design of NES.

ACKNOWLEDGMENTS

We thank the anonymous reviewers for their careful reading of our manuscript and their many relevant
comments and suggestions.
*International Journal of Non-Linear Mechanics 150 (2023) : 104351* [a3]

# EFFECT OF STOCHASTIC FORCING ON THE DYNAMIC BEHAVIOR OF A SELF-SUSTAINED OSCILLATOR COUPLED TO A NON-LINEAR ENERGY SINK

Baptiste BERGEOT

## ABSTRACT

In this paper the influence of stochasticity (i.e. a Gaussian white noise forcing) on the dynamic behavior of a self-sustained oscillator coupled to a non-linear energy sink is investigated. To this end, the standard stochastic averaging is used to compute the slow flow dynamics of the system. Preliminary results show that the reasoning which allows to predict the system behavior in the deterministic case can be contradicted in presence of stochasticity. Then, by means of the Monte Carlo method, the stochastic averaging procedure is validated. Finally, two quantities are introduced to highlight more precisely the special features of the stochastic system behavior compared to that of the deterministic system. These are the probability of being in a harmless regime and the First-Passage Time to reach a harmful regime which are computed and investigated combining again the Monte Carlo approach with numerical integrations of the slow flow dynamics. The results obtained show afresh that the stochastic forcing can modify significantly the dynamic behavior of the corresponding deterministic system. Indeed, when they are computed on the latter, the two quantities aforementioned have a discontinuity at the mitigation limit (i.e. the value of the bifurcation parameter under consideration below which the NES acts and above which it no longer acts) revealing an abrupt change of behavior of the coupled system. The paper shows that this typical characteristic of the deterministic system is lost in the presence of stochasticity, the stochastic system becoming smooth at the mitigation limit.

**Keywords** : Passive vibration control, Non-linear energy sink, Self-sustained oscillations, Stochastic averaging, Stochastic forcing

## E.1 INTRODUCTION

The non-linear vibration absorbers known as non-linear energy sinks (NESs) are these days well-known devices used for passive mitigation of unwanted oscillations caused by either external, parametric or self-excitations of a mechanical or acoustical primary structure. In general an NES is defined as a non-linear attachment consisting of a light mass (compared to the total mass of the primary structure), an essentially non-linear spring (most of the time purely cubic) and a viscous linear damper. The dependence between the vibratory amplitude and the oscillating frequency of the NES (because of its strongly non-linear nature) makes it able to resonate at any frequency. After tuning to the primary structure and having absorbed and then dissipated its energy, the NES can detune from the primary structure to avoid the return of energy. This irreversible transfer of vibrational energy from the primary system to the NES is called targeted energy transfert (TET) or energy pumping. In their seminal papers [76, 203] Gendelman, Vakakis and al. explain the TET phenomenon by the interaction between two non-linear modes of vibrations of the system producing a 1 : 1 resonance capture. Reviews of these concepts can be found in [206] and more recently in [129].

Using an NES to alleviate or even suppress limit cycle oscillations (LCOs) resulting from dynamic instabilities has been extensively studied in the literature in a deterministic framework. Mitigation of LCOs of the Van der Pol oscillator has been studied numerically by Lee and al. in their seminal paper [122] and then theoretically by Gendelman and Bar [75]. This latter work has been extended to a Van der Pol-Duffing oscillator coupled to one NES by Domany and Gendelman [59]. A number of works focused on the problem of mitigation, using one or several NESs, of LCOs due to flutter instabilities in aircraft wings. This problem has been first studied both numerically and experimentally, again by Lee. and al. [121, 123, 124]. The theoretical prediction of observed operating regimes has been performed by means of multiple time scales approaches [78, 131] and improved more recently using the center manifold reduction technique [a2]. Also in the context of
self-oscillations induced by fluid-structure interactions, the mitigation of vortex-induced vibrations caused by the non-linear interaction of a laminar flow and a rigid circular cylinder has been first investigated by Tumkur and al. [201] constructing a two-DOFs reduced-order model and validating it by means of a comparison with a finite-element model. An improved and experimentally validated reduced-order model has been proposed by Dai and al. [47]. On the same issue, a recent theoretical work [39] (i) explains the global dynamics of the system by means a three-DOFs reduced-order model coupling fluid-structure interaction framework with the NES attachment and (ii) identifies, through a parametric study, the optimal operational parameter ranges for an efficient NES design. Concerning aeronautics, the use of NESs to control helicopter ground resonance instability has been studied in [a10, a12]. The use of several NESs to alleviate self-sustained oscillations has been analytically studied by Bergeot and Bellizzi in two papers. First, the use of several parallel NESs attached to a Van der Pol oscillator has been analyzed in [a7]. Then, the prediction of the dynamic behavior of a multi-DOFs mechanical system with only one unstable mode and coupled to a set of NESs has been proposed in [a8] with an application to flutter instability of an airplane wing model. Finally, the possibility of mitigating self-sustained oscillations of a linear friction system having two unstable modes has been studied by Bergeot and al. [a9] by means of a sophisticated multiple time scales analysis allowing the understanding of the phenomena underlying the appearance of the many possible regimes of the system.

Many engineering applications, especially in the context of fluid-structure interaction and aeronautics, involve stochasticity which can have an important influence on the deterministic dynamics. Although that often leads to a better understanding of the dynamics of a mechanical system under more realistic conditions, there are very few papers that study the TET mechanisms taking into account stochasticity. The polynomial chaos approach has been used by Gourdon and Lamarque [84] to analyze the NESs behavior during instationary regimes in a two-DOFs academic system. The goal was to verify the robustness of the TET mechanism when the model has uncertain parameters. In the same vein, Cataldo and al. [35] used the Monte Carlo method to study the energy pumping robustness considering the uncertainties of the parameters of a linear primary structure coupled to an NES. The design optimization of parallel NESs taking into account uncertainties has been performed in [30]. In the context of limit cycles mitigation, Pidaparthi and Missoum [162] performed optimization under uncertainties of an NES used to alleviate limit cycles created by aeroelastic instability in an aircraft wing. Snoun and al. [a17] used a multi-element generalized polynomial chaos based method to predict the dynamic behavior of an uncertain friction system coupled to two non-linear energy sinks. A robust method was proposed to optimize non-linear energy sinks (NES) in the same model in [a17].

In the previous cited works, the stochasticity comes from parametric uncertainties considering that para-
meters values cannot be known exactly. This approach, however, assumes that for a given realization of the
process the parameters are constant. Each realization is therefore deterministic. Another approach consists
in investigating the TET taking into account a stochastic time dependent forcing. Even fewer works consider
this approach. Schmidt and Lamarque [176] studied a linear one-DOF primary structure coupled to one NES
in presence of a white noise forcing by solving numerically the Fokker-Planck equation associated with the
stochastic equations of motion of the system. Starosvetsky and Gendelman [189] showed, using numerical
simulations, that an NES can produce complete elimination of undesired response of a linear one-DOF system
subject to a randomly modulated and narrow-band excitation. A TET problem from a linear medium to a
non-linear attachment is studied in the presence of stochasticity by Sapsis and al. [173]. Using a stochastic
averaging method, the authors determine the equation of the stochastic slow flow dynamics of the system. Nu-
merical integration of the associated Fokker-Planck equation reveals that the optimal TET regimes, predicted
in the deterministic case, are preserved and even enhanced because of the interaction between non-linearity
and stochasticity.

In the present paper we consider the problem of mitigation of self-sustained oscillations by means of an
NES in presence of stochasticity (in the form of a Gaussian white noise forcing) which, to the knowledge
of the author, has not been treated in the literature. The problem is studied using the standard stochastic
averaging method [114, 169, 191] to compute the slow flow dynamics of a Van der Pol oscillator undergoing
a Gaussian white noise forcing and coupled to an NES. Then two quantities are introduced to highlight the
special features of the stochastic system, in comparison with the behavior of the deterministic system. These
are the *probability of being in a harmless regime* and the *First-Passage Time to reach a harmful regime* which
are studied using the Monte Carlo method [64] combined with numerical simulations of the stochastic slow
flow.

The paper is organized as follows. In Section E.2 the equations of motion of the full-order system under study - i.e. a Van der Pol oscillator forced by a stochastic forcing (Gaussian white noise) and coupled to a purely cubic non-linear energy sink - are derived. The standard stochastic averaging method is used in Section G.3 to compute the slow flow dynamics of the full-order system. Section E.4 recalls the typical behavior of the deterministic system (in Section E.4.1) with the aim of highlighting some special features caused by stochasticity (in Section E.4.2). The probability of being in a harmless regime and the First-Passage Time to
<figure><img src="image_106.png" /><figcaption>FIGURE E.1. A Van der Pol oscillator (primary system) with stochastic forcing (Gaussian white noise) $\hat{\nu}\xi$ and coupled to an NES.</figcaption></figure>

reach a harmful regime are defined and then investigated by means of the Monte Carlo method in Sections E.5 and E.6 respectively. Finally, concluding remarks and perspectives are formulated in Section F.8.

## E.2 EQUATIONS OF THE MODEL

In this section we derive the equations of motion of the full-order system under study, i.e. a Van der Pol (VdP) oscillator - used as an archetypal self-sustained oscillator - forced by a Gaussian white noise and coupled to an ungrounded purely cubic non-linear energy sink (NES). A schematic representation of this system is depicted in Fig. E.1. The equations of motion of this system are as follows

$$ m_1 \frac{d^2 y_1}{dt^2} + c_1^{\text{NL}} \frac{dy_1}{dt} (ry_1^2 - R^2) + k_1 y_1 + c_2 \left( \frac{dy_1}{dt} - \frac{dy_2}{dt} \right) + k_2^{\text{NL}} (y_1 - y_2)^3 = \hat{\nu}\xi(t) \quad (\text{E.1a}) $$

$$ m_2 \frac{d^2 y_2}{dt^2} + c_2 \left( \frac{dy_2}{dt} - \frac{dy_1}{dt} \right) + k_2^{\text{NL}} (y_2 - y_1)^3 = 0 \quad (\text{E.1b}) $$

where $m_1$ and $m_2$ are the masses of the primary VdP oscillator and of the NES respectively. The parameter $k_1$ is the linear stiffness of the VdP oscillator and $c_1^{\text{NL}}$ and $rc_1^{\text{NL}}$ characterize its negative and non-linear dampings respectively. The NES is a purely cubic oscillator characterized by its linear damping coefficient $c_2$ and non-linear stiffness $k_2^{\text{NL}}$. The term $\xi(t)$ is a unitary idealized white noise process with zero mean, i.e.

$$ \mathbb{E}[\xi(t)] = 0 \quad \text{and} \quad \mathbb{E}[\xi(t)\xi(t + \tau)] = \delta(\tau) \quad (\text{E.2}) $$

where $\hat{\nu}$ is the noise level.

For convenience Eq. (G.8) is rescaled leading to the following dimensionless system of differential equations

$$ \ddot{x}_1 + \epsilon \rho \dot{x} (rx_1^2 - 1) + x_1 + \epsilon \mu (\dot{x}_1 - \dot{x}_2) + \epsilon \alpha (x_1 - x_2)^3 = \epsilon \nu \xi(t') \quad (\text{E.3a}) $$

$$ \epsilon \ddot{x}_2 + \epsilon \mu (\dot{x}_2 - \dot{x}_1) + \epsilon \alpha (x_2 - x_1)^3 = 0 \quad (\text{E.3b}) $$

where $x_i = y_i/R$ ($i = 1, 2$), $\omega_1 = \sqrt{k_1/m_1}$, $t' = \omega_1 t$, $\{\}$ = $d\{\}/dt'$, $\epsilon = m_2/m_1$ is the mass ratio between the NES and the VdP oscillator, $\rho = c_1^{\text{NL}} R^2/(m_2\omega_1)$, $\mu = c_2/(m_2\omega_1)$, $\alpha = k_2^{\text{NL}} R^2/(m_2\omega_1^2)$ and $\nu = \hat{\nu}/(m_2 R \omega_1^{3/2})$. In the latter, the exponent 3/2 is obtained recalling that a normalized white noise $\xi(t)$ is defined as the time derivative of the normalized Wiener process $W(t')$ and using the scaling property of Wiener process (see [153], Chap. 2) which states that $W(t'/\omega_1)$ and $1/\sqrt{\omega_1}W(t')$ are the same stochastic processes. Therefore, we have

$$ \begin{aligned} \xi(t) &= \frac{dW(t)}{dt} \sim \omega_1 \frac{dW(\frac{t'}{\omega_1})}{dt'} \sim \sqrt{\omega_1} \frac{dW(t')}{dt'} \\ &= \sqrt{\omega_1} \xi(t') \end{aligned} \quad (\text{E.4}) $$

that explains the expression of parameter $\nu$.
Finally, using the change of variable $u_1 = x_1 + \epsilon x_2$ and $u_2 = x_1 - x_2$ the equations of motion (E.3) become

$$ \ddot{u}_1 + u_1 - \epsilon [\rho \dot{u}_1 (ru_1^2 - 1) + u_1 - u_2] = \epsilon \nu \xi(t') \quad (E.5a) $$

$$ \epsilon [\rho \dot{u}_1 (ru_1^2 - 1) + u_1 - \mu \dot{u}_2 - u_2 - \alpha u_2^3] = \epsilon \nu \xi(t') \quad (E.5b) $$

in which a first order Taylor expansion around $\epsilon = 0$ has been performed assuming a small mass ratio between the NES and the primary oscillator (i.e. $0 < \epsilon \ll 1$).

In the next sections, the time $t'$ is denoted by $t$ for the sake of conciseness.

## E.3 THE STOCHASTIC SLOW FLOW DYNAMICS

The standard stochastic averaging method [169, 191], whose general formulation is recalled in Appendix G.8, is used in this section to obtain the equations governing the stochastic slow flow dynamics of Eq. (E.5).

TET is due to the interaction between two non-linear modes of the coupled structure [76, 203]. This phenomenon, called a 1 : 1 resonance capture, occurs at a frequency close to the natural frequency of the primary structure, here the VdP oscillator. It is customary to study the dynamic behavior of the system in the neighborhood of this 1 : 1 resonance capture. Then, the system is simplified by averaging it over a natural period of the primary structure using the Krylov-Bogolyubov averaged approximation (see e.g. [145]). The resulting averaged dynamics is called slow flow. In general, in the context of NES studies, the complexification-averaging method [135, 206] is used. Here a real amplitude-angle representation (completely equivalent to the latter) is preferred by stating

$$ u_1 = a_1 \cos (t + \varphi_1) \quad (E.6a) $$

$$ \dot{u}_1 = -a_1 \sin (t + \varphi_1) \quad (E.6b) $$

and

$$ u_2 = a_2 \cos (t + \varphi_2) \quad (E.7a) $$

$$ \dot{u}_2 = -a_2 \sin (t + \varphi_2). \quad (E.7b) $$

The imposed forms of Eqs. (E.6b) and (E.7b) require that

$$ \dot{a}_1 \cos \phi_1 - a_1 \dot{\varphi}_1 \sin \phi_1 = 0, \quad (E.8) $$

and

$$ \dot{a}_2 \cos \phi_2 - a_2 \dot{\varphi}_2 \sin \phi_2 = 0, \quad (E.9) $$

with

$$ \phi_1 = t + \varphi_1 \quad \text{and} \quad \phi_2 = t + \varphi_2 \quad (E.10) $$

that leads to

$$ \dot{\varphi}_1 = \frac{\dot{a}_1 \cos \phi_1}{a_1 \sin \phi_1}, \quad \dot{a}_1 = a_1 \dot{\varphi}_1 \frac{\sin \phi_1}{\cos \phi_1}, \quad (E.11a) $$

$$ \dot{\varphi}_2 = \frac{\dot{a}_2 \cos \phi_2}{a_2 \sin \phi_2} \quad \text{and} \quad \dot{a}_2 = a_2 \dot{\varphi}_2 \frac{\sin \phi_2}{\cos \phi_2}. \quad (E.11b) $$

Then differentiation of Eqs. (E.6b) and (E.7b) yields

$$ \ddot{u}_1 = -a_1 \cos \phi_1 - \dot{a}_1 \sin \phi_1 - a_1 \dot{\varphi}_1 \cos \phi_1. \quad (E.12) $$

and

$$ \ddot{u}_2 = -a_2 \cos \phi_2 - \dot{a}_2 \sin \phi_2 - a_2 \dot{\varphi}_2 \cos \phi_2. \quad (E.13) $$

The substitution of Eqs. (E.6), (E.7), (E.12) and (E.13) into Eq. (E.5) and the use of Eq. (E.11) leads to the following (2, 2)-fast-slow non-autonomous system of stochastic differential equations

$$ \dot{a}_1 = \epsilon f_1(a_1, \varphi_1, a_2, \varphi_2, t) - \epsilon \nu \xi \sin \phi_1 \quad (E.14a) $$

$$ \dot{\varphi}_1 = \epsilon f_2(a_1, \varphi_1, a_2, \varphi_2, t) - \epsilon \nu \xi \frac{\cos \phi_1}{a_1} \quad (E.14b) $$
$$ \dot{a}_2 = f_3(a_1, \varphi_1, a_2, \varphi_2, t, \epsilon) - \epsilon\nu\xi \sin\phi_2 \quad (E.14c) $$

$$ \dot{\varphi}_2 = f_4(a_1, \varphi_1, a_2, \varphi_2, t, \epsilon) - \epsilon\nu\xi \frac{\cos\phi_2}{a_2} \quad (E.14d) $$

where the expressions of the functions $f_1, f_2, f_3$ and $f_4$ are given in Appendix E.8.

The terms of order $O(\epsilon)$ are neglected in the equations of the fast variables (i.e. Eqs. (E.14c) and (E.14d)), that leads to

$$ \dot{a}_1 = \epsilon f_1(a_1, \varphi_1, a_2, \varphi_2, t) - \epsilon\nu\xi \sin\phi_1 \quad (E.15a) $$

$$ \dot{\varphi}_1 = \epsilon f_2(a_1, \varphi_1, a_2, \varphi_2, t) - \epsilon\nu\xi \frac{\cos\phi_1}{a_1} \quad (E.15b) $$

$$ \dot{a}_2 = f_3(a_1, \varphi_1, a_2, \varphi_2, t, 0) \quad (E.15c) $$

$$ \dot{\varphi}_2 = f_4(a_1, \varphi_1, a_2, \varphi_2, t, 0) \quad (E.15d) $$

which has a similar form as (G.70) with

$$ \begin{aligned} \mathbf{x} &= (a_1, \varphi_1, a_2, \varphi_2)^T \\ \mathbf{f}(\mathbf{x}, t) &= \begin{pmatrix} \epsilon f_1(a_1, \varphi_1, a_2, \varphi_2, t) \\ \epsilon f_2(a_1, \varphi_1, a_2, \varphi_2, t) \\ f_3(a_1, \varphi_1, a_2, \varphi_2, t, 0) \\ f_4(a_1, \varphi_1, a_2, \varphi_2, t, 0) \end{pmatrix}, \quad (E.16) \\ \mathbf{g}(\mathbf{x}, t) &= \epsilon\nu \begin{pmatrix} -\sin\phi_1 & 0 \\ 0 & -\frac{\cos\phi_1}{a_1} \end{pmatrix}, \\ \mathbf{\eta} &= (\xi, \xi, \xi, \xi)^T \end{aligned} $$

where $()^T$ denotes the transpose operator.

Following the standard stochastic averaging method described in Appendix G.8, the drift vector **m** (see Eq. (G.72)) is first computed. The first term corresponds to classical (deterministic) Krylov-Bogolyubov averaging of the vector function **f**, here over one period equal to $2\pi$. This term is written as follows

$$ M^{\text{av}} \{\mathbf{f}\} = \begin{pmatrix} \epsilon \hat{f}_1(a_1, a_2, \varphi) \\ \epsilon \hat{f}_2(a_1, a_2, \varphi) \\ \hat{f}_3(a_1, a_2, \varphi) \\ \hat{f}_4(a_1, a_2, \varphi) \end{pmatrix} \quad (E.17) $$

where $\varphi = \varphi_1 - \varphi_2$ and

$$ \hat{f}_1(a_1, a_2, \varphi) = \frac{1}{8} (4a_2 \sin(\varphi) - a_1 (ra_1^2 - 4) \rho) \quad (E.18) $$

$$ \hat{f}_2(a_1, a_2, \varphi) = \frac{1}{2} \left( \frac{a_2 \cos(\varphi)}{a_1} - 1 \right) \quad (E.19) $$

$$ \hat{f}_3(a_1, a_2, \varphi) = \frac{1}{8} (-4a_2\mu - 4a_1 \sin(\varphi)) \quad (E.20) $$

$$ \hat{f}_4(a_1, a_2, \varphi) = \frac{3\alpha a_2^3 + 4a_1 \cos(\varphi) - 4a_2}{8a_2} \quad (E.21) $$

Then the second part of the drift vector **m** (see again Eq. (G.72)) is determined as

$$ \begin{aligned} & M^{\text{av}} \left\{ \int_{-\infty}^{0} \mathbb{E} \left[ \left( \frac{\partial(\mathbf{g}\eta)}{\partial \mathbf{x}} \right)_t (\mathbf{g}\eta)_{t+s} \right] ds \right\} \\ &= \frac{1}{2\pi} \int_{t_0}^{t_0+2\pi} \int_{-\infty}^{0} \mathbb{E} \left[ \left( \frac{\partial(\mathbf{g}\eta)}{\partial \mathbf{x}} \right)_t (\mathbf{g}\eta)_{t+s} \right] ds dt \\ &= \begin{pmatrix} \frac{\epsilon^2 \nu^2}{2a_1} \int_{-\infty}^{0} \delta(s) \cos(s) ds \\ 0 \\ 0 \\ 0 \end{pmatrix} \quad (E.22) \end{aligned} $$
in which the symbol $\mathbb{E}$ has been removed by means of Eq. (G.8). Because

$$ \int_{-\infty}^{0} \delta(s) \cos(s) ds = \frac{1}{2} \int_{-\infty}^{+\infty} \delta(s) \cos(s) ds = \frac{1}{2} \qquad (E.23) $$

the final expression of the drift vector is obtained as

$$ \mathbf{m} = \begin{pmatrix} \epsilon \hat{f}_1(a_1, a_2, \varphi) + \frac{\epsilon^2 \nu^2}{4a_1} \\ \epsilon \hat{f}_2(a_1, a_2, \varphi) \\ \hat{f}_3(a_1, a_2, \varphi) \\ \hat{f}_4(a_1, a_2, \varphi) \end{pmatrix}. \qquad (E.24) $$

The expression of the diffusion matrix $\sigma$ is now determined from (G.73). First, $\sigma\sigma^T$ is computed as follows

$$ \begin{aligned} \sigma\sigma^T &= M^{\text{av}} \left\{ \int_{-\infty}^{+\infty} \mathbb{E} \left[ (g\eta)_t (g\eta)_{t+s}^T \right] ds \right\} \\ &= \frac{1}{2\pi} \int_{t_0}^{t_0+2\pi} \int_{-\infty}^{+\infty} \mathbb{E} \left[ (g\eta)_t (g\eta)_{t+s}^T \right] ds dt \\ &= \epsilon^2 \nu^2 \int_{-\infty}^{+\infty} \left( \frac{\frac{\delta(s) \cos(s)}{2}}{2a_1} - \frac{\frac{\delta(s) \sin(s)}{2a_1}}{2a_1^2} \right) ds \\ &= \epsilon^2 \nu^2 \begin{pmatrix} \frac{1}{2} & 0 & 0 & 0 \\ 0 & \frac{1}{2a_1^2} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}. \end{aligned} \qquad (E.25) $$

Therefore, from (G.24) a possible solution of (G.73) is

$$ \sigma = \epsilon\nu \begin{pmatrix} \frac{1}{\sqrt{2}} & 0 & 0 & 0 \\ 0 & \frac{1}{a_1\sqrt{2}} & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}. \qquad (E.26) $$

According to (G.71) the equations governing the stochastic slow dynamics of (E.5) are

$$ \dot{a}_1 = \epsilon \hat{f}_1(a_1, a_2, \varphi) + \frac{\epsilon^2 \nu^2}{4a_1} + \frac{\epsilon\nu}{\sqrt{2}} \xi(t) \qquad (E.27a) $$

$$ \dot{\varphi}_1 = \epsilon \hat{f}_2(a_1, a_2, \varphi) + \frac{\epsilon\nu}{a_1\sqrt{2}} \xi(t) \qquad (E.27b) $$

$$ \dot{a}_2 = \hat{f}_3(a_1, a_2, \varphi) \qquad (E.27c) $$

$$ \dot{\varphi}_2 = \hat{f}_4(a_1, a_2, \varphi) \qquad (E.27d) $$

The final form of the stochastic slow flow is obtained by combining Eqs. (E.27b) and (E.27d), that leads to

$$ \dot{a}_1 = \epsilon f(a_1, a_2, \varphi) + \epsilon \sigma \xi(t) \qquad (E.28a) $$

$$ \dot{a}_2 = g_1(a_1, a_2, \varphi) \qquad (E.28b) $$

$$ \dot{\varphi} = g_2(a_1, a_2, \varphi) \qquad (E.28c) $$

where $f = \hat{f}_1$, $g_1 = \hat{f}_3$, $\sigma = \frac{\nu}{\sqrt{2}}$ and the term $\frac{\epsilon^2 \nu^2}{4a_1}$ (of order $O(\epsilon^2)$) has been neglected in Eq. (E.28a). The difference between Eqs. (E.27b) and (E.27d) should yield $\epsilon \hat{f}_2(a_1, a_2, \varphi) + \frac{\epsilon\nu}{a_1\sqrt{2}}\xi - \hat{f}_4(a_1, a_2, \varphi)$ as the right-hand side of Eq. (E.28c). However, the terms of order $O(\epsilon)$ are neglected, one has therefore only $g_2 = -\hat{f}_4(a_1, a_2, \varphi)$ remains. Eq. (E.28) appears as a stochastic (2, 1)-fast-slow system where $a_2$ and $\varphi$ are the fast variables and $a_1$ is the slow variable on which the noise acts.
A key mathematical tool for the description of a fast-slow system such as Eq. (E.28) is its critical manifold. The latter is defined in this section. First, the slow flow (E.28) is written with respect to the slow time $\tau = \epsilon t$ as follows

$$a'_{1} = f(a_{1}, a_{2}, \varphi) + \sqrt{\epsilon\sigma\xi(\tau)} \qquad (E.29a)$$

$$\epsilon a'_{2} = g_{1}(a_{1}, a_{2}, \varphi) \qquad (E.29b)$$

$$\epsilon\varphi' = g_{2}(a_{1}, a_{2}, \varphi) \qquad (E.29c)$$

where $(\cdot)'$ denotes the derivative with respect to the slow time $\tau$. Again the scaling property of the Wiener process has been used, i.e. $\xi(t) = \sqrt{\epsilon\xi(\tau)}$. Considering $\epsilon = 0$ respectively in Eqs. (E.28) and (E.29) yields the slow subsystem

$$a'_{1} = f(a_{1}, a_{2}, \varphi) \qquad (E.30a)$$

$$0 = g_{1}(a_{1}, a_{2}, \varphi) \qquad (E.30b)$$

$$0 = g_{2}(a_{1}, a_{2}, \varphi), \qquad (E.30c)$$

which is a differential-algebraic equation, and the fast subsystem

$$\dot{a}_{1} = 0 \qquad (E.31a)$$

$$\dot{a}_{2} = g_{1}(a_{1}, a_{2}, \varphi) \qquad (E.31b)$$

$$\dot{\varphi} = g_{2}(a_{1}, a_{2}, \varphi). \qquad (E.31c)$$

The critical manifold of the slow flow is the solution of the algebraic part of Eq. (E.30) and it is expressed as follows

$$M_{0} = \left\{ (a_{1}, a_{2}, \varphi) \in \mathbb{R}^{+2} \times [-\pi, \pi] \mid g_{1}(a_{1}, a_{2}, \varphi) = 0 \text{ and } g_{2}(a_{1}, a_{2}, \varphi) = 0 \right\}. \qquad (E.32)$$

Note that because of the assumptions made (i.e. a noise which acts only on the primary structure with a level of the order $O(\epsilon)$), the effect of the noise on the fast variables $a_2$ and $\varphi$ has been neglected. Consequently both slow and fast subsystems are deterministic. This has two major consequences. First, the critical manifold is deterministic too and it can therefore be obtained with the usual approach (the classical analysis of the deterministic slow flow including the computation of the critical manifold is recalled in Appendix E.10). Then, it means that the role of the noise cannot be investigated within the zeroth-order approximation, i.e. within the limit case in which the perturbation parameter $\epsilon$ is equal to zero, as it is usually done in the context of NES studies.

## E.4 PROBLEM STATEMENT AND PRELIMINARY RESULTS

In this section we first recall in Section E.4.1 how the deterministic system behaves using both numerical simulations of the deterministic slow flow (i.e. Eq. (E.29) with $\sigma = 0$) and the classical analytical results presented in Appendix E.10. Then, in Section E.4.2, the comparison with numerical integration of the stochastic slow flow (E.29) (using the function ItoProcess of Mathematica software [216]) allows us to highlight the influence of noise on the system behavior.

### E.4.1 Reminder of the deterministic behavior of the system

The slow-fast nature of the slow flow implies that it evolves at two time scales : the slow time scale $\tau$ in which the slow flow is on the critical manifold and approximatively described by the slow subsystem (E.30) and the fast time scale $t$ in which the slow flow is outside the critical manifold and described by the fast subsystem (E.31). The particular S-shape of the critical manifold (see Fig. E.9) together with the stability analysis of the fixed points of the slow flow allows to explain and predict its different responses and consequently those of the full-order system. Four scenarios are possible and in previous works by the authors [a7, a8] these responses are classified into two categories called *harmless situations* and *harmful situations*. In harmless situations, the NES acts, resulting in three possible responses. The first scenario is called *Complete suppression*. In this case the trivial fixed point of the slow flow is stable and then reached. This is a linear effect of the NES and, in the context of non-linear vibrations absorption, the complete suppression is not the desired effect. The second possible regime is called *mitigation through periodic responses (PRs)*. This corresponds to the situation in which a nontrivial stable fixed point of the slow flow is reached leading to a periodic regime for the initial full-order system. The last harmless situations correspond to *mitigation through Strongly Modulated Responses (SMRs)*. SMRs are quasi-periodic regimes (amplitude and phase modulated) of the full-order system
<figure><img src="image_107.png" /><figcaption>**FIGURE E.2.** Transition from relaxation oscillations scenario (Figs. E.2(a) and E.2(c)) to no mitigation scenario (Figs. E.2(b) and E.2(d)). The top figures show the times series (obtained from numerical integration of the deterministic slow flow, i.e. Eq. (E.29) with $\sigma = 0$) of the amplitudes $a_1$ (solid black lines) and $a_2$ (solid gray lines) at the slow time scale $\tau$ and the bottom figures show the comparison between the trajectory of the slow flow (solid red lines) and the critical manifold (E.32) in the ($a_2, a_1$)-plane. The following set of parameters is used: $r = 0.1$, $\epsilon = 0.01$, $\mu = 0.4$, $\alpha = 3$ and $\rho = 1.8$ in (a) and (c) and $\rho = 2$ in (b) and (d).</figcaption></figure>

corresponding to relaxations oscillations of the slow flow. There is only one type of harmful situations for which the NES is not able to produce small amplitude responses. In this case, a limit cycle with an amplitude close to that of the Van der Pol oscillator alone is observed.

In general, the transition from harmless to harmful situations corresponds to the transition from an SMR to a no mitigation regime. This transition is illustrated in Fig. E.2: Figs. E.2(a) and E.2(c) show a relaxation oscillation scenario and Figs. E.2(b) and E.2(d) show a no mitigation scenario. The top figures depict the times series of the amplitude $a_1$ and $a_2$ at the slow time scale $\tau$ and the bottom figures depict the comparison between the trajectory of the slow flow and the critical manifold (E.32) in the ($a_2, a_1$)-plane. In the latter (see Figs. E.2(c) and E.2(d)), in both relaxation oscillations and no mitigation scenarios, from an initial condition near zero and outside the critical manifold, the trajectory of the slow flow evolves rapidly and almost horizontally to the left attracting branch of $M_0$. During this fast epoch the slow flow dynamics is approximatively described by the fast subsystem (E.31). Then, the slow flow evolves slowly close to this branch. During this slow epoch the slow flow dynamics is now approximatively described by the slow subsystem (E.30). If the trajectory does not meet a stable fixed point<fnref n="1" />, the trajectory goes toward the left fold point ($a_2^{LF}, a_1^{LF}$) at which $M_0$ becomes repelling (see Eqs. (E.49), (E.51) and (E.52)).

**Remark.** When the full-order system passes from a PR to an SMR (i.e. when the slow flow passes from a reached stable fixed point to stable relaxation oscillations) *folded singularities* can appear. The latter are points for which fold and fixed points of the slow flow coincide and are hints of particular solutions of fast-slow systems, called *canards* (see e.g. [54] for rigorous mathematical definitions of these concepts). In general, canard solutions appear for a very small range of the bifurcation parameter and are not investigated in this

<footnote n="1">1. If this is the case the trajectory stops and a mitigation through a PR is observed.</footnote>
paper.

The zeroth-order approximation assumes that the left fold point is actually reached. However numerical simulations show that when approaching the left fold point the trajectory deviates slightly from the critical manifold (this behavior is analytically described in [a2]). Then the slow flow undergoes a fast and nearly horizontal jump to the right attracting part of $M_0$ on a point called arrival point$^2$. This is where the two scenarios differ. In the case of relaxation oscillations the largest unstable fixed point of the slow flow (depicted by a green bullet on the figures) is above the arrival point and prevents the trajectory from reaching the stable fixed point located higher on the critical manifold. This fixed point is not depicted on Figs. E.2(c) and E.2(d), its coordinates in the $(a_2, a_1)$-plane are (1.51, 6.28). If this fixed point is reached the system undergoes a no mitigation regime. This is what happens when the largest unstable fixed point of the slow flow is below the arrival point (see Fig. E.2(d)). The end of the scenario of relaxation oscillations is as follows (see Fig. E.2(c)): from the arrival point a second slow epoch occurs in which the trajectory goes toward the right fold point $(a_2^{\text{RF}}, a_1^{\text{RF}})$. A second nearly horizontal jump occurs and the trajectory returns to the left attracting part of $M_0$. A third slow epoch on $M_0$ towards $(a_2^{\text{LF}}, a_1^{\text{LF}})$ is observed, and so on.

Then the **mitigation limit**, defined below, is introduced to quantify the NES efficiency.

**DEFINITION E.1 (MITIGATION LIMIT).** Considering a set of initial conditions (for the slow flow) as a small perturbation of the trivial solution, the **mitigation limit** is defined as the value of the bifurcation parameter $\rho$ which separates harmless situations from the harmful situation.

For the set of parameters used in Fig. E.2, using numerical simulations of the slow flow we find a mitigation limit equal to $\rho_{\text{ml}} = 1.9$. The numerical simulations of the deterministic full-order system (i.e. Eq. (E.5) with $\nu = 0$) give $\rho_{\text{ml}} = 1.86$.

In the deterministic case the mitigation limit can be predicted analytically by means of a slow-fast partition of the slow flow dynamics [a2, 75].

## E.4.2 Influence of noise on the system behavior

In this section we highlight how noise affects the system behavior and more precisely how it influences the mitigation limit. For this purpose, we show in Fig. E.3 the result of the numerical simulations of the stochastic slow flow (E.29) using the same parameters as in Section E.4.1 with, in addition, $\sigma = 0.5$.

In Figs. E.3(a), E.3(b), E.3(e) and E.3(f), two samples are shown for $\rho = 1.8$, a value of the bifurcation parameter smaller than the mitigation limit $\rho_{\text{ml}} = 1.9$. Therefore, in the deterministic case, we observe persistent relaxation oscillations (see Figs. E.2(a) and E.2(c)). In the deterministic case all cycles of relaxation oscillations are equivalent, whereas in the stochastic case the trajectory of the slow flow is, due to the presence of noise, slightly different at each cycle. The two samples depicted in the figures show that persistent relaxation oscillations can still occur. We can see in Fig. E.3(e) that over the 8 cycles observed (see Fig. E.3(a)), 7 correspond to a scenario similar to the deterministic case (i.e. the arrival point on the right attracting part of $M_0$ is below the largest fixed point of the deterministic slow flow). There is however one cycle where the arrival point is above the fixed point and where the trajectory still descends towards the right fold point instead of ascending towards the stable fixed point. For the second sample (see Figs. E.3(b) and E.3(f)) it is the contrary. Indeed, we observe one and a half cycle, and halfway through the second we see that even if the arrival point is below the unstable fixed point the trajectory ascends towards the stable fixed point and a no mitigation regime finally occurs.

In Figs. E.3(c), E.3(d), E.3(g) and E.3(h) two samples are depicted for $\rho = 2$, a value of the bifurcation parameter larger than the mitigation limit $\rho_{\text{ml}} = 1.9$. Consequently, in the deterministic case, we observe a no mitigation regime (see Figs. E.2(b) and E.2(d)). The first sample (see Figs. E.3(c) and E.3(g)) shows that during the first cycle, the noise prevents the trajectory from reaching the stable fixed point corresponding to the no mitigation regime, which is finally reached in the second cycle. In this example, at each cycle the regime actually reached can be deduced in a similar way as in deterministic case comparing the relative position of the arrival point and the largest unstable fixed point. However, as in the previous example for $\rho = 1.8$, this reasoning can be contradicted in a positive way (favoring the relaxation oscillations) or in a negative way (favoring the no mitigation regimes). For example, in the second sample (see Figs. E.3(d) and E.3(h)) we can see that during the first cycle, the arrival point is above the largest fixed point but the trajectory still descends towards the right fold point.

<footnote n="2">2. The zeroth-order approximation supposes that the arrival point is $(a_2^U, a_1^{\text{LF}})$ (see Eqs. (E.53) and (E.9)). In fact, the arrival point is slightly higher on the critical manifold [a2].</footnote>
<figure><img src="image_108.png" /></figure>

**FIGURE E.3.** Numerical simulations of the stochastic slow flow (E.29) using the same parameters as in Section E.4.1 with, in addition, $\sigma = 0.5$. In Figs. E.3(a), E.3(b), E.3(e) and E.3(f) (resp. Figs. E.3(c), E.3(d), E.3(g) and E.3(h)) two samples are shown for $\rho = 1.8$ (resp. $\rho = 2$). Figs. E.3(a), E.3(b), E.3(c) and E.3(d) show the time series of $a_1$ (black lines) and $a_2$ (gray lines) and Figs. E.3(e), E.3(f), E.3(g) and E.3(h) show the trajectory (red lines) of the slow flow and the critical manifold $M_0$ (gray and black lines) in the $(a_1, a_2)$-plane.

The previous preliminary results confirm that the effect of the noise cannot be only studied with the zeroth-order approximation. Indeed, as in the deterministic case the trajectory of the slow flow dynamics moves away from the critical manifold when it passes near the left fold point [a2, cn1]. The zeroth-order approximation [75] assumed that the trajectory goes along the critical manifold and reaches the left jump before jumping to the right attractive part of the critical manifold. However, in the deterministic case the zeroth-order approximation is sufficient to describe the slow epochs of the slow flow except in the neighborhood of the left fold points. Combining the center manifold reduction technique to describe the slow flow dynamics near the left fold point and the zeroth-order approximation elsewhere, and following the reasoning described in Section E.4.1, allows an accurate prediction of the mitigation limit [a2]. Here, in the stochastic case, the latter reasoning can be wrong to predict the actual regime and therefore the mitigation limit. That means that the noise acts at the slow time scale in contradiction with the zeroth-order approximation (see Eqs. (E.30) E.31 and the comments below).

In the next sections the Monte Carlo method is used to study the effect of the noise on the dynamic behavior of the system in a more systematic way. To achieve that, two quantities are introduced : (i) the *probability of being in a harmless regime* and (ii) the *First-Passage Time to reach a harmful regimes*.

## E.5 PROBABILITY OF BEING IN A HARMLESS REGIME

**DEFINITION E.2 (THE PROBABILITY OF BEING IN A HARMLESS REGIME).** The probability of being in a harmless regime (PBHR), denoted by $p_{h,n}$, is, as its name suggests, the probability for the system of being in a harmless regime after a given number $n$ of full cycles of relaxation oscillations.

For instance, in the case shown in Figs. E.3(b) and E.3(f), the trajectory of the slow flow ascends towards the stable fixed point (to reach a no mitigation regime) halfway through the second cycle of relaxation oscillations, therefore after one full cycle. In practice, $p_{h,n}$ is computed as the proportion of samples for which we observe at least $n + 1$ consecutive full cycles of relaxation oscillations from the beginning of the sample. For example, if 4 samples are computed and we observe : 1 full cycle in the first sample (as in Figs. E.3(b), E.3(c) a,d E.3(d)), 0 full cycle in the second sample (as in the deterministic case shown in Fig. E.2(b)) and 3 and 2
<figure><img src="image_109.png" /><figcaption>FIGURE E.4. Probabilities $p_{h,0}$ as functions of $\rho$ obtained from 1000 samples of the slow flow (E.29) (blue squares) and 1000 samples of the full order system (E.14) (red bullets) and for several values of the bifurcation parameter $\rho$. The deterministic cases (green dashed line for the slow flow and black line for the full-order system) are also represented. The same parameters as in Section E.4.1 are used with, in addition, $\sigma = 0.5$.</figcaption></figure>

<figure><img src="image_110.png" /><figcaption>FIGURE E.5. Probabilities $p_{h,0}$ as functions of $\rho$ for three different values of the noise level, i.e. $\sigma = 0.2$ (blue bullets), $\sigma = 0.5$ (red squares) and $\sigma = 0.8$ (green diamonds) and for the same parameters as in Section E.4.1.</figcaption></figure>

full cycles in the third and fourth samples respectively, then the corresponding PBHR are $p_{h,0} = \frac{3}{4}$, $p_{h,1} = \frac{2}{4}$, $p_{h,2} = \frac{1}{4}$ and $p_{h,n} = 0$ for $n \ge 3$.

Fig. E.4 shows the PBHR $p_{h,0}$ obtained from 1000 samples of the slow flow (E.29) and 1000 samples of the full-order system (E.14) and for several values of the bifurcation parameter $\rho$. The deterministic cases for which $p_{h,0} = 1$ if $\rho < \rho_{ml}$ and $p_{h,0} = 0$ if $\rho > \rho_{ml}$ are also represented (remember that $\rho_{ml} = 1.86$ for the full order system and $\rho_{ml} = 1.9$ for the slow flow). The same parameters as in Section E.4.1 are used with, in addition, $\sigma = 0.5$.

First, the figure validates the stochastic averaging procedure. Therefore, in the remaining of the paper, only the slow flow will be used, which allows us to perform the Monte Carlo method with a reasonable computational cost. Indeed with the full-order system, the convergence of the numerical results is obtained with a very small time step, which is highly expensive in terms of computing time.

Then we can see that the step shape observed in the deterministic case, turns into an error function shape in the stochastic case, with $p_{h,0}$ tending to 1 when $\rho \ll \rho_{ml}$ and to 0 when $\rho \gg \rho_{ml}$.

Fig. E.5 plots $p_{h,0}$ for three different values of the noise level, i.e. $\sigma = 0.2$, $\sigma = 0.5$ and $\sigma = 0.8$. As expected, the lower the noise level is, the closer of the deterministic case the system behavior is.

Finally, in Fig. E.6 we can see the comparison between $p_{h,0}$, $p_{h,1}$ and $p_{h,2}$. First, as expected, we have $p_{h,2} < p_{h,1} < p_{h,0}$ (see Fig. E.6(a)). Then, we observe that $p_{h,1} \approx p_{h,0}^2$ and $p_{h,2} \approx p_{h,0}^3 \approx p_{h,0}p_{h,1}$ (see Fig. E.6(b)). From these numerical results (and others not shown here) the following generalization is proposed:

$$p_{h,n} \approx \prod_{i=0}^{n-1} p_{h,i} \approx p_{h,0}^{n+1}. \qquad (E.33)$$
<figure><img src="image_111.png" /><figcaption>FIGURE E.6. (a) Comparison between $p_{h,0}$ (blue bullets), $p_{h,1}$ (red squares) and $p_{h,2}$ (green diamonds) as functions of the bifurcation parameter $\rho$. (b) Comparison between $p_{h,0}^2$ (blue bullets), $p_{h,1}$ (red squares), and comparison between $p_{h,0}^3$ (green diamonds) and $p_{h,2}$ (orange triangles). In each figure the deterministic cases (black line) is also represented. The same parameters as in Section E.4.1 are used with, in addition, $\sigma = 0.5$.</figcaption></figure>

This result highlights that at each cycle of relaxation oscillations the probability of being in a harmless regime is independent of the past trajectory of the system (we verified that this also true for the full-order system) and it is equal to $p_{h,0}$.

The PBHR quantifies, in terms of probability, one of the results previously obtained in Section E.4.2, namely the fact that no mitigation (resp. mitigation) regimes may occur even if the $\rho < \rho_{ml}$ (resp. $\rho > \rho_{ml}$). It is interesting to know the probability to be in a harmless regime after a given number of cycles of relaxation oscillations as the value of $p_{h,n}$ tells us. However, the time spent in a harmless regime, i.e. the time the system is safe, is also interesting from an engineering point of view. This is quantified in the next section by computing the First-Passage Time Density to reach a harmful regime.

## E.6 FIRST-PASSAGE TIME TO REACH A HARMFUL REGIME

Let us denote by $(a_1^e, a_2^e, \varphi^e)$ the coordinates of the largest stable fixed point of the deterministic slow flow. Remember that if the trajectory of the slow flow reaches the fixed point the system is in a no mitigation (i.e. harmful) regime (see Section E.4.1). Then, we choose that the system is safe as long as $a_1 < a_1^{th} = 0.5a_1^e$. This choice is arbitrary, it may change depending of the unstable primary structure studied. Next, the First-Passage Time to reach a harmful regime is defined below.

**DEFINITION E.3 (FIRST-PASSAGE TIME TO REACH A HARMFUL REGIME).** If the process starts at $a_1(0) = a_{1,0}$, the first time it reaches the value of the threshold $a_1^{th}$ is called the First-Passage Time to reach a Harmful Regime (FPTHR) and it is denoted by $T$.

At the slow time scale $\tau$ and from the previous definition, the probability $Pr(\tau < T)$ that $\tau < T$ is the probability $Pr(a_1 < a_1^{th})$ that $a_1 < a_1^{th}$. Therefore, one also has

$$ \mathrm{Pr}(T < \tau) = \mathrm{Pr}(a_1 > a_1^{\mathrm{th}}) := P(\tau). \qquad (\mathrm{E.34}) $$

Consequently, the probability that the FPTHR $T$ is between the actual time $\tau$ and $\tau + d\tau$ is $R(\tau)d\tau$ where

$$ R(\tau) = \frac{dP(\tau)}{d\tau} \qquad (\mathrm{E.35}) $$

is the First-Passage Time to reach a Harmful Regime Density (FPTHRD).

The Mean First-Passage Time to reach a Harmful Regime (MFTHR) denoted by $\mathbb{E}[T]$ is also defined as

$$ \mathbb{E}[T] = \int_{0}^{\infty} \tau R(\tau) d\tau. \qquad (\mathrm{E.36}) $$
<figure><img src="image_112.png" /><figcaption>**FIGURE E.7.** Histograms representing $R(\tau)$ (Figs. E.7(a) and E.7(c)) and $P(\tau)$ (Figs. E.7(b) and E.7(d)). Vertical lines of equations $\tau = T_{\text{det}}$ (blue lines), $\tau = T_{\text{LF}}$ (red lines) and $\tau = \mathbb{E}[T]$ (green lines) are also represented. The same parameters as in Section E.4.2 are used with $\rho = 1.8 < \rho_{\text{ml}}$ (Figs. E.7(a) and E.7(b)) and $\rho = 2 > \rho_{\text{ml}}$ (Figs. E.7(c) and E.7(d)).</figcaption></figure>

In practice, 2000 samples of Eq. (E.29) are computed. On each sample, the FPTHR is measured as the first value of the time $\tau$ for which $a_1(\tau) > a_1^{\text{th}}$. We obtain therefore a list of 2000 values the FPTHR and from it $R(\tau)$, $P(\tau)$ (as histograms) and $\mathbb{E}[T]$ are evaluated.

The Deterministic Passage Time to reach a Harmful Regime (DPTHR), denoted by $T_{\text{det}}$, is defined as the first time $a_1(\tau)$, obtained from the numerical integration of Eq. (E.29) for $\sigma = 0$ and $a_1(0) = a_{1,0}$, reaches the threshold $a_1^{\text{th}}$. For given initial conditions, $T_{\text{det}}$ has a unique finite value if $\rho > \rho_{\text{ml}}$ and tends to infinity if $\rho < \rho_{\text{ml}}$ (i.e. the deterministic system never reaches the threshold $a_1^{\text{th}}$).

Finally, for comparison purposes, a characteristic time is needed. We choose the time, denoted by $T_{\text{LF}}$, at which the deterministic system reaches the neighborhood of the left fold point of the critical manifold (i.e. when $a_1(\tau) = a_1^{\text{LF}}$, obtained from the numerical integration of Eq. (E.29) for $\sigma = 0$ and $a_1(0) = a_{1,0}$, reaches $a_1^{\text{LF}}$).

Examples of histograms representing $R(\tau)$ and $P(\tau)$ are shown in Fig. E.8 (on the left and the right respectively). The same parameters are used as in Section E.4.2 with $\rho = 1.8 < \rho_{\text{ml}}$ (see Figs. E.7(a) and E.7(b)) and $\rho = 2 > \rho_{\text{ml}}$ (see Figs. E.7(c) and E.7(d)). Vertical lines of equations $\tau = T_{\text{det}}$, $\tau = T_{\text{LF}}$ and $\tau = \mathbb{E}[T]$ are also represented. In Figs. E.7(a) and E.7(b) the deterministic system is in a mitigation regime and in theory $T_{\text{det}} \to +\infty$. In practice it is equal to the final time of the numerical simulation, here the final time is $\tau = 140$ which is outside the frame of the figures. The numerical values of the important times are $T_{\text{det}} = 60$, $T_{\text{LF}} = 3.85$ and $\mathbb{E}[T] = 12$ for $\rho = 1.8$ and $T_{\text{det}} = 7.86$, $T_{\text{LF}} = 3.37$ and $\mathbb{E}[T] = 8.11$ for $\rho = 2$.

We can see in Figs. E.7(a) and E.7(b) that even if $\rho < \rho_{\text{ml}}$ (i.e. the corresponding deterministic system undergoes a harmless regime), it is almost certain that if we wait for a sufficiently long time the threshold $a_1^{\text{th}}$ will be reached. Even for relatively short times (approximately 2 to 4 times $T_{\text{LF}}$) the probability of reaching the threshold is not negligible. To complete the information given in Fig. E.7, the values of $P(6)$, $P(10)$ and $P(14)$ for $\rho = 1.7, 1.6, \dots, 2.4$ are given in Table. E.1. For example (see Fig. E.7(b) and Table. E.1), we have $P(10) = 0.433$, i.e. if $\tau = 10$ we have a 43.3% chance that the threshold has already been reached. Not surprisingly, for a given value of $\rho$, the larger the time $\tau$ is, the larger $P(\tau)$ becomes; for a given value of $\tau$, the larger $\rho$ is, the larger $P(\tau)$ becomes again. Moreover, these results (those shown in Figs. E.7(a) and E.7(c)) show that the most likely value is not the deterministic value $\tau = T_{\text{det}}$ which is very close to the mean time
<figure><img src="image_113.png" /><figcaption>**FIGURE E.8.** Comparison between $T_{det}$ (blue line), $\mathbb{E}[T]$ (green dashed line), $\mathbb{E}[T] - \sqrt{\mathbb{V}[T]}$ (red line) and $\mathbb{E}[T] + \sqrt{\mathbb{V}[T]}$ (orange line) as function of the bifurcation parameter $\rho$. The same parameters as in Section E.4.2 are used.</figcaption></figure>

$\mathbb{E}[T]$ (see Fig. E.7(c)). Contrary to the deterministic case, in the stochastic case we do not see any qualitative difference between the results obtained for $\rho = 1.8$ and those obtained for $\rho = 2$.

<table><caption>**TABLEAU E.1.** Values of $P(6)$, $P(10)$ and $P(14)$ for $\rho = 1.7, 1.6, \dots, 2.4$. The same parameters are used as in Fig. E.7.</caption><thead><tr><th>ρ</th><th>1.7</th><th>1.8</th><th>1.9</th><th>2</th><th>2.1</th><th>2.2</th><th>2.3</th><th>2.4</th></tr></thead><tbody><tr><th>P(6)</th><td>0.003</td><td>0.013</td><td>0.051</td><td>0.144</td><td>0.276</td><td>0.425</td><td>0.563</td><td>0.706</td></tr><tr><th>P(10)</th><td>0.219</td><td>0.433</td><td>0.652</td><td>0.828</td><td>0.932</td><td>0.973</td><td>0.985</td><td>0.997</td></tr><tr><th>P(14)</th><td>0.5</td><td>0.745</td><td>0.913</td><td>0.977</td><td>0.998</td><td>1</td><td>1</td><td>1</td></tr></tbody></table>

The latter observation is even more highlighted in Fig. E.8 which compares $T_{det}$, $\mathbb{E}[T]$ and $\mathbb{E}[T] \pm \sqrt{\mathbb{V}[T]}$ where $\sqrt{\mathbb{V}[T]}$ is the standard deviation of the FPTHR with

$$ \mathbb{V}[T] = \mathbb{E}[T^2] - \mathbb{E}[T]^2 \qquad (E.37) $$

the variance. Indeed, we can see a discontinuity in the graph of $T_{det}$ at $\rho = \rho_{ml}$ which does not appear in the graphs of $\mathbb{E}[T]$ and $\mathbb{E}[T] \pm \sqrt{\mathbb{V}[T]}$. Again, note that the constant value of the DPTHR, i.e. $T_{det} = 140$ for $\rho < \rho_{ml}$, is a numerical artefact. Indeed, when the deterministic system is in a mitigation regimes (here relaxation oscillations of the slow flow) we have $T_{det} \to +\infty$ in theory. In practice, $T_{det}$ is detected as the final time of the numerical integration which is here $\tau = 140$. Both mean and standard deviation increase when $\rho$ decreases, moving away from the deterministic mitigation limit $\rho_{ml}$. Finally, we can see again that when $\rho > \rho_{ml}$ one has $\mathbb{E}[T] \approx T_{det}$.

In my opinion, what should be retained from the results of this section is the fact that even if the system is supposed to be in an mitigation regime from a deterministic point of view, the amplitude of the self-sustained oscillator can reach a harmful value in a relatively short time with a non-negligible probability.

## E.7 CONCLUSION

In this paper a Van der Pol oscillator (used as an archetypal self-sustained oscillator) undergoing a white noise forcing and coupled to a non-linear energy sink (NES) has been studied. It was highlighted that the stochastic forcing can modify significantly the dynamic behavior of the corresponding deterministic system. Indeed we showed first, by means of numerical integration of the slow flow dynamics of the system - derived
using the standard stochastic averaging method - that the reasoning which allows the theoretical prediction of the system behavior in the deterministic case can be contradicted in presence of stochasticity.

Then the system behavior has been studied in a more systematic way by introducing two quantities : (i) the probability of being in a harmless regime (PBHR), i.e. the probability of obtaining a regime of small amplitude compared to the case without NES; and (ii) the first-passage time to reach a harmful regime (FPTHR), i.e. the time for the system to reach an arbitrary amplitude threshold considered as dangerous if reached by the primary structure. Both quantities have been investigated combining the Monte Carlo approach with numerical integration of the stochastic slow flow dynamics. When they are computed on the corresponding deterministic regime, the PBHR and FPTHR reveal themselves to be discontinuous functions of the bifurcation parameter under consideration. Indeed, at the mitigation limit (i.e. the value of the bifurcation parameter below which the NES acts and above which it no longer acts) the PBHR goes abruptly from 1 (we are sure to be in a mitigation regime) to 0 (we are sure to be in a no mitigation regime), and the FPTHR from infinity (i.e. the deterministic system never reaches the no mitigation regime) to a finite value. We showed that in the stochastic case, both PBHR and FPTHR become smooth at the mitigation limit. This work highlights that taking into account stochasticity allows a better understanding of targeted energy transfer mechanisms in the context of limit cycles mitigation by means of an NES under more realistic conditions than in a deterministic study.

An interesting sequel to this research would be the theoretical analysis of stochastic slow flow inspired by methodologies dedicated to the analysis of noisy slow-fast systems as developed by Berglund and Gentz [22, 24]. This would help to understand the mechanisms underlying the phenomena highlighted in this article.

## E.8 APPENDIX A. EXPRESSIONS OF FUNCTIONS

The expressions of the functions $f_1$, $f_2$, $f_3$ and $f_4$ which appear in Eq. (E.14) are

$$f_1(a_1, \varphi_1, a_2, \varphi_2, t) = \frac{1}{2} \left( a_1 - a_1 \rho \cos(2t + 2\varphi_1) - a_1 \sin(2t + 2\varphi_1) \right. \\ \left. + a_2 \sin(2t + \varphi_1 + \varphi_2) + a_2 \sin(\varphi_1 - \varphi_2) \right. \\ \left. - a_1^3 \rho r + a_1^3 \rho r \cos(4t + 4\varphi_1) \right), \qquad (E.38)$$

$$f_2(a_1, \varphi_1, a_2, \varphi_2, t) = \frac{1}{2a_1} \left( -a_1 + 4a_1 \rho \sin(2t + 2\varphi_1) - a_1 \cos(2t + 2\varphi_1) \right. \\ \left. + a_2 \cos(2t + \varphi_1 + \varphi_2) + a_2 \cos(\varphi_1 - \varphi_2) \right. \\ \left. - 2a_1^3 \rho r \sin(2t + 2\varphi_1) + a_1^3 \rho r \sin(4t + 4\varphi_1) \right), \qquad (E.39)$$

$$f_3(a_1, \varphi_1, a_2, \varphi_2, t, \epsilon) = \frac{1}{8} \left( -4a_2\mu + 4a_1\sin(2t + \varphi_1 + \varphi_2) - 4a_1\sin(\varphi_1 - \varphi_2) \right. \\ + 4a_2\mu\cos(2t + 2\varphi_2) - 4a_2\sin(2t + 2\varphi_2) \\ + 2\alpha a_2^3 \sin(2t + 2\varphi_2) + \alpha a_2^3 \sin(4t + 4\varphi_2) \\ + \frac{\epsilon}{8} \left[ 4a_1\rho\cos(\varphi_1 - \varphi_2) - 4a_1\rho\cos(2t + \varphi_1 + \varphi_2) \right. \\ + 4a_1\sin(\varphi_1 - \varphi_2) - 4a_1\sin(2t + \varphi_1 + \varphi_2) \\ - 4a_2\mu + 4a_2\mu\cos(2t + 2\varphi_2) + 4a_2\sin(2t + 2\varphi_2) \\ + 2\alpha a_2^3 \sin(2t + 2\varphi_2) + \alpha a_2^3 \sin(4t + 4\varphi_2) \\ - a_1^3\rho r \cos(\varphi_1 - \varphi_2) - a_1^3\rho r \cos(2t + 3\phi_1 - \phi_2) \\ + a_1^3\rho r \cos(2t + \phi_1 + \phi_2) + a_1^3\rho r \cos(4t + 3\phi_1 + \phi_2) \right] \qquad (E.40)$$

and

$$f_4(a_1, \varphi_1, a_2, \varphi_2, t, \epsilon) = \frac{1}{8a_2} \left( -4a_2 + 3\alpha a_2^3 + 4\alpha a_2^3 \cos(2t + 2\varphi_2) + \alpha a_2^3 \cos(4t + 4\varphi_2) \right. \\ + 4a_1 \cos(2t + \varphi_1 + \varphi_2) + 4a_1 \cos(\varphi_1 - \varphi_2) \\ \left. - 4a_2\mu \sin(2t + 2\varphi_2) - 4a_2\cos(2t + 2\varphi_2) \right)$$
$$ + \frac{\epsilon}{8a_2} \left[ \begin{aligned} & 4a_2 + 3\alpha a_2^3 - 4a_1 \cos (2t + \varphi_1 + \varphi_2) - 4a_1 \cos (\varphi_1 - \varphi_2) \\ & + 4a_1 \rho \sin (\varphi_1 - \varphi_2) + 4a_1 \rho \sin (2t + \varphi_1 + \varphi_2) \\ & + 4a_2 \cos (2t + 2\varphi_2) - 4a_2 \mu \sin (2t + 2\varphi_2) \\ & + 4\alpha a_2^3 \cos (2t + 2\varphi_2) + \alpha a_2^3 \cos (4t + 4\varphi_2) \\ & - a_1^3 \rho r \sin (\varphi_1 - \varphi_2) - a_1^3 \rho r \sin (2t + 3\phi_1 - \phi_2) \\ & a_1^3 \rho r \sin (2t + \phi_1 + \phi_2) a_1^3 \rho r \sin (4t + 3\phi_1 + \phi_2) \end{aligned} \right]. \quad (E.41) $$

## E.9 APPENDIX B. GENERAL FORMULATION OF THE STANDARD STOCHASTIC AVERAGING METHOD

In this appendix the stochastic averaging method [114, 191] is briefly described. For that we consider the following system of differential equations in standard form

$$ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}, t) + \mathbf{g}(\mathbf{x}, t)\boldsymbol{\eta} \quad (E.42) $$

where $\mathbf{x} \in \mathbb{R}^n$. If the deterministic vector function $\mathbf{f}(\mathbf{x}, t) \in \mathbb{R}^n$ and matrix function $\mathbf{g}(\mathbf{x}, t) \in \mathbb{R}^n \times \mathbb{R}^n$ satisfy certain requirements [114] and if the elements of the vector $\boldsymbol{\eta}$ are broadband processes, with zero means, then the slow (or averaged) dynamics of Eq. (G.70) may be approximated by the following stochastic differential equations

$$ \dot{\mathbf{x}} = \mathbf{m}(\mathbf{x}, t) + \boldsymbol{\sigma}(\mathbf{x}, t)\boldsymbol{\xi}, \quad (E.43) $$

where $\boldsymbol{\xi} \in \mathbb{R}^n$ is a vector of $n$ white noise processes with zero means. The vector $\mathbf{m}$ and the matrix $\boldsymbol{\sigma}$ are called drift vector and diffusion matrix respectively and are defined by

$$ \mathbf{m} = M^{\text{av}} \left\{ \mathbf{f} + \int_{-\infty}^{0} \mathbb{E} \left[ \left( \frac{\partial(\mathbf{g}\boldsymbol{\eta})}{\partial \mathbf{x}} \right)_t (\mathbf{g}\boldsymbol{\eta})_{t+s} \right] ds \right\}, \quad (E.44) $$

where $\frac{\partial(\mathbf{g}\boldsymbol{\eta})}{\partial \mathbf{x}}$ is the Jacobian matrix of $\mathbf{g}\boldsymbol{\eta}$, and

$$ \boldsymbol{\sigma}\boldsymbol{\sigma}^T = M^{\text{av}} \left\{ \int_{-\infty}^{+\infty} \mathbb{E} \left[ (\mathbf{g}\boldsymbol{\eta})_t (\mathbf{g}\boldsymbol{\eta})_{t+s}^T \right] ds \right\}, \quad (E.45) $$

where $\{\cdot\}^T$ and $\mathbb{E}[\{\cdot\}]$ denote respectively the transpose and the expected value of $\{\cdot\}$. $M^{\text{av}}$ is an averaging operator defined as follows

$$ M^{\text{av}} \{\cdot\} = \lim_{T \to +\infty} \frac{1}{T} \int_{t_0}^{t_0+T} \{\cdot\} dt. \quad (E.46) $$

It should be noted that in the case of periodic variables with period $T_0$ (which is the case in this paper), the operator $M^{\text{av}}$ becomes a classical Krylov-Bogolyubov time averaging over one period $T_0$, i.e.

$$ M^{\text{av}} \{\cdot\} = \frac{1}{T_0} \int_{t_0}^{t_0+T_0} \{\cdot\} dt \quad (E.47) $$

and the result is independent of $t_0$.

## E.10 APPENDIX C. ELEMENTS OF THE CLASSICAL DETERMINISTIC ANALYSIS OF THE SLOW FLOW DYNAMICS

Combining Eqs. (E.30b) and (E.30c) leads to the following amplitude and angle equations

$$ a_1 = a_2 \sqrt{\mu^2 + \left(1 - \frac{3\alpha a_2^2}{4}\right)^2} = H(a_2) \quad (E.48a) $$

$$ \tan \varphi = \frac{4\mu}{4 - 3\alpha a_2^2}. \quad (E.48b) $$
<figure><img src="image_114.png" /><figcaption>FIGURE E.9. Typical example of the critical manifold in the $(a_2, a_1)$-plane given by Eq. (E.48a) for $\mu = 0.25$ and $\alpha = 3$.</figcaption></figure>

As usual, the stability of the critical manifold is carried out studying the fast subsystem (E.31). A fixed point of Eq. (E.31) is shown to be stable if $d_{a_2}H(a_2) > 0$ and unstable if $d_{a_2}H(a_2) < 0$. Exploiting the polynomial properties of $H$, it can be shown that the local extrema of $H$ (i.e $d_{a_2}H(a_2) = 0$) occur at

$$ 
\begin{aligned}
a_2^{\text{LF}} &= \frac{2}{3\sqrt{\alpha}} \sqrt{2 - \sqrt{1 - 3\mu^2}} \\
a_2^{\text{RF}} &= \frac{2}{3\sqrt{\alpha}} \sqrt{2 + \sqrt{1 - 3\mu^2}}
\end{aligned} 
\qquad (\text{E.49}) $$

if the following relation holds

$$ \mu < \frac{1}{\sqrt{3}}. \qquad (\text{E.50}) $$

If the condition (E.50) is not satisfied, the $H$ function no longer has local extrema. In the rest of the paper, one considers that (E.50) always holds.

In Eq. (E.49), $a_2^{\text{LF}}$ and $a_2^{\text{RF}}$ are the abscissa values in the $(a_2, a_1)$-plane of the maximum and the minimum of the $H$ function respectively. The superscripts (')$^{\text{LF}}$ and (')$^{\text{RF}}$ refer to left fold point and right fold point respectively. Indeed, in the $(a_2, a_1, \varphi)$-space, the two points $(a_2^{\text{RF}}, a_1^{\text{RF}}, \varphi^{\text{RF}})$ and $(a_2^{\text{LF}}, a_1^{\text{LF}}, \varphi^{\text{LF}})$ (where $a_1^{\text{LF}}, a_1^{\text{RF}}, \varphi^{\text{LF}}$ and $\varphi^{\text{RF}}$ are obtained from $a_2^{\text{LF}}$ and $a_2^{\text{RF}}$ using Eq. (E.48)) are generally called fold points.

Each point of the critical manifold is a fixed point of the fast subsystem (E.31). It can be shown that the stability of these fixed points changes at the fold points: they are stable (resp. unstable) for $a_2 < a_2^{\text{LF}}$ and $a_2 > a_2^{\text{RF}}$ (resp. $a_2^{\text{LF}} < a_2 < a_2^{\text{RF}}$).

Given the previous results, the critical manifold $M_0$ is hyperbolic$^3$ except at the fold points and it consists of an attracting branch

$$ M_{0,a} = \left\{ (a_1, a_2, \varphi) \in \mathbb{R}^{+2} \times [-\pi, \pi] \mid a_2 < a_2^{\text{LF}} \text{ and } a_2 > a_2^{\text{RF}} \right\} \qquad (\text{E.51}) $$

and a repelling branch

$$ M_{0,r} = \left\{ (a_1, a_2, \varphi) \in \mathbb{R}^{+2} \times [-\pi, \pi] \mid a_2^{\text{LF}} < a_2 < a_2^{\text{RF}} \right\}. \qquad (\text{E.52}) $$

A typical amplitude part of the critical manifold (see Eq. (E.48a)) is shown in Fig. E.9. One can see that the fold points connect attracting parts to the repelling part of the critical manifold $M_0$.

Solving $H(a_2^{\text{RF}}) = H(a_2^{\text{D}})$ and $H(a_2^{\text{LF}}) = H(a_2^{\text{U}})$ the two following scalars $a_2^{\text{D}}$ and $a_2^{\text{U}}$ are obtained

$$ 
\begin{aligned}
a_2^{\text{D}} &= \frac{2\sqrt{2}}{3\sqrt{\alpha}} \sqrt{1 - \sqrt{1 - 3\mu^2}} \\
a_2^{\text{U}} &= \frac{2\sqrt{2}}{3\sqrt{\alpha}} \sqrt{1 + \sqrt{1 - 3\mu^2}}
\end{aligned} 
\qquad (\text{E.53}) $$

which are the horizontal projection of the fold points on the critical manifold.

<footnote n="3">3. $M_0$ is hyperbolic when all eigenvalues of the Jacobian of the fast subsystem on $M_0$ have nonzero real parts, here that means $d_{a_2}H(a_2) \neq 0$.</footnote>
# ROBUST OPTIMIZATION OF NONLINEAR ENERGY SINKS USED FOR MITIGATION OF FRICTION-INDUCED LIMIT CYCLE OSCILLATIONS

Cherif SNOUN, BAPTISTE BERGEOT AND CHRISTOPHE VERGEZ

## ABSTRACT

This paper aims at proposing robust methods to optimize nonlinear energy sinks (NES) used for the mitigation of friction-induced vibrations due to mode coupling instabilities. The study is based on a mechanical system composed of two NES coupled to the well-known two-degrees-of-freedom Hultèn's model. In such an unstable system coupled with NES, it is usual to observe a discontinuity in the steady-state amplitude profile which separates the parameter space into two parts corresponding respectively to the mitigated and unmitigated regimes. The discontinuity is predicted by a methodology previously developed by the authors and based on Multi-Element generalized Polynomial Chaos. The method allows to determine the Propensity of the system to undergo a Harmless Steady-State Regime (PHSSR). The objective of the present work is therefore to maximize the value of the PHSSR to obtain a robust optimal design of the NES. To this end, several stochastic optimization problems are presented that take into account the dispersion of the uncertain parameters using two approaches; in the first one, the parameters of the NES are considered as deterministic, and in the second one they are also supposed uncertain but with a known probability law.

**Keywords:** Friction-induced vibration, Nonlinear Energy Sink, Uncertainty, Robust optimisation, Polynomial Chaos

## F.1 INTRODUCTION

Friction-induced vibrations rank among the major issues that friction systems may face. This phenomenon can result from the generation of Limit Cycle Oscillations (LCO) induced by dry friction, which are explained by the coupling of two modes of the system [72, 97, 152]. Given the parametric uncertainties present in friction systems, it is very difficult to design them so that they always have stable behaviors. It is therefore necessary to use strategies to attenuate vibrations when they appear. As a way to model this phenomenon, the well-known two-degrees-of-freedom Hultèn's model [101, 102] has been widely used in the past.

Nonlinear Energy Sinks (NES) are nonlinear passive vibration absorbers which consist in a spring mass damper with a strong nonlinear stiffness (usually cubic as is the case in this work). The NES can adapt themselves to the Primary Structure (PS) without being tuned to a specific frequency. Their operation is based on the concept of Targeted Energy Transfer (TET) which has become an important passive control technique for reducing or eliminating unwanted vibrations [203, 205]. The design and optimization of NES constitute a very dynamic field of research. Recently, a lot of work has been produced concerning the study of new configurations and technologies allowing/enhancing TET [197, 219], the analysis of the underlying dynamics of the coupled system [a2, a9, 96], the deterministic optimization of NES [80, 115], the stochastic optimization of NES [161, 162], the coupling of the NES to an energy harvester [60, 108, 109] or the use of a NES as an energy harvester [125, 226], etc. Among the many research subfields about NES, the present work takes place within the framework of robust optimization of NES when the latter have uncertain parameters and are used to mitigate Limit Cycle Oscillations (LCO).

The mitigation of LCO by NES has also been widely studied in the past. The first work reported in this framework concerns the LCO mitigation of the Van der Pol oscillator [122]. In the field of friction-induced vibrations mitigation, [a14] analyzed the behavior of a Hultèn's model coupled to two ungrounded NES; the same model in considered in this work.

Lately, special attention has been paid to the optimal design of NES in a deterministic context. Nguyen and Pernot [150] have studied the effectiveness of a NES to attenuate free oscillations of a single degree-of-freedom oscillator. The optimization procedure consisted in studying the influence of the NES damping on the efficiency
of the energy transfer from the primary system to the NES. Wang et al. [210] presented multi-objective designs of a track NES and a single-sided vibro-impact track NES used as effective control devices to mitigate the seismic response of high-rise buildings. In [165] a design criteria for optimally tuned Vibro-Impact NES to control vibration under periodic and transient excitation has been proposed. Oliva et al. [154] developed an approximate design approach based on the use of the Statistical Linearization Technique, and an accurate empirical formulation linking the NES optimal parameters to the characteristics of the main structure and the random excitation.

More recently, the optimization of NES under uncertainties has also been considered. Boroson [31] has introduced the “discontinuity in the effectiveness” of the NES, i.e., the fact that the effectiveness can change from a high value to a low value following a small perturbation of the design parameters or of the initial conditions. In [30], the optimal design of several parallel NES configurations for maximum mean efficiency has been investigated. It has been shown again that the variation of some system parameters or of the initial conditions generates a discontinuity on the energy shape of the system. An optimization algorithm based on the Support Vector Machine (SVM) technique and a kriging meta-model have been used to predict this discontinuity according to the dispersion of some uncertain parameters. Khazaee et al. [115] studied the stochastic optimization of multiple NES (configured in parallel and in series) attached to a simply supported pipe conveying fluid. The results show that the optima obtained from stochastic optimization are much more robust than those obtained from deterministic optimization. Pidaparthi and Missoum [161, 162] performed optimization of NES under uncertainties in the context of attenuating self-sustained vibrations of an aircraft wing created by aeroelastic instability (which is also due to mode coupling). Depending on the nonlinear properties of the system, the latter may undergo LCO that are either supercritical or subcritical. In this case the bifurcation diagram of the system (which represents the amplitude of the LCO as a function of the airflow velocity over the wing) may exhibit two types of discontinuities, respectively due to the operation of the NES, and to the generation of subcritical LCO. Using again the SVM technique and a kriging meta-model, several stochastic optimization problems have been presented in order to maximize the average reduction of the LCO amplitude.

Various probabilistic approaches may be used to study the propagation of uncertainties in the deterministic model (DM) of a mechanical system. A review of the numerical methods dedicated to the stochastic prediction can be found in [151]. Efficient methods such as the generalized Polynomial Chaos (gPC) or the Multi-Element generalized Polynomial Chaos (ME-gPC) have been developed as a less costly alternative to the common Monte Carlo approach. These methods have been for example used to perform a stochastic analysis of the dynamic behavior of friction systems [146-149, 173, 175, 177].

In this paper, two original optimization methods under uncertainties, based on polynomial chaos theory, are developed for the robust design of two ungrounded NES used to mitigate LCO due to mode coupling instability in the Hultèn's model, which is the primary structure (PS). [a14] classified the steady-state response regimes in two main categories related to the dispersion of some uncertain parameters : the mitigated regimes (harmless situations) and the unmitigated regimes (harmful situations). In the context of instability mitigation by means of NES, the LCO amplitude profile as a function of a given bifurcation parameter usually presents a discontinuity between these two regimes, which makes the NES potentially very sensitive to uncertainties. In [a17], this discontinuity is predicted by means of a ME-gPC based method. This allows to determine the Propensity of the system to undergo a Harmless Steady-State Regime (PHSSR), which is defined as the probability for the system to undergo a mitigated regime. The aim of this paper is therefore to maximize the PHSSR value to obtain a robust optimal design of the NES with a low computational cost. For this purpose, several stochastic optimization problems are presented for a PS having uncertain parameters, in which the NES parameters are either considered as deterministic, or are also supposed uncertain but with a known probability law.

The article is organized as follows. In Section F.2, the two degrees-of-freedom Hultèn's model coupled to two NES is presented. In Section F.3, the possible steady-state regimes of the coupled system are studied, and the PHSSR is defined. Section F.4 describes the polynomial chaos theory and recalls the stochastic optimization algorithm used to compute the PHSSR. The optimization formulation under uncertainties is described in Section F.5, and the results are presented in Section F.6. Finally, conclusions are given in Section F.8.

## F.2 THE FRICTION SYSTEM UNDER STUDY

The studied mechanical system is similar to that previously used by the authors in [a14, a17]. It is recalled hereafter for the sake of clarity.
<figure><img src="image_115.png" /><figcaption>FIGURE F.1. The mechanical model. (a) Mechanical system with NES; (b) Zoom on NES1; (c) Zoom on NES2.</figcaption></figure>

## F.2.1 The primary system

The primary system (PS) consists in the two degree-of-freedom (DOF) Hultèn's phenomenological model ([101, 102]). This simple model can reproduce the mode-coupling phenomenon and is composed of a mass $m$ held against a moving strip. The contact between the mass and the strip is modeled by two plates supported by two different springs, with linear components $k_1$ and $k_2$ and cubic components $k_1^{NL}$ and $k_2^{NL}$, and two dampers with damping coefficients $c_1$ and $c_2$ (see Fig. F.1). We assume that the friction coefficient is constant and that the strip moves at a constant velocity. The relative velocity between the strip velocity and $dx_1/dt$ or $dx_2/dt$ is supposed to be positive, and the friction contact is modeled by the Coulomb's law, i.e. $F_T = \mu F_N$ where $F_T$ is the tangential friction force, $F_N$ is the normal force and $\mu$ is the friction coefficient. Using Newton's second law the equations of motion, that constitute the primary system in this work, are derived as follows:

$$m\ddot{x}_1 + c_1\dot{x}_1 + k_1x_1 - \mu k_2x_2 + k_1^{NL}x_1^3 - \mu k_2^{NL}x_2^3 = 0 \quad (F.1a)$$

$$m\ddot{x}_2 + c_2\dot{x}_2 + k_2x_2 + \mu k_1x_1 + \mu k_1^{NL}x_1^3 + k_2^{NL}x_2^3 = 0. \quad (F.1b)$$

## F.2.2 Mechanical model with Nonlinear Energy Sinks

Two strongly cubic and identical ungrounded NES with masses $m_h$, damping coefficients $c_h$ and cubic stiffnesses $k_h^{NL}$, are attached on the primary system (see Fig. F.1). Introducing the NES displacements $h_1(t)$ and $h_2(t)$ in Eq. (F.1) leads to the following equations of motion of the coupled system

$$\ddot{x}_1 + \eta_1\omega_1\dot{x}_1 + \omega_1^2x_1 - \mu\omega_2^2x_2 + \varphi_1x_1^3 - \mu\varphi_2x_2^3 + \eta_h\omega_1(\dot{x}_1 - \dot{h}_1) + \xi_h(x_1 - h_1) + \varphi_h(x_1 - h_1)^3 = 0 \quad (F.2a)$$

$$\ddot{h}_1 + \eta_h\omega_1(\dot{h}_1 - \dot{x}_1) + \xi_h(h_1 - x_1) + \varphi_h(h_1 - x_1)^3 = 0 \quad (F.2b)$$

$$\ddot{x}_2 + \eta_2\omega_2\dot{x}_2 + \omega_2^2x_2 + \mu\omega_1^2x_1 + \mu\varphi_1x_1^3 + \varphi_2x_2^3 + \eta_h\omega_1(\dot{x}_2 - \dot{h}_2) + \xi_h(x_2 - h_2) + \varphi_h(x_2 - h_2)^3 = 0 \quad (F.2c)$$

$$\ddot{h}_2 + \eta_h\omega_1(\dot{h}_2 - \dot{x}_2) + \xi_h(h_2 - x_2) + \varphi_h(h_2 - x_2)^3 = 0, \quad (F.2d)$$

where $\eta_i = c_i/\sqrt{mk_i}$ (with $i = 1, 2$) are the relative damping coefficients, $\omega_i = \sqrt{k_i/m}$ are the natural angular frequencies, $\varphi_i = k_i^{NL}/m$; as for the NES parameters, $\epsilon = m_h/m$ is the mass ratio between the NES and the primary structure, $\xi_h = k_h/m$, $\eta_h = c_h/\sqrt{mk_1}$ and $\varphi_h = k_h^{NL}/m$. Because strongly cubic NES are considered, one has $\xi_h \ll \varphi_h$. In general, the mass ratio $\epsilon$ is in the range $0.01 - 0.1$ as in the present paper.
<figure><img src="image_116.png" /><figcaption>FIGURE F.2. Comparison between time series $x_1(t)$ resulting from the numerical integration of the friction system with and without NES. (a) Complete suppression, $\mu = 0.16$; (b) Mitigation : PR, $\mu = 0.18$; (c) Mitigation : SMR, $\mu = 0.2$; (d) No mitigation, $\mu = 0.22$. The set of parameters Eq. (F.4) has been used.</figcaption></figure>

## F.3 POSSIBLE STEADY-STATE REGIMES

Four main types of steady-state regimes (SSR) may be generated when one or several NES are attached on a primary system which undergoes an LCO : complete suppression of the instability, mitigation through Periodic Response (PR), mitigation through Strongly Modulated Response (SMR) or no mitigation. Those four regimes can also be observed on the studied system (F.2) ([a14]), and are presented in Fig. F.2 which shows for each case the displacements $x_1(t)$ - with and without NES attachments - with respect to time. The regimes described above are classified into two categories depending on whether the LCO is considered as mitigated (Figs. F.2(a), F.2(b) and F.2(c)) or not (Fig. F.2(d)).

In this work the Quantity of Interest (QoI) under consideration is the amplitude $A_{x_1}^{wNES}$ of the variable $x_1$ in the system with NES (F.2) and within a steady-state regime. It is defined as follows

$$A_{1}^{wNES} = \frac{\max[x_{1}^{SSR}(t)] - \min[x_{1}^{SSR}(t)]}{2}, \qquad (F.3)$$

where $x_1^{SSR}(t)$ is the time series of the variables $x_1$ obtained from the numerical integration of the coupled system (F.2) within the steady-state regime. For comparison purposes, the amplitude $A_{x_1}^{wNES}$ of the system without NES (F.1) is also computed.

The time integration of system (F.2) (resp. (F.1)) is performed between $t_b = 0$ and $t_e = 4$ seconds and the amplitude $A_{x_1}^{wNES}$ (resp. $A_{x_1}^{wNES}$) is computed in the interval $[0.9t_e, t_e]$ in which we assume that the system has reached its steady-state regime.
<figure><img src="image_117.png" /><figcaption>FIGURE F.3. Amplitudes $A_1^{wNES}$ and $A_1^{wNES}$ as functions of the friction coefficient $\mu$. The set of parameters Eq. (F.4) has been used.</figcaption></figure>

Using the following set of parameters

$$ 
\begin{aligned}
\omega_1 &= 2\pi 100 \text{ (rad}\cdot\text{s}^{-1}), & \omega_2 &= 2\pi 85 \text{ (rad}\cdot\text{s}^{-1}), \\
\eta_1 &= 0.02, & \eta_2 &= 0.06, \\
\varphi_1 &= 10^5 \text{ (N}\cdot\text{kg}^{-1}\cdot\text{m}^{-3}), & \varphi_2 &= 0 \text{ (N}\cdot\text{kg}^{-1}\cdot\text{m}^{-3}), \\
\epsilon &= 0.05, & \xi_h &= 0.001 \text{ (N}\cdot\text{kg}^{-1}\cdot\text{m}^{-1}), \\
\eta_h &= 0.02, & \varphi_h &= 1.4 \cdot 10^5 \text{ (N}\cdot\text{kg}^{-1}\cdot\text{m}^{-3}),
\end{aligned}
\qquad (F.4) $$

the amplitudes $A_1^{wNES}$ and $A_1^{wNES}$ are plotted as functions of the friction coefficient $\mu$ in Fig. F.3. The specific values of $\mu$ denoted as $\mu_b^w$ and $\mu_b^w$ correspond to the Hopf bifurcation points without and with the NES attachments, respectively. The four steady-state regimes described above are visible on the Figure, which highlights a discontinuity (or jump) in the profile of $A_1^{wNES}$. For increasing values of $\mu$, this discontinuity corresponds in general to the transition from the SMR to the no suppression regime and separates mitigated regimes from unmitigated ones. The value of $\mu$ at the discontinuity is called mitigation limit and is denoted as $\mu_{ml}$.

The Propensity to be in a Harmless Steady-State Regime (PHSSR) ([a17]) can be computed from the mitigation limit $\mu_{ml}$. The PHSSR is the probability of obtaining a mitigated regime given the probability laws governing the dispersion of the primary system parameters. Here, only $\mu$ is considered as uncertain and it is assumed to follow a uniform probability law within the interval $[0, \mu_e = 0.4]$; the PHSSR is then simply defined as

$$ \text{PHSSR} = \int_{0}^{\mu_{ml}} \frac{1}{\mu_{e}} d\mu = \frac{\mu_{ml}}{\mu_{e}}. \qquad (F.5) $$

## F.4 MULTI-ELEMENT POLYNOMIAL CHAOS BASED ALGORITHM TO DETECT A DISCONTINUITY

The ME-gPC method previously developed by the authors in [a17] to locate the discontinuity in the QoI profile is recalled in Section F.4.2. Elements of the Polynomial Chaos theory are first presented in Section F.4.1.

### F.4.1 Elements of the Polynomial Chaos theory

#### F.4.1.1 Generalized Polynomial Chaos (gPC)

Let us consider a dynamical system with $r$ uncertain physical parameters $\beta_j$ ($j = 1, \dots, r$) assumed to be uniformly distributed within a given space $\prod_{i=1}^r [a_i, b_i]$, and let $\xi_j$ ($j = 1, \dots, r$) be the independent random variables within the space $[-1]^r$, linked to $\beta_j$ by

$$ \beta_j(\xi_j) = \frac{a_j + b_j}{2} + \frac{b_j - a_j}{2} \xi_j, \quad (j = 1, \dots, r). \qquad (F.6) $$
The gPC theory ([34, 213, 220]) states that the system output variable X(ξ₁, ..., ξₜ), called the quantity of Interest (QoI), may be approximated by a truncated orthogonal polynomial function series as

$$X(\xi_1, ..., \xi_r) \approx \sum_{j=0}^{N_p} \bar{x}_j \phi_j(\xi_1, ..., \xi_r), \qquad (F.7)$$

where φⱼ(ξ) are orthogonal polynomials which represent the stochastic part of the process, whereas x̄ⱼ are the gPC coefficients that take into account the deterministic part of the process. In this paper, according to the Askey-scheme ([8]), the orthogonal polynomials φⱼ(ξ₁, ..., ξₜ) are built from Legendre polynomials because the uncertain parameters follow a uniform probabilistic law. From [220], the number of terms Nₚ + 1 of the expansion is given by

$$N_p + 1 = \frac{(p+r)!}{p!r!}, \qquad (F.8)$$

where r is the number of uncertain parameters and p is the order of the gPC. To obtain the approximated values of the QoI X, the coefficients x̄ⱼ of the truncated series (F.7) must be computed. To that end, the non-intrusive regression method ([28]) is used in this paper. In this method, the gPC coefficients are built from Q values of the QoI X. The corresponding Q simulations may be performed at points chosen with the Latin Hypercube Samples (LHS) method [138] that will be referred to as the Numerical Experimental Design (NED) in this paper. A minimum of Q = Nₚ simulations is required but in practice, Q = kNₚ simulations are used, with k a small integer usually equal to 2, 3 or 4.

### F.4.1.2 Multi-Element Generalized Polynomial Chaos (ME-gPC)

When the QoI is nonlinear, high polynomial orders may be required to reduce the approximation error. Thus, when the number of uncertain parameters r is high, the number of simulations required to compute the gPC coefficients may be excessive, leading to a prohibitive computational cost. An alternative is to split the vector ξ = (ξ₁, ..., ξₜ) into a collection of m non-intersecting elements and to use a low order polynomial approximation on each element. That is the basic idea behind the ME-gPC [209].

The local physical variables βⱼᵏ, in the kᵗʰ element ∏ᵢ=₁ᵏ [aᵢᵏ, bᵢᵏ], are expressed in terms of independent uniform random variables ξⱼᵏ in [-1, 1]ᵣ through

$$\beta_j^k = \frac{b_j^k + a_j^k}{2} + \frac{b_j^k - a_j^k}{2} \xi_j^k. \qquad (F.9)$$

Consequently, Legendre polynomials can be used locally and a gPC expansion developed in each element as follows:

$$X_k(\bar{\xi}^k) \approx \sum_{j=0}^{N_p} \bar{x}_{k,j} \phi_j(\bar{\xi}^k), \qquad (F.10)$$

where Xₖ(ξ̄ᵏ) is the random process corresponding to the kᵗʰ element.

From the gPC theory, the approximated local mean Ẑₚ,ₖ and variance σ²ₚ,ₖ in the kᵗʰ element can be expressed analytically from the gPC coefficients as

$$\left\{ \begin{array}{l} \hat{X}_{p,k} = \bar{x}_{k,0} \\ \sigma_{p,k}^2 = \frac{1}{2r} \sum_{j=1}^{N_p} \bar{x}_{k,j}^2 \langle \phi_j^2 \rangle. \end{array} \right. \qquad (F.11)$$

Several criteria have been proposed in the literature for the convergence of the ME-gPC algorithm. They depend on the aim of the ME-gPC. The criteria chosen in this study are presented in the next section.

### F.4.2 Detection of a discontinuity in the QoI with the ME-gPC

This section recalls the method based on the ME-gPC and proposed in [a17] to locate a discontinuity in a Quantity of Interest (QoI) derived from the DM of a mechanical system. The aim of the method is not necessarily to obtain an accurate representation of the QoI, but to locate precisely the discontinuity, i.e. to know the values of the uncertain parameters for which the discontinuity in the QoI profile appears.

The underlying principle of the method is that the variance of a QoI that contains a discontinuity is higher than that of a QoI that does not contain any discontinuity. Therefore, if the discontinuity lies in a given element of the stochastic parameter space, the variance of the ME-gPC expansion in that element will
be high. From these considerations a ME-gPC algorithm has been built, in which at each step and for each element the local variance $\sigma_{p,k}^2$, computed directly from the gPC coefficients using Eq. (F.11), is compared to a threshold $\theta_1$. If for a given element $k$

$$ \sigma_{p,k}^2 \geq \theta_1, \qquad (F.12) $$

it is assumed that the discontinuity lies within this element. The latter is then divided by two in each direction of the stochastic space in order to locate more accurately the discontinuity. On the contrary, if the condition Eq. (F.12) is not satisfied, the element $k$ is supposed not to contain the discontinuity and is removed from the algorithm.

The threshold $\theta_1$ is chosen so that the variance of the element which contains the discontinuity always satisfies the condition Eq. (F.12). Therefore, two other criteria are introduced to stop the algorithm. The first one is the minimum element size $J_{min}$ defined as

$$ J_{min} = \theta_2 J_0, \qquad (F.13) $$

where $J_0$ is the size of the initial element corresponding to the whole stochastic parameter space and $\theta_2$ is a percentage. At a given iteration $i$, the retained elements of size $J_i$ are actually divided into two equal parts in each direction of the stochastic space if the following condition

$$ J_i \geq J_{min} \qquad (F.14) $$

holds. It is worth recalling that all the elements at iteration $i$ have the same size $J_i$ because at each step the division is performed in each direction of the stochastic. If the condition (F.14) does not hold, the algorithm is stopped and we assume that the discontinuity in the QoI profile lies within the remaining elements.

In Section F.6, a comparison is performed with a reference method that needs $N$ simulations of the DM to locate a discontinuity in the QoI profile. Therefore, the proposed method is only effective if the number of numerical simulations $N'$ required to obtain the gPC coefficients in all the elements is smaller than $N$. Consequently, the second criterion is that $N'$ should verify

$$ N' \leq N, \qquad (F.15) $$

otherwise the algorithm is stopped.

The algorithm is summarized in Fig. F.4. At the first step, only one element is present (i.e. the whole stochastic parameter space) and for a given iteration the parameter space has been divided into $K$ elements during the previous iterations. At each iteration, the gPC coefficients are computed for the $K$ elements and we check if the corresponding computational cost is larger than the reference one (see Eq. (F.15)). If that is the case, the algorithm is stopped and in general the discontinuity in the QoI profile has not been located. If Eq. (F.15) holds, we check the variance in each of the $K$ elements (see Eq. (F.13)). If it is smaller than the threshold $\theta_1$, the element is removed from the algorithm. If it is larger, it is kept and $K'$ elements remain at the exit of the loop in $k$. If Eq. (F.14) does not hold, the algorithm is stopped and we assume that the discontinuity in the QoI profile lies within the $K'$ remaining elements. If Eq. (F.14) holds, then each of the $K'$ elements is divided into two parts in each direction of the stochastic parameter space, then a new iteration begins with $K = 2^r K'$.

As for the reference method, the detection of the discontinuity can be performed by the so-called Jump Criterion (JC). To be considered as an effective jump, the difference between the amplitude at the right of the jump, denoted as $A(\mu_{ml}^+)$, and the amplitude at the left, denoted as $A(\mu_{ml}^-)$, must be greater than the half of $A(\mu_{ml}^-)$, i.e.

$$ JC = \frac{A(\mu_{ml}^{+}) - A(\mu_{ml}^{-})}{A(\mu_{ml}^{-})} \geq \frac{1}{2}. \qquad (F.16) $$

As an example, Fig. F.5 shows the last iteration of the algorithm applied to the system (F.2) for $\theta_1 = 0.7 \cdot 10^{-3}$, $\theta_2 = 1\%$ and $p = 1$. The Reference curve is similar to that displayed in Fig. F.3. One can see that the last interval obtained with the ME-gPC method contains the discontinuity observed on the reference.

## F.5 OPTIMIZATION STRATEGY UNDER UNCERTAINTIES

This section focuses on the optimization methodology proposed to maximize the PHSSR defined in Section F.3. Two NES parameters are retained as the design parameters to optimize : the damping coefficient $\eta_h$ and the nonlinear stiffness $\phi_h$. The uncertain parameter is again the friction coefficient $\mu$. Following [30, 162], the problem is treated in two different ways depending on whether the physical design parameters $\eta_h$ and $\phi_h$ are deterministic or stochastic.
<figure><img src="image_118.png" /><figcaption>FIGURE F.4. Algorithm of the proposed method based on the ME-gPC.</figcaption></figure>

F.5.1 Deterministic design variables

In this case, the parameters to optimize are directly the physical design parameters which are assumed to be known with certainty. The objective of the method is to find the optimal values of the NES parameters $x_d$ that maximize the PHSSR value. For this purpose, the PHSSR is determined in the space of uncertain parameters $x_u$ for each value of $x_d \in [x_d^{(\min)} \ x_d^{(\max)}]$. The maximum of the PHSSR is then detected, which
<figure><img src="image_119.png" /><figcaption>FIGURE F.5. Last iteration of the detection algorithm when the friction coefficient μ is uncertain, with θ₁ = 0.7 · 10⁻³, θ₂ = 1% and p = 1; reference curve obtained from Fig. F.3.</figcaption></figure>

corresponds to the optimal value of $x_d$. The method can be summarized as follows :

$$
\begin{array}{ll}
\text{Maximize} & \text{PHSSR } (A(x_d, x_u)) \\
\text{Subject to} & x_d^{(\min)} \le x_d \le x_d^{(\max)} \\
\text{where} & x_d = \{\eta_h, \varphi_h\}, \\
& x_u = \{\mu\}
\end{array}
\qquad (\text{F.17})
$$

where

* A(x_d, x_u) is the amplitude of the displacement of the coupled system (F.2) defined by Eq. (F.3) (from now on the amplitude is simply written A instead of A_1^WNES to lighten the notations); it depends on both the uncertain parameters x_u and the design variables x_d;

* $x_d^{(\min)}$ and $x_d^{(\max)}$ are the minimum and maximum values of the design variables $x_d$ respectively; the optimal values of the design parameters are searched within the intervals $[x_d^{(\min)}, x_d^{(\max)}]$, called *design intervals*.

F.5.2 Stochastic design variables

The physical design parameters are considered here to be random with a known probability law. This
case is closer to a real situation, as it is known in particular that obtaining an exact value of a damping
coefficient is almost impossible. The parameters to be optimized are then statistical characteristics of these
physical parameters (for instance their means or variances) and are called *hyperparameters*. In this study the
hyperparameters considered are the means $\hat{\eta}_h$ and $\hat{\varphi}_h$ of the damping coefficient and of the cubic stiffness of
the NES, respectively.

The optimization process consists in dividing the design space into subspaces. In each subspace a draw of N values of the design parameters is performed according to their probability laws (in this work we use uniform laws). For each of these N values the PHSSR is calculated and the optimal value of x_d (here of $\hat{\eta}_h$ or $\hat{\varphi}_h$) corresponds to the subspace where the mean of the N PHSSR values is maximum. This can be summarized as follows :

$$
\begin{array}{ll}
\text{Maximize} & \mathbb{E}[\text{PHSSR}(A(x_d, x_u))] \\
\text{Subject to} & x_d^{(\min)} \le x_d \le x_d^{(\max)} \\
\text{where} & x_d = \{\hat{\eta}_h, \hat{\varphi}_h\}, \\
& x_u = \{\mu\}
\end{array}
\qquad (\text{F.18})
$$

where

* $\mathbb{E}$ is the value of the mathematical expectation (the mean) of all PHSSR evaluations;
* $\hat{\eta}_h$ is the mean of the NES damping coefficient values;

* $\hat{\phi}_h$ is the mean of the NES cubic stiffness values.

## F.6 RESULTS FOR DETERMINISTIC DESIGN VARIABLES

As explained before, the first method of NES optimization considers that the NES parameters are deterministic and constitute therefore directly the design variables to optimize.

### F.6.1 Optimization of the NES damping coefficient

In this section, we consider first that the design space is a one-dimensional space according to the variable $\eta_h$, represented by the interval $[0, 0.04]$. The optimization problem (F.17) is therefore solved with $x_d = \eta_h \in [0, 0.04]$ with the other parameters given Eq. (F.4). The results are expressed in terms of computational cost, on the one hand, and of precision compared to the reference method, on the other hand. For this purpose, the following quantities are introduced:

* $\mu_{ml}^{\max}$ is the maximum value of the mitigation limit $\mu_{ml}$, and corresponds to the maximum value of the PHSSR, denoted as $\text{PHSSR}^{\max}$. The mitigation limit is here searched within the interval $[0, 0.4]$. The PHSSR value estimated by the ME-gPC method is determined using the upper bounds of the last intervals found by the algorithm presented in Fig. F.5. The precision of the corresponding $\mu_{ml}^{\max}$ value is therefore given by:

$$ \mu_{ml}^{\max} \pm \frac{\Delta\mu}{\mu_{ml}^{\max}}. \qquad (F.19) $$

For the reference method, $\Delta\mu$ is given by the discretization step of the design space, whereas $\Delta\mu = 0.4\theta_2$, with $\theta_2$ defined by Eq. (F.13), for the ME-gPC method.

* $E_{\text{PHSSR}^{\max}}$ is the relative error of the maximum PHSSR obtained with the ME-gPC method ($\text{PHSSR}_{\text{ME-gPC}}^{\max}$) compared to the maximum PHSSR value obtained with the reference method ($\text{PHSSR}_{\text{ref}}^{\max}$), and is defined as follows:

$$ E_{\text{PHSSR}^{\max}} = \frac{|\text{PHSSR}_{\text{ME-gPC}}^{\max} - \text{PHSSR}_{\text{ref}}^{\max}|}{\text{PHSSR}_{\text{ref}}^{\max}}. \qquad (F.20) $$

* $\eta_h^{\text{opt}}$ is the optimal value of the NES damping coefficient, i.e. the value for which $\mu_{ml}^{\max}$, and thus $\text{PHSSR}^{\max}$, are obtained. Since the design space is discretized with a step $\Delta\eta_h$, the optimal value is expressed as follows:

$$ \eta_h^{\text{opt}} \pm \frac{\Delta\eta_h}{\eta_h^{\text{opt}}}. \qquad (F.21) $$

* $Err$ is the relative error of the optimal design value $\eta_h^{\text{opt}}$ obtained with the ME-gPC method, compared to that obtained with the reference method:

$$ Err = \frac{|(\eta_h^{\text{opt}})_{\text{ME-gPC}} - (\eta_h^{\text{opt}})_{\text{ref}}|}{(\eta_h^{\text{opt}})_{\text{ref}}}. \qquad (F.22) $$

### F.6.1.1 Reference optimization

The reference value of the optimal NES damping coefficient, denoted as $(\eta_h^{\text{opt}})_{\text{ref}}$, is obtained by solving (F.17) where the amplitude $A(x_d, x_u)$, and thus the mitigation limit, directly result from the numerical simulation of (F.2). The design space for $\eta_h$ $[0, 0.04]$ is uniformly discretized into 51 values. To reduce the computational cost, the mitigation limit $\mu_{ml}$ and the PHSSR are first roughly estimated for each $\eta_h$ value using only 1000 values of the friction coefficient $\mu \in [0, 0.4]$ and applying the criterion (F.16). Fig. F.6 shows the resulting evolution of the PHSSR as a function of $\eta_h$. This first step enables the locating of the interval - here $[0.0208, 0.028]$ - in which the optimal value of the design parameter $(\eta_h^{\text{opt}})_{\text{ref}}$ lies. In this interval the PHSSR is maximum ($\text{PHSSR}^{\max} = 50.10\% \pm 0.19\%$) but constant due to the too coarse discretization in $\mu$. To refine this reference result, the PHSSR is then calculated in the design interval $\eta_h \in [0.02, 0.03]$ with the same discretization step for $\eta_h$ as previously, but now with 10 000 values of $\mu$ to determine the mitigation limit for each value of $\eta_h$. The total number of simulations needed to obtain the reference value is therefore equal to $50 \times 1000 + 10 \times 10000 = 150000$.

Fig. F.7 shows the new PHSSR as a function of $\eta_h$. In this case, $(\mu_{ml}^{\max})_{\text{ref}} = 0.2004 \pm 0.02\%$ which gives ($\text{PHSSR}^{\max}$)$_{\text{ref}} = 50.1\% \pm 0.02\%$. The corresponding optimal design value of the damping coefficient is $(\eta_h^{\text{opt}})_{\text{ref}} = 0.0240 \pm 3.33\%$.
<figure><img src="image_120.png" /><figcaption>**FIGURE F.6.** Result of the first step of the reference optimization procedure: PHSSR as a function of $\eta_h$, for $\phi_h = 140\,000 \text{ N}\cdot\text{kg}^{-1}\cdot\text{m}^{-3}$.</figcaption></figure>

<figure><img src="image_121.png" /><figcaption>**FIGURE F.7.** Result of the second step of Reference optimization: PHSSR as a function of $\eta_h$, for $\phi_h = 140\,000 \text{ N}\cdot\text{kg}^{-1}\cdot\text{m}^{-3}$.</figcaption></figure>

### F.6.1.2 Optimization with the ME-gPC method

In the proposed approach, the optimal value of the NES damping coefficient, denoted as $(\eta_h^{\text{opt}})_{\text{ME-gPC}}$, is still obtained by solving (F.17), but the amplitude $A(x_d, x_u)$ and the mitigation limit (F.17) are now determined by the method based on the ME-gPC.

Three optimization procedures are applied with a gPC order $p=1$, a variance threshold $\theta_1 = 2 \cdot 10^{-3}$ and three values of the element size threshold: $\theta_2 = 1\%$, 0.1% and 0.001%. The results of this ME-gPC based optimization are presented in Tab. F.1, along with the reference results for comparison.

Fig. F.8 displays the variations of the PHSSR as functions of $\eta_h$ respectively obtained with the reference method and the three ME-gPC methods with $\theta_2 = 1\%$, 0.1% and 0.001%. For $\theta_2 = 1\%$ the number of simulations is equal to 4324, the optimal value is $(\eta_h^{\text{opt}})_{\text{ME-gPC}} = 0.0241 \pm 15.14\%$ and the maximum value of the corresponding PHSSR is 50.78%, implying an error compared to the reference solution $E_{\text{PHSSRmax}} = 1.35\%$. In the case of $\theta_2 = 0.1\%$, although the new value of $\theta_2$ has no influence on $\eta_h^{\text{opt}}$, it changes the maximum value of the PHSSR to $\text{PHSSR}_{\text{ME-gPC}}^{\text{max}} = 50.09\%$ which induces a relative error $E_{\text{PHSSRmax}} = 0.02\%$. This error is lower than the previous one (with $\theta_2 = 1\%$) but the number of simulations has increased and is now equal to 5165. For $\theta_2 = 0.001\%$ the optimal value becomes $(\eta_h^{\text{opt}})_{\text{ME-gPC}} = 0.0237 \pm 3.33\%$, implying a relative error $Err = 1.25\%$ compared to the reference value, which is higher than the errors in the previous cases. This is not illogical since the error is calculated according to the midpoints of the intervals. Nevertheless, the interval obtained with $\theta_2 = 0.001\%$ is much narrower ($\pm 3.33\%$) than that obtained with $\theta_2 = 0.01\%$ ($\pm 15.14\%$), and contains the value resulting from the reference method. The total number of simulations is in this last case equal to 7026, which remains 22 times smaller than the number needed with the reference method.

The plateaus observed in Fig. F.8 are due to the fact that for two close values of $\eta_h$ the mitigation limits
<table><caption>**TABLEAU F.1.** Comparison between the reference optimization and the ME-gPC optimization considering deterministic design variables. The design parameter is $\eta_h$. The reference results are obtained as described in Sect. F.6.1.1.</caption><thead><tr><th></th><th>Nb of simulations</th><th>Computational time</th><th>PHSSR<sup>max</sup></th><th>E<sub>PHSSR<sup>max</sup></sub></th><th>η<sub>h</sub><sup>opt</sup></th><th>Err</th></tr></thead><tbody><tr><td>Reference</td><td>150 000</td><td>2.15 days</td><td>50.10% ± 0.02%</td><td>—</td><td>0.0240 ± 3.33%</td><td>—</td></tr><tr><td>ME-gPC<br>(θ<sub>1</sub> = 2 · 10<sup>-3</sup><br>and θ<sub>2</sub> = 1%)</td><td>4324</td><td>1.70 hours</td><td>50.78% ± 1.97%</td><td>1.35%</td><td>0.0241 ± 15.14%</td><td>0.41%</td></tr><tr><td>ME-gPC<br>(θ<sub>1</sub> = 2 · 10<sup>-3</sup><br>and<br>θ<sub>2</sub> = 0.1%)</td><td>5165</td><td>2.06 hours</td><td>50.09% ± 0.2%</td><td>0.02%</td><td>0.0241 ± 15.14%</td><td>0.41%</td></tr><tr><td>ME-gPC<br>(θ<sub>1</sub> = 2 · 10<sup>-3</sup><br>and<br>θ<sub>2</sub> = 0.001%)</td><td>7026</td><td>3.23 hours</td><td>50.09% ± 0.02%</td><td>0.02%</td><td>0.0237 ± 3.33%</td><td>1.25%</td></tr></tbody></table>

<figure><img src="image_122.png" /><figcaption>**FIGURE F.8.** Comparison between the optimization strategies based on the reference method and the ME-gPC for $\eta_h \in [0 \ 0.04]$ with $\phi_h = 140 000 \ N \cdot kg^{-1} \cdot m^{-3}$. (a) PHSSR and (b) Zoom of the PHSSR. The reference results are obtained as described in Sect. F.6.1.1.</figcaption></figure>

are found in the same element. They are related to the value of the threshold $\theta_2$: the lower it is, the more accurate the localization of the mitigation limit is, with in counterpart an increase of the computational cost. The precision is therefore driven by the value of the threshold $\theta_2$.

Finally, with a computational cost much lower than the reference one, the proposed approach proves efficient to identify the value of the optimal NES damping coefficient while taking into account the dispersion of the friction coefficient.

## F.6.2 Optimization of the NES nonlinear stiffness

The aim of this section is to find the optimal value of the cubic stiffness of the NES $\phi_h$ assumed here to be deterministic. The damping coefficient value $\eta_h$ is here set to 0.024. Similarly to the previous section, the following quantities are defined:

* $\phi_h^{\text{opt}}$ is the optimal value of the NES cubic stiffness for which $\mu_{\text{ml}}^{\text{max}}$, and thus $\text{PHSSR}^{\text{max}}$, are obtained. It may be again expressed with respect to the discretization step $\Delta\phi_h$ of the design space as

$$ \phi_h^{\text{opt}} \pm \frac{\Delta\phi_h}{\phi_h^{\text{opt}}}. \qquad (\text{F.23}) $$