Detecting singular weak-dissipation limit for flutter onset in reversible systems

TL;DR

This paper introduces an experimental setup to study the transition between classical and reversible Hopf bifurcations in nonconservative systems, confirming the Whitney umbrella singularity's role in dissipation-induced instabilities.

Contribution

It provides the first laboratory demonstration of the singular interface between classical and reversible bifurcations in a viscoelastic system, validating theoretical predictions.

Findings

High sensitivity of bifurcation onset to damping ratios

Experimental confirmation of Whitney umbrella singularity

Discontinuities in dissipation-induced instabilities

Abstract

A `flutter machine' is introduced for the investigation of a singular interface between the classical and reversible Hopf bifurcations that is theoretically predicted to be generic in nonconservative reversible systems with vanishing dissipation. In particular, such a singular interface exists for the Pfl\"uger viscoelastic column moving in a resistive medium, which is proven by means of the perturbation theory of multiple eigenvalues with the Jordan block. The laboratory setup, consisting of a cantilevered viscoelastic rod loaded by a positional force with non-zero curl produced by dry friction, demonstrates high sensitivity of the classical Hopf bifurcation onset {to the ratio between} the weak air drag and Kelvin-Voigt damping in the Pfl\"uger column. Thus, the Whitney umbrella singularity is experimentally confirmed, responsible for discontinuities accompanying dissipation-induced instabilities in a broad range of physical contexts.

Equations

$$\frac{d\mathbf{x}}{dt}=\mathbf{g}(\mathbf{x}),\mathbf{x}\in {\mathbb{R}}^n$$

$$\det (\mathbf{A}-\lambda \mathbf{I})=\det (-\mathbf{R}\mathbf{A}\mathbf{R}-\mathbf{R}\lambda \mathbf{R})=(-1{)}^n\det (\mathbf{A}+\lambda \mathbf{I}),$$

$$\frac{d^2\mathbf{x}}{d{t}^2}=\mathbf{f}(\mathbf{x}),$$

$$\ddot{\mathbf{x}}(t)+(\delta \mathbf{D}+\Omega \mathbf{G})\dot{\mathbf{x}}(t)+(\mathbf{K}+\nu \mathbf{N})\mathbf{x}(t)=0,\mathbf{x}\in {\mathbb{R}}^2$$

$${\bf G}={\bf N}=\left(\begin{array}[]{rr}0&-1\\ 1&0\\ \end{array}\right),$$

$${\nu }_f=\sqrt{{\omega }_f^4-\det \mathbf{K}},{\omega }_f^2=\frac{\mathrm{tr}\mathbf{K}}{2},$$

$${\nu }_H(\Omega ,\delta )\approx {\nu }_f-\frac{2{\nu }_f}{(\mathrm{tr}\mathbf{D}{)}^2}{\left(\frac{\Omega }{\delta }-\frac{\mathrm{tr}(\mathbf{KD}-{\omega }_f^2\mathbf{D})}{2{\nu }_f}\right)}^2.$$

$${\Omega }_H\approx {\Omega }_f+\frac{2{\Omega }_f}{({\omega }_f\mathrm{tr}\mathbf{D}{)}^2}{\left(\frac{\nu }{\delta }-\frac{\mathrm{tr}(\mathbf{KD}+({\Omega }_f^2-{\omega }_f^2)\mathbf{D})}{2{\Omega }_f}\right)}^2,$$

$$\mathcal{M}(x,t)=-EJ{v}^{\prime \prime }(x,t)-E^{\ast }J{\dot{v}}^{\prime \prime }(x,t),$$

$$\xi =\frac{x}{l},\tau =\frac{t}{l^2}\sqrt{\frac{EJ}{m}},p=\frac{P{l}^2}{EJ},\alpha ={\tan }^{-1}\left(\frac{M}{ml}\right),$$

$$\eta =\frac{E^{\ast }{l}^2}{\sqrt{mEJ}}\frac{J}{l^4},\gamma =\frac{K{l}^2}{\sqrt{mEJ}},\beta =\frac{\gamma }{\eta },\chi =1-\bar{\chi },$$

$$v^{\prime \prime \prime \prime }(\xi ,\tau )+\eta {\dot{v}}^{\prime \prime \prime \prime }(\xi ,\tau )+p{v}^{\prime \prime }(\xi ,\tau )+\gamma \dot{v}(\xi ,\tau )+\ddot{v}(\xi ,\tau )=0,$$

$$(1+\eta \omega ){\tilde{v}}^{\prime \prime \prime \prime }+p{\tilde{v}}^{\prime \prime }+(\gamma \omega +{\omega }^2)\tilde{v}=0,$$

$$(1+\eta \omega ){\tilde{v}}^{\prime \prime \prime }(1)-(\chi -1){\tilde{v}}^{\prime }(1)p-{\omega }^2\tan (\alpha )\tilde{v}(1)=0,$$

$$\tilde{v}(0)={\tilde{v}}^{\prime }(0)=0,{\tilde{v}}^{\prime \prime }(1)=0.$$

$$\tilde{v}(\xi )$$

$${\lambda }_{1,2}^2=\frac{\sqrt{p^2-4(1+\eta \omega )(\gamma \omega +{\omega }^2)}\mp p}{2(1+\eta \omega )}$$

$$0={\lambda }_1{\lambda }_2(1+\eta \omega )({\lambda }_1^4+{\lambda }_2^4)+{\lambda }_1{\lambda }_{2}p(\chi -1)({\lambda }_2^2-{\lambda }_1^2)$$

$$+{\lambda }_1{\lambda }_2[2(1+\eta \omega ){\lambda }_1^2{\lambda }_2^2-p(\chi -1)({\lambda }_2^2-{\lambda }_1^2)]\cosh {\lambda }_1\cos {\lambda }_2$$

$$-{\omega }^2\tan \alpha ({\lambda }_1^2+{\lambda }_2^2)[{\lambda }_2\sinh {\lambda }_1\cos {\lambda }_2-{\lambda }_1\cosh {\lambda }_1\sin {\lambda }_2]$$

$$+{\lambda }_1^2{\lambda }_2^2[2p(\chi -1)+(1+\eta \omega )({\lambda }_2^2-{\lambda }_1^2)]\sinh {\lambda }_1\sin {\lambda }_2.$$

$$({\omega }^2[\mathbf{I}+4{\mathbf{M}}_1\tan \alpha ]+\omega [\gamma \mathbf{I}+\eta {\mathbf{D}}_i]+\left[{\mathbf{K}}_1-p{\mathbf{K}}_2+\chi p\mathbf{N}\right])\mathbf{a}=0,$$

$$(\mathbf{M}(\alpha ){\omega }^2+\mathbf{D}(\gamma ,\eta )\omega +\mathbf{A}(p,\chi ))\mathbf{a}=0,$$

$$\frac{\partial \det \mathbf{X}}{\partial p}=\mathrm{tr}\left({\mathbf{X}}^{\ast }\frac{\partial \mathbf{X}}{\partial p}\right).$$

$$q(\omega ,\alpha ,\chi ,p,\gamma ,\eta )=\det \mathbf{M}{\omega }^4+\mathrm{tr}({\mathbf{D}}^{\ast }\mathbf{M}){\omega }^3$$

$$+(\mathrm{tr}({\mathbf{A}}^{\ast }\mathbf{M})+\det \mathbf{D}){\omega }^2+\mathrm{tr}({\mathbf{A}}^{\ast }\mathbf{D})\omega +\det \mathbf{A}.$$

$${\sigma }_0^2=\frac{\mathrm{tr}({\mathbf{A}}_0^{\ast }{\mathbf{M}}_0)}{2\det {\mathbf{M}}_0}=\sqrt{\frac{\det {\mathbf{A}}_0}{\det {\mathbf{M}}_0}},$$

$${\mathbf{M}}_0=\mathbf{M}({\alpha }_0),{\mathbf{A}}_0=\mathbf{A}(p_0,{\chi }_0)$$

$$(\mathrm{tr}({\mathbf{A}}_0^{\ast }{\mathbf{M}}_0){)}^2=4\det {\mathbf{A}}_0\det {\mathbf{M}}_0.$$

$$q_0(\omega ,p(\epsilon ))=\sum _{r=0}^{2N}\frac{(\omega (\epsilon )-i{\sigma }_0{)}^r}{r!}\left(\frac{{\partial }^r{q}_0}{\partial {\omega }^r}+\epsilon \frac{{\partial }^r{q}_1}{\partial {\omega }^r}+o(\epsilon )\right),$$

$$\frac{{\partial }^r{q}_1}{\partial {\omega }^r}=\frac{{\partial }^{r+1}{q}_0}{\partial {\omega }^r\partial p}\frac{dp}{d\epsilon },$$

$$\omega (\epsilon )=i{\sigma }_0+{\epsilon }^{1/2}{\sigma }_1+\epsilon {\sigma }_2+\dots ,$$

$$q_0(i{\sigma }_0,p_0)=0,{{\sigma }_1\frac{\partial q_0}{\partial \omega }|}_{\omega =i{\sigma }_0,p=p_0}=0,$$

$${\left(q_1+\frac{1}{2}{\sigma }_1^2\frac{{\partial }^2{q}_0}{\partial {\omega }^2}+{\sigma }_2\frac{\partial q_0}{\partial \omega }\right)|}_{\omega =i{\sigma }_0,p=p_0}=0.$$

$${\sigma }_1^2=-q_1{\left(\frac{1}{2}\frac{{\partial }^2{q}_0}{\partial {\omega }^2}\right)}^{-1}=-{\left(\frac{1}{2}\frac{{\partial }^2{q}_0}{\partial {\omega }^2}\right)}^{-1}\frac{\partial q_0}{\partial p}\frac{dp}{d\epsilon }.$$

Full text

I Introduction

In a dissipative system oscillatory flutter instability, an example of a classical Hopf bifurcation, shifts a complex-conjugate pair of eigenvalues to the right in the complex plane. This instability mechanism is modified for a non-dissipative system possessing a reversible symmetry, defined with reference to the differential equation

[TABLE]

which is said to be $\bf R$-reversible (${\bf R}^{-1}={\bf R}$) if it is invariant with respect to the transformation $({\bf x},t)\mapsto({\bf R}{\bf x},-t)$, implying that the right hand side must satisfy ${\bf R}{\bf g}({\bf x})=-{\bf g}({\bf R}{\bf x})$.

If ${\bf x}={\bf x}_{0}$ is a reversible equilibrium such that ${\bf R}{\bf x}_{0}={\bf x}_{0}$, and ${\bf A}=\nabla{\bf g}$ is the linearization matrix about ${\bf x}_{0}$, then ${\bf A}=-{\bf R}{\bf A}{\bf R}$, and the characteristic polynomial

[TABLE]

implies that $\pm\lambda,\pm\overline{\lambda}$ are eigenvalues of ${\bf A}$ OR1996 ; W2003 ; Clerc1999 ; Clerc2001 . Due to the spectrum’s symmetry with respect to both the real and imaginary axes of the complex plane, the reversible-Hopf bifurcation requires the generation of a non-semi-simple double pair of imaginary eigenvalues and its subsequent separation into a complex quadruplet OR1996 ; W2003 ; Clerc1999 ; Clerc2001 .

All equations of second order

[TABLE]

are reversible OR1996 ; W2003 , including the case when the positional force ${\bf f}({\bf x})$ has a non-trivial curl, $\nabla\times{\bf f}({\bf x})\neq 0$, which makes the reversible system nonconservative.

Such nonconservative curl forces BS2012 appear in modern opto-mechanical applications, including optical tweezers G2003 ; Wu2009 ; SH2010 . In mechanics, they are known as circulatory forces for producing non-zero work along a closed circuit. Circulatory forces are common in the models of friction-induced vibrations HG2003 , rotordynamics Clerc2001 , biomechanics BD2016 and fluid-structure interactions MM2010 ; G2017 , to name a few. A circulatory force acting on an elastic structure and remaining directed along the tangent line to the structure at the point of its application during deformation is known as follower Z1952 ; B1963 ; KS1998 .

Since the dynamics of an elastic structure under a follower load is described by reversible equations OR1996 , flutter instability may occur via the reversible-Hopf bifurcation mechanism OR1996 ; Clerc2001 . In these conditions, Ziegler Z1952 discovered that, when viscosity is present, the location of the curve for the onset of the classical Hopf bifurcation is displaced by an order-one distance in the parameter space, with respect to the curve for the onset of the reversible-Hopf bifurcation in the elastic structure. This occurs even if the viscous damping in the structure is infinitesimally small Z1952 . Other velocity-dependent forces, such as air drag (or even gyroscopic forces), can also destabilize an elastic structure under a follower load OR1996 ; SW1975 ; T2016 ; K2013dg ; KV2010 . However, acting together, the velocity-dependent forces, e.g., the air drag and the material (Kelvin-Voigt) viscous damping, can inhibit the destabilizing effect of each other at a particular ratio of their magnitudes due to the singular interface between the classical Hopf and reversible-Hopf bifurcations B1956 ; Z1952 ; KV2010 ; T2016 .

For instance, the system

[TABLE]

where $\delta$, $\Omega$, $\nu$ are scalar coefficients and matrices ${\bf D}>0$, ${\bf K}>0$ are real and symmetric, while matrices G and N are skew-symmetric as follows

[TABLE]

is nonconservative and reversible for $\delta=\Omega=0$.

The reversible-Hopf bifurcation in the system (1) with $\delta=\Omega=0$ occurs at

[TABLE]

where “$\rm tr$” denotes the trace operator, which yields flutter instability when $\nu>\nu_{f}$. However, when $\delta>0$, $\Omega>0$ the classical-Hopf bifurcation occurs at a different value of $\nu$ K2013dg

[TABLE]

The expression for $\nu_{H}(\Omega,\delta)$ defines a surface in the $(\delta,\Omega,\nu)$-space that has a Whitney’s umbrella singular point at $(0,0,\nu_{f})$ L2003 ; footnote . Near that singular point, the neutral stability surface is a ruled surface, with a self-intersection degenerating at the singularity, so that a unique value of the ratio $\Omega/\delta$ is produced, for which the onsets of the classical and reversible Hopf bifurcations tend to coincide B1956 ; KV2010 ; footnote .

For a dissipative nearly-reversible system, the singular dependence of the classical Hopf bifurcation onset on the parameters of velocity-dependent forces has a general character B1956 , which follows from the codimension 3 (for dissipative systems) and 1 (for reversible vector fields) of non-semi-simple double imaginary eigenvalues B1956 ; A1972 ; KV2010 ; KH2010 ; K2013dg .

Since the singularity is related to a double imaginary eigenvalue arising from a Jordan block A1972 , it can be found in other dissipative systems that are close to undamped systems with the ‘reversible’ symmetry of spectrum KV2010 .

Indeed, the system (1) with $\delta=0$, $\Omega=0$, and $\nu=0$ is a conservative Hamiltonian system, which is statically unstable for ${\bf K}<0$. Adding gyroscopic forces with $\Omega>0$, keeps this system Hamiltonian and yields its stabilization if $\Omega>\Omega_{f}=\sqrt{-\kappa_{1}}+\sqrt{-\kappa_{2}}$, where $\kappa_{1,2}<0$ are eigenvalues of $\bf K$. Owing to the ‘reversible’ symmetry of its spectrum OR1996 ; MK1991 ; BKMR94 , the Hamiltonian system displays flutter instability via the collision of imaginary eigenvalues at $\Omega=\Omega_{f}$ and their subsequent splitting into a complex quadruplet as soon as $\Omega$ decreases below $\Omega_{f}$. This is the so-called linear Hamilton-Hopf bifurcation L2003 ; K2013dg ; K2007 .

If $\delta>0$, $\nu>0$ the gyroscopic stability is destroyed at the threshold of the classical-Hopf bifurcation K2013dg ; K2007

[TABLE]

where $\omega_{f}^{2}=\sqrt{\kappa_{1}\kappa_{2}}$ and ${\bf D}>0$. The dependency of the new gyroscopic stabilization threshold just on the ratio $\nu/\delta$ implies that the limit of $\Omega_{H}$ as both $\nu$ and $\delta\rightarrow 0$ is higher than $\Omega_{f}$ for all ratios except a unique one. Similarly to the case of nonconservative reversible systems, this happens because the classical Hopf and the Hamilton-Hopf bifurcations meet in the Whitney umbrella singularity that exists on the stability boundary of a nearly-Hamiltonian dissipative system and corresponds to the onset of the Hamilton-Hopf bifurcation L2003 ; KH2010 ; K2007 ; KV2010 ; K2013dg ; BKMR94 ; KM2007 .

The singular weak-dissipation limit for the flutter onset in nearly-Hamiltonian systems in the presence of two different damping mechanisms has been discovered first in the problem of secular instability of equilibria of rotating and self-gravitating masses of fluid, when dissipation due to both fluid viscosity TT1879 ; RS1963 ; BO2014 and emission of gravitational waves C1970 ; C1984 is taken into account LD1977 ; A2003 . Later on this phenomenon manifested itself as the ‘Holopäinen instability mechanism’ for a baroclinic flow H1961 ; WE2012 and as an enhancement of modulation instability with dissipation BD2007 . Analysis of this effect based on the method of normal forms and perturbation of multiple eigenvalues has been developed, among others by B1956 ; OR1996 ; AY1974 ; A1972 ; KH2010 ; KI2005 ; K2013dg ; K2007 ; BKMR94 ; Clerc2001 ; KI2017 ; LF2016 ; KS2005 ; MK1991 ; KV2010 ; FS1978 ; S1980 .

Although the destabilizing effect of damping for equilibria of Hamiltonian and reversible systems has been discussed for decades, no experimental evidence is known for the singular limit of the classical Hopf bifurcation in a nearly-Hamiltonian, or a nearly-reversible system, when the dissipation tends to zero. The main difficulty for such experiments is the accurate identification and control of at least two different damping mechanisms. For reversible elastic structures an additional challenge lies in the realization of circulatory follower loads, acting for a sufficiently long time. Previous attempts are reported to create a follower load through the thrust produced either by water flowing through a nozzle WSS1969 , or by a solid rocket motor mounted at the end of an elastic rod in a cantilever configuration SKK1995 ; SMR1995 ; SKK2000 ; LS2000 . In the former realization hydrodynamical effects enter into play and in the latter the duration of the experiments is limited to a few seconds. In contrast, the frictional follower force acting on a wheel mounted at the free end of the double-link Ziegler pendulum allowed Bigoni and Noselli to significantly relax the limitation on time BN2011 .

In the present article, an experimental realization is reported for the Pflüger column B1952 ; P1955 ; SKK1976 ; RS2003 , a viscoelastic cantilevered rod carrying a point mass at the free end and loaded with a follower force (Fig. 1) obtained via friction, similarly to BN2011 . Two dissipation mechanisms—the air resistance and the internal Kelvin-Voigt damping—are identified and controlled by changing the geometrical characteristics of the sample rods. The measured critical flutter loads demonstrate a high sensitivity to the ratio between the two damping coefficients, being almost insensitive to each of the damping coefficients that both are very close to zero, in agreement with both numerical modeling of T2016 and perturbation theory developed for the Pflüger column in the present work.

II Pflüger’s column and its Galerkin discretization

Consider a rod of length $l$, mass density per unit length $m$ and end mass $M$, its deflection $v$, function of the $x$ coordinate, obeys the Bernoulli law that the rotation of the cross-section $\phi$ is given by $\phi(x)=-v^{\prime}(x)$, where a prime denotes derivative with respect to $x$. A moment-curvature viscoelastic constitutive relation of the Kelvin-Voigt type is assumed in the form

[TABLE]

where a superimposed dot denotes the time derivative, $E$ and $E^{*}$ are respectively the elastic and the viscous moduli of the rod, which has a cross section with moment of inertia $J$. The rod is clamped at one end and is loaded through the force $P$ that is inclined with respect to the tangent to the rod at its free end such that $v^{\prime}(l)\bar{\chi}=const.$, Fig. 1.

Assuming that a distributed external damping $K$ caused by the air drag is acting on the rod, and introducing the dimensionless quantities

[TABLE]

the linearized partial differential equation of motion governing the dynamics of the rod can be written as

[TABLE]

where now a prime and a dot denote partial differentiation with respect to $\xi$ and $\tau$, respectively. Separating time in (II) with $v(\xi,\tau)=\tilde{v}(\xi)\exp(\omega\tau)$ yields a non-self-adjoint boundary eigenvalue problem T2016

[TABLE]

Assuming that $\tilde{v}(\xi)$ has the form

[TABLE]

with $A_{i}$ ($i=1,..,4$) arbitrary constants and

[TABLE]

and substituting (5) into (II) yields an algebraic system of equations which admits non-trivial solutions if T2016

[TABLE]

Results from experiments are compared with the eigenvalues, eigenfunctions and critical parameters of the boundary eigenvalue problem (II) which are directly found by numerical solution of the transcendental characteristic equation (II).

For theoretical purposes, the $N$-dimensional Galerkin discretization of the continuous problem (II) is also considered:

[TABLE]

where $\bf a$ is an $N$-vector and $\bf I$ is the $N\times N$ identity matrix. The entries of the $N\times N$ mass matrix ${\bf M}_{1}$ are $M_{1,ij}=(-1)^{i+j}$, the matrix of internal damping ${\bf D}_{i}$ is ${\bf D}_{i}={\rm diag}(\omega_{1}^{2},\omega_{2}^{2},\ldots,\omega_{N}^{2})$, and the stiffness matrix ${\bf K}_{1}$ is ${\bf K}_{1}={\rm diag}(\omega_{1}^{2},\omega_{2}^{2},\ldots,\omega_{N}^{2})$. The values of the frequencies $\omega_{1},\ldots,\omega_{N}$ as well as the entries of the symmetric stiffness matrix ${\bf K}_{2}$ and the non-symmetric matrix of circulatory forces $\bf N$ are given in the Appendix A.

III Theory of dissipation-induced flutter instability

For a Galerkin-discretized model of the Pflüger column (II) a perturbation theory is developed of the singular weak-dissipation limit for the onset of flutter.

III.1 The $N=2$ modes approximation and its stability analysis

The eigenvalue problem (II) has the form

[TABLE]

where ${\bf M}={\bf M}^{T}$, ${\bf D}={\bf D}^{T}$, ${\bf D}(0,0)=0$, and ${\bf A}\neq{\bf A}^{T}$, with the superscript $T$ denoting transposition.

Recall that the adjugate ${\bf X}^{*}$ of a $N\times N$ matrix $\bf X$ is defined as ${\bf X}^{*}={\bf X}^{-1}\det{\bf X}$ and, in particular,

[TABLE]

Since ${\rm tr}({\bf X}^{*}{\bf Y})={\rm tr}({\bf Y}^{*}{\bf X})$ for $N=2$, the characteristic polynomial of (9) in the case of $N=2$ can be written by means of the Leverrier algorithm in a compact form:

[TABLE]

Assuming that for $\eta=0$, $\gamma=0$, $\alpha=\alpha_{0}$, $\chi=\chi_{0}$, and $p=p_{0}$ the undamped system with $N=2$ degrees of freedom be on the flutter boundary, on this boundary its eigenvalues are imaginary and form a double complex-conjugate pair $\omega=\pm i\sigma_{0}$ of a Jordan block. In these conditions, the real critical frequency $\sigma_{0}$ at the onset of flutter follows from the characteristic polynomial in the closed form

[TABLE]

and the flutter boundary is described by the equation

[TABLE]

Since ${\bf M}_{0}={\bf I}+4{\bf M}_{1}\tan\alpha_{0}$ and ${\bf A}_{0}={\bf K}_{1}-p_{0}{\bf K}_{2}+\chi_{0}p_{0}{\bf N}$ is a linear function of $p_{0}$, equation (13) is quadratic with respect to $p_{0}$, which can thus be easily solved. The red solid curve in Fig. 2(a) shows the flutter boundary (13) of the undamped discretized model (II) of the Pflüger column with $N=2$ modes for $\chi_{0}=1$ in the $(\alpha_{0},p_{0})$-plane. The red solid curves in Fig. 2(b) demonstrate the movement of the eigenvalues of the undamped system at given $\chi=\chi_{0}=1$ and $\alpha=\alpha_{0}=0.1$ when the load parameter $0\leq p\leq 70$. The equilibrium is stable for $0\leq p<p_{0}$ where the critical flutter load is $p_{0}\approx 17.83368$, corresponding to a double pair of imaginary eigenvalues with the imaginary part $\sigma_{0}\approx 9.366049$ (see Eq. (III.1)). The value $p=p_{0}$ corresponds to the linear reversible-Hopf bifurcation, yielding the splitting of the double eigenvalues into a complex quadruplet causing flutter instability.

III.2 Reversible-Hopf bifurcation in the undamped model

A perturbation formula is now derived for the splitting of a double eigenvalue $\omega=i\sigma_{0}$, when $\gamma=\gamma_{0}$ and $\alpha=\alpha_{0}$ are fixed and $p$ is left to vary. Introducing a small parameter $0\leq\varepsilon\ll 1$ and assuming in the polynomial $q_{0}(\omega,p)=q(\omega,\alpha_{0},\chi_{0},p,\gamma=0,\eta=0)$ that $p(\varepsilon)=p_{0}+\varepsilon\frac{dp}{d\varepsilon}+\ldots$ (where the derivative is taken at $\varepsilon=0$) yields

[TABLE]

where the partial derivatives are evaluated at $p=p_{0}$ and $\omega=i\sigma_{0}$.

Assuming for the perturbed double non-semisimple eigenvalue the Newton-Puiseux series

[TABLE]

substituting equations (III.2) and (15) into the equation $q_{0}(\omega,p)$ and collecting the terms of the same powers of $\varepsilon$, leads to

[TABLE]

and

[TABLE]

Conditions (16) are satisfied for the double eigenvalue $\omega=i\sigma_{0}$, so that an account of this into (17) yields

[TABLE]

Hence, the splitting of the double non-semisimple eigenvalue due to the variation of $p$ is governed by the formula

[TABLE]

With the help of Eq. (10), Eq. (III.1), and the relations

[TABLE]

the following result is finally obtained

[TABLE]

For instance, for $\alpha_{0}=0.1$, $\chi_{0}=1$, $p_{0}\approx 17.83368$, and $\sigma_{0}\approx 9.366049$, the expression (19) becomes

[TABLE]

confirming the splitting of the double $i\sigma_{0}$ into two complex eigenvalues with opposite real parts (flutter) at $p>p_{0}$.

III.3 Dissipative perturbation of simple imaginary eigenvalues

At $p<p_{0}$ the eigenvalues of the undamped system $\omega=\omega(p)$ remain simple and imaginary. To investigate how they are affected by dissipation, it is assumed that $\eta(\varepsilon)=\frac{d\eta}{d\varepsilon}\varepsilon+o(\varepsilon)$, and $\gamma(\varepsilon)=\frac{d\gamma}{d\varepsilon}\varepsilon+o(\varepsilon)$ in the polynomial (III.1), where $\alpha=\alpha_{0}$, $\gamma=\gamma_{0}$ and $0\leq p<p_{0}$ are also fixed. Then, $\omega=\omega(p)+\frac{d\omega}{d\varepsilon}\varepsilon+o(\varepsilon)$, with

[TABLE]

The following approximation is therefore obtained

[TABLE]

where the partial derivatives are evaluated at $p<p_{0}$ and $\omega=\omega(p)$. An account of the following derivatives

[TABLE]

leads to

[TABLE]

III.4 Linear approximation to the stability boundary and the exact zero-dissipation limit of the critical flutter load

The correction, linear in $\eta$ and $\gamma$, to the simple imaginary eigenvalue in (22) due to damping is real and therefore it determines whether the dissipative perturbation is stabilizing or destabilizing. Equating this linear term to zero and taking into account that ${\bf A}={\bf K}_{1}-p({\bf K}_{2}-\chi_{0}{\bf N})$ and ${\bf D}_{i}={\bf K}_{1}={\rm diag}(\omega_{1}^{2},\omega_{2}^{2})$ yields the following approximation to the flutter boundary, which represents the onset of the classical Hopf bifurcation

[TABLE]

where ${\bf M}_{0}={\bf I}+4{\bf M}_{1}\tan\alpha_{0}$ and $\omega(p)$ is a root of the polynomial $q_{0}(\omega,p)$ in equation (III.2)1 at $p<p_{0}$. In the $(\eta,\gamma)$-plane the equation (III.4) defines a straight line, Fig. 3(a). In fact, at every $p<p_{0}$ there exist two lines (III.4) corresponding to two different eigenvalues $\omega(p)$ that participate in the reversible-Hopf bifurcation at $p=p_{0}$. However, as $p$ tends to $p_{0}$, the angle between the two lines decreases and completely vanishes in the limit $p\rightarrow p_{0}$, Fig. 3(a). This suggests that the approximation (III.4) defines a ruled surface in the $(\eta,\gamma,p)$-space. As a consequence, every fixed damping ratio $\beta=\gamma/\eta$ corresponds to a ruler at some $p<p_{0}$. Therefore, the condition for which the damping tends to zero at fixed damping ratio will occur along this ruler for the corresponding constant value of $p<p_{0}$ and will result in the limiting value of the critical flutter load that is lower than the critical load at the onset of the reversible-Hopf bifurcation, $p_{0}$, see Fig. 3(b). Note that equation (III.4) gives the exact dependency of the limit of the critical flutter load at vanishing dissipation as a function of the damping ratio, $\beta$, if the exact solution $\omega(p)$ of the polynomial $q_{0}(\omega,p)$ is used, see K2013dg ; KV2010 ; AY1974 ; K2007 ; LF2016 .

III.5 Quadratic approximation in $\beta$ to the exact zero-dissipation limit of the critical flutter load

In the vicinity of $p=p_{0}$ the two roots participating in the reversible-Hopf bifurcation are approximated by the expression (19). Using this expression in equation (III.4), the limit of zero dissipation can be found for the critical flutter load as a function of the damping ratio, $p(\beta)$, in the form of a series

[TABLE]

where

[TABLE]

From the quadratic approximation (24) it is evident that $p(\beta)\leq p_{0}$ for all $\beta$ except for the specific case of $\beta=\beta_{0}$, at which it exactly coincides with the critical flutter load of the undamped system: $p(\beta_{0})=p_{0}$. For instance, for $\alpha_{0}=0.1$ and $\chi_{0}=1$, the approximation (24) is

[TABLE]

as shown in Fig. 3(b) with a red solid curve.

III.6 The Whitney umbrella singularity

Truncating the series (24) and substituting $\beta=\gamma/\eta$ into the result, yields an expression for the ruled surface in the $(\eta,\gamma,p)$-space

[TABLE]

This expression is in the form $Z=X^{2}/Y^{2}$, which is the well-known normal form for the Whitney umbrella surface [12-15]. The surface (27) has a singular point at $p=p_{0}$, corresponding to the onset of the reversible-Hopf bifurcation, and a self-intersection at $p<p_{0}$.

In Fig. 4 the cross-sections are plotted in the $(\eta,\gamma)$-plane for different values of $p$ of the exact stability boundary calculated with the use of the Routh-Hurwitz criterion applied directly to the polynomial (III.1). Physically relevant is the first quadrant of the $(\eta,\gamma)$-plane.

For every $\alpha_{0}\in[0,\pi/2]$ the cross-sections look qualitatively similar. For $p>p_{0}$ the stability domain is bounded by a smooth curve departing from the origin, Fig. 4. For $p=p_{0}(\alpha_{0})$ the stability boundary has a cuspidal point at the origin with the single tangent line to the boundary specified by the ratio $\beta_{0}$ given by the equation (25); the stability region is inside the cusp. For $p<p_{0}(\alpha_{0})$ the stability boundary has a point of intersection at the origin in the $(\eta,\gamma)$-plane; the stability region is inside the narrow angle-shaped domain, which becomes wider as $p$ decreases and for $p=0$ spreads over the first quadrant of the plane for every possible mass distribution.

A comparison between Fig. 3(a) and Fig. 4(b) shows that equation (III.4) gives a correct linear approximation to the stability domain provided by the Routh-Hurwith criterion in the $(\eta,\gamma)$-plane and, therefore, to the singular interface between the classical-Hopf and reversible-Hopf bifurcations in the $(\eta,\gamma,p)$-space.

III.7 Stabilizing damping ratio $\beta_{0}$ for different mass distributions $\alpha_{0}$

Fig. 4 demonstrates that the contour plot patterns of the stability boundary in the $(\eta,\gamma)$-plane remain qualitatively the same for different values of $\alpha_{0}$, but differ in the orientation of the cusp, which is determined by the stabilizing damping ratio $\beta_{0}$. Evaluating (25) at the points of the stability boundary of the undamped system, provides the plot of the function $\beta_{0}(\alpha_{0})$ reported in Fig. 3(c). One can see that two intervals of $\alpha_{0}$ exist with opposite signs of $\beta_{0}$. The intervals are bounded by the values $\alpha_{0}=0$ and $\alpha_{0}\approx 0.342716$, at which the graph $\beta_{0}(\alpha_{0})$ displays a vertical asymptote, Fig. 3(c). Positive values of $\beta_{0}$ correspond to sufficiently small $\alpha_{0}\leq 0.342716$, cf. Fig. 4(b); negative values of $\beta_{0}$ are characteristic for $0.342716\leq\alpha_{0}\leq\pi/2$.

The above critical values of $\alpha_{0}$ are determined by the zeros of the denominator of equation (25). Indeed, taking into account that

[TABLE]

the denominator can be obtained in the form

[TABLE]

Evidently, one of the roots is $\alpha_{0}=0$, corresponding to the case of the Beck column (which is the Pflüger column without the end mass). In this case, the cusp in the $(\eta,\gamma)$-plane is oriented vertically, see Fig. 4(a). This confirms the well-known fact that for the Beck column the internal Kelvin-Voigt damping $(\eta)$ is destabilizing, and the external air drag damping $(\gamma)$ is stabilizing KS2005 ; B1963 ; T2016 . As soon as $\alpha_{0}$ departs from zero, the external damping becomes a destabilizing factor due to the change in the orientation of the cusp in Fig. 4. Nevertheless, at a specific mass distribution $\alpha_{0}\approx 0.342716$, which is given by the root of the equation

[TABLE]

the cusp restores its vertical orientation, as is visible in Fig. 4(c). For this specific mass ratio the external damping is stabilizing again.

The revealed behaviour of the stabilizing damping ratio as a function of the mass distribution is reflected in Fig. 2(a) that shows the red solid curve of the onset of the reversible-Hopf bifurcation in the undamped system together with the onset of the classical-Hopf bifurcation in the limit of vanishing (the green dash-dotted curve) internal damping and (the blue dashed curve) external damping. The latter curve has two common points with the stability boundary of the undamped system exactly at $\alpha_{0}=0$ and $\alpha_{0}\approx 0.342716$.

Remarkably, $\beta_{0}$ and its sign determine which mode will be destabilized by either of the two damping mechanisms or by their combination. For instance, in the case of $\beta_{0}>0$ the cusp of the stability boundary in the $(\eta,\gamma)$-plane is directed to the first quadrant, Fig. 4(b). Therefore, a dominating external damping will destabilize the mode with the higher frequency, whereas a dominating internal damping will destabilize the mode with the lower frequency, see Fig. 2(b). In the case of $\beta_{0}<0$ the cusp is oriented towards the second quadrant, Fig. 4(d), so that for every choice of internal and external damping with $\eta>0$ and $\gamma>0$, the mode with the lower frequency will be the destabilizing one.

Finally, using equation (III.4), the critical flutter load in the limit of vanishing dissipation is plotted in Fig. 5 as a function of the damping ratio $\beta$, for different mass ratios $\alpha\in[0,\pi/2]$. It is worth noting that in the range $0.7\leq\alpha\leq\pi/2$ the curves form a dense family. According to Fig. 3(c), for $0.342716\leq\alpha_{0}\leq\pi/2$ the stabilizing damping ratio $\beta_{0}$ is negative and tends to infinity as $\alpha_{0}\rightarrow+0.342716\ldots$, which corresponds to the vertically oriented cusp in Fig. 4(c).

III.8 Agreement with the solution of the boundary eigenvalue problem (II)

When $N$ is increased, the eigenvalues, eigenvectors, and stability boundary based on the finite-dimensional approximation (II) converge to those solutions of the eigenvalue problem (II). However, already the $N=2$ approximation is in an excellent

qualitative agreement and in a very reasonable quantitative agreement with the solution of (II).

For completeness, Appendix B reports the perturbation formulas for the singular flutter boundary, which are valid for arbitrary dimension $N$ of the discretized model.

IV Experimental detection of the singular flutter limit

IV.1 Experimental realization of the Pflüger column

Inspired by the Ziegler set-up BN2011 , a new mechanical device (Fig. 6) has been designed and realized to induce a follower force at the end of a Pflüger column. The force (whose magnitude is continuously acquired with a miniaturized load cell) is produced by friction generated through sliding of a freely rotating wheel against a conveyor belt and can be calibrated as proportional (through the Coulomb friction rule) to a vertical load (provided via frictionless contact with a glass plate, loaded through a pulley system) pressing the wheel against the conveyor belt (which was running at a constant speed of $0.1$ m/s in all experiments) sm .

IV.2 Identification of internal and external damping

During vibration of a rod two types of dissipations arise: an external (due to the air drag) and an internal (due to the viscosity of the constitutive material of the rod) damping. Often external and internal damping are condensed in a single coefficient, but it was shown SW1975 ; T2016 that for problems of flutter a careful distinction has to be maintained between the different sources of damping, as both strongly influence results. Therefore, experiments were performed to identify the two damping parameters introduced in the model, namely, a viscous modulus $E^{*}$ (modelling the internal damping) and an air drag coefficient $K$ (corresponding to a distributed external damping). To this purpose, the viscoelastic rod used for the flutter experiments was mounted on a shaker in a cantilever configuration

and the acceleration of its free end measured when the basis was imposed a sinusoidal displacement of a frequency corresponding to the first two modes of resonance. Results from these experiments were used with a modified logarithmic decrement approach detailed in Appendix C, to obtain the following values of the internal and external damping coefficients: $E^{*}=2.139796\cdot 10^{6}\,kg\,m^{-1}s^{-1}$ and $K=1.75239\cdot 10^{-5}\,kg\,m^{-1}s^{-1}$.

IV.3 Detection of the singular limit for the flutter onset

Our experiments are compared with the numerical solution of the boundary eigenvalue problem (II). The roots of the characteristic equation (II) are the eigenvalues $\omega$ governing the vibrations of the Pflüger column. The first two eigenvalues with their conjugates are plotted in Fig. 7 versus the load $p$, with all the other parameters kept fixed. In the absence of both the Kelvin-Voigt damping ($\eta$) and the air drag ($\gamma$), the Pflüger column is a reversible system and loses stability by flutter via collision of imaginary eigenvalues in a linear reversible-Hopf bifurcation, Fig. 7(a). In the presence of the two dissipation mechanisms, the merging of modes is imperfect, thus yielding flutter through the classical Hopf bifurcation at a value of $p$ significantly lower than in the case when the dissipation source is absent, Fig. 7(b). The theory of the previous section predicts that when the damping coefficients tend to zero while their ratio is kept constant, a limiting value of the flutter onset is reached, which generically differs from the flutter onset of the undamped column, thus justifying the numerical results of T2016 .

The critical flutter load for the Pflüger column was experimentally investigated covering a wide range of values of the mass ratio $\alpha$, Table 1. Note that, since $E^{*}$ and $K$ are constant, the geometry of the tested rods parameterizes the dimensionless damping coefficients $\eta$ and $\gamma$ according to Eqs. (II), so that different values of $\gamma$ and $\eta$ are obtained for rods of different length $(l)$ and thickness $(b)$.

The results of the measurements, together with the numerical calculations T2016 ; P1955 , are shown in Fig. 8 for eleven samples (see Table 1) in the plane $p$ versus $\alpha$. Theoretical critical curves, pertaining to samples of different lengths and thicknesses, are plotted and highlighted for the relevant intervals of $\alpha$. These boundaries are well-separated from the flutter boundary of the undamped system, represented by the upper dashed curve. In cases when either $\eta=0$ (the dot-dashed curves) or $\gamma=0$ (the lower dashed curves) the difference between the flutter boundaries corresponding to samples of various geometry is hardly visible, as it should be, in agreement with the theory, when the damping coefficients are very small T2016 ; P1955 ; B1956 ; KV2010 ; Z1952 . In contrast, when both damping mechanisms are taken into account, the critical curves dramatically differ for samples of different length and thickness. This is because the ratio $\beta=\gamma/\eta=(K/E^{*})(l^{4}/J)$ between the two damping coefficients increases almost 25 times from the first sample to the eleventh (see Table 1), although the damping coefficients $\gamma$ and $\eta$ vary weakly with the sample geometry.

Assuming $\gamma=\beta\eta$ in equation (II) and fixing $\eta$ to be one of the values reported in Table 1, the flutter boundary is plotted in Fig. 9 in the $p$ versus $\beta$ representation. Since for every length and thickness the critical flutter load depends weakly on $\alpha$, see Fig. 8, the flutter boundaries in Fig. 9, inset (a), are situated very close to each other (cf. Fig. 5). If the results of the measurements are superimposed, the experimental points perfectly fit this family of boundaries, within the error bands. Both the theoretical curves and the experimental points lie below the critical values of the undamped system for all values of $\alpha$. Nevertheless, the critical flutter load of the weakly damped Pflüger column is very sensitive to the damping ratio and increases as $\beta$ increases with the tendency to touch the lowest of the ideal flutter boundaries at $\beta>1000$, where the critical loads of the damped and undamped system tend to coincide (within the error bands), Fig. 9.

IV.4 The flutter modes

The analysis of the experiments is complemented by the determination of the flutter modes, which can be pursued by calculating the eigenvectors associated to the eigenvalues solutions of equation (7). The knowledge of the flutter modes is in fact useful to identify the shape of the vibrating rod during experiments. The analysis of the eigenvectors is reported in Fig. 10, relative to the first (lower frequency) vibration branch for the sample n. 5 of Table I, with dimensionless dampings $\eta=0.557\cdot 10^{-3}$ and $\gamma=48.368\cdot 10^{-3}$. All modes 1-3 in the figures refer to stable vibrations, while the onset of flutter corresponds to the mode numbered 4 and the onset of divergence to the mode numbered 9.

It is evident from Fig. 10(1) that the shape of the vibration mode corresponds (as it should be) at null $p$ to the free vibrations of a cantilever rod with a concentrated mass on its tip, vibrating at first resonance frequency. When the load $p$ increases beyond the threshold of the classical-Hopf bifurcation and approaches the higher value of the load corresponding to the threshold of the reversible-Hopf bifurcation in the undamped case, the vibrations become more and more similar to the second vibration mode of the free cantilever rod. This is not surprising in view of the fact that in the undamped case the eigenvectors of the first and the second mode merge at the flutter threshold because of the formation of a double imaginary eigenvalue with the Jordan block. In all the performed experiments the modes sketched in Fig. 10 have been observed.

V Conclusion

The theoretically predicted singular limiting behavior for the onset of the classical Hopf bifurcation has been detected and can now be considered as experimentally confirmed for a nearly-reversible system in the limit of vanishing dissipation.

This effect has been both theoretically and experimentally analyzed on a classical paradigmatic model of a nearly-reversible system, namely, the Pflüger viscoelastic column moving in a resistive medium under the action of a tangential follower force. For the theoretical treatment the continuous non-self-adjoint boundary eigenvalue problem has been Galerkin-discretized and reduced to a finite-dimensional matrix eigenvalue problem. With the use of perturbation theory of multiple eigenvalues, explicit expressions for the critical flutter load with and without dissipation have been derived thus proving the Whitney umbrella singularity at the interface between the classical Hopf bifurcation of the dissipative Pflüger system and the reversible-Hopf bifurcation of its undamped version. The conducted experiments with the laboratory realization of the Pflüger column confirmed the high sensitivity of the flutter onset to the damping ratio and accurately fitted both the theoretically and numerically predicted laws.

The designed, manufactured and tested ‘flutter machine’ opens a way to dedicated experiments on dissipation-induced instabilities with multiple damping mechanisms in a controlled laboratory environment.

Acknowledgments– We thank I. Hoveijn for useful discussions. The authors gratefully acknowledge financial support from the ERC Advanced Grant ERC-2013-ADG-340561-INSTABILITIES.

Appendix A Discretization

A.1 Adjoint boundary eigenvalue problems

The boundary eigenvalue problem for the Pflüger column with partial follower load is given by equation (II). The problem is self-adjoint only for $\chi=0$ and non-self-adjoint otherwise. Indeed, integration by parts of the differential equation (II) together with the boundary conditions lead to the following adjoint boundary eigenvalue problem

[TABLE]

The problem (A.1) coincides with (II) only for $\chi=0$. Otherwise, the boundary conditions of the two problems differ.

A.2 Variational principle

Let us consider now the functional

[TABLE]

Integrating by parts the first two terms in equation (A.2) and accounting for the boundary conditions for the problems (II) and (A.1), leads to

[TABLE]

On the other hand, the last of the boundary conditions (II) provides

[TABLE]

Hence,

[TABLE]

Stationarity of this functional with respect to arbitrary smooth variations $\delta\tilde{v}$, $\delta\tilde{w}$, which satisfy kinematic boundary conditions, is equivalent to the boundary value problems (II), (A.1).

A.3 Discretization and reduced finite-dimensional model

Let us consider solutions to the self-adjoint problems (II) and (A.1), with $\chi=0$, $p=0$, $\eta=0$, $\gamma=0$, and $\alpha=0$

[TABLE]

where $\omega_{j}$ is a root of the characteristic equation

[TABLE]

which provides for instance,

[TABLE]

The functions (A.3) are orthogonal and normalized as follows:

[TABLE]

Therefore, the eigenmodes $\tilde{v}$ and $\tilde{w}$ can be represented in the form of the expansions

[TABLE]

where $\tilde{w}_{j}=\tilde{v}_{j}$.

Substituting the expansions (36) into the functional (A.2) yields the discretized version of the functional (A.2)

[TABLE]

The gradient of the discretized functional, $I_{N}$, calculated with respect to the vector of coefficients ${\bf b}=(b_{1},b_{2},\ldots,b_{N})$, and equated to zero, provides the discretized eigenvalue problem for the Pflüger column

[TABLE]

where ${\bf a}=(a_{1},a_{2},\ldots,a_{N})$ and the elements of the matrices are

[TABLE]

with $\delta_{ij}$ denoting the Kronecker symbol. The entries of the matrix ${\bf K}_{2}$ in the explicit form are

[TABLE]

where $A=\frac{4\sqrt{\omega_{i}\omega_{j}}}{(-1)^{i}\omega_{i}-(-1)^{j}\omega_{j}}$. All the matrices are real. In addition, the matrices of mass, $\bf M$, external damping, ${\bf D}_{e}$, internal damping, ${\bf D}_{i}$, and stiffness, ${\bf K}_{1}$ and ${\bf K}_{2}$, are symmetric. The matrix of nonconservative positional forces with non-zero curl, $\bf N$, is real and non-symmetric. Note that $\det{\bf M}=1+4N\tan\alpha>0$.

Appendix B Perturbation formulas for arbitrary $N$

The eigenvalue problem (9) can be formulated as the eigenvalue problem

[TABLE]

for the matrix polynomial

[TABLE]

where ${\bf k}=(p,\chi,\gamma,\eta,\alpha)$ is a vector of parameters. The adjoint matrix polynomial ${\mathbf{L}}^{\dagger}={\mathbf{A}}^{T}+{\bf D}{\overline{\omega}}+{\mathbf{M}}{\overline{\omega}}^{2}$ is introduced, so that $({\mathbf{L}}{\mathbf{a}},{\mathbf{b}})=({\mathbf{a}},{\mathbf{L}}^{\dagger}{\mathbf{b}})$, where the inner product is defined as $({\bf a},{\bf b})=\overline{\bf b}^{T}{\bf a}$. With this definition, the adjoint eigenvalue problem can be rewritten as

[TABLE]

Let us assume that, for the values of the parameters $\chi=\chi_{0}$, $\alpha=\alpha_{0}$, $\gamma=0$, $\eta=0$, and $p=p_{0}$, an algebraically double imaginary eigenvalue $\omega_{0}=i\sigma_{0}$ exists with the Jordan block which satisfies the following equations

[TABLE]

where ${\mathbf{a}}_{0}$ is an eigenvector and ${\mathbf{a}}_{1}$ is an associated vector at $\omega_{0}$. Then, an eigenfunction ${\mathbf{b}}_{0}$ and an associated function ${\mathbf{b}}_{1}$ at the complex-conjugate eigenvalue $\overline{\omega}_{0}=-i\sigma_{0}$ are governed by the adjoint equations

[TABLE]

Note the orthogonality between the eigenvectors, that is

[TABLE]

When the parameter $p$ is perturbed in the vicinity of $p_{0}$ as $p=p_{0}+\Delta p$, an approach similar to that used for $N=2$ yields

[TABLE]

where ${\mathbf{A}}^{\prime}_{p}=\left.\frac{\partial{\mathbf{A}}}{\partial p}\right|_{p=p_{0}}$. Therefore, the eigenvalues and eigenvectors of the undamped reversible system can be approximated in the vicinity of $p=p_{0}$, i.e. in the vicinity of the flutter boundary corresponding to the reversible-Hopf bifurcation.

Assume that at $p<p_{0}$ the eigenvalues of the undamped reversible system are imaginary, $\omega(p)=i\sigma(p)$, with an eigenvector ${\mathbf{a}}(p)$ and the eigenvector of the adjoint problem ${\mathbf{b}}(p)$. Then, at $p>p_{0}$ the eigenvalues (B) are complex-conjugate (denoting instability). A dissipative perturbation with the matrix ${\mathbf{D}}(\eta,\gamma)$ where ${\mathbf{D}}(0,0)=0$ changes the eigenvalue $\omega(p)=i\sigma(p)$ as follows

[TABLE]

The following condition for the imaginary eigenvalue is assumed to hold

[TABLE]

so that the eigenvalue remains imaginary after a dissipative perturbation. This means that the neutral stability surface is not abandoned after the dissipative perturbation. Using the perturbation formulas (B) for ${\mathbf{a}}(p)$ and ${\mathbf{b}}(p)$ in (46), introducing the damping ratio $\beta=\gamma/\eta$, and defining

[TABLE]

the following quadratic approximation in $\beta$ can be found to the critical flutter load in the limit of vanishing dissipation

[TABLE]

From the orthogonality of eigenvectors (43) and the expression for the mass matrix ${\bf M}_{0}={\bf I}+4{\bf M}_{1}\tan\alpha_{0}$ it follows immediately that the denominator in (47) vanishes at $\alpha_{0}=0$, thus confirming that in the case of the Beck column the external air drag damping is stabilizing. Now this result has been established for the discretized model of the Pflüger column of arbitrary dimension $N$.

In the case of $N=2$, $\chi_{0}=1$, $\alpha_{0}=0.1$, $p_{0}\approx 17.83368$, $\sigma_{0}\approx 9.366049$, the following vectors are obtained

[TABLE]

With these vectors the formula (B) exactly reproduces equation (20). The formula (47) provides $\beta_{0}\approx 1478.074$ in full accordance with equation (25) in the case of $N=2$. Finally, equation (48) exactly reproduces equation (26).

For $N>2$ the procedure is the same: one only needs to find the vectors ${\bf a}_{0}$, ${\bf a}_{1}$, ${\bf b}_{0}$, ${\bf b}_{1}$ solving (B) and (B) with the corresponding $N\times N$ matrices which entries are given by equations (A.3) and (A.3).

Appendix C Modified logarithmic decrement approach

Equations of motion

A viscoelastic rod is considered, made up of a material which follows the Kelvin-Voigt model

[TABLE]

where $\sigma_{z}$ and $\varepsilon_{z}$ are the longitudinal stress and strain, respectively, and $E$ and $E^{*}$ are the elastic and the viscous moduli. In an Euler rod the strain is defined as

[TABLE]

where $\phi^{\prime}$ is the curvature and $y$ the coordinate orthogonal to the rod’s axis $x$, so that the bending moment can be computed as

[TABLE]

and rewritten in terms of displacement $v(x,t)$ as

[TABLE]

The equation governing the dynamics of a straight rod is

[TABLE]

where $m$ is the mass density per unit length of the rod and $p$ is the transversal load per unit length, which can be identified with the sum of an applied load $f(t)$ and a force proportional (through a coefficient $K$) to the velocity $\dot{v}$, to model external damping. A substitution of Eq. (62) into Eq. (63) yields

[TABLE]

A sinusoidal excitation at the clamped end of a rod in a cantilever configuration can be modeled with a specific form of external load, namely

[TABLE]

where $U_{0}$ is the amplitude of the displacement imposed at the clamp, which varies in sinusoidally in time with pulsation $\bar{\omega}$.

C.1 Free vibration of a cantilever rod

The solution of Eq. (64) with an imposed sinusoidal displacement in terms of $v(x,t)$ can be found exploiting the separation of variables

[TABLE]

where the function $Y_{i}(x)$ and $y_{i}(t)$ are mode functions, respectively, in space $x$ and in time $t$. The force $f(t)$ acting on the rod plays a role only in the definition of the $y_{i}(t)$ modes. Assuming a function of time $y(t)=\exp(-i\omega t)$ yields the characteristic equation

[TABLE]

where $\Lambda_{n}^{4}$ is a real quantity (dimensionally equal to [length]-4)

[TABLE]

The solution to Eq. (C.1) is a sum of periodic and hyperbolic functions

[TABLE]

where the constants $C_{i,n}$ depend on the boundary conditions. For a cantilever rod, the boundary conditions are

[TABLE]

A substitution of the boundary conditions in Eq. (69) yields in a matrix form

[TABLE]

The first two equations yield

[TABLE]

so that Eq. (71) reduces to

[TABLE]

Imposing the determinant of the matrix of the system (C.1) to vanish provides

[TABLE]

Eq. (73) defines the $\Lambda_{n}$ values as

[TABLE]

Now, the solution for the functions $Y_{n}(x)$ can be expressed in terms of one arbitrary constant $C_{1,n}$, so that

[TABLE]

which leads to the general solution for the free vibrations of a cantilever rod expressed as an infinite sum of the following mode functions

[TABLE]

C.2 Properties of the function $Y_{n}(x)$

The free vibration shape equations $Y_{n}(x)$ satisfy the orthogonality relations

[TABLE]

Morover, equation (C.1) allows to write

[TABLE]

It is expedient now to define the quantity

[TABLE]

so that equation (76) yields

[TABLE]

C.3 Expression of $y(t)$ for a cantilever rod with a base motion excitation

The differential equations governing the sinusoidal motion of the clamped rod subject to the force $f(t)$, Eq. (65), are

[TABLE]

In order to exploit the orthogonality property of the shape functions $Y_{n}(x)$, each term of the previous equation is multiplied by $Y_{k}(x)$ and integrated over the length of the rod $l$, which provides the expression

[TABLE]

where $F_{n}=\int_{0}^{l}Y_{n}(x)dx$. Eq. (C.3) reminds the equation of motion which governs a single-degree-of-freedom system with a mass $m_{n}$, a damper with constant $c_{n}$ and a spring with stiffness $k_{n}$

[TABLE]

where

[TABLE]

Another form of Eq. (81) is

[TABLE]

where

[TABLE]

The solution of the differential equation (83) is expressed as the sum of the solution of the associated homogeneous equation and of a particular integral. The latter can be found in the form

[TABLE]

where the coefficients $A_{n}$ and $B_{n}$ satisfy Eq. (83) and assume the form

[TABLE]

in which $N_{n}$ is the so-called ‘dynamic amplification factor’

[TABLE]

The solution of the homogeneous equation is

[TABLE]

where $\alpha_{n,d}=\alpha_{n}\sqrt{1-\zeta_{n}^{2}}$ are the damped pulsations of the system.

The coefficients $C_{n}$ and $D_{n}$ can be found by imposing the initial conditions

[TABLE]

in the complete solution of

[TABLE]

which leads to the expressions

[TABLE]

C.4 Relation between $\zeta_{n}$, $E^{*}$, $K$

The relation between the damping ratio $\zeta_{n}$, the internal ($E^{*}$), and the external ($K$) damping is described by Eq. (C.3)2, which can be rewritten as

[TABLE]

The problem of the identification of the two damping coefficients thus reduces to the quantification of the damping ratio $\zeta_{n}$ relative to two different modes. The logarithmic decay over $k$ cycles can be written as

[TABLE]

where $\delta_{k}=\log\left(y_{1}/y_{k+1}\right)$.

The dimensionless internal and external damping coefficients can be finally expressed through the relations

[TABLE]