On the effects of linear damping on the nonlinear Ziegler’s column

In this paper, the effects of linear damping on the post-critical behavior of the Ziegler’s column are discussed. To this end, the well-known double-pendulum, loaded at the free-end by a follower force, firstly introduced by Ziegler, is considered in regime of finite displacements. The multiple scale method is applied to the equations of motion expanded up to the cubic terms, to analyze the nonlinear behavior of a generically damped column, close to the simple-Hopf bifurcation triggered by the follower force. The obtained bifurcation equations are shown to be useful in providing qualitative information about the nonlinear mechanical response of the column in the whole damping plane. Validation of the asymptotic solution, carried out via numerical analyses of the exact equations of motion, points out the effectiveness of the proposed analysis also on the quantitative side.


Introduction
The 'destabilizing effect of damping' or the 'Ziegler's Paradox' is a mechanical phenomenon, which has been in depth studied by many researchers in the last decades . The paradox consists of a finite reduction in the circulatory Hopf critical load, occurring when a linear circulatory system, i.e., an undamped system loaded by follower forces (see, e.g., [22,[24][25][26][27][28][29]), is perturbed by a small and positive-definite damping. In this framework, the well-known Ziegler's column, namely a double-pendulum, visco-elastically restrained at the hinges, and loaded at the free end by a follower force, can be regarded as a paradigmatic system of such a paradoxical behavior.
Literature on the linear stability analysis of the column, aiming to explain the destabilization phenomenon, is wide and rich; a review of the subject can be found in [16]. Theories sided on geometrical, numerical and asymptotic methods have been proposed, among the others, in [12,14,18,30,31]. The existence of a singular surface in the space of damping and load parameters, called the 'Whitney's umbrella surface' [32], is shown in [30,31]. Sensitivity analyses of the eigenvalues are applied in [14,18] through perturbation methods, by starting the asymptotic expansions from the critically loaded circulatory system in [14], and from a marginally stable, sub-critically loaded undamped system in [18]. In all these papers, it is shown that the effect of adding a small damping on the undamped Ziegler's column is mostly detrimental, since, in except for a small region in the plane of damping parameters, the critical load of the damped system is lower than that of the circulatory one.
Very few papers have been focused to investigate the nonlinear behavior of the Ziegler's column, i.e., to characterize the dynamics of the system once the critical load is exceeded and finite displacements are considered. Among them, in [33,34] the Normal Form Theory and second method of Lyapunov are applied to investigate the stability of the Ziegler's column, also in the presence of different damping forms: no information is given about the amplitude of the limit-cycle, arising in the post-critical regime. In [35], the nonlinear dynamics of the Ziegler's double-pendulum, loaded by a partially follower force, is investigated by the multiple scale method [36][37][38]: the supercritical nature of the limit-cycle is shown in the case of equal damping coefficients. In [39], a first-order multiple scale analysis is developed in the presence of nonlinear damping of Van der Pol-type: it is shown that, when a proper nonlinear damping form is added to a linearly damped Ziegler's column, it can reduce the, otherwise large, limit-cycle amplitude. Moreover, in [39] it is also found that nonlinear damping triggers the so-called hard loss of stability phenomenon, according to which a periodic subcritical branch turns back, regaining stability. In [40], a second-order multiple scale analysis is developed to shed light on the hard loss of stability of a nonlinearly damped Ziegler's column. It is worth to notice that, both in [39,40] numerical proofs of the supercritical nature of the limit-cycle are found for some linearly damped columns, taken as case studies.
The present work aims to enhance the research so far carried out. In particular, attention is focused on important aspects of the nonlinear behavior of the linearly damped Ziegler's double-pendulum which have been left open by the above-mentioned papers, namely: how does the linear damping affect the dynamics close to the simple-Hopf bifurcation? more precisely, is the limit-cycle supercritical in the whole damping plane? is it always stable? how is the limit-cycle amplitude changed by damping?
To give an answer to these questions, the multiple scale method is suitably applied in this work, in order to analyze the post-critical behavior of the generically damped Ziegler's column, undergoing simple-Hopf bifurcation. The relevant bifurcation equations are shown to be useful in providing qualitative and quantitative information on the nonlinear mechanical behavior of the system. The asymptotic results are validated against benchmark solutions, obtained via a direct numerical approach, applied on the exact equations of motion.
The paper is organized as follows. In Sect. 2, the equations of motion of the Ziegler's column are presented. In Sect. 3, a background on the linear stability analysis is given. In Sect. 4, the nonlinear problem is tackled via the multiple scale method, while, in Sect. 5, numerical results are presented. In Sect. 6, the main conclusions are drawn. Finally, three appendixes give some definitions relevant to the model, to the linear stability analysis and to the multiple scale algorithm.

Model
The system depicted in Fig. 1, known as the Ziegler's column, consists in two hinged weightless rigid bars of equal length , carrying the lumped masses, m 1 := 2m and m 2 := m. The bars are visco-elastically restrained at the hinges by two linear springs of stiffness k 1 := k and k 2 := k, and two linear dashpots, of viscosity coefficients c 1 and c 2 . Moreover, a follower force of intensity F, which maintains its direction parallel to the upper bar, is applied at the free end.
By taking the rotations of the two bars, viz. ϑ 1 (t) and ϑ 2 (t) as Lagrangian coordinates, t being the time, the nondimensional form of the equations of motion reads (see [34,40]): where dot denotes t-differentiation and, moreover, the following quantities have been introduced: and tilde omitted for simplifying the notation. When Eq. (1) is expanded up to the cubic terms and recast in state-space form, it reads: where x := (x 1 , x 2 , x 3 , x 4 ) T is the vector collecting the state variables, having defined x 1 = ϑ 1 , x 2 = ϑ 2 , x 3 =θ 1 , x 4 =θ 2 , μ is the non-dimensional load parameter and ξ 1 ≥ 0, ξ 2 ≥ 0 are the non-dimensional damping coefficients. Matrix A is the so-called state matrix, depending on mass, stiffness, damping and circulatory contributions, namely: and F (x, x, x), F j (x, x, x) ( j = μ, ξ 1 , ξ 2 ) are vectors collecting the nonlinear terms. All the quantities, appearing in Eqs. (3) and (4), are defined in "Appendix A".

Background on linear stability analysis
In this section, the background on linear stability analysis is briefly illustrated. The aim of such an analysis is to build up the plot of the organized parameter-space, i.e., the linear stability diagram, which gives a comprehensive scenario of the system behavior near to the equilibrium. Linear stability of the column, around the equilibrium x = 0, is governed by the eigenvalue problem associated to the linear part of Eq. (3). By taking x = ue λt in the linear part of Eq. (3), the following eigenvalue problem is obtained: (λ, u) being the eigenpairs and I the identity matrix. When ξ 1 = ξ 2 = 0, i.e., the system is undamped (also referred to as circulatory), and μ = 0, Eq. (5) admits two pairs of purely imaginary eigenvalues, which lie on the imaginary axis. In this condition, the system is (marginally) stable. If μ increases from zero, the eigenvalues move on the imaginary axis, still remaining distinct (see Fig. 2a), until the load reaches a critical value μ = μ c , at which they collide and a circulatory (or reversible) Hopf bifurcation takes place; if an infinitesimal increment δμ > 0 is given to the critical value μ c (subscript c means circulatory), they separate and instability occurs. The exact solution of the eigenvalue problem (5) furnishes (see, e.g., [17]): Therefore, at the critical state, the eigenvalue problem (5) admits a double defective eigenvalue ±i ω c , and, consequently, the basis of the proper eigenvectors must be completed by generalized eigenvectors (see, e.g., [12]). When the system is damped, and μ = 0, the complex conjugate eigenvalues are on the left of the imaginary axis, i.e., the system is asymptotically stable. When μ increases and reaches a threshold value, namely μ = μ d (subscript d means damped), a complex conjugate eigenvalue ±i ω d crosses from the left to the imaginary axis, (see Fig. 2b). In this condition, a generic Simple-Hopf bifurcation takes place, for which the exact analysis of problem (5) gives (see, e.g., [17]): At the critical state, the eigenvalue problem (5) admits proper right and left eigenvectors, u d and v d associated with ±iω d , which are the solutions of the following linear problems: where Proper (left and right) eigenvectors have been reported in "Appendix B". The 'destabilizing effect of damping' or the Ziegler's paradox consists in the fact that, when small and positive-defined damping is added to a circulatory system, a finite reduction in critical load is experienced, i.e., μ d < μ c (see, e.g., [2,3,5,12,14,16,[18][19][20]). The exact boundary of the linear stability diagram of the Ziegler's column is given by Eq. (7a), which defines the critical locus in the (μ, ξ 1 , ξ 2 )-space, at which incipient instability manifests itself. The locus is displayed in Fig. 3a (see [18]) and separates the stable region (below the surface, label S in the figure) from the unstable one (above the surface, label U in the figure). The surface is known in the literature as the 'Whitney's umbrella' surface [16,32]. The points of the surface represent generic Hopf bifurcation states of damped systems, while at μ = μ c and ξ 1 = ξ 2 = 0, instability is triggered by a circulatory Hopf bifurcation. In Fig. 3b, the contour lines μ = const of the surface are displayed: each of them divides the damping plane in a stable region (label S, on the right of the contour line) and in an unstable one (label U , on the left), and, below the value of μ = 0.33, the column is stable for any (positive) damping. It is apparent that a generic small damping, of modulus ρ := ξ 2 1 + ξ 2 2 , and polar angle ϕ := arctan (ξ 2 /ξ 1 ), entails a detrimental effect on the stability of the circulatory system: indeed, in except for a small region in the (ξ 1 , ξ 2 )-plane (shaded in blue in Fig. 3b), which is close to an optimal direction, i.e., ϕ opt := arctan 4 Fig. 3b), the critical load μ d = μ d (ϕ, ρ) is less than μ c . In the limit for ρ → 0, i.e., for an evanescently damped system, and when ϕ ≡ ϕ opt , it is μ d → μ c , i.e., the optimal direction is the only one for which the critical load of a damped system recovers that of the circulatory one: this aspect is emphasized in Fig. 3c, where the lim ρ→0 μ d , versus the polar angle ϕ, is represented.
It is important to remark that, a point on the surface, when ρ is large, or for small ρ and ϕ = ϕ opt , represents a bifurcation state of the column at which the eigenvalues are well-separated, one critical and the other stable; in this case instability is, therefore, triggered by a simple-Hopf bifurcation. On the contrary, for small ρ and ϕ ϕ opt , the eigenvalues are close to each other and an interaction between can occur; here the system is nearly-defective, and its loss of stability may be governed by the defective-Hopf bifurcation of the (close) circulatory system. However, the region of the space of damping parameters at which the nearly defectiveness of the eigenvalues plays a role is not a-priori known and, therefore, deserves to be investigated.
Finally, mechanical interpretations of the Ziegler's paradox have been given in the literature also through asymptotic approaches: these are based on the eigenvalue sensitivities of the circulatory system, perturbed by an evanescent damping (see, e.g., [14,18,19,53]. Both perturbations of defective double eigenvalue [14], and simple eigenvalue [18,19,53], the latter relevant to a (marginally stable) and sub-critically loaded circulatory system (μ < μ c ), have been developed to build up an asymptotic expression of the locus (7a), to shed the light on the true essence this fascinating paradoxical phenomenon.
Notwithstanding the different bifurcation mechanisms described in the previous section, i.e., the Simple-Hopf (SH) and the Defective-Hopf (DH) bifurcations, a MSM algorithm is applied around the SH bifurcation, of a generically damped system, which takes place at μ = μ d . It is worth to notice that, the algorithm is expected to give reliable results for almost all the family of damped systems, i.e., in except for evanescently damped columns, whose damping ratio is close to the optimal one, namely for nearly-defective systems. Contrarily to this conjecture, it will be shown ahead that, the developed MSM algorithm is very effective, when (heuristically) applied close to the circulatory Hopf bifurcation.
The starting point of the asymptotic expansion is the SH bifurcation point, namely: μ = μ d and (i ω d , u d ) the critical eigenpairs. The variable and the bifurcation parameter are rescaled as x → ε 1/2 x, μ → μ d + εδμ, where δμ denotes the increment of load measured from the bifurcation point and 0 < ε 1 is a bookkeeping parameter to be reabsorbed at the end of the procedure. Several independent integer time-scales, namely t k = ε k t, k = 0, 1, 2, . . ., are introduced, such that d/dt = d 0 + ε 1 d 1 + ε 2 d 2 + · · · , in which d k = ∂ k /∂t k . By expanding the Lagrangian coordinates in series of integer powers of ε, namely x = x 0 +ε 1 x 1 +ε 2 x 2 +· · · , and substituting in Eq. (3), the following chain of perturbation equations is obtained The ε 0 -order problem admits the solution: where A is an unknown complex amplitude, depending on slow time-scales, and c.c. stands for complexconjugate terms. By substituting Eq. (10) in the εorder problem (9b), a non-homogeneous problem in the unknown x 1 is obtained. The removal of secular terms of this latter, i.e., imposing the so-called solvability condition, leads to a first-order amplitude modulation equation (AME). After reconstitution (see "Appendix C") and coming back to the true time and quantities, the AME reads: where the following definitions hold: c = c (ξ 1 , ξ 2 ) being a complex function of the damping parameters. Expressions of b and c are reported in "Appendix C". Equation (11) can be written in real form by the way of A = 1 2 ae iθ , from which: and (·) R = Re (·), (·) I = Im (·) denotes real and imaginary parts of the involved quantities, respectively. Equations (13) govern the dynamics around the bifurcation. Once they have been solved, the state reads: where h.o.t. stands for higher-order terms.

Qualitative analysis of the post-critical scenario
The coefficient α in Eq. (16a) governs, for a given increment of load δμ = μ−μ d , the amplitude a s : the higher α, the smaller a s . In Fig. 4a, the plot of the coefficient α, as function of damping parameters ξ 1 , ξ 2 , is shown. Remarkably, it resembles the shape of the Whitney's umbrella surface of the critical load (see Fig. 3a). In Fig. 4b, the contour lines α = const of the surface are shown. It is seen that for a generic small damping modulus ρ, there exists a polar angle ϕ nl , which maximizes α,  i.e., α = α max 0.24, when ϕ = ϕ nl 0.21 = ϕ opt : therefore, in except for a small region in the (ξ 1 , ξ 2 )plane (shaded in blue in Fig. 4b), which is close to ϕ nl (dashed red line in Fig. 4b), the coefficient α < α max . As a consequence, a generic damping produces a detrimental effect also in the nonlinear field since α significantly decreases when ϕ moves away from ϕ nl and, hence, a s increases.
Moreover, Fig. 4c shows the plot of the coefficient α as a function of the polar angle ϕ and for different ρ values. It is apparent that, there exist two intervals of the polar angle, for which the coefficient α can increase or decrease with an increasing damping modulus. Indeed, for any ϕ ∈ I 1 := (0, 1.464), the larger is the damping modulus, the larger is α, i.e., a s decreases when ρ increases; on the contrary, for ϕ ∈ I 2 := (1.464, π/2], α decreases when ρ increases, i.e., a s increases with ρ. Consequently, when ϕ = 0 or ϕ 1.464, α and a s are not influenced by damping modulus. Remarkably, damping has a detrimental effect when ϕ ∈ I 2 , since it enlargers the amplitude of the stationary solution, counterintuitively from what it is expected from a physical point of view. It is worth to notice that, the amplitude of the limitcycle depends on a s , as described here, but also on u d = u d (ξ 1 , ξ 2 ), see Eq. (15). Therefore, for a given δμ, the effect of an increasing ρ on the maximum value attained by the components of motion relies also on how u d is changed by damping. Due to the chosen normalization (see "Appendix B"), the first two components of u d are u d1 = u d1 (ξ 1 , ξ 2 ) and u d2 = 1, so that only max (ϑ 2 ) has the same qualitative behavior of a s when damping changes.
As a final remark, it is believed that such a qualitative analysis can be useful for a damping design, based not only on the linear analysis, but also on the optimization of the post-critical response of the system, i.e., of the limit-cycle amplitude.

Numerical results
Numerical analyses, carried out in this section, aims at confirming, on a quantitative side, the previously dis- To these purposes, five damped systems, taken as case studies, are defined as reported in Table 1, namely by varying the polar angle ϕ, and for different damping moduli. Consequently, a sufficiently wide region of the damping plane parameters is explored. In particular, the polar angle in cases from I to III is chosen to be close or coincident with ϕ opt , while in cases IV and V it is far from ϕ opt but, respectively, inside and outside the interval I 1 , introduced in the qualitative analysis of Sect. 4. Figures 5, 6, 7, 8 and 9 show the bifurcation diagrams relevant to the case studies from I to V, respectively; here the values of max ϑ j , j = 1, 2, are plotted versus the load parameter μ. The black curves represent limit-cycles obtained through the numerical continua- The bifurcation diagrams relevant to the case study I, II and III, (see Figs. 5, 6 and 7) show that an increase of ρ entails: (i) for a given δμ, an increase of max (ϑ 1 ) and a decrease of max (ϑ 2 ); (ii) a forward shifting of the curves, due to the increasing of μ d with ρ. A very good accordance between numerical and asymptotic results, also for small damping modulus, is obtained. This occurrence points out that the MSM algorithm, starting from the SH bifurcation, gives reliable results also close to the DH bifurcation. Such an unexpected result is probably due to the fact that, the region of the damping plane in which the eigenvalues are nearly defective and can interact is very small; however, an indepth investigation of these aspects calls for the analysis of the defective double-Hopf, which is beyond the scope of the present paper. The results relevant to the case study IV, displayed in Fig. 8, show a qualitative behavior similar to that of the previous cases. However, it is worth to notice that, in this case, the good agreement between numerical and asymptotic solutions is expected, since linear damping has a medium detrimental effect on the linear stability, and the spectrum of eigenvalues is well-separated. It is important to remark that the case of uniform damping was also discussed in [35].
In Fig. 9, the bifurcation diagrams of case study V are reported. Similarly to the other case studies, the MSM algorithm is in excellent agreement with the numerical solution. However, differently from them, here an increase of ρ entails: (i) for a given δμ, a decrease in both max (ϑ 1 ) and max (ϑ 2 ); (ii) no forward shifting since when ϕ = π/2, μ d is constant with respect to ρ (remember Fig. 3b). As a final remark, it is worth to notice that in all the analyzed case studies, large amplitude limit-cycles are detected: these are very well captured by the asymptotic solution even if some discrepancy arises for large δμ.
Finally, Figs. 10 and 11 report the results of numerical integrations, relevant to case studies I (at μ = 2.2) and IV (at μ = 1.7), respectively. These are performed on both the exact equations of motion (1), black curves, and on Eq. (13), blue curves, which have been used to build-up the state (15). In particular, in subfigures (a) and (b) the time-histories of ϑ 1 (t) and ϑ 2 (t) are displayed, while in subfigures from (c) to (f), projections of the trajectories onto phase planes, in steady-states, are plotted. In both the case studies, the effectiveness of the asymptotic procedure is confirmed. In particular, time-histories show that, after a transient in which the effect of the passive mode is exhausted, the asymptotic steady response is in good agreement with the numerically exact one. Small differences between the two solutions can be detected in the phase portraits, as due to the fact that the critical mode is not corrected by the asymptotic procedure at the first order; however, extending the algorithm to higher orders, as done in [40], could provide the needed modal corrections.

Conclusions
The effects of linear damping on the nonlinear dynamics of the Ziegler's column have been discussed in this paper. The work has been motivated by former papers by the same authors, in which the presence of supercritical limit-cycles, arising from the Hopf bifurcation, was detected for some damped columns, taken as case studies. A first-order asymptotic analysis has been developed through the multiple scale method, applied to the equations of motion expanded up to the cubic terms. The bifurcation equations, ruling the asymptotic dynamics of a generally damped system, around a simple-Hopf bifurcation, have been obtained. These have been shown to provide valuable qualitative information on how linear damping changes the post-critical behavior of the Ziegler's column. Asymptotic results have been validated on the quantitative side, with numerical solutions of the exact equations of motion. The following conclusions can be drawn.
1. The limit-cycle is supercritical and stable for any combination of (positive) damping parameters. 2. For a given load increment overcoming the threshold of linear stability, the amplitude of the asymptotic stationary solution is governed by a positive coefficient, whose representation, as a function of damping parameters, resembles the Whitney's umbrella surface of the critical load. 3. A very good agreement between the asymptotic and numerical solutions is detected in several case studies, which are representative of the behavior of the Zielger's column in a wide region of the damping plane. 4. The effectiveness of the asymptotic results has been shown also close to the defective double-Hopf bifurcation of the circulatory system. 5. In all the considered case studies, large amplitude limit-cycles have been detected. sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

A Details on the equations of motion
Matrices involved in Eq. (4) are defined as: The nonlinear (cubic) vector functions of Eq. (3) read: