# Low-order coupled model for vortex-induced vibrations mitigation by resonant piezoelectric shunt

Arthur Haudeville, Xavier Amandolese, Boris Lossouarn, Christophe

Giraud-Audine, Olivier Thomas

► To cite this version:

Arthur Haudeville, Xavier Amandolese, Boris Lossouarn, Christophe Giraud-Audine, Olivier Thomas. Low-order coupled model for vortex-induced vibrations mitigation by resonant piezoelectric shunt. Journal of Fluids and Structures, 2025, 139, pp.104442. 10.1016/j.jfluidstructs.2025.104442. hal-05331960

HAL Id: hal-05331960

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

Submitted on 26 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Low-order coupled model for vortex-induced vibrations mitigation by
resonant piezoelectric shunt

Arthur Haudeville<sup>a,b,c,*</sup>, Xavier Amandolese<sup>a</sup>, Boris Lossouarn<sup>a</sup>, Christophe Giraud-Audine<sup>b</sup>, Olivier Thomas<sup>c</sup>

<sup>a</sup>Laboratoire de Mécanique des Structures et des Systèmes Couplés (LMSSC), Conservatoire national des arts et métiers
(Cnam), 292 Rue Saint-Martin, Paris, 75003, France

<sup>b</sup>Laboratoire d'Electrotechnique et d'Electronique de Puissance (L2EP), Univ. Lille, Arts et Metiers Institute of Technology,
Centrale Lille, Junia, ULR 2697, France

<sup>c</sup>Laboratoire d'Ingénierie des Systèmes Physiques et Numériques (LISPEN), Arts et Métiers Institute of Technology
(Ensam), F-59000 Lille, France

**Abstract**

The present work investigates the ability of a reduced-order fluid-structure model to estimate the vortex-induced vibrations (VIV) of marine lifting surface under hydrodynamic flow, such as hydrofoils. A particular VIV area is scrutinized, for which a hydrodynamic excitation mechanism due to a Kármán-type vortex wake organization successively locks the first twisting and second bending mode of a cantilever truncated hydrofoil. Coupling two structural oscillators with a Van der Pol wake oscillator satisfactorily reproduces the amplitude response and the lock-in frequency. This work also investigates the ability of a piezo-elasto-hydrodynamic model to anticipate the vibration amplitude of the hydrofoil when mitigated thanks to a resonant piezoelectric shunt. Composed of an inductance in series with a resistance connected to a piezoelectric patch, the passive shunt was tuned to minimize the vibration amplitude in the frequency lock-in range. The proposed semi-empirical models are fitted to experimental results in order to reproduce the coupled system's dynamic.

**Keywords:** Vortex-induced vibrations, Vibration mitigation, Resonant piezoelectric shunt, Passive control, Low-order model, Cantilever hydrofoil

# 1. Introduction

Flow-induced vibrations (FIV) of hydrodynamic lifting structures such as hydrofoils have to be considered in order to avoid destructive or fatigue behaviors, to improve hydrodynamic performances and to reduce vibration-induced noises. Mitigation of these vibrations is an important issue in marine engineering and FIV analysis must be performed in order to determine the optimal damping solutions.

According to Blake [1], vibrations of a hydrofoil operating at low angles of attack are due to a random excitation mechanism (generated by the turbulence of the incoming flow or the turbulent boundary layer on the hydrofoil surface), in addition to a harmonic flow excitation generated by a periodic vortex shedding from the trailing edge. In the latter case, referred to as vortex-induced vibration (VIV), tonal noise emission can be observed. Both hydrofoils with either a blunt trailing-edge or a sharp trailing edge are concerned with VIV and tonal noise emissions. For the former, it is the consequence of a Kármán-type vortex wake organization [2]. For the latter, it can be due to a hydro-elastic coupling between a boundary layer transition mechanism occurring close to the trailing edge and a structural mode of the hydrofoil involving significant displacement in the trailing edge area [3]. For both scenarios, vibrations of significant amplitude can be observed in a limited range of the flow velocity, which makes this fluid-structure behavior classified as vortex-induced vibrations.

Regarding structures of circular (or low aspect ratio) cross-section, of great importance for wind and off-shore engineering, numerous phenomenological models have been developed in the past decades [4]. Regarding hydrodynamic lifting surfaces, Toebes & Eagleson [2] have proposed a so-called qualitative fluid-structure model to capture the rotational motion of a rigid flat plate suspended by a twisting spring through its leading edge. This section model considers fluid added mass and damping, based on Theodorsen [5] formulation, in addition to a vortex shedding lifting term introducing a “magnification factor” that depends on the amplitude of motion and a “geometry factor” to take account of the influence of the trailing edge

*Corresponding author
Email addresses: arthur.haudeville@ensam.eu (Arthur Haudeville), xavier.amandolese@lecnam.net (Xavier Amandolese), boris.lossouarn@lecnam.net (Boris Lossouarn), christophe.giraud-audine@ensam.eu (Christophe Giraud-Audine), olivier.thomas@ensam.eu (Olivier Thomas)
geometry. Recently, a hydro-elastic model derived from the Facchinetti, de Langre and Biolley wake oscillator model [6] has been proposed by Amandolese et al. [7], to fit the bending or twisting VIV observed on a cantilever rectangular plate in a hydrodynamic tunnel. Assuming a single mode shape in either bending or twisting (Rayleigh's method) and neglecting spanwise fluid effects, reduced order models coupling a structural oscillator and a wake oscillator have been proposed for either bending or twisting VIV and a set of empirical parameters that fit the experimental data has been identified.

As mentioned previously, vibrations of hydrodynamic lifting surfaces are a matter of great concern. Numerous techniques for vibrations mitigation exist, such as tuned mass dampers (TMD) [8, 9] or electromagnetic and piezoelectric shunts [10, 11]. Vibration reduction applications of submerged structures using piezoelectric devices, by active or passive control, have already demonstrated strong performances [12, 13, 14, 15, 16, 17] and analytic models for linear vibration reduction [18, 19, 20] enable optimization of the mitigation performances. Passive piezoelectric damping is traditionally resistive with simply a resistance (R-shunt) or resonant with a resistance and an inductance in series with the piezoelectric patch (RL-shunt) [21]. Greater vibration reduction is observed with resonant piezoelectric shunts since the shunt will also counteract the oscillations thanks to the occurrence of an electrical resonance [19]. This resonant shunt solution is considered in this study, in line with Watine et al. [12] which provides experimental data obtained in the hydrodynamic tunnel of the French Naval Academy Research Institute (IRENav), on a cantilever hydrofoil of truncated NACA 66-306 cross-section, equipped with a resonant piezoelectric shunt. Following the work of Amandolese et al. [7] on a flat plate without piezoelectric damping, the aim of this paper is to validate an enhanced low-order fluid-structure model based on the experiments. Then, the piezoelectric patch coupling to the structure is added to the system in order to also validate a low-order piezo-elasto-hydrodynamic model. A short range of Reynolds number $2.55 \times 10^5 \le Re \le 4.68 \times 10^5$ is scrutinized for which a Kármán-type vortex wake organization successively locks the first twisting and second bending mode of the cantilever hydrofoil, with and without the resonant shunt control solution.

The paper is organized as follows. The enhanced low-order piezo-elasto-hydrodynamic model, considering hydrofoil mode shapes and a single piezoelectric patch for the mitigation of both the first twisting mode and second bending mode, is developed in section 2. In section 3, two methods for computing the dynamics of the resulting non-linear system are compared before analyzing the effect of the different parameters on the numerical results. Then, parameters that best fit experimental results in both amplitude response and frequency lock-in, without and with control, are identified in section 4, prior to a conclusion.

## 2. Piezo-elasto-hydrodynamic model

This section introduces the piezo-elasto-hydrodynamic model developed to compute the amplitude and frequency responses of a cantilever hydrofoil for which a Kármán-type vortex wake organization successively locks its first two structural modes.

### 2.1. Structural equation of motion

#### 2.1.1. Elastic model of the hydrofoil

<figure><img src="image_2.png" /><figcaption>Figure 1 – (a) Representation of the hydrofoil as a variable thickness plate model of chord length *B*, span length *S* and thickness *h*(*x*), under hydrodynamic flow of velocity *U*, with a uniform fluid forcing per unit length *f*<sub>L</sub>(*t*) induced by the vortex wake organization and acting at the distance *B*<sub>VIV</sub> from the leading edge; the middle plane of the hydrofoil is shown in gray; (b) 3D view of the foil model used for the experiments in [12] and in section 4.</figcaption></figure>

The hydrofoil structure has the form of a cylinder, its base with in the form of a hydrodynamic profile and with its generatrix along the *y*-axis. As shown in Fig. 1(a), it is modeled as a rectangular plate, of span
length S (along the y-axis) and chord length B (along the x-axis). Any point of the middle surface of the plate has coordinate (x, y). The particular 3D geometry of the hydrofoil is taken into account with a variable thickness h(x), that depends only on x because of the cylinder geometry. Assuming that the hydrofoil only deforms in the vertical direction, the motion of the structure is described thanks to a Kirchhoff-Love plate theory by the transverse displacement field w(x, y, t) of the middle surface (i.e. along the z-axis), with t the time, that verifies the following equation [22, 23, 24, 25]:

$$ \mathcal{M}(\ddot{w}) + \mathcal{K}(w) = \delta p, \qquad (1) $$

where $\dot{o} = \partial o / \partial t$, $\delta p(x, y, t)$ [Pa] is the external distributed fluid force per unit area applied to the structure (it is the pressure difference between the upper and lower sides of the hydrofoil, respectively named as suction side and pressure side) and $\mathcal{M}(o)$, $\mathcal{K}(o)$ are mass and stiffness operators, that takes into account the variable thickness h(x) [24] as well as a possible non-homogeneous material because of the piezoelectric patches. In particular, in the case of a constant stiffness plate made in a homogeneous material of mass density $\rho$, Young's modulus E and Poisson's ratio $\nu$, $\mathcal{M}(\ddot{w}) = \rho h \ddot{w}$ and $\mathcal{K}(w) = D \Delta \Delta w$, with $D = Eh^3 / [12(1 - \nu^2)]$ the bending stiffness of the plate and $\Delta$ the Laplacian operator [26]. Figure 1(b) shows and example of hydrofoil, that includes a thin piezoelectric patch embedded in the structure, close to the top face.

In order to introduce a reduced order model, the transverse displacement w is expanded on the natural mode basis of the hydrofoil, truncated to N modes, as follows:

$$ w(x, y, t) \approx \sum_{i=1}^{N} \Phi_i(x, y) q_i(t), \qquad (2) $$

where $\Phi_i$ [-] and $q_i$ [m] are the i-th mode shape and generalized coordinate of the hydrofoil's transverse displacement. Mode shapes are normalized to unit displacement at the tip of the trailing edge of the hydrofoil such that $\Phi_i(B, S) = 1$. Note that hydrofoil vibrations were measured in [12] with a vibrometer focused on a single point located close to the (B, S) point (see Fig. 1(a)). Using the expression for w(x, y, t) from Eq. (2) into Eq. (1), multiplying by $\Phi_j(x, y)$, integrating along the middle surface domain $\mathcal{D} = \{(x, y) \in [0, B] \times [0, S]\}$ (see Fig. 1(a)) and adding structural damping yields, thanks to the orthogonality property of the modes [25], to:

$$ \ddot{q}_j + 2\zeta_j^s \omega_j^s \dot{q}_j + \omega_j^s q_j = \frac{1}{m_j^s} F_j^v, \quad \forall j \in \{1, \dots, N\}, \qquad (3) $$

where $\omega_j^s$ [rad.s$^{-1}$] is the j-th natural frequency, $\zeta_j^s$ is the associated structural damping ratio and $m_j^s$ [kg] is the j-th modal mass of the hydrofoil, defined as:

$$ m_j^s = \iint_{\mathcal{D}} \rho h(x) \Phi_j^2(x, y) d\mathcal{D}. \qquad (4) $$

Eq. (3) also introduces $F_j^v$, the resulting modal forcing term applied on the hydrofoil:

$$ F_j^v = \iint_{\mathcal{D}} \delta p(x, y, t) \Phi_j(x, y) d\mathcal{D}. \qquad (5) $$

### 2.1.2. Fluid-added inertia and damping

A vibrating structure submerged into a hydrodynamic flow exhibits motion-induced loading that can be taken into account as a fluid-added inertia, damping and stiffness. Because the present study concerns a quasi-symmetric hydrofoil at zero angle of attack, the fluid-added stiffness can be neglected. The equations of motions, Eq. (3), can then be expressed as follows:

$$ \ddot{q}_j + 2\zeta_j^{\text{tot}} \omega_j \dot{q}_j + \omega_j^2 q_j = \frac{1}{m_j^{\text{tot}}} F_j^{\text{v}} \quad \forall j \in \{1, \dots, N\}, \qquad (6) $$

where $\omega_j$ is the j-th natural frequency of the hydrofoil in still water (note that $\omega_j \neq \omega_j^s$ because of the fluid-added mass), $\zeta_j^{\text{tot}} = \zeta_j^s + \zeta_j^f$ is the sum of the structural and fluid-added damping ratios and $m_j^{\text{tot}} = m_j^s + m_j^f$ is the j-th total modal mass of the hydrofoil including the fluid-added inertia. Because the motion-induced fluid loading is taken into account as added mass and damping, the fluid forcing on the hydrofoil is now only due to the vortex shedding at the truncated trailing edge of the hydrofoil. The associated modal forcing term, $F_j^v$, will be developed in section 2.2.

Assuming that the mode shapes are the same in still air and water, air to water frequency ratios (identified in still fluid) allow to compute the fluid-added inertia:

$$ m_j^f = m_j^s \left[ \left( \frac{\omega_j^s}{\omega_j} \right)^2 - 1 \right]. \qquad (7) $$
Regarding the fluid-added damping, Theodorsen [5] proposed an analytic solution for a thin airfoil in harmonic plunging and pitching motions. In the present study, mode shapes of the hydrofoil are combinations of twisting and bending deformations, but a dominant deformation in twisting or in bending is observed for the first three modes (refer to [12] and Fig. 2). Thus, fluid-added damping evaluation can be based on the approximation of a cantilever thin airfoil under hydrodynamic flow with pure bending or twisting mode shape. Following the models proposed by Amandolese et al. [7], which are derived from Theodorsen's analytic solution [5], the velocity-dependent fluid-added damping ratios relative to twisting $\zeta_{j,\alpha}^f$ and bending $\zeta_{j,z}^f$ modes are simplified in comparison to [7] and are given linearly dependent of the flow velocity $U$ (see Appendix A for details):

$$ \zeta_{j,\alpha}^f = P_{\alpha} \frac{\pi \rho_f B^3}{64 \omega_j I_j^{\text{tot}}} U, \quad \zeta_{j,z}^f = P_z \frac{\pi \rho_f B}{4 \omega_j M_j^{\text{tot}}} U, \qquad (8) $$

where $\rho_f$ [kg.m$^{-3}$] is the fluid density, $U$ [m.s$^{-1}$] is the mean flow velocity and $0 < P_{\alpha}, P_z < 1$ are weighting coefficients. These weighting coefficients were introduced by Amandolese et al. [7] to moderate the super-position principle of the motion-induced and vortex-induced fluid loading on the structure. Fitted values of these weighting coefficients are given in [7]. They will nevertheless be adjusted in the present study for the case of the hydrofoil's geometry.

<figure><img src="image_3.png" /><figcaption>Figure 2 – (a) Photograph of the aluminum truncated hydrofoil with its piezoelectric patch; deformed shape of (b) the second (first twisting mode) and (c) third (second bending) modes of the cantilever hydrofoil in still water, obtained from a finite element hydro-elastic model [12]. The transparent frame represents the structure without deflection.</figcaption></figure>

Expressions of the fluid added damping ratios, Eq. (8), also introduce $M_j^{\text{tot}}$ [kg.m$^{-1}$] and $I_j^{\text{tot}}$ [kg.m] that are total modal mass and moment of inertia, per unit length, for the cases of pure bending and twisting mode shapes. More details are given in Appendix A. Note that depending on whether the j-th mode shape is dominated by twisting or bending, the fluid-added damping ratio $\zeta_j^f$ that will be used to compute the total damping ratio in Eq. (6) will be either $\zeta_{j,\alpha}^f$ or $\zeta_{j,z}^f$.

## 2.2. Vortex wake forcing model

### 2.2.1. Hydrodynamic pressure field equivalence

The fluid forcing on the hydrofoil is due to the vortex shedding around the truncated trailing edge of the hydrofoil. Some classical VIV models (see [27]) consider only the two-dimensional plunging and pitching motion of the section of the hydrofoil in the $(x, z)$ plane, assumed rigid, whose loading can be considered equivalent to a force applied at a particular point shifted of a distance $B_{\text{viv}}$ from the leading edge.

In our case of a foil with a 3D geometry, we first assume that the fluid loading is uniform in the $y$ (span-wise) direction, such that the load distribution $\delta p(x, t)$ is independent of $y$, a relevant hypothesis in the context of vortex induced vibration. Then, we assume that the mode shapes of the hydrofoil are combinations of bending and twisting motion, with a simplified expression such that each cross section of the foil has a rigid body motion. The j-th mode shape can be then modeled by coefficients $a_j(y)$ and $b_j(y)$ that depend only on the coordinate along the span, such that:

$$ \Phi_j(x, y) = a_j(y) + b_j(y) x. \qquad (9) $$
In that way, the associated modal force can be expressed by introducing Eq. (9) into the modal forcing expression Eq. (5):

$$F_j^v = \int_0^B \delta p(x,t) \, dx \cdot \int_0^S a_j(y) \, dy + \int_0^B x \, \delta p(x,t) \, dx \cdot \int_0^S b_j(y) \, dy, \quad (10a)$$

$$= f_L(t) \int_0^S \Phi_j(B_{\text{viv}}, y) \, dy, \quad (10b)$$

where $f_L(t) = \int_0^B \delta p(x,t) \, dx$ and $B_{\text{viv}} f_L(t) = \int_0^B x \, \delta p(x,t) \, dx$ are the lift force and pitching moment resulting from the pressure distribution. As a consequence, the fluid loading is equivalent to a force per unit length distribution, applied on a line parallel to the y-axis at a distance $B_{\text{viv}}$ from the leading edge, as illustrated in Fig. 1(a). In the following, the lift force distribution in $[\text{N.m}^{-1}]$ is classically written as:

$$f_L(t) = \frac{1}{2} \rho_f U^2 B C_L(t), \quad (11)$$

introducing $C_L(t)$, the time-dependent lift coefficient due to the vortex shedding.

For the real foil shown in Fig. 2, the mode shapes $\Phi_i(x, y)$ do not exactly fulfill Eq. (9) since the cross sections slightly deform. However, this deformation appears negligible at first sight, such that considering the fluid forcing as the line loading of Eq. (11) appears reasonable.

### 2.2.2. Van der Pol oscillator

Numerous VIV phenomenological models have been developed in the last decades and, according to de Langre [4], these models can be classified into three types: the forced system models, the fluid-elastic system models and the coupled system models. This last type of model, referring to a non-linear oscillator of the lift force coupled to a structural equation of motion, is selected here since it is able to take into account both the frequency lock-in and the motion-dependent vortex wake effects. Introduced by Hartlen and Currie [28], it has been improved by Facchinetti et al. [6] through the introduction of a van der Pol damping term in the wake oscillator. More recently, the model has been extended and validated for a cantilever flat plate under hydrodynamic cross-flow by Amandolese et al. [7]. Following this work, a dimensionless time-varying reduced lift coefficient is used,

$$q_v(t) = 2 \frac{C_L(t)}{C_L^0}, \quad (12)$$

where $C_L^0$ is a reference fluctuating lift coefficient. This reduced lift coefficient is the consequence of the Kármán vortex wake organization at the trailing edge of the hydrofoil, whose motion-dependent behavior is well characterized by the following non-linear Van der Pol's equation:

$$\ddot{q}_v + \omega_v \varepsilon (q_v^2 - 1) \dot{q}_v + \omega_v^2 q_v = \sum_{k=1}^{N} A_k \frac{\ddot{q}_k}{H}, \quad (13)$$

where

$$\omega_v = \frac{2\pi St U}{H}, \quad (14)$$

is the wake oscillator frequency, in $[\text{rad.s}^{-1}]$, proportional to the Strouhal number $St$ associated with the Kármán vortex wake organization, $H$ [m] is the maximal thickness of the hydrofoil section, $\varepsilon$ is a growth rate coefficient (also called Van der Pol parameter) and $A_k$ [–] is a coupling parameter scaling the structure acceleration to the wake dynamics of the k-th mode. As a reminder, $N$ is the number of modes in the modal expansion of Eq. (2).

Using expressions (11), (10b) in Eq. (6) and the definition of the reduced lift coefficient, the coupled fluid-structure system describing the vortex-induced vibration of the hydrofoil can be expressed as follows:

$$\ddot{q}_j + 2\zeta_j^{\text{tot}} \omega_j \dot{q}_j + \omega_j^2 q_j = \Psi_j U^2 q_v \quad \forall j \in \{1, \dots, N\}, \quad (15a)$$

$$\ddot{q}_v + \omega_v \varepsilon (q_v^2 - 1) \dot{q}_v + \omega_v^2 q_v = \sum_{k=1}^{N} A_k \frac{\ddot{q}_k}{H}, \quad (15b)$$

where we introduce $\Psi_j$ the coupling parameter scaling the wake dynamics to the structure motion. Thanks to the expression of the fluid forcing in Eq. (10b) and of the reduced lift coefficient of the Van der Pol equation in Eq. (12), the development of the forcing term of the structural oscillator in Eq. (6) yields to:

$$\Psi_j = \frac{\rho_f B C_L^0}{4 m_j^{\text{tot}}} S_{\Phi_j} \quad \text{with} \quad S_{\Phi_j} = \int_0^S \Phi_j(B_{\text{viv}}, y) \, dy, \quad (16)$$

where the variable of the flow velocity $U$ has been sidelined and with the proposed notation of $S_{\Phi_j}$ to reduce this expression of the coupling parameter $\Psi_j$ and the following.
## 2.3. Resonant piezoelectric shunt damping

<figure><img src="image_4.png" /><figcaption>Figure 3 – (a) Vibrating elastic structure coupled to a piezoelectric shunt composed of an inductance L and a resistance R in series; (b) A piezoelectric patch with Q, $\dot{Q}$ and V the electric charge, the current and the difference of electric potentials, respectively.</figcaption></figure>

In line with the experimental work of Watine et al. [12], a passive resonant piezoelectric shunt damping is considered in this model. Using the piezoelectric effect, mechanical energy of a vibrating structure is converted into electrical energy and vice-versa. To do so, as in Fig. 3, a piezoelectric patch is coupled to the structure and the patch is connected to a resistance in series with an inductance. We denote $V(t)$ [V] the electrical potential difference (voltage) between the electrodes of the piezoelectric patch and $Q(t)$ [C] the electrical charge contained in one of the patch electrodes.

This piezoelectric patch is connected to a dedicated electronic circuit which is able to damp the vibration of the vibrating structure. In the case of a resistive shunt (a R-shunt, the circuit includes only a resistor), the electrical energy is simply dissipated in the resistor through Joule effect. Then, adding an inductance in series with the resistor in the electric circuit, one obtains a resonant shunt (RL-shunt), equivalent to a mechanical tuned mass damper since the inductance resonates with the capacitance of the piezoelectric patch and counteracts the structural oscillations. If the values of R and L are correctly tuned on a given resonance of the primary structure, its frequency response function has the form of a blunted resonance, thus reducing amplitude of vibrations in its vicinity [21, 19, 10].

### 2.3.1. Electromechanical equations

Let us consider a single piezoelectric patch applied on the hydrofoil (as in Fig. 2) where the modal model of a linear electromechanical coupling between the j-th mode of the structure and the piezoelectric patch corresponds to Eq. (6) in which the piezoelectric forcing is introduced [18, 29, 30],

$$ \ddot{q}_j + 2\zeta_j^{\text{tot}}\omega_j \dot{q}_j + \omega_j^2 q_j = \frac{1}{m_j^{\text{tot}}} F_j^{\text{v}} + \frac{\chi_j}{m_j^{\text{tot}}} V \quad \forall j \in \{1, \dots, N\}, \quad (17a) $$

$$ CV - Q + \sum_{k=1}^{N} \chi_k q_k = 0, \quad (17b) $$

where the modal piezoelectric forcing term is defined with $\chi_j$ [C.m$^{-1}$] the piezoelectric coupling coefficient relative to the j-th mode and C is the electric capacitance of the patch. For a refined model, C must be taken as the average value of the dynamic capacitance in the vicinity of the k-th resonance, as defined in [31]. One can note that, in the case of a resonant piezoelectric shunt, the electrical potential applied to the piezoelectric transducer is $V = -L\ddot{Q} - R\dot{Q}$. When introduced into Eq. (17b), the electromechanical equation becomes:

$$ \frac{1}{\omega_e^2} \ddot{Q} + \frac{2\zeta_e}{\omega_e} \dot{Q} + Q = \sum_{k=1}^{N} \chi_k q_k, \quad (18) $$

where $\omega_e = 1/\sqrt{LC}$ is the electrical frequency and $\zeta_e = R/(2\sqrt{C/L})$ is the electrical damping ratio. Moreover, for better convenience, the electrical charge can be set as the electrical unknown in Eq. (17a). Then, by introducing Eq. (17b) into Eq. (17a) to eliminate $V(t)$ and assuming that the electrical potential is truncated at the j-th mode in open-circuit, one obtains:

$$ \ddot{q}_j + 2\zeta_j^{\text{tot}}\omega_j \dot{q}_j + \dot{\omega}_j^2 q_j = \frac{1}{m_j^{\text{tot}}} F_j^{\text{v}} + \frac{\omega_j k_j}{\sqrt{m_j^{\text{tot}}C}} Q \quad \forall j \in \{1, \dots, N\}, \quad (19) $$
with $\hat{\omega}_j$ the open-circuit natural frequency and $k_j$ the modal electromechanical coupling factor, both relative to the j-th mode and defined as:

$$ \hat{\omega}_j^2 = \omega_j^2 + \frac{\chi_j^2}{m_j^{\text{tot}} C} = \omega_j^2 (1 + k_j^2), \quad k_j = \frac{1}{\omega_j} \frac{\chi_j}{\sqrt{m_j^{\text{tot}} C}} = \sqrt{\frac{\hat{\omega}_j^2 - \omega_j^2}{\omega_j^2}}. \qquad (20) $$

Therefore, the VIV piezo-elasto-hydrodynamic model of the hydrofoil, in Eq. (15), with a single resonant piezoelectric shunt, truncated to N modes, is:

$$ \left\{ \begin{array}{l} \ddot{q}_j + 2\zeta_j^{\text{tot}}\omega_j \dot{q}_j + \hat{\omega}_j^2 q_j = \Psi_j U^2 q_v + \frac{\omega_j k_j}{\sqrt{m_j^{\text{tot}} C}} Q \quad \forall j \in \{1, \dots, N\}, \\ \ddot{q}_v + \omega_v \varepsilon(q_v^2 - 1) \dot{q}_v + \omega_v^2 q_v = \sum_{k=1}^{N} A_k \frac{\ddot{q}_k}{H}, \end{array} \right. \qquad (21a) $$

$$ \ddot{q}_v + \omega_v \varepsilon(q_v^2 - 1) \dot{q}_v + \omega_v^2 q_v = \sum_{k=1}^{N} A_k \frac{\ddot{q}_k}{H}, \qquad (21b) $$

$$ \frac{1}{\omega_e^2} \ddot{Q} + \frac{2\zeta_e}{\omega_e} \dot{Q} + Q = \sum_{k=1}^{N} \omega_k k_k \sqrt{m_k^{\text{tot}} C} q_k . \qquad (21c) $$

### 2.3.2. Dimensionless equations

The piezo-elasto-hydrodynamic model can be written in dimensionless form by introducing a dimensionless structural generalized coordinate based on the maximal thickness of the hydrofoil's section H and setting a dimensionless time based on the i-th ($i \in \{1, \dots, N\}$) structural frequency, as:

$$ \bar{q}_j = \frac{q_j}{H}, \quad \tau = \omega_i t. \qquad (22) $$

Resulting dimensionless angular frequencies are denoted by an over-bar $\bar{\omega}_n = \omega_n/\omega_i$ with $n = \{j, v, e\}$. Moreover, one should note that the fluid-to-structure frequency ratio can be expressed as $\omega_v/\omega_i = St U_{r_i}$, where $U_{r_i} = 2\pi U/(\omega_i H)$ is the reduced flow velocity. Thus, the previously defined coupling parameters that scale the wake dynamics to the structure motion in Eq. (16) can be rewritten thanks to the introduction of a dimensionless mass number parameters $\mu_j$:

$$ \bar{\Psi}_j = \frac{C_L^0}{16\pi^2 St^2 \mu_j} \quad \text{with} \quad \mu_j = \frac{m_j^{\text{tot}}}{\rho_f B H S_{\Phi_j}}, \qquad (23) $$

In the same way, the fluid-added damping ratios in Eq. (8) can be rewritten with the total moment of inertia $I_{\alpha}^{\text{tot}}$ and mass $M_{j,z}^{\text{tot}}$ per unit length introduced in Appendix A:

$$ \zeta_{j,\alpha}^{f} = \begin{cases} C_{\alpha} \bar{\omega}_v & \text{if } j = i, \\ C_{\alpha} \frac{\bar{\omega}_v}{\bar{\omega}_j} & \text{if } j \neq i, \end{cases} \quad \text{with} \quad C_{\alpha} = \frac{P_{\alpha}}{128 St \mu_{\alpha}} \quad \text{and} \quad \mu_{\alpha} = \frac{I_{\alpha}^{\text{tot}}}{\rho_{\text{f}} B^3 H}, \qquad (24a) $$

$$ \zeta_{j,z}^{f} = \begin{cases} C_{j,z} \bar{\omega}_v & \text{if } j = i, \\ C_{j,z} \frac{\bar{\omega}_v}{\bar{\omega}_j} & \text{if } j \neq i, \end{cases} \quad \text{with} \quad C_{j,z} = \frac{P_z}{8 St \mu_{j,z}} \quad \text{and} \quad \mu_{j,z} = \frac{M_{j,z}^{\text{tot}}}{\rho_f B H}. \qquad (24b) $$

For the electrical oscillator, the following change of variable is proposed to render the electrical charge variable dimensionless,

$$ \bar{Q} = \frac{1}{H} \frac{1}{\omega_i \sqrt{m_i^{\text{tot}} C}} Q. \qquad (25) $$

Then, by introducing the modal mass ratios $\bar{m}_{i,j} = \sqrt{m_i^{\text{tot}}/m_j^{\text{tot}}}$ and $\bar{m}_{k,i} = \sqrt{m_k^{\text{tot}}/m_i^{\text{tot}}}$, the dimensionless form of the piezo-elasto-hydrodynamic model can be written as follows:

$$ \ddot{\bar{q}}_j + 2\zeta_j^{\text{tot}}\bar{\omega}_j \dot{\bar{q}}_j + \hat{\omega}_j^2 \bar{q}_j = \bar{\Psi}_j \bar{\omega}_v^2 q_v + k_j \bar{\omega}_j \bar{m}_{i,j} \bar{Q} \quad \forall j \in \{1, \dots, N\}, \qquad (26a) $$

$$ \ddot{q}_v + \bar{\omega}_v \varepsilon(q_v^2 - 1) \dot{q}_v + \bar{\omega}_v^2 q_v = \sum_{k=1}^{N} A_k \ddot{q}_k, \qquad (26b) $$

$$ \frac{1}{\bar{\omega}_e^2} \ddot{\bar{Q}} + \frac{2\zeta_e}{\bar{\omega}_e} \dot{\bar{Q}} + \bar{Q} = \sum_{k=1}^{N} k_k \bar{\omega}_k \bar{m}_{k,i} \bar{q}_k. \qquad (26c) $$
Experimental study in [12] focused on a narrow range of Reynolds number, $2.55 \times 10^5 \le Re \le 4.68 \times 10^5$, for which a Kármán-type vortex wake organization successively locks the first twisting mode (i.e. the second natural mode of the hydrofoil) and second bending mode (i.e. the third natural mode of the hydrofoil). If the lock-in with the first twisting mode has the most significant amplitude response and lock-in range, lock-in with the second bending mode is also observed, although not clearly separated from the first twisting VIV. Both the structural modes must be taken into account simultaneously. The experimental study in [12] also focuses on a mitigation of the first twisting mode (the second natural mode). In agreement with this, the dimensionless time in Eq. (22) and the electrical charge change of variable in Eq. (25) have been defined with respect to the second natural mode of the hydrofoil, i.e. $i = 2$. Thus, the generalized dimensionless piezo-elasto-hydrodynamic system of Eqs. (26a, 26b, 26c) becomes:

$$ \ddot{q}_2 + 2(\zeta_2^s + C_\alpha \bar{\omega}_v) \dot{q}_2 + (1 + k_2^2) \bar{q}_2 = \bar{\Psi}_2 \bar{\omega}_v^2 q_v + k_2 \bar{Q}, \quad (27a) $$

$$ \ddot{q}_3 + 2(\zeta_3^s \bar{\omega}_3 + C_{3,z} \bar{\omega}_v) \dot{q}_3 + \bar{\omega}_3^2 \bar{q}_3 = \bar{\Psi}_3 \bar{\omega}_v^2 q_v + k_3 \bar{\omega}_3 \bar{m}_{2,3} \bar{Q}, \quad (27b) $$

$$ \ddot{q}_v + \bar{\omega}_v \varepsilon(q_v^2 - 1) \dot{q}_v + \bar{\omega}_v^2 q_v = A_2 \ddot{q}_2 + A_3 \ddot{q}_3, \quad (27c) $$

$$ \frac{1}{\bar{\omega}_e^2} \ddot{Q} + \frac{2\zeta_e}{\bar{\omega}_e} \dot{Q} + \bar{Q} = k_2 \bar{q}_2 + k_3 \bar{\omega}_3 \bar{m}_{3,2} \bar{q}_3. \quad (27d) $$

Note that in the following, all numerical and experimental results will consider the Root Mean Square (RMS) values of the transverse dimensionless displacement $\bar{w}_{\text{tip}}(\tau) = \sum_{j=1}^{N} \Phi_j(B, S) q_j(\tau)/H$ signal, computed over several periods in the periodic steady state. Since at most two structural modes are considered, the modal expansion (2) is truncated to the two modes involved, and the transverse dimensionless displacement at the tip of the trailing edge of the hydrofoil (considering the normalization $\Phi_j(B, S) = 1$) is:

$$ \bar{w}_{\text{tip}}(\tau) = \Phi_2(B, S) \bar{q}_2(\tau) + \Phi_3(B, S) \bar{q}_3(\tau) = \bar{q}_2(\tau) + \bar{q}_3(\tau). \quad (28) $$

### 3. Numerical analysis

In the present section, we investigate the behavior of the VIV model using two complementary methods for nonlinear systems solving. We are interested in the vibratory response of the system when the mean flow velocity $U$ is varied around a range of values for which VIV appear, around the resonances of the structural modes. Instead of using $U$ as the control parameter, we equivalently use the wake oscillation frequency $\omega_v$, thanks to the proportionality relation (14).

The first numerical solving relies on standard time integration, thanks to a Runge-Kutta fourth-order scheme [32]. This method necessitates formalizing equations in the form $\dot{X} = f(X, \tau)$. It is used with the time-step fixed to $d\tau = 1/16$ to ensure the convergence of the solution. For each mean flow velocity $U$ or wake oscillation frequency $\omega_v$ (increasing or decreasing by step), and after a transient regime, a periodic solution is observed for all the variables of the system. In the following, we will focus on the RMS value of the tip amplitude of the foil $\bar{w}_{\text{tip}}(t)$. The frequency $\omega$ of the periodic solution, associated to the dominant harmonic solution, is also of great interest and extracted from the power spectral density (PSD) of the solution.

The second numerical method is the numerical continuation of periodic solutions, based on the Harmonic Balance Method (HBM) and an asymptotic numerical method (ANM). Continuation method also permits computing the solutions of non-linear ODEs, but targets directly the periodic steady state and compute it as a function of a control parameter. We used here the Manlab software, which computes the numerical continuation of periodic solutions of the non-linear ODEs given in the form $\dot{X} = f(X, \tau, \lambda)$, where $\lambda \in \mathbb{R}$ is the control parameter (also called continuation parameter), taken here equal to $U$ or $\omega_v$. The solver is based on the harmonic balance method (HBM) and the asymptotic numerical method [33, 34, 35]. It also uses a particular method to compute the bifurcation points and the stability of the branches [36, 37]. It naturally targets the periodic solutions of the system, enabling to compute the values of all the harmonics of all variables (and thus their RMS amplitude), as well as the frequency $\omega$ of the periodic solution.

#### 3.1. Comparison of numerical methods: time-integration and continuation

Let us consider the VIV model of Eqs. (27a, 27c), with a single twisting mode and the piezoelectric patch in short-circuit. The patch in short-circuit corresponds to the electromechanical system of Eqs. (17a, 17b) with $V(t) = 0$ and, to obtain an equivalent behavior, one can set the coupling factors to $k_2 = k_3 = 0$. The two numerical methods are compared with a set of values for the main parameters of the system: $\zeta_2^s, C_\alpha, \bar{\Psi}_2, \varepsilon$ and $A_2$. Values are chosen so that we can observe nonlinear effects in the amplitude and frequency responses. Dimensionless amplitudes of vibrations are plotted as a function of the dimensionless fluid frequency (i.e. the reduced flow velocity $U_r$). Fig. 4(a) shows amplitude and frequency hysteresis results obtained using the numerical integration method by sweeping the fluid frequency $\bar{\omega}_v$ forward and
<figure><img src="image_5.png" /><figcaption>Figure 4 – Typical RMS amplitude of the hydrofoil tip transverse oscillations $\bar{w}_{\text{tip}}(t)$ and frequency $\omega$ of the periodic solution as a function of the control parameter $\bar{\omega}_v$, for the coupled fluid-structure model in Eqs. (27a, 27c), using $k_2 = k_3 = 0$ (without control and with a single structural degree of freedom). The values of the parameters are $\zeta_2^s = 3 \times 10^{-3}$, $C_\alpha = 2 \times 10^{-2}$, $\bar{\Psi}_2 = 3 \times 10^{-4}$, $\varepsilon = 0.3$ and $A_2 = 120$; (a) Results obtained using the time-integration method (arrows indicate the direction of the jumps); (b) Results obtained using the continuation method.</figcaption></figure>

backward. It shows that jumps with specific directions between the solutions occur for different values of the driven frequency. These regions correspond to unstable solutions for which amplitude and system frequency depend on the direction of sweeping of the fluid frequency. These instabilities in the solutions are confirmed by the continuation method, in Fig. 4(b), where unstable regions are distinguished from stable ones with dotted lines. For practical reasons related to a better observation of the unstable solutions, the continuation method is preferred in this study and following results are obtained using the Manlab software. Such amplitude hysteresis from VIV modeled by a Van der Pol oscillator coupled to a structure is mathematically detailed in [38].

Regarding the frequency responses, we clearly observe the linear evolution of the system frequency within the flow velocity (for $\bar{\omega}_v < 0.8$ with $\bar{\omega}_v = St U_{r2}$) until the lock-in between the fluid and the structural frequencies happens. The lock-in effect is characterized by the synchronization between the fluid and the structural frequencies [39], resulting to a normalized system frequency as observed in the plotted frequency responses. With the fluid frequency increasing, once the lock-in range passed (for $\bar{\omega}_v > 1.3$), the system frequency resumes its linear evolution with the flow velocity. For better readability, Figs. 4(a) and (b) are plotted separately, but it has been verified that merging the two responses obtained by the two numerical methods leads to the same response away from the unstable parts, for which the continuation method gives the right hysteretic behavior.

## 3.2. Influences of the VIV model parameters

Again, we consider the same VIV model as previously of Eqs. (27a, 27c) with a single twisting mode and the piezoelectric patch in short-circuit (which is equivalent to setting $k_j = 0$). This model is composed of unknown parameters, namely $C_\alpha$, $\bar{\Psi}_2$, $\varepsilon$ and $A_2$, and their effects on the amplitude and frequency of vibrations are scrutinized here.

Fig. 5 shows the effects of $C_\alpha$, a parameter proportional to the weighting coefficient of the fluid-added damping expressed in Eq. (24a), and $\bar{\Psi}_2$ which is the coupling parameter scaling the structure acceleration to the wake dynamics introduced in Eq. (23). Both are fluid-related structural parameters affecting the structural equation of motion Eq. (27a). Fig. 5(a) shows that when decreasing the $C_\alpha$ parameter, it increases the amplitude of vibration and the lock-in range as a function of the fluid frequency. The results show that
unstable solution segments are greatly impacted by this damping parameter and thus the non-linear behavior of the system. Fig. 5(b) shows that the maximal amplitude response and the lock-in range increase quasi-linearly when increasing the $\bar{\Psi}_2$ parameter. Unlike the damping parameter $C_\alpha$ that the coupling parameter has less effect on the unstable solution segments. In both cases, by varying $C_\alpha$ or $\bar{\Psi}_2$, one can also observe a right-side shift of the maximal amplitude reaching point with the growing amplitude of vibration. This is a characteristic amplitude behavior of a system coupled to the non-linear Van der Pol's oscillator.

<figure><img src="image_6.png" /><figcaption>Figure 5 – Amplitude and frequency response of the coupled fluid-structure model in Eqs. (27a, 27c) using $k_2 = k_3 = 0$ (without control and with a single structural degree of freedom): $\zeta_3^2 = 3 \times 10^{-3}$, $\varepsilon = 0.3$, $A_2 = 45$; Solid lines (—) represent a stable solution and dotted lines (····) represent unstable ones; (a) Effect of the fluid damping parameter $C_\alpha$ (three values are selected, see the legend) on the system behavior, with $\bar{\Psi}_2 = 3 \times 10^{-4}$; (b) Effect of the coupling parameters scaling the structure acceleration to the wake dynamics $\bar{\Psi}_2$ (three values are selected, see the legend) on the system behavior, with $C_\alpha = 2 \times 10^{-2}$.</figcaption></figure>

Fig. 6 shows the effects of the growth rate coefficient $\varepsilon$ from the Van der Pol equation in Eq. (13), and the coupling parameter $A_2$ scaling the structure acceleration to the wake dynamics introduced from the same Van der Pol equation. Both of these parameters affect the Van der Pol wake oscillator equation Eq. (27c) of the system. Fig. 6(a) shows that the increase of the growth rate parameter $\varepsilon$ slightly impacts the system's lock-in frequency range while shifting it to higher values of the flow velocity (i.e. the fluid frequency $\bar{\omega}_v$). Regarding the maximum amplitude of vibrations, the behavior is complex and not developed here. Similarly, there exists a relation between the growth rate parameter $\varepsilon$ and the fluid frequency which could allow to fix its value [40]. Investigations will be performed in a future study. Fig. 6(b) shows that both the amplitude response and the lock-in range increase within the coupling parameter $A_2$. In comparison with the effect of the coupling parameter $\bar{\Psi}_2$ exposes in Fig. 5, one can note that $A_2$ has a stronger influence on the lock-in range with less impact on the maximal amplitude of vibrations.

Tuning effects of the electrical resistance $R$ and the inductance $L$, on the dynamical behaviors and the shunt damping performances, are left for the upcoming study in the subsection 4.4.

## 4. Experimental results and model identifications

In this section, the piezo-elasto-hydrodynamic model is updated based on the experimental results of Wamine et al. [12]. Relevant experimental parameters are recalled and an evaluation of mode-shape-dependent parameters is proposed before fitting the VIV model (i.e. piezoelectric patch in short-circuit) and the piezo-elasto-hydrodynamic model (i.e. piezoelectric patch with RL shunt) with the experimental data.
<figure><img src="image_7.png" /><figcaption>Figure 6 – Amplitude and frequency response of the coupled fluid-structure model in Eqs. (27a, 27c) using $k_2 = k_3 = 0$ (without control and with a single structural degree of freedom): $\zeta_2^s = 3 \times 10^{-3}$, $C_\alpha = 2 \times 10^{-2}$, $\bar{\Psi}_2 = 3 \times 10^{-4}$; Solid lines (—) represent a stable solution and dotted lines (····) represent unstable ones; (a) Effect of the fluid damping parameter $\varepsilon$ (three values are selected, see the legend) on the system behavior, with $A_2 = 45$; (b) Effect of the coupling parameters scaling the structure acceleration to the wake dynamics $A_2$ (three values are selected, see the legend) on the system behavior, with $\varepsilon = 0.3$.</figcaption></figure>

### 4.1. Experimental data

The cantilever hydrofoil studied in [12] is made of aluminum (see picture in Fig. 2(a)). It has a span length $S = 191$ mm and a NACA 66-306 section truncated at 85% of the original chord, resulting in a chord length of $B = 85$ mm. The maximum thickness, located at 45.46% of the original chord, is chosen as the characteristic length $H = 6.13$ mm that includes an extra thickness of 0.13 mm due to the surface treatment applied to reduce corrosion. The truncated cantilever hydrofoil with chord-to-thickness ratio $B/D = 13.9$ is set at zero angle of attack (maximal uncertainty of 0.1 deg). Natural modes of the hydrofoil have been measured in still water (shapes of the second and third modes in Fig. 2(b, c)). In the present study we will focus on the VIV of the first twisting mode (i.e. the second natural mode) for which $f_2 = 281.2$ [Hz] with $\zeta_2^s = 0.32\%$ and of the second bending mode (i.e. the third natural mode) for which. $f_3 = 302.8$ [Hz] with $\zeta_3^s = 0.67\%$. Relevant structural parameters are summarized in Tab. 1.

Measured vibrations of the hydrofoil by Watine et al. [12] are given as velocity time series for which the PSD of the derived displacement signal is computed. The RMS values of the vibration amplitude are obtained by a local integration around the dominant peak of the PSD so that the impact of residual turbulence-induced vibrations is reduced. The RMS dimensionless vibration amplitudes of the hydrofoil and associated dominant frequencies are plotted in Fig. 7 as a function of the flow velocity for both with and without control.

Experimental determination of the Strouhal number in [12] reflects two fluid excitation mechanisms which are characterized by a Strouhal number $St_2 = 0.34$ (due to a Kármán vortex shedding at the trailing edge) and $St_1 = 0.49$ (most likely due to an impinging shear layer instability mechanism from the leading edge). To simplify the VIV model, only one fluid excitation mechanism has been considered in the present study (i.e. a single wake oscillator equation in Eqs. (27a, 27b, 27c, 27d)) and an intermediate Strouhal number is fixed to $St = 0.41$. As mentioned in the previous section, a relation may exist between the Strouhal number and the growth rate coefficient of the Van der Pol equation, as $\varepsilon$ modifies the fluid natural frequencies and thus the Strouhal number [41, 42]. In an objective to estimate vibration amplitudes, it is nevertheless sufficient to fix the Strouhal number and adapt the growth rate coefficient to fit back lock-in frequencies. Note that if the Strouhal number $St_1 = 0.49$, due to the Kármán vortex shedding at the trailing edge, would
<table><caption>Table 1 – Main structural parameters of the hydrofoil in still water and air, as provided in [12].</caption><thead><tr><th>Quantity</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Chord length</td><td><i>B</i></td><td>85</td><td>mm</td></tr><tr><td>Maximal thickness</td><td><i>H</i></td><td>6.13</td><td>mm</td></tr><tr><td>Span length</td><td><i>S</i></td><td>191</td><td>mm</td></tr><tr><td>Aluminum density</td><td>ρ</td><td>2650</td><td>kg.m<sup>-3</sup></td></tr><tr><td>Water density</td><td>ρ<sub>f</sub></td><td>1000</td><td>kg.m<sup>-3</sup></td></tr><tr><td>Second mode frequency in water</td><td><i>f</i><sub>2</sub></td><td>281.2</td><td>Hz</td></tr><tr><td>Second mode frequency in air</td><td><i>f</i><sub>2</sub><sup>air</sup></td><td>518.4</td><td>Hz</td></tr><tr><td>Second mode damping in still water</td><td>ζ<sub>2</sub><sup>s</sup></td><td>0.32</td><td>%</td></tr><tr><td>Third mode frequency in water</td><td><i>f</i><sub>3</sub></td><td>302.8</td><td>Hz</td></tr><tr><td>Third mode frequency in air</td><td><i>f</i><sub>3</sub><sup>air</sup></td><td>637.2</td><td>Hz</td></tr><tr><td>Third mode damping in still water</td><td>ζ<sub>3</sub><sup>s</sup></td><td>0.67</td><td>%</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 7 – Experimental data provided by Watine et al. [12] in the hydrodynamic tunnel for the hydrofoil without and with piezoelectric shunt; (a) RMS dimensionless amplitudes of vibration of the hydrofoil as a function of the flow velocity; (b) Dominant vibration frequencies of the hydrofoil as a function of the flow velocity; Natural frequency $f_2$ and $f_3$ are plotted in dashed lines (---).</figcaption></figure>

have been kept, the growth rate coefficient $\varepsilon$ enabling to place the lock-in range at the right frequencies in comparison to the experimental data would have been fitted to values greater than one. Yet, applications of the Van der Pol equation in VIV models [6, 7, 43] tend to fix the growth rate coefficient $\varepsilon < 1$ and so the same limit is applied here.

Electromechanical coupling factors of the piezoelectric patch have also been determined in [12]. $k_2 = 8.5\%$ will be used for the second natural mode and $k_3 = 4.7\%$ for the third natural mode. Values of the electrical inductance $L$, resistance $R$ and equivalent capacitance of the piezoelectric patch $C$ are however considered as unknowns because of the numerous uncertainties encountered when evaluating the equivalent electrical properties (non-linearities, leakage...).

## 4.2. Evaluation of mode shapes dependent parameters

Mode shape-dependent parameters have been computed with the help of the finite element solver of Comsol Multiphysics [44] allowing a precise definition of the mode shapes. Considering a mass matrix $M$ and an eigenvector $\Phi_j$, we are interested in computing the modal mass of the $j$-th mode in kg, which is given by the relation:

$$m_j^s = \Phi_j^T M \Phi_j. \qquad (29)$$

The solver of finite element softwares usually provide eigenvectors $\tilde{\Phi}_j$ scaled by the mass matrix, such that:

$$\tilde{\Phi}_j^T M \tilde{\Phi}_j = 1. \qquad (30)$$

Considering that $\tilde{\Phi}_j = \alpha \Phi_j$, with $\alpha \in \mathbb{R}$ and that we choose the scaling of the mode such that $\Phi_j|_{x=B,y=S} = \Phi_j(B, S) = 1$, one has $\alpha = \tilde{\Phi}_j|_{x=B,y=S} = \tilde{\Phi}_j(B, S)$. Then, Eqs. (29) and (30) simply give the $j$-th modal mass, defined in Eq. (4), as:

$$m_j^s = \frac{1}{(\tilde{\Phi}_j|_{x=B,y=S})^2}, \qquad (31)$$
where $\tilde{\Phi}_j|_{x=B,y=S}$ is obtained as the value, at the edge of the hydrofoil, in $(x, y) = (B, S)$, of the j-th eigenvector computed by the finite element code and mass normalized taken.

Using the same approach, one can also calculate the integral along the span of the j-th mode shape evaluated at the application line of the lifting force $S_{\Phi_j}$ [m] (previously introduced in Eq. (16)), as follows:

$$S_{\Phi_j} = \frac{\int_0^S \tilde{\Phi}_j(B_{\text{viv}}, y) \, dy}{\tilde{\Phi}_j(B_{\text{viv}}, S)} \quad (32)$$

Regarding the coupling parameters $\Psi_2$ and $\Psi_3$ which depend on $S_{\Phi_j}$, so on the location of $B_{\text{viv}}$, values have been computed by numerical interpolation of the mode shape $\tilde{\Phi}_j$. Results from previous calculation and those from Appendix A are reported in Tab. 2.

<table><caption>Table 2 – Evaluated inertial parameters of the hydrofoil using finite element model.</caption><thead><tr><th>Quantity</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Structural mass per unit length</td><td>$M_2^s$</td><td>1.10</td><td>kg.m<sup>-1</sup></td></tr><tr><td>Fluid-added mass per unit length</td><td>$M_{3,z}^f$</td><td>3.78</td><td>kg.m<sup>-1</sup></td></tr><tr><td>Structural moment of inertia per unit length</td><td>$I_\alpha^s$</td><td>4.73 &times; 10<sup>-3</sup></td><td>kg.m</td></tr><tr><td>Fluid-added moment of inertia per unit length</td><td>$I_\alpha^f$</td><td>3.40 &times; 10<sup>-3</sup></td><td>kg.m</td></tr><tr><td>Second modal mass</td><td>$m_2^s$</td><td>2.88 &times; 10<sup>-2</sup></td><td>kg</td></tr><tr><td>Second mode fluid-added mass</td><td>$m_2^f$</td><td>6.90 &times; 10<sup>-2</sup></td><td>kg</td></tr><tr><td>Third modal mass</td><td>$m_3^s$</td><td>5.61 &times; 10<sup>-2</sup></td><td>kg</td></tr><tr><td>Third mode fluid-added mass</td><td>$m_3^f$</td><td>1.92 &times; 10<sup>-1</sup></td><td>kg</td></tr></tbody></table>

### 4.3. Identification of the vortex-induced vibrations unknowns parameters

Regarding the developed piezo-elasto-hydrodynamic system in Eqs. (27a, 27b, 27c, 27d)), with the first two equations related to the structural oscillators, the finite elements model and experimental tests allow to fix the structural damping $\zeta_2^s$ and $\zeta_2^s$, as well as the dimensionless frequency $\bar{\omega}_3$ and mass $\bar{m}_{2,3}$, and also the modal coupling factors $k_2$ and $k_3$. The remaining parameters to fit are therefore the fluid-added damping ratios $C_\alpha$ and $C_{3,z}$ because of the weighting parameters $P_\alpha$ and $P_z$, introduced in Eq. (8), but also the coupling parameters $\bar{\Psi}_2$ and $\bar{\Psi}_3$ because of the static lift coefficient $C_L^0$ and the distance from the leading edge to the pressure center $B_{\text{viv}}$, as shown in Eq. (23). About the von der Pol equation, each of the three parameters, namely the growth rate parameter $\varepsilon$ and the coupling parameters scaling the structure acceleration to the wake dynamics $A_2$ and $A_3$ are unknown. Considering the electrical equation, details are given in the following section about the electrical frequency $\bar{\omega}_e$ and damping ratio $\zeta_e$. To summarize, remaining fluid-structure unknown parameters to fit include $P_\alpha$, $P_z$, $C_L^0$, $B_{\text{viv}}$, $A_2$, $A_3$ and $\varepsilon$.

We initialized numerical calculations with the value of the Van der Pol parameter set to $\varepsilon = 0.3$, which is known to be suitable for a structure with a circular cross section in water [6, 43]. The fluctuating lift coefficient is initialized to $C_L^0 = 0.008$, in accordance with the best fitted value for a cantilever flat plate [7] of geometrical dimensions close to those of this studied hydrofoil. The fluid-added damping weighting parameters are also initialized to the values fitted for a flat plate: $P_\alpha = 0.25$ and $P_z = 0.15$. The coupling parameters scaling the structure acceleration to the wake dynamics are initialized to $A_2 = 10$ and $A_3 = 10$, in the same order of magnitude as in [6]. Finally, the center of pressure of the lift force per unit length is first located at the trailing edge of the hydrofoil, such that $B_{\text{viv}} = B$, where it is expected to be applied. As mentioned previously, the Strouhal number has been fixed to $St = 0.41$ and one should start to fit the Van der Pol parameter $\varepsilon$ in order to conform to the experimental data lock-in range.

Results reported in Fig. 8(a) indicate that a better lock-in range location is reached for $\varepsilon = 0.7$ by anticipating the incoming lock-in range expansion because of the other parameters. The fitting is done by trial and error. Then, maximum vibration amplitudes reached by the second and third modes are balanced thanks to the fluid-added damping weighting parameters. Fig. 8(b) shows that $P_\alpha = 0.5$ and $P_z = 0.075$ best adjust the difference in maximal amplitude between the two modes. Compared to the weighting parameters fitted by Amandolese et al. [7] for a cantilever flat plate, $P_z$ is lower in our case and the difference can be attributed to the fact that we are working with the second bending mode of the structure, and not with the first one as for the flat plate. For the weighting parameters relative to the twisting mode, $P_\alpha$ is however higher than the value identified for the flat plate.

Next, one should fit the coupling parameters scaling the structure acceleration to the wake dynamics $A_2$ and $A_3$ in order to adjust the extent of the lock-in ranges without impacting too strongly the maximal vibration amplitudes. Fig. 9(a) shows that best fitted coupling parameters values are for $A_2 = 20$ and
<figure><img src="image_9.png" /><figcaption>Figure 8 – Amplitude and frequency responses of the coupled fluid-structure model in Eqs. (27a, 27b, 27c) using $k_2 = k_3 = 0$, in comparison with experimental data for the short-circuited shunt configuration, with $C_L^0 = 0.008$, $B_{viv} = B$, $A_2 = 10$ and $A_3 = 10$; (a) Influence of the Van der Pol parameter $\varepsilon$ (three values are selected, see the legend) on the system behavior, with $P_\alpha = 0.25$, $P_z = 0.15$; (b) Influence of the fluid-added damping weighting parameters $P_\alpha$ and $P_z$ (three values are selected, see the legend) on the system behavior, with $\varepsilon = 0.7$.</figcaption></figure>

$A_3 = 25$. Finally, the fluctuating lift coefficient is adjusted to fit maximal vibration amplitudes without impacting too strongly the lock-in ranges of the two modes. As shown in Fig. 9(b), the best fitted value is $C_L^0 = 0.0097$. All previously fitted parameters are summarized in Tab. 3. One should note that the pressure center of the unsteady lift force due to the vortex shedding, initially placed at the trailing edge of the hydrofoil $B_{viv} = B$, does not require modification in this fitting.

<table><caption>Table 3 – Best fitted parameters of the coupled VIV model for the second and third modes of the hydrofoil.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Unit</th></tr></thead><tbody><tr><td>Van der Pol parameter</td><td>ε</td><td>0.7</td><td>—</td></tr><tr><td>Static lifting coefficient</td><td>C<sub>L</sub><sup>0</sup></td><td>0.0097</td><td>—</td></tr><tr><td>Second mode fluid-added damping weighting coefficient</td><td>P<sub>α</sub></td><td>0.5</td><td>—</td></tr><tr><td>Third mode fluid-added damping weighting coefficient</td><td>P<sub>z</sub></td><td>0.075</td><td>—</td></tr><tr><td>Second mode acceleration coupling parameter</td><td>A<sub>2</sub></td><td>20</td><td>—</td></tr><tr><td>Third mode acceleration coupling parameter</td><td>A<sub>3</sub></td><td>25</td><td>—</td></tr><tr><td>Pressure center distance from the leading edge</td><td>B<sub>viv</sub></td><td>85</td><td>mm</td></tr></tbody></table>

Numerical results show that the proposed VIV model, using the parameters reported in Tab. 3, accurately reproduces the experimental data in both amplitude and frequency within the lock-in range corresponding to the second bending mode of the hydrofoil. Regarding the first twisting mode, the model fails to fit the experimental data for $\bar{\omega}_v \in [0.8, 0.95]$ (i.e. $U \in [3.4, 4]$ m.s$^{-1}$). It can be explained by the fact that we considered a single fluid excitation on the structure while a second one is expected to reach the first twisting mode in this range of flow velocity, according to the associated Strouhal number, $St_1 = 0.49$, reported in [12]. Moreover, the turbulence-induced vibration model is not considered in the current fluid excitation. Yet, it impacts the structural dynamics and could contribute to extending the lock-in range. Also, one should note that there is no uniqueness of the solution. By changing some of the fitted parameters and modifying the others adequately, similar results shown in the last plot of Fig. 9(b) can be obtained. The authors tried
<figure><img src="image_10.png" /><figcaption>Figure 9 – Amplitude and frequency responses of the coupled fluid-structure model in Eqs. (27a, 27b, 27c) using $k_2 = k_3 = 0$, in comparison with experimental data for the short-circuited shunt configuration, with $\epsilon = 0.7$, $B_{viv} = B$, $P_{\alpha} = 0.5$, $P_z = 0.075$; (a) Influence of the coupling parameters scaling the structure acceleration to the wake dynamic $A_2$ and $A_3$ (three values are selected, see the legend) on the system behavior, with $C_L^0 = 0.008$; (b) Influence of the static lifting coefficient $C_L^0$ (three values are selected, see the legend) on the system behavior, with $A_2 = 20$, $A_3 = 25$.</figcaption></figure>

to keep parameters in the same order of magnitude as in the cited VIV studies. Still, the proposed VIV model is validated in the lock-in range area where vibration amplitudes are highest and where the effect of piezoelectric control will be scrutinized.

### 4.4. Confrontation of the piezo-elasto-hydrodynamic model to experimental data

Thanks to the VIV model fitted to experimental data, we now consider the piezo-elasto-hydrodynamic system of Eqs. (27a, 27b, 27c, 27d) with the piezoelectric shunt activated. As mentioned previously, the electromechanical coupling factors were experimentally measured by Watin et al. [12]: $k_2 = 8.5\%$ and $k_3 = 4.7\%$, for the second and third modes respectively.

To the best knowledge of the authors, no optimal electrical parameters are available in the literature for the piezo-elasto-hydrodynamic coupled system under lock-in conditions (i.e. for non-linear system). Despite that, optimal values of the electrical frequency $\omega_e^{\text{lin}}$ and the electrical damping ratio $\zeta_e^{\text{lin}}$ which maximize vibration amplitude reduction are given for reducing the resonance of a linear system subjected to external forcing in [19]:

$$ \omega_e^{\text{lin}} \approx \omega_j \sqrt{1 + k_j^2} \quad \text{and} \quad \zeta_e^{\text{lin}} \approx \sqrt{\frac{3}{8}} k_j, \qquad (33) $$

where, in our case, $j = 2$ since the piezoelectric patch is initially configured with respect to the second natural mode of the hydrofoil. In addition, experimental values of the electrical inductance $L$ and resistance $R$ used for vibration control of the hydrofoil in [12] are not known precisely, and the aim is, therefore, to find an equivalent level of vibration mitigation. These values of the electrical inductance and resistance should be measured in future work. Thus, as initial electrical parameters, we use the ones in Eq. (33) since they are optimal for weakly coupled VIV, for which the external excitation can be seen as a sinusoidal forcing. In Fig. 10(a) is shown the evolution of the vibration control with the electrical frequency parameter. By modifying its value, the anti-resonance created by the electrical inductance is moved in frequency and, in doing so, one can adjust the mitigation of the two modes as we can observe in Fig. 10(b).

Even if the values of the electrical components $L$ and $R$ were not precisely known to the authors, numerical results show that the piezo-elasto-hydrodynamic model in Eqs. (27a, 27b, 27c, 27d) is able to
<figure><img src="image_11.png" /><figcaption>Figure 10 – Amplitude and frequency responses of the coupled piezo-elasto-hydrodynamic model in Eqs. (27a, 27b, 27c, 27d), in comparison with experimental data for the shunt activated, with fitted parameters in Tab. 3; (a) Influence of the electrical frequency $\omega_e$ (three values are selected, see the legend) on the system behavior, with $\zeta_e = \zeta_e^{\text{lin}}$; (b) Influence of the electrical damping ratio $\zeta_e$ (three values are selected, see the legend) on the system behavior, with $\omega_e = 1.08 \times \omega_e^{\text{lin}}$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11 – Amplitude response of the coupled piezo-elasto-hydrodynamic model, with fitted parameters in Tab. 3 and $\omega_e = 1.08 \times \omega_e^{\text{lin}}$ and $\zeta_e = 2 \times \zeta_e^{\text{lin}}$, comparing experimental data for the hydrofoil without and with piezoelectric shunt; Numerical results obtained by the continuation method where solid lines (—) represent stable solutions and dotted lines (····) represent unstable one.</figcaption></figure>

reproduce the vibration mitigation of the resonant piezoelectric shunt used in [12]. Fig. 11 summarizes the results obtained with the VIV model and the piezo-elasto-hydrodynamic one, in comparison with the experimental data. Aside from the limitation pointed out previously for $\bar{\omega}_v \in [0.8, 0.95]$, the model allows a satisfactory computation of the overall lock-in range and the maximum vibration amplitude without and with control around two structural modes. One can, however, notice that the transition from the lock-in
with the twisting mode to the lock-in with the bending mode occurs sooner (i.e for lower $\bar{\omega}_v$) for the VIV model without control.

## 5. Conclusion

A fluid-structure-piezoelectric reduced order model of a hydrofoil has been proposed in order to analyze mitigation of vortex-induced vibrations. The first twisting mode and the second bending mode of the cantilever hydrofoil have been scrutinized in a range of Reynolds number $2.55 \times 10^5 \le Re \le 4.68 \times 10^5$ where strong vibrations, including lock-in, occur with the considered structures. The proposed hydro-elastic VIV model has been expressed with a set of seven unknown parameters. Considering a single fluid excitation mechanism (i.e. the Kármán vortex shedding at the trailing edge), parameters have been identified by confrontation with experimental data to best fit the amplitude and frequency response at both lock-in with the structural modes considered. The model could be improved in order to better capture (1) the impact of a secondary fluid excitation mechanism (probably due to an impinging shear layer instability mechanism from the leading edge) and, (2) the transition between the lock-in with the twisting mode and the lock-in with the bending mode. Values of both the electrical frequency $\omega_e$ and damping ratio $\zeta_e$, leading to the best mitigation results in the water tunnel, could not be precisely identified during the experiments (the equivalent electrical inductance $L$ and resistance $R$ used in the experiments were unknown). Then, a parametric study has been done on both these electrical parameters in order to fit the experimental data. Using $\omega_e = 1.08 \times \omega_e^{\text{lin}}$ and $\zeta_e = 2 \times \zeta_e^{\text{lin}}$, where $\omega_e^{\text{lin}}$ and $\zeta_e^{\text{lin}}$ are the optimal parameters for a linear system subjected to a harmonic forcing, the model shows mitigation results, in both twisting and bending, in accordance with the experimental data.

Future work includes an analytic or parametric optimization of the electrical inductance $L$ and resistance $R$ (i.e. the electrical frequency $\omega_e$ and damping ratio $\zeta_e$) for strong VIV coupling, in order to maximize the vibration mitigation performed by the shunt. Adaptive piezoelectric shunts will also be investigated in order to keep the value of the inductance $L$ optimized while the natural frequency of the structure or the electrical capacitance is changing. New experiments in the water tunnel on a composite hydrofoil equipped with a piezoelectric shunt are also planned which will be an opportunity to experimentally measure some parameters initially fitted in this study (such as the fluid-added damping ratios $C_\alpha$, $C_{3,z}$ and the coupling parameters $\bar{\Psi}_2, \bar{\Psi}_3$).

## Acknowledgments

The authors would like to express sincere gratitude to the members of the Institut de Recherche et d'Études Navales (IRENav), more particularly to Yann Watine and Céline Gabillet, for their assistance and data sharing. The authors also acknowledge the Agence de l'Innovation de Défense (AID) and the Agence Nationale de la Recherche (ANR) for supporting the HYDRAVIB project (ANR-22-ASTR-0012).

## Appendix A. Fluid-added damping approximation

Regarding the fluid-added damping ratios, we consider the model of Amandolese et al. [7], with the assessment that the mode shapes are mainly bending or twisting deformations. This leads to the following expression:

$$ \Phi_{j,z}(x, y) = \Phi_{j,z}(y), \quad \Phi_{j,\alpha}(x, y) = \bar{\Phi}_{j,\alpha}(y) \left(x - \frac{B}{2}\right), \qquad (A.1) $$

where $\Phi_{j,z}$ and $\Phi_{j,\alpha}$ are the $j$-th bending mode and twisting mode shapes of the structure, respectively, with the elastic center along the span approximately located at $x = B/2$. Moreover, $\bar{\Phi}_{j,\alpha}(y)$ is the dimensionless twisting mode shape along the span. Therefore, for a bending mode, the modal mass in Eq. (4) can be rewritten as:

$$ m_{j,z}^{s} = M_{z}^{s} \int_{0}^{S} \Phi_{j,z}^{2}(y) \, dy \quad \text{with} \quad M_{z}^{s} = \int_{0}^{B} \rho h(x) \, dx, \qquad (A.2) $$

where $M_z^s$ [kg.m$^{-1}$] is the mass per unit length span-wise. The total modal mass per unit length is therefore $M_j^{\text{tot}} = M_z^s + M_{j,z}^f$ with $M_{j,z}^f = M_z^s [(\omega_j^s/\omega_j)^2 - 1]$ the fluid-added modal mass per unit length.

For a twisting mode, the modal mass becomes:

$$ m_{j,\alpha} = \frac{2}{B} I_{\alpha}^{s} \int_{0}^{S} \bar{\Phi}_{j,\alpha}^{2}(y) \, dy \quad \text{with} \quad I_{\alpha}^{s} = \int_{0}^{B} \rho h(x) \left(x - \frac{B}{2}\right)^{2} \, dx, \qquad (A.3) $$

where $I_\alpha^s$ [kg.m] is the moment of inertia per unit length around the elastic center. The total moment of inertia per unit length is therefore $I_j^{\text{tot}} = I_\alpha^s + I_{j,\alpha}^f$ with $I_{j,\alpha}^f = I_\alpha^s [(\omega_j^s/\omega_j)^2 - 1]$. Note that the fluid-added
modal mass per unit length can also be approximated using theoretical expressions for flat plate of chord B in an inviscid fluid [45]: $M_z^f \approx \rho \pi B^2/4$ and $I_\alpha^f \approx \rho \pi B^4/128$.

Therefore, derived from Theodorsen's analytic solution [5] and following the model proposed by Amandolese et al. [7], the velocity-dependent fluid-added damping ratios relative to bending or twisting modes, can be expressed as:

$$ \zeta_{j,z}^f = P_z \frac{\rho_f U B}{4 \omega_j M_j^{\text{tot}}} \frac{dC_L}{d\alpha} F(\kappa) \quad (A.4a) $$

$$ \zeta_{j,\alpha}^f = P_{\alpha} \frac{\rho_f U B^3}{4 \omega_j I_j^{\text{tot}}} \left[ \frac{\pi}{8} - \left. \frac{dC_L}{d\alpha} \right|_{\alpha_0} \left( \frac{F(\kappa)}{16} + \frac{G(\kappa)}{4\kappa} \right) \right], \quad (A.4b) $$

where $dC_L/d\alpha \approx 2\pi$ is the slope of the lift coefficient, $F(\kappa)$ and $G(\kappa)$ are respectively the real and imaginary parts of the Theodorsen's function [5] with $\kappa = B\omega/U$ the reduced frequency (and $\omega$ the system frequency). Thanks to the short range of the fluid velocity for which the hydrofoil is studied [12], Theodorsen's functions vary only slightly, and the reduced frequency $\kappa$ is far greater than one. Then, approximations can be made:

$$ \lim_{\kappa \gg 1} F(\kappa) \approx \frac{1}{2}, \quad \lim_{\kappa \gg 1} \frac{F(\kappa)}{16} + \frac{G(\kappa)}{4\kappa} \approx \frac{1}{32}. \quad (A.5) $$

Thus, the fluid-added damping ratios in Eqs. (A.4b, A.4a) can be rewritten as in Eq. (8).

## References

[1] W. K. Blake, Excitation of plates and hydrofoils by trailing edge flows, Journal of Vibration, Acoustics, Stress, and Reliability in Design 106 (3) (1984) 351-363. doi:10.1115/1.3269201.

[2] G. H. Toebes, P. S. Eagleson, Hydroelastic vibrations of flat plates related to trailing edge geometry, Journal of Basic Engineering 83 (4) (1961) 671-678. doi:10.1115/1.3662292.

[3] P. Francois, J. A. Astolfi, X. Amandolese, Experimental analysis of hydrofoil hydroelastic trailing edge vibrations, Journal of Fluids and Structures 125 (2024) 104078. doi:https://doi.org/10.1016/j.jfluidstructs.2024.104078.

[4] M. P. Paidoussis, S. J. Price, E. de Langre, Fluid-Structure Interactions: Cross-Flow-Induced Instabilities, Cambridge University Press, 2011. doi:10.1017/CBO9780511760792.

[5] T. Theodorsen, General theory of aerodynamic instability and the mechanism of flutter, Tech. rep., Aerodynamic Flutter, American Institute of Aeronautics and Astronautics (1935).

[6] M. L. Facchinetti, E. De Langre, F. Biolley, Coupling of structure and wake oscillators in vortex-induced vibrations, Journal of Fluids and Structures 19 (2) (2004) 123-140. doi:10.1016/j.jfluidstructs.2003.12.004.

[7] X. Amandolese, C. Gabillet, Y. Watine, Adaptation of the Facchinetti, De Langre and Biolley model for the hydroelastic vortex-induced vibrations of a cantilevered flat plate, Journal of Fluids and Structures 124 (2024) 104017. doi:10.1016/j.jfluidstructs.2023.104017.

[8] S. Krenk, Frequency analysis of the tuned mass damper, Journal of Applied Mechanics-transactions of The ASME 72 (2005) 936-942. doi:10.1115/1.2062867.

[9] A. M. S. Baz, Active and Passive Vibration Damping, John Wiley & Sons, Ltd, 2019. doi:https://doi.org/10.1002/ 9781118537619.ch1.

[10] M. Auleley, O. Thomas, C. Giraud-Audine, H. Mahé, Enhancement of a dynamic vibration absorber by means of an electromagnetic shunt, Journal of Intelligent Material Systems and Structures 32 (3) (2021) 331-354. doi:10.1177/ 1045389X20957097.

[11] K. Marakakis, G. K. Tairidis, P. Koutsianitis, G. E. Stavroulakis, Shunt piezoelectric systems for noise and vibration control: A review, Frontiers in Built Environment 5 (2019) 1-17. doi:10.3389/fbuil.2019.00064.

[12] Y. Watine, B. Lossouarn, C. Gabillet, J. A. Astolfi, J. F. Deü, Mitigation of hydrofoil torsional flow induced vibrations by resonant piezoelectric shunt, Ocean Engineering 313 (2024) 119598. doi:https://doi.org/10.1016/j.oceaneng.2024. 119598.

[13] L. Pernod, B. Lossouarn, J. A. Astolfi, J. F. Deü, Vibration damping of marine lifting surfaces with resonant piezoelectric shunts, Journal of Sound and Vibration 496 (2021) 115921. doi:10.1016/j.jsv.2020.115921.

[14] Z. A. Shami, C. Giraud-Audine, O. Thomas, A nonlinear piezoelectric shunt absorber with a 2:1 internal resonance: Theory, Mechanical Systems and Signal Processing 170 (2022) 108768. doi:https://doi.org/10.1016/j.ymssp.2021.108768.

[15] K. K. Moon, D. H. Yang, Dynamic modelling and active vibration control of a submerged rectangular plate equipped with piezoelectric sensors and actuators, Journal of Fluids and Structures 54 (2015) 848-867. doi:https://doi.org/10.1016/j.jfluidstructs.2015.02.001.

[16] A. Presas, Y. Luo, Z. Wang, D. Valentin, M. Egusquiza, A review of PZT patches applications in submerged systems, Sensors 18 (7) (2018). doi:10.3390/s18072251.

[17] L. Sheng, Active modal control simulation of vibro-acoustic response of a fluid-loaded plate, Journal of Sound and Vibration 330 (23) (2011) 5545-5557. doi:https://doi.org/10.1016/j.jsv.2011.07.001.

[18] J. Ducarne, O. Thomas, J. F. Deü, Placement and dimension optimization of shunted piezoelectric patches for vibration reduction, Journal of Sound and Vibration 331 (14) (2012) 3286-3303. doi:10.1016/j.jsv.2012.03.002.

[19] O. Thomas, J. Ducarne, J. F. Deü, Performance of piezoelectric shunts for vibration reduction, Smart Materials and Structures 21 (1) (2012) 015008. doi:10.1088/0964-1726/21/1/015008.

[20] B. Lossouarn, M. Aucejo, J. F. Deü, B. Multon, Design of inductors with high inductance values for resonant piezoelectric damping, Sensors and Actuators A: Physical 259 (2017) 68-76. doi:10.1016/j.sna.2017.03.030.

[21] N. W. Hagood, A. von Flotow, Damping of structural vibrations with piezoelectric materials and passive electrical networks, Journal of Sound and Vibration 146 (2) (1991) 243-268. doi:https://doi.org/10.1016/0022-460X(91)90762-9.

[22] G. J. Efstathiades, A new approach to the large-deflection vibrations of imperfect circular disks using Galerkin's procedure, Journal of Sound and Vibration 16 (2) (1971) 231-253. doi:https://doi.org/10.1016/0022-460X(71)90485-8.
[23] O. Thomas, S. Bilbao, Geometrically nonlinear flexural vibrations of plates: In-plane boundary conditions and some symmetry properties, Journal of Sound and Vibration 315 (3) (2008) 569-590. doi:https://doi.org/10.1016/j.jsv. 2008.04.014.

[24] V. Denis, A. Pelat, C. Touzé, F. Gautier, Improvement of the acoustic black hole effect by using energy transfer due to geometric nonlinearity, International Journal of Non-Linear Mechanics 94 (2017) 134-145. doi:https://doi.org/10. 1016/j.ijnonlinmec.2016.11.012.

[25] L. Meirovitch, Analytical methods in vibrations, Macmillian Publishing Co., 1967.

[26] M. Géradin, D. Rixen, Mechanical Vibrations: Theory and Applications to Structural Dynamics, 3rd Edition, J. Wiley & Sons, 2015.

[27] T. Sarpkaya, Vortex-induced oscillations: A selective review, Journal of Applied Mechanics 46 (2) (1979) 241-258. doi: 10.1115/1.3424537.

[28] R. T. Hartlen, I. G. Currie, Lift-oscillator model of vortex-induced vibration, Journal of the Engineering Mechanics Division. Proceedings of the American Society of Civil Engineers 96 (5) (1970) 577-591.

[29] O. Thomas, J. F. Deü, J. Ducarne, Vibrations of an elastic structure with shunted piezoelectric patches: efficient finite element formulation and electromechanical coupling coefficients, International Journal for Numerical Methods in Engineering 80 (2) (2009) 235-268. doi:10.1002/nme.2632.

[30] A. Lazarus, O. Thomas, J. F. Deü, Finite element reduced order models for nonlinear vibrations of piezoelectric layered beams with applications to NEMS, Finite Elements in Analysis and Design 49 (1) (2012) 35-51. doi:10.1016/j.finel. 2011.08.019.

[31] M. Berardengo, O. Thomas, C. Giraud-Audine, S. Manzoni, Improved resistive shunt by means of negative capacitance: new circuit, performances and multi-mode control, Smart Materials and Structures 25 (7) (2016) 075033. doi:10.1088/ 0964-1726/25/7/075033.

[32] E. Isaacson, H. Keller, Analysis of numerical methods, John Wiley and Sons, New York, 1994.

[33] L. Guillot, B. Cochelin, C. Vergez, A Taylor series-based continuation method for solutions of dynamical systems, Nonlinear Dynamics 98 (2019) 2827-2845. doi:https://doi.org/10.1007/s11071-019-04989-5.

[34] L. Guillot, B. Cochelin, C. Vergez, A generic and efficient Taylor series-based continuation method using a quadratic recast of smooth nonlinear systems, International Journal for Numerical Methods in Engineering 119 (4) (2019) 261-280. doi:https://doi.org/10.1002/nme.6049.

[35] O. Thomas, B. Cochelin, Numerical continuation in nonlinear dynamics with the asymptotic numerical method, in: O. Gendelman, G. Rega (Eds.), Handbook on Nonlinear Dynamics, Vibrations and Acoustics. Volume 2: Numerical Methods, World Scientific, to appear in 2025.

[36] L. Guillot, A. Lazarus, O. Thomas, C. Vergez, B. Cochelin, A purely frequency based Floquet-Hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems, Journal of Computational Physics 416 (2020) 109477. doi:https://doi.org/10.1016/j.jcp.2020.109477.

[37] B. Cochelin, M. Medale, Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations, Journal of Computational Physics 236 (2013) 594-607. doi:https://doi.org/10. 1016/j.jcp.2012.11.016.

[38] O. Ellingsen, X. Amandolese, Amplitude hysteresis and the synchronization region: Prediction of vortex-induced vibration using a freely forced van der pol oscillator (03 2020). doi:10.48550/arXiv.2003.03838.

[39] A. Balanov, N. Janson, D. Postnov, O. Sosnovtseva, Synchronization: From Simple to Complex, 1st Edition, Springer Berlin, Heidelberg, 2009. doi:https://doi.org/10.1007/978-3-540-72128-4.

[40] F. Verhulst, Nonlinear Differential Equations and Dynamic Systems. 1996. doi:10.1007/978-3-642-61453-8.

[41] L. Liu, E. H. Dowell, K. C. Hall, A novel harmonic balance analysis for the Van der Pol oscillator, International Journal of Non-Linear Mechanics 42 (1) (2007) 2-12. doi:https://doi.org/10.1016/j.ijnonlinmec.2006.09.004.

[42] C. M. Andersen, J. F. Geer, Power series expansions for the frequency and period of the limit cycle of the van der pol equation, SIAM Journal on Applied Mathematics 42 (3) (1982) 678-693. doi:10.1137/0142047.

[43] A. Farshidianfar, H. Zanganeh, A modified wake oscillator model for vortex-induced vibration of circular cylinders for a wide range of mass-damping ratio, Journal of Fluids and Structures 26 (3) (2010) 430-441. doi:https://doi.org/10. 1016/j.jfluidstructs.2009.11.005.

[44] COMSOL AB, COMSOL multiphysics® v. 6.2., stockholm, sweden, www.comsol.com (2024).

[45] R. D. Blevins, Flow-induced vibration, 2nd Edition, Van Nostrand Reinhold, New York, 1990.