Nonlinear Dynamics of an Internally Resonant Base-Isolated Beam under Turbulent Wind Flow

: A base isolation system, aimed to passively control the nonlinear dynamics of an internally resonant tower, exposed to turbulent wind ﬂow, is studied. A continuous visco-elastic beam, constrained at the bottom end by a nonlinear visco-elastic device and free at the top end, is considered. All the nonlinearities, structural, inertial and aeroelastic, these latter computed via the quasi-static theory, are accounted in the model. The interaction between self- and parametric excitations, trig-gered by the mean wind velocity and the turbulent component, respectively, are analyzed. The Multiple Scale Method is applied to the partial differential equations of motion, to investigate critical and post-critical behaviors, when two modes in internal 1:3 resonance are involved in the response. The ﬁrst mode is found to lead the phenomenon, while the second mode is marginally involved. The effectiveness of the visco-elastic nonlinear isolation system is assessed, both in increasing the mean wind bifurcation value and in reducing the limit-cycle amplitude. The contribution of structural nonlinearities is found to weakly affect the response.

To passively suppress the aeroelastic instabilities, different kinds of added control devices have been proposed. In [27] a new aeroelastic device for the wind turbines was studied, yielding an attenuation of loads. In [28], a single-loop control system was applied to long-span suspension bridges, by controlling the leading-and trailing-edge flaps by sensing the main deck pitch angle. In [29], a passive control tool to increase the supersonic flutter boundary of composite panels, using a multimodal shunted piezoceramic in parallel topology, was identified. A wide use of the nonlinear energy sinks (NES), applied to aeroelastic problems, was also made in the literature (see, e.g., [30][31][32][33][34][35]). Here, the passive devices, consisting of essentially nonlinear oscillators, were designed to absorb energy from the main structure. They are able to shift forward the bifurcation point and, moreover, to reduce the amplitude of the nonlinear vibrations.
Referring to tall buildings, a visco-elastic device, made of an elastic spring in parallel with a nonlinear dashpot, was applied to the base of a tower, to mitigate the effects of a

Nonlinear Aeroelastic Model
The aeroelastic behavior of a base-isolated tall building, with a square cross-section and subjected to a turbulent wind flow, is investigated. The system is similar to the one proposed in [37], but here the structural model, besides the aerodynamic one, is assumed to be nonlinear as well. The building is modeled by a planar visco-elastic Euler-Bernoulli beam (Figure 1), obeying the Kelvin-Voigt constitutive law. The beam is free at the top, and constrained at ground by a visco-elastic device, consisting of a rheological model made of a linear elastic spring of stiffness κ in-parallel with a nonlinear dashpot. This last is characterized by a linear viscosity coefficient ζ 1 and a nonlinear (van der Pol-like) viscosity coefficient ζ 3 . The wind velocity U(t) is assumed to be uniform in space and equal to: in whichŪ is the (leading) steady-state wind velocity, andû ( Ū ) the amplitude and Ω the fundamental frequency of the turbulent component (the higher harmonics being neglected). The flow induces aerodynamic loads p a (s, t) nonlinearly depending on the structural velocityv(s, t), described by the quasi-steady theory [44]. The relevant equations of motion and boundary conditions are derived in Appendix A, extending those obtained in [45] to the case of Kelvin-Voigt material and base isolation. Introduction of the following positions is done: where v(s, t) is the transverse displacement at the abscissa s ∈ [0, l] and time t ∈ [0, ∞), with l the length of the beam; EI is the flexural stiffness of the beam; m is the mass per unit length; η is the internal viscous damping coefficient; c e is the external damping coefficient, accounting for dissipation of the beam in motionless air; b i := 1 2 a DA i (for i = 1, 3) where a is the air mass density and A i are dimensionless aerodynamic coefficients, depending on the cross-section shape, of which D is a characteristic length; ω r is a reference frequency; prime and dot denote differentiation with respect to the non-dimensional abscissa and time, respectively.
Accounting for Equation (2) and omitting the asterisks, the nondimensional form of the partial differential equation is found to be: with boundary conditions:

Modal Proprieties
First, free dynamics of the undamped and unloaded beam, constrained at ground by the elastic spring κ, are addressed, to evaluate the natural frequencies and modes. To this aim, by letting v(s, t) = ϕ(s)e iωt and substituting it in Equations (3) and (4), the following boundary value problem is derived, when retaining linear terms and neglecting structural and aerodynamic damping: with (ω, ϕ(s)) a real eigenpair. The solution reads: ϕ(s) = c 1 sin(αs) + c 2 cos(αs) + c 3 sinh(αs) + c 4 cosh(αs) (6) with α := √ ω a root of the transcendental characteristic equation: The system natural frequencies, therefore, depend on the stiffness κ of the base elastic spring. Aimed to investigate possible internal resonances, the ratio ω 2 /ω 1 between the second and the first root of Equation (7) is evaluated by varying the stiffness κ, see Figure 2. It is observed that, for κ = 6.5, a 1:3 internal resonance occurs.  The relevant frequencies ω 1 and ω 2 of the two first natural modes, and the coefficients of the associated modal shapes ϕ 1 (s) and ϕ 2 (s), normalized according to ϕ(1) = 1, are reported in Table 1, for both the uncontrolled (κ → ∞) and controlled (κ = 6.5) systems. The modal shapes are also shown in Figure 3. They show that the chosen spring stiffness allows quite large (but again acceptable) displacements at the base, according which the isolation system is effective (in the linear field). Therefore, this particular internally resonant case deserves to be analyzed.

Perturbation Analysis
The perturbation analysis is carried out here to evaluate the bifurcation equations and analyze in detail the post-critical solutions, into a small region of the parameter space. It is performed using the Multiple Scale Method in direct approach [46], i.e., directly applying the method to the PDE of motion (3) and (4), on the same line followed in [37], where a 1:2 resonance between the first mode of the beam and the turbulence excitation term is considered. However, here the procedure is extended to the concurrent 1:3 internal resonance condition between the first two modes of the beam, and taking into account the structural nonlinear contributions as well.
Introducing a small and positive bookkeeping parameter , the following scaling of variable is chosen: v = 1/2v (8) consistently with systems with cubic nonlinear terms; for the parameters, the scaling is: which brings linear damping and wind effect to the highest perturbation order here considered. Time scales are introduced, as t 0 = t, t 1 = t, and time-derivatives consistently become: where ∂ j = ∂ ∂t j , j = 0, 1. Moreover, series expansion of the dependent variable is applied: Resonance conditions read: where σ e , σ i are the external and internal detuning parameters, respectively. Substituting Equations (8)- (11) in Equations (3) and (4), dividing by 1/2 and collecting terms multiplying 0 and 1 , respectively, the following perturbation equations are obtained (omitting hat): and Equation (13) provides the generating solution: where, due to the presence of damping, only internally resonant modes are retained, with A 1 (t 1 ), A 2 (t 1 ) slow time varying complex amplitudes to be determined, c.c. standing for complex conjugate and i for the imaginary unit. Substitution of Equation (15) in Equation (14), with use of Equation (12), provides: where q j , Q j , C j , j = 1, 2, having the meaning of distributed loads, point-forces and couples, respectively, collect resonant terms and N.R.T. stands for non-resonant terms. For the sake of brevity, definition of q j , Q j , C j is omitted here.
Solvability conditions must be imposed in order to remove resonant terms at the right hand side of Equation (16). It assumes the form of the Virtual Power Equation, namely: Coming back to the time t and reabsorbing , Equation (17) produces the following Amplitude Modulation Equations (AME): where the coefficients are given in Appendix B. It is worth noting that the structural geometric and inertial nonlinear terms give rise to the coefficients indicated asd j , j = 2, 5, 7, 13, 16, 20, which contribute as imaginary parts of the relevant terms. The real form of the AME is obtained using the polar transformation, namely A 1 (t) = a 1 (t) 2 e iϑ 1 (t) and A 2 (t) = a 2 (t) 2 e iϑ 2 (t) , to be substituted in Equation (18). Then, separating real and imaginary parts and defining the phase differences ψ 1 (t) = 2ϑ 1 (t) − σ e t and ψ 2 (t) = ϑ 2 (t) − 3ϑ 1 (t) + σ i t, the real AME become: a 1 = d 10 +Ūd 10 a 1 +ûd 12 cos ψ 1 a 1 + 1 4û Equilibrium solutions are sought imposingȧ 1 =ȧ 2 =ψ 1 =ψ 2 = 0 in Equation (19), and their stability is evaluated through the sign of the real part of the eigenvalues of the Jacobian matrix of the same system. More specifically, to evaluate the stability of the trivial solution, the Cartesian form of Equation (18) is evaluated as well, after substituting A 1 (t) = (p 1 (t) + iq 1 (t)) exp(itσ e /2), A 2 (t) = (p 2 (t) + iq 2 (t)) exp(it(3σ e /2 − σ i )) and separating real and imaginary parts. Consequently, the Jacobian matrix evaluated at the trivial solution reads: − σ e 2 d 10 +Ūd 10 −ûd 12 0ûd 11 ud 21 0 d 22 and the changes in sign of the real parts of its eigenvalues are sought applying the Hurwitz criterion [47] to the characteristic polynomial det(J c − λI) = 0, where I is the identity matrix of dimension 4. It is worth noting that the number of bifurcation parameters to perform a complete unfolding of system (19) is four (e.g., (Ū,û, σ e , σ i )), which corresponds to its codimension. The expression of the beam displacement v(s, t) is reconstituted through Equation (15), which becomes: v(s, t) = a 1 (t) cos (21) where Equation (12) is used. From Equation (21), equilibra of Equation (19) correspond to periodic evolution of the displacement of the beam points, whereas periodic solutions in terms of amplitudes a 1 , a 2 give rise to quasi-periodic evolution of the beam.
Numerical calculations and parameter-continuation of the solutions of Equation (19), as well as bifurcation and stability analysis, are carried out within the software AUTO [48], combined with integration of the same equations in Mathematica [49].

Numerical Results
The complete unfolding of the codimension-four bifurcation is out of the scope of this paper. Here, however, the turbulence amplitude is kept fixed to the valueû = 0.1, whereas, for different values of the internal detuning σ i , the critical and post-critical behavior is analyzed on the plane (σ e ,Ū). More specifically, the values σ i = 0, i.e., perfect internal resonance, and σ i = −0.06, i.e., slightly far from the internal resonance, are alternatively considered, in order to assess the transition from a resonant to non-resonant case [50]. The other system parameters assume the following values: κ = 6.5, η = 0.00026, c e = 0.025, Different configurations of the viscous-elastic device are also addressed, namely the case of linear controlled beam (LC), i.e., ζ 3 = 0, to be compared to the outcomes of the nonlinear controlled system (NLC), where ζ 3 = 50. A further comparison with the uncontrolled (UC) case, realized when the base is fixed (κ → +∞), is carried out as well. In the above mentioned cases, geometric and inertia nonlinear contributions are not considered (d j set to zero), in order to be consistent with the model proposed in [37], where nonlinearities come from aerodynamic and viscous device only and where, however, exclusively the first mode was considered (NLC-1M). Finally, the effect of structural geometric and inertial nonlinear contributions is addressed (NLC-NL), through the modification produced by the coefficientsd j to the NLC case. The considered cases are summarized in Table 2. First, the stability domain of the trivial solution (a 1 , a 2 ) = (0, 0) is shown in Figure 4 for the uncontrolled and controlled cases, where only the first bifurcations are reported. The plots are in good agreement with those reported in [37] since, due to the nature of the considered damping, always the first mode is directly involved in the bifurcation, being the second one pulled along (as it will be confirmed later). For the same reason, the domains turn out to be insensitive to σ i , however, in dependence on σ e , the first bifurcation is of Hopf type (H, for |σ e | < 0.0022) or of Neimark type (N, for |σ e | ≥ 0.0022). For |σ e | < 0.0022, second Hopf bifurcations are found as well (SH, upper semi-perimeter of the ellipses). The control through the viscous-elastic device (black lines) generally induces an increasing of the critical values ofŪ, as expected, whereas the nonlinear terms in the equations do not affect the domains, they being ruled only by the linear ones.
In Figure 5, the bifurcation diagrams for amplitudes a 1 , a 2 as functions ofŪ are shown, obtained as vertical path in Figure 4 for σ e = 0. In particular, solid and dashed lines represent stable and unstable equilibria, respectively, whereas the filled regions indicate limit cycles in the amplitudes a 1 , a 2 , the former corresponding to (stable and unstable) periodic motions of the beam and the latter to quasi-periodic evolution. Comparison between controlled and uncontrolled (grey line, UC) cases shows that the control system pushes forward the first and second Hopf bifurcation points. Specifically, from the first point a stable branch arises, whereas from the second one the branch is unstable. Moreover, the NLC case (red line) shows smaller amplitudes in the post-critical evolution than the LC case (black line). If only one mode is considered (NLC-1M, blue line), a similar behavior to NLC is found for a 1 . Considering the nonlinear geometric and inertia terms (NLC-NL, orange line) modifies very slightly the amplitudes of the post-critical branches compared to the NLC case, except for their nature: periodic motions occur up toŪ = 1.27; then, for increasingŪ, quasi-periodic solutions are found, even if their modulation amplitudes are small. As a final comments, amplitude a 1 is one order of magnitude larger than a 2 and the bifurcated branch arises perpendicularly to the trivial one for a 1 , whereas it presents horizontal tangent for a 2 , as a confirmation of the pulling along of the second mode by the first one. The combined 3D bifurcation diagram, relevant to the case NCL-NL only, is shown in Figure 6, where it is evident the initial tangency of the bifurcated branch to the (Ū, a 2 ) plane, as it is found also in other internally resonant problems [51].
For the same values of parameters, and forŪ = 1.25, the phase plots of the reconstituted displacement (Equation (21)) of the beam, at cross-section s = 0.5 and s = 1, respectively, are shown in Figure 7. It is confirmed the reduction of the limit cycle amplitude due to the nonlinear control device. Moreover, slight differences appear between NLC-NL and: (i) NLC, where structural geometric and inertial nonlinearities are neglected; (ii) NLC-1M, where the second mode is ignored.    When σ e = 0.003 is considered, the bifurcation diagram for amplitudes a 1 , a 2 as functions ofŪ are shown in Figure 8. The major difference with respect to the perfect resonant case of Figure 5 is related to the different kind of bifurcation point, which is on Neimark type. Indeed, filled regions indicate that quasi-periodic solutions arise, however with modulation amplitudes which appear quite small. Furthermore, periodic motions are found for the UC and LC cases for larger values ofŪ, in some regions even coexisting with quasi-periodic evolution (LC). Quite far from internal resonance, Figure 9 shows the bifurcation diagram for amplitudes a 1 , a 2 as functions ofŪ at σ i = −0.06, and for σ e = 0. There, for better visibility, quasi-periodic solutions are omitted, whereas only periodic motion branches are shown. The behavior is qualitatively similar to the corresponding internally resonant case ( Figure 5), even if smaller amplitudes a 2 occur here, due to the weaker nonlinear interaction between the two modes. As final comments, the efficiency of the nonlinear control device is confirmed, both in pushing forward the bifurcation and in reduction of the post-critical amplitudes. The internal resonance for the analyzed cases appears as a weak phenomenon compared to the galloping induced on the first mode, slightly modifying the post-critical evolution with the presence of the minor component a 2 . Nonlinear structural geometry and inertia give minor contributions as well, confirming the almost linear features of cantilevers [52], even in the case of internal resonance.

Conclusions
The aeroelastic behavior of a tall building modeled as an Euler-Bernoulli beam is considered here. Base isolation with a nonlinear viscous-elastic device is assumed, for vibration mitigation purpose. The steady component of the wind triggers galloping, whereas the turbulent component, with harmonic evolution, is in principal parametric resonance with the first mode of the beam. Internal resonance in ratio 1:3 between the first and second beam mode, reached after suitable calibration of the base device, is considered as well.
After the application of the Multiple Scale Method in direct approach, a bifurcation analysis is carried out to evaluate the effects of the nonlinear viscous term of the device, the internal resonance and the structural nonlinear geometric and inertial terms.
The main conclusions can be summarized as follows: • The base device, through its linear viscous component, pushes forward the bifurcation point, as expected; • The nonlinear viscous term of the device reduces the amplitude of the post-critical evolution, as compared to the use of a purely linear device; • The galloping induced by the mean wind on the first mode represents the leading phenomenon, being the second mode marginally involved in the motion due to its smaller amplitudes; • The structural nonlinearities provide weak modification to the response, confirming the almost linear behavior of the cantilever, even when internal resonance occurs. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Nonlinear Equations of Motion
The nonlinear equations of motion of the cantilever, realized with viscous-elastic material and equipped with base isolation device, are evaluated here as an extension of what was done in [45] to the case under analysis. Being u(s, t), v(s, t) the longitudinal and transversal displacements of the axis-line points, respectively, the kinematic problem of the cantilever reads: where ε is the axial strain and k the bending curvature, Taylor-expanded up to the quadratic and cubic terms, respectively. Essential boundary conditions are written at cross-section A, namely: As typical for the case of longitudinally movable end (here at cross-section B), the internal constraint ε = 0 is applied and it is implemented introducing the Lagrange multiplier λ(s, t), which represents the reactive axial force along the beam. Moreover, a response function relevant to linear Kelvin-Voigt material is assumed, namely: being EI the bending stiffness and η the internal viscous factor. The extended Hamiltonian Principle, consequently, reads: where m is the mass per unit length, c e the external damping coefficient, p a the aerodynamic force, κ the base device stiffness, ζ 1 , ζ 3 the linear and cubic damping coefficients of the device, respectively, δ is the variation operator and t 1 , t 2 two generic time instants, at which δu, δv, δλ all vanish. After substitution of Equations (A1) and (A2) in (A6) and application of the part integration formula, the equations of motion are obtained. Among them, the one relevant to δu reads: which can be integrated to evaluate the expression of the Lagrange multiplier: where the boundary condition λ(1, t) = 0 is used. The equation relevant to δλ reproduces the constraint ε = 0. Then, using Equation (A1) for ε and after integration, the expression for u is obtained: where the boundary condition (A3) is used. Substitution of Equations (A8) and (A9) in the equation relevant to δv, provides the following condensed equation of motion, where only the variable v appears and terms up to the cubic order are retained: with boundary conditions:

Appendix B. Coefficients of Equations
Defining then the coefficients of the AME (18) are: