Acta Polytechnica Vol. 51 No. 4/2011 Singularities in Structural Optimization of the Ziegler Pendulum O. N. Kirillov Abstract Structural optimization of non-conservative systems with respect to stability criteria is a research area with important applications influid-structure interactions, friction-induced instabilities, andcivil engineering. Incontrast tooptimization of conservative systems where rigorously proven optimal solutions in buckling problems have been found, for non- conservative optimization problems only numerically optimized designs have been reported. The proof of optimality in non-conservative optimization problems is a mathematical challenge related to multiple eigenvalues, singularities in the stability domain, and non-convexity of the merit functional. We present here a study of optimal mass distribution in a classical Ziegler pendulum where local and global extrema can be found explicitly. In particular, for the undamped case, the two maxima of the critical flutter load correspond to a vanishing mass either in a joint or at the free end of the pendulum; in the minimum, the ratio of the masses is equal to the ratio of the stiffness coefficients. The role of the singularities on the stability boundary in the optimization is highlighted, and an extension to the damped case as well as to the case of higher degrees of freedom is discussed. Keywords: circulatory system, structural optimization, Ziegler pendulum, Beck column, flutter, divergence, damping, Whitney umbrella. 1 Introduction Structural optimization of conservative and non-conservative systems with respect to stability criteria is a rapidly growing research area with important applications in industry [1–3]. Optimization of conservative elastic systems such as the problem of the optimal shape of a column against buckling is already non-trivial, because some optimal solutions could be multi-modal and thus correspond to a multiple semi-simple eigenvalue which creates a conical singularity of the merit functional [3]. Despite these complications, a number of rigorous optimal solutions are known in conservative structural optimization. Nevertheless, the increase in the critical divergence load given by the optimal design in such problems is usually not very large in comparison with the initial design [2, 3]. In contrast to conservative systems, non-conservative systems can loose stability both by divergence and by flutter. It is known that mass and stiffness modification can increase the critical flutter load by hundreds percent, which is an order of magnitude higher than typical gains achieved in optimization of conservative systems [4–18]. For example, Ringertz [9] reported an 838 % increase of the critical flutter load for the Beck column [19], from 20.05 for a uniform design to 188.1 for an optimized shape. Recently, Temis and Fedorov [17] found for a free-free beam moving under the follower thrust an optimal design with a critical flutter load that exceeds the load for a uniform beam by 823 %. We note that although the very notion of the follower forces was debated in [20–22], the Beck column [19] as well as its discrete analogues [23–25] including the Ziegler pendulum [26], remain popular models for investigating mode-coupling instabilities in non-conservative systems and related optimization problems. In both conservative and non-conservative problems of structural optimization of slender structures, their optimal or optimized shapes often possess places with small or even vanishing cross-sections. The known optimized shapes of the Beck column or of a free-free rod moving under follower thrust have an almost van- ishing cross-section, e.g., at the free end, which means vanishing mass of a finite element in the corresponding discretization [9–12, 17, 18]. Another intriguing feature of optimizing non-conservative systems is the ‘wandering’ critical frequency at the optimal critical load. During optimization the eigenvalue branches experience numerous mutual overlappings and veerings [4, 6, 7, 27–29] with a tendency for the critical frequency to increase and to correspond to higher modes [9, 14–17, 30]. This puzzling behavior of the critical frequency still awaits an explanation. In some problems, such as the optimal placement of the point mass along a uniform free-free rod moving under the follower thrust [5,8], the local maxima were found to correspond to singularities of the flutter boundary such as cuspidal points where multiple eigenvalues with the Jordan block exist [11, 12]. In order for the last phenomenon to happen, we need at least three modes [11, 12], meaning that two-mode approximations [5, 8] 32 Acta Polytechnica Vol. 51 No. 4/2011 are unable to detect such optima. This reflects a general question on model reduction and the validity of low- dimensional approximations in non-conservative problems, already discussed by Bolotin [30] and Gasparini et al. [23] and recently raised again in the context of friction-induced vibrations by Butlin and Woodhouse [31]. The above-mentioned phenomena make rigorous proofs of optimality in non-conservative optimization prob- lems substantially more difficult than in conservative ones. To the best of our knowledge, no rigorously proven optimal solutions in optimization problems for distributed circulatory systems have been found. Although the situation is not much better in the finite-dimensional case, it seems reasonable to try to understand the nature of the observed difficulties of optimization on final dimensional non-conservative systems that depend on a finite number of control parameters. Let us consider a circulatory system Mẍ + Kx = 0, (1) where dot indicates time differentiation, M is a real symmetric m×m mass matrix and K is a real non-symmetric m × m matrix of positional forces that include both potential and non-potential (circulatory) forces. Equation (1) typically originates after linearization and discretization of stability problems for structures under follower loads, in problems of friction-induced vibrations, and even in rotor dynamics when the damping is not taken into account [30, 28, 32]. The characteristic equation for the circulatory system (1) is given by the modified Leverrier algorithm [33]. In the case when m = 2, it reads detMλ4 + (trMtrK− tr(MK))λ2 + detK = 0, (2) where λ is an eigenvalue that determines the stability of the trivial solution. The coefficients of matrix K usually contain the loading parameter, say p, that we need to increase by varying the coefficients of, e.g., the mass matrix. During this optimization process some masses can come close to zero, so that the mass matrix can degenerate yielding detM = 0. As a consequence, some eigenvalues λ can increase substantially. However, such a singular perturbation of the characteristic polynomial may cause large values of the gradient of the critical load with respect to the mass distribution. We see that a problem of optimal mass distribution in a finite-dimensional circulatory system (1) in order to increase the critical flutter load, looks promising for explaining the peculiarities of optimizing distributed non-conservative structures. However, it makes sense to tackle first not the most general system (1) but rather a particular finite-dimensional non- conservative system with m degrees of freedom. In this paper, we propose to take an m-link Ziegler pendulum [23–26] as a toy model for an investigation of the optimal mass and stiffness distributions that give an extremum to the critical flutter load. It appears that even the classical two-link Ziegler pendulum has rarely been studied from this point of view, in contrast to its continuous analogue — the Beck column. The only example of such a study known to the author is contained in the book by Gajewski and Zyczkowski [1]. This paper is organized as follows. In the next Section, we first consider optimization of an undamped two-link Ziegler pendulum. We find an explicit expression for the critical flutter load as a function of the mass or stiffness distribution, and demonstrate that in the space of the two mass coefficients and the flutter load as well as in the space of the two stiffness coefficients and the flutter load, the critical flutter load forms a surface with a self-intersection and with the Whitney umbrella singular point. We consider the problem of optimal mass redistribution and find that the only two maxima of the critical flutter load correspond to a vanishing mass either in a joint or at the free end of the pendulum; in the only minimum, the ratio of the masses is equal to the ratio of the stiffness coefficients. The maxima are attained at the singular cuspidal points of the stability domain and are characterized by a dramatic increase in the critical frequency of the vibrations. Then, we write down the equations of motion of an m-link undamped Ziegler pendulum and consider the case m = 3, in which we again find that the optimal mass distributions maximizing the critical flutter load correspond to vanishing of some of the three point masses. Other types of local extrema are also found that correspond to points where the flutter boundary has a vertical tangent, such as cuspidal points with triple eigenvalues, cf. [11, 12], or points where the flutter boundary experiences a sharp turn, cf. [27–29]. Finally we consider the problem of optimal mass distribution for a two-link Ziegler pendulum with dissipation. In conclusion, we formulate an optimization problem for an m-link Ziegler pendulum and discuss some hypotheses on plausible optimal solutions and their expected properties. 2 Structural optimization of the Ziegler pendulum Let us consider the classical Ziegler pendulum consisting of two light and rigid rods of equal length l. The pendulum is attached to a firm basement by a viscoelastic revolute joint with stiffness coefficient c1 and damping 33 Acta Polytechnica Vol. 51 No. 4/2011 Fig. 1: Undamped 2-link Ziegler pendulum. (a) The critical load P(c1, c2) as a function of the stiffness coefficients forms a conical surface in the (c1, c2, P)-space (the casewhen m1 = m2 =1and l =1 is shown); (b)The critical load p(m1, m2) as a function of the masses, forms a self-intersecting surface with the Whitney umbrella singularity at point (0,0,2) of the (m1, m2, p)-space (the case when c1 = c2 =1 is shown) coefficient d1. Another viscoelastic revolute joint with stiffness coefficient c2 and damping coefficient d2 connects the two rods [26, 32]. The point masses m1 and m2 are located at the second revolute joint and at the free end of the second rod, respectively. The second rod is subjected to a tangential follower load P [26, 32]. 2.1 Undamped case Small deviations from the vertical equilibrium for an undamped Ziegler pendulum are described by equation (1) with mass and stiffness matrices that have the following form [26, 32] M = l2 ( m1 + m2 m2 m2 m2 ) , K = ( c1 + c2 − P l P l − c2 −c2 c2 ) , (3) where x = (θ1, θ2) T is the vector consisting of small angle deviations from the vertical equilibrium position. Calculating the characteristic equation det(Mλ2 +K) = 0 for the Ziegler pendulum without dissipation, we find m1m2l 4λ4 + (m1c2 + 4m2c2 + c1m2 − 2P lm2)l2λ2 + c1c2 = 0. (4) By direct calculation of the roots of the characteristic equation (4) or by using the Gallina criterion [24], we find a critical surface that separates the flutter instability and the marginal stability domains (2lm2P − 4m2c2 − m1c1 − m2c1)2 + (m1c1 − m2c2)2 = (m1c1 + m2c2)2. (5) Equation (5) defines a conical surface in the (c1, c2, P )-space when m1,2 and l are fixed: flutter inside the cone, Figure 1(a). However, in the (m1, m2, P )-space the critical surface (5) looks different and has the form of a self-intersecting surface with the Whitney umbrella singularity at the (0, 0, 2c2/l)-point. Indeed, expressing the critical load P from (5) we get P = 4m2c2 + (√ m1c2 ± √ m2c1 )2 2lm2 ≥ 2c2 l . (6) Defining the new critical flutter load as p := P l/c2 we come to the more symmetric expression p = 2 + 1 2 (√ m1 m2 ± √ c1 c2 )2 ≥ 2. (7) 34 Acta Polytechnica Vol. 51 No. 4/2011 Fig. 2: Undamped2-link Ziegler pendulumwith c1 =1, c2 =1: The critical flutter load p(α) as a function of the azimuth angle α indicating the direction in the (m1, m2)-plane. Point A is an absolute minimum of the flutter load: pA = 2, point B corresponds to the local maximum: pB =2+ c1/(2c2) with m1 =0, and the absolute maximum corresponds to point C (not shown) with pC =+∞ and m2 =0. Point Z corresponds to Ziegler’s original design: m1/m2 =2 The case when c1 = c2 = 1, m1 = 2 and m2 = 1 corresponds to the classical result of Ziegler [26] p = 7 2 ± √ 2. (8) The lower value of the critical load corresponds to the boundary between marginal stability and flutter, while the higher critical load corresponds to the conventional transition from flutter to divergence through the double zero eigenvalue with the Jordan block. The critical load (7) as a function of the masses p = p(m1, m2) is plotted in Figure 1(b). It is seen that the stability boundary has a self-intersection along a ray of the p-axis that starts at the Whitney umbrella singularity with the coordinates (0, 0, 2) in the (m1, m2, p)-space. Indeed, for small absolute values of m1c2 − m2c1 we can expand the critical load in a series p = 2 + (m1c2 − m2c1)2 8c1c2m22 + O((m1c2 − m2c1)3), (9) which gives an approximation to the flutter boundary having a canonical form for the Whitney umbrella Z = X2/Y 2. Due to the symmetry of the expression (7), the critical load as a function of the stiffness coefficients, p = p(c1, c2), forms the identical surface in the (c1, c2, p)-space. In the following we will study the function p = p(m1, m2). According to inequality (7), the critical load is always not less than p0 = 2. The minimum is reached when the masses satisfy the constraint m1c2 = m2c1. (10) Note that the equal stiffness coefficients c1 = c2 imply equal masses m1 = m2. This situation corresponds to uniformly distributed mass and stiffness in continuous systems such as the Beck column [19]. In structural optimization problems the uniformly distributed stiffness and mass are usually considered as the initial design that is the starting point in optimization procedures. The critical load of the optimized structure is conventionally compared to that of the same structure with uniform mass and stiffness distributions [4–10,13–18]. Since p(m1, m2) is a ruled surface and thus p effectively depends on the mass ratio only, it is convenient to introduce the azimuth angle α by assuming m1 = cos α and m2 = sin α and to plot the critical load as a function of α. In Figure 2, the curves p = p(α) bound the flutter domain shown in light gray. When α tends to zero, which corresponds to the vanishing mass m2, the critical load increases to infinity. When α tends to π 2 35 Acta Polytechnica Vol. 51 No. 4/2011 and, correspondingly, mass m1 is vanishing, then the critical flutter load increases to the value pB = 2 + 1 2 c1 c2 . (11) At point B in the stability diagram of Figure 2 the flutter boundary has a vertical tangent, which is a typical phenomenon in non-conservative optimization [27–29]. The lower part of the flutter boundary corresponds to designs with a complex conjugate pair of pure imagi- nary double eigenvalues with the Jordan block; the upper part is the designs with two real double eigenvalues of the same magnitude and different sign, each with the Jordan block. Above the upper flutter boundary lies the domain of statical instability or divergence, with two unstable modes corresponding to two different positive real eigenvalues. At point B, the flutter boundary is tangent to the vertical part of the divergence boundary. To the right of this vertical line, there is a pair of pure imaginary eigenvalues and a pair of real eigenvalues with the same magnitude and different signs. The transition from the stability boundary to the divergence boundary below point B occurs when a pair of pure imaginary eigenvalues goes out of the origin in the complex plane to infinity, merges there and returns back along the real axis. This happens because at the boundary m1 = 0, i.e. detM = 0. A similar divergence boundary exists at α = 0, which corresponds to m2 = 0. Transition through the vertical line above point B is accompanied by another eigenvalue bifurcation at infinity: two real eigenvalues of the same magnitude and different signs go out of the origin in the complex plane in order to merge at infinity and then come back along the imaginary axis. The very point B corresponds to an antagonist of a quadruple zero eigenvalue with the Jordan block, i.e. to a quadruplet of complex eigenvalues that merge at infinity in the complex plane. To summarize, the popular initial design corresponding to uniformly distributed mass and stiffness turns out to give an absolute minimum to the critical flutter load of the Ziegler pendulum. The critical flutter load attains its local maximum, pB , for m1 = 0 at the singular cuspidal point B of the stability boundary where the flutter domain has a vertical tangent and touches the boundary of the divergence domain. Note that in [11] a local extremum of the flutter load for the free-free beam carrying a point mass was also found to be at the cuspidal point on the flutter boundary. The global maximum of the critical flutter load for the undamped Ziegler pendulum is at infinity when m2 = 0. The global maximum corresponds to a vanishing mass at the free end of the column, which is qualitatively in agreement with the numerically found optimized designs of the Beck column available in the literature [9,13–18]. Indeed, all known optimized designs of the Beck column are characterized by the vanishing cross-sections at the free end. Moreover, the gradients of the critical flutter load with respect to the mass or stiffness distribution of the Beck column are large, which is, again, in qualitative agreement with our stability diagram of Figure 2. Most interestingly, with the increase in the critical flutter load the higher and higher modes were reported to be involved into the coupling, which indicates the onset of flutter [9, 13–18]. Our simple model shows that this phenomenon seems to be natural for the optimal design that causes the degeneracy in the mass matrix that gives rise to the critical frequency that increases without bounds. 2.2 The m-link Ziegler pendulum It would be very desirable to extend our study to the case of the multiple-degrees-of-freedom Ziegler pendu- lum. The corresponding models and recursive schemes for deriving the equations of motion were proposed by Gasparini et al. [23] and by Lobas [25]. A three-link Ziegler pendulum was considered by Gallina [24]. We take the linearized equations of Lobas [25], since their model deals with the arbitrary masses and stiffnesses in the joints of an m-link Ziegler pendulum, in contrast to the models of Gasparini and Gallina [23,24]. The mass and stiffness matrices of the Ziegler-Lobas model look like M = l2 ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ m∑ i=1 mi m∑ i=2 mi · · · m∑ i=m−1 mi m∑ i=m mi m∑ i=2 mi m∑ i=2 mi · · · m∑ i=m−1 mi m∑ i=m mi · · · · · · · · · · · · · · · m∑ i=m−1 mi m∑ i=m−1 mi · · · m∑ i=m−1 mi m∑ i=m mi m∑ i=m mi m∑ i=m mi · · · m∑ i=m mi m∑ i=m mi ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ , (12) 36 Acta Polytechnica Vol. 51 No. 4/2011 Fig. 3: Undamped 3-link Ziegler pendulum with l = 1, c1 = c2 = c3 = 1: The critical flutter load P(α) as a function of the azimuth angle α indicating the direction in the (m2, m3)-plane at the fixed radial distance r =1 for (a) m1 =0, (b) m1 =10, and (c) m1 =200 and for the fixed m1 =5 and (d) r =1, (e) r =0.65 and (f) r =0.1 K = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ c1 + c2 − P l −c2 · · · 0 P l −c2 c2 + c3 − P l · · · 0 P l · · · · · · · · · · · · · · · 0 0 · · · cm−1 + cm − P l −cm−1 + P l 0 0 · · · −cm−1 cm ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ . (13) For m = 2, matrices (12) and (13) are reduced to (3). Note that detM = l2m m∏ i=1 mi. (14) 37 Acta Polytechnica Vol. 51 No. 4/2011 Let us consider the Ziegler-Lobas pendulum with m = 3 links. Now the mass at the free end is m3, while the masses m2 and m1 are located at the third and second joints, respectively. The first joint connects the first rod with the basement. The length of each of the three rods is equal to l. The stiffness coefficients of the joints are c3, c2, and c1, respectively. For simplicity we assume that c1 = c2 = c3 = c. Then, the characteristic polynomial has the form a0λ 6 + a1λ 4 + a2λ 2 + a3 = 0, (15) with the coefficients a0 = l 6m1m2m3, a1 = cl 4(6m2m3 + 5m1m3 + m1m2) − 2l5P m3(m1 + m2), a2 = 3P 2l4m3 − 2(7m3 + m2)P l3c + (m1 + 5m2 + 14m3)l2c2, a3 = c 3. (16) Composing the discriminant matrix [24] S = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ a0 a1 a2 a3 0 0 0 3a0 2a1 a2 0 0 0 a0 a1 a2 a3 0 0 0 3a0 2a1 a2 0 0 0 a0 a1 a2 a3 0 0 0 3a0 2a1 a2 ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ (17) and calculating the discriminant sequence consisting of the determinants of the three main minors of even order, we find that Δ1 = 3l 12m21m 2 2m 2 3 > 0. The expressions for determinants Δ2 and Δ3 = detS are rather involved, and for this reason we omit them here. However, numerical experiments evidence that the stability boundary is given by the equation Δ3 = 0 for the stability condition Δ3 > 0 implies Δ2 > 0. In Figure 3 using the inequality Δ3 > 0 we present the stability diagrams in the (α, P )-plane, where the azimuth angle α in the (m2, m3)-plane is introduced by assuming m2 = r cos α and m3 = r sin α. We assume equal lengths of the links l = 1 and equal stiffness coefficients c1 = c2 = c3 = 1 and vary the radial distance r in the (m2, m3)-plane and the mass m1. Since for m = 3 the critical surface P (m2, m3) is no longer a ruled surface as it was in the case m = 2, the pictures in the (α, P )-plane change with the variation of the radial distance r, which complicates the optimization problem. Nevertheless, such diagrams are convenient for analyzing the geometry of the stability boundary and thus for identifying the potential extrema. Moreover, the critical surface P (m2, m3) has a self-intersection along a ray of the P -axis that starts from the singularity Whitney umbrella, as in the case of the 2-link pendulum. Therefore, at small values of m2 and m3 the critical load can be locally approximated by a ruled surface. In the left column of Figure 3 the radial distance r in the (m2, m3)-plane is fixed to r = 1 while the mass m1 is increasing. As in the case m = 2, (marginal) stability is possible for α ∈ [0, π/2]. For m1 = 0, two finite maximal values PA and PB are identified at α = 0 (m3 = 0) and α = π/2 (m2 = 0), respectively, Figure 3(a). Both maxima are attained at the cuspidal points of the stability boundary, where it has vertical tangents. However, the stability diagram changes when m1 = 10, Figure 3(b). Again, local extrema exist at the boundary points α = 0 (m3 = 0) and α = π/2 (m2 = 0), while the global maximum is at the point of the sharp turn of the flutter boundary with the vertical tangent near the cuspidal point C1 � (0.040 347 7, 11.961 144), which corresponds to triple pure imaginary eigenvalues λ � ±i1.163 524 3 with the Jordan block of the third order, cf. [11, 12] where at such a singular point a maximum of the critical flutter load was found for a free-free beam under the follower thrust. With a further increase in the first mass up to m1 = 200, the stability diagram converges to that similar to the diagram of the two-link pendulum, cf. Figure 2 and Figure 3(c). This is not surprising because big inertia of the first joint makes the three-link pendulum effectively a two-link pendulum. The right column in Figure 3 corresponds to the fixed first mass m1 = 5 and varied radial distance r in the (m2, m3)-plane. Small values of r correspond effectively to a two-link pendulum. That is why Figure 3(f) with r = 0.1 looks similar to Figure 3(c) and Figure 2. An increase in r is accompanied by a complication of the stability diagram. In particular, two cusp point singularities originate corresponding to triple pure imaginary eigenvalues with the Jordan block, Figure 3(d,e). 38 Acta Polytechnica Vol. 51 No. 4/2011 Near these singularities, the stability boundary experiences a sharp turn at points B1 and B2, where the tangent to the boundary is vertical. At such points the eigencurves Imλ(P ) form a crossing that can change either to avoided crossing or to overlapping with the origination of a bubble of complex eigenvalues in dependence on the direction of variation of azimuth angle α. This phenomenon has been observed in many numerical studies of non-conservative optimization problems [4, 6, 7, 9, 10, 13, 14, 17, 18] and was described analytically by Kirillov and Seyranian in [27, 28]. Moreover, the critical load at these points can jump to a higher value corresponding to the merging of other (often higher) modes, and is thus undefined [4, 6, 7]. Nevertheless, these points could be local extrema of the merit functional, see [29] where the necessary conditions for that were derived. We stress that due to the finite number of degrees of freedom and the finite number of control parameters the stability boundary of an m-link Ziegler pendulum can be thoroughly analyzed both analytically and numerically. In particular, the coordinates of the singular points can be calculated with arbitrary precision and thus the issues of high sensitivity of the critical flutter load at the optima could be successfully resolved in this model, unlike the optimization problems for distributed systems. The possibility to work with singularities related to coalescence of more than two eigenvalues allows us to make a qualitative investigation of the question ‘Should low-order models be believed?’ [31]. Although two-mode approximations work well far from such singularities, in their vicinity we have to take into account higher modes. 2.3 Damped case In the presence of dissipation, the equation of the Ziegler pendulum is Mẍ + Dẋ + Kx = 0 (18) with the damping matrix [25, 26] D = ⎛ ⎜⎜⎜⎜⎜⎜⎜⎝ d1 + d2 −d2 · · · 0 0 −d2 d2 + d3 · · · 0 0 · · · · · · · · · · · · · · · 0 0 · · · dm−1 + dm −dm−1 0 0 · · · −dm−1 dm ⎞ ⎟⎟⎟⎟⎟⎟⎟⎠ . (19) For the two-link damped Ziegler pendulum with m = 2, we find the characteristic equation a0λ 4 + a1λ 3 + a2λ 2 + a3λ + a4 = 0. (20) with the coefficients a0 = l 4m1m2, a1 = l 2(m1d2 + d1m2 + 4m2d2), a2 = d1d2 + m1l 2c2 + 4m2l 2c2 + c1m2l 2 − 2P l3m2, a3 = d1c2 + c1d2, a4 = c1c2. (21) Applying the Routh-Hurwitz criterion, we find the critical flutter load P = 4m22(d 2 2c 2 1 + d 2 1c 2 2) + d1d2(8m2(m1 + 2m2)c 2 2 + (c1m2 − m1c2)2) 2(m2l(4m2d2 + d1m2 + m1d2)(c1d2 + d1c2) + 1 2 d1d2 m2l3 . (22) For c1 = c2 = 1, l = 1, m1 = 2 and m2 = 1 it was found to be [34] P = 4d21 + 33d1d2 + 4d 2 2 2(6d2 + d1)(d2 + d1) + 1 2 d1d2. (23) Equation (23) defines a surface with the Whitney umbrella singularity in the (d1, d2, P )-space which explains Ziegler’s destabilization paradox by vanishing dissipation [26], as was first demonstrated by Bottema already in 1956 [32, 35]. In contrast to earlier studies, e.g. [32, 34–36], we consider the critical flutter load (22) as a function of the masses P = P (m1, m2) for the fixed damping distribution. 39 Acta Polytechnica Vol. 51 No. 4/2011 Fig. 4: Damped 2-link Ziegler pendulum with l = 1, c1 = c2 = 1: The critical flutter load P(α) as a function of the azimuth angle α indicating the direction in the (m1, m2)-plane for (a) d1 = d2 = 1 and r = 1, (b) d1 = d2 = 1 and r =0.4, (c) d1 =1, d2 =0.1 and r =1, (d) d1 =1, d2 =0.1 and r =0.1 In Figure 4, the stability diagrams for the damped two-link Ziegler pendulum are shown in the assumption of c1 = c2 = 1, l = 1, m1 = r cos α and m2 = r sin α. The critical load P (m1, m2) does not constitute a ruled surface and thus the function P (α) depends on the radial distance r in the (m1, m2)-plane. We see that the extrema again correspond to the boundary points m1 = 0 and m2 = 0, although the singularity at α = π/2 is an intersection point. Nevertheless, with an increase in the number of degrees of freedom and control parameters new types of singularities will appear. The planar stability diagrams of Figure 4 depend on the damping distribution but do not tend to that of the undamped pendulum when damping goes to zero (destabilization paradox [32, 35, 36]). 2.4 ‘Problema novum ad cuius solutionem Mathematici invitantur’ ‘A new problem that mathematicians are invited to solve’ is a translation of the Latin title of the work by J. Bernoulli published in 1696 where he proposed the famous Brachistochrone problem [37]. Supporting this good old tradition, we would like to formulate the following optimization problem: Given a circulatory system (1) with matrices M and K defined as in (12) and (13). Find all local ex- trema, the absolute maximum of the critical flutter load P , and the corresponding extremal mass distributions {m1, m2, . . . , mm}. We can also consider the problem of optimal stiffness distribution or even vary both stiffness and mass. The same problems can be formulated for the damped pendulum with the damping matrix (19). We expect that both in the undamped case and in the damped case there exists a class of extrema corre- sponding to the distributions {m1, m2, . . . , mm} with some masses mi = 0. It should be possible to find these 40 Acta Polytechnica Vol. 51 No. 4/2011 optimal mass distributions explicitly, perhaps with the use of the Pontryagin’s maximum principle. It would be interesting to identify the singularities of the stability boundary that correspond to these extrema; some of them could be related to infinite eigenvalues λ. On the other hand, some local extrema should exist with mass distributions that do not contain vanishing masses mi. It would be interesting to understand at which points — smooth or non-smooth — of the stability boundaries they are attained. Since the system is finite-dimensional and contains a finite number of control parameters with a clear physical meaning, the locations of the singularities corresponding to multiple eigenvalues can easily be found numerically with high accuracy. In the vicinity of such points where at least three pure imaginary eigenvalues couple, the question ‘Should low-order models be believed’ [31] makes sense, because here one more degree of freedom is crucial for the correct solution. Knowledge of rigorously established optimal solutions at every m should shed light on the behavior of the optimal mass distributions and critical frequencies with an increase in the number of degrees of freedom. As a by-product, such a study will give an insight into the problem of dimension reduction, and will serve as a nice playground both for applying the existing methods of nonsmooth analysis and eigenvalue optimization [38– 40] and for developing them further in the direction of a tighter relation both with singularity theory and with the needs of applications. We expect that the proposed model optimization problem will yield useful recommendations for improving the algorithms for optimizing real non-conservative structures in industry. References [1] Gajewski, A., Zyczkowski, M.: Optimal structural design under stability constraints. Dordrecht : Kluwer, 1988, in particular, p. 137. [2] Zyczkowski, M. (ed.): Structural optimization under stability and vibration constraints. Wien–New York : Springer-Verlag, 1989. [3] Seyranian, A. P., Lund, E., Olhoff, N.: Multiple eigenvalues in structural optimization problems. Journal of Structural and Multidisciplinary Optimization. 8 (1994), 207–227. [4] Claudon, J. L.: Characteristic curves and optimum design of two structures subjected to circulatory loads, Journal de Mecanique. 14(3) (1975), 531–543. [5] Sundararajan, C.: Optimization of a nonconservative elastic system with stability constraint. Journal of Optimization Theory and Applications. 16(3/4) (1975), 355–378. [6] Hanaoka, M., Washizu, K.: Optimum design of Beck’s column. Computers and Structures. 11(6) (1980) 473–480. [7] Kounadis, A. N., Katsikadelis, J. T.: On the discontinuity of the flutter load for various types of cantilevers, International Journal of Solids and Structures. 16 (1980), 375–383. [8] Park, Y. P., Mote, C. D.: The maximum controlled follower force on a free-free beam carrying a concentrated mass, Journal of Sound and Vibration. 98 (1985), 247–256. [9] Ringertz, U.: On the design of Beck’s column, Structural Optimization. 8(2–3) (1994), 120–124. [10] Kirillov, O. N., Seyranian, A. P.: Optimization of stability of a flexible missile under follower thrust, AIAA Paper 98-4969, (1998), 2 063–2 073. [11] Kirillov, O. N.: Optimization of stability of the flying bar. Young Scientists Bulletin. Applied Mathematics and Mechanics 1(1) (1999), 64–78. [12] Kirillov, O. N., Seyranian, A. P.: Optimization of stability of a flying column. 3rd World Congress of Structural and Multidisciplinary Optimization. Buffalo. New York (USA) : May 17–21, 1999. Short paper proceedings, 2 (1999) 355–357. [13] Langthjem, M. A., Sugiyama, Y.: Optimum shape design against flutter of a cantilevered column with an end-mass of finite size subjected to a non-conservative load, Journal of Sound and Vibration. 226(1) (1999), 1–23. 41 Acta Polytechnica Vol. 51 No. 4/2011 [14] Langthjem, M. A., Sugiyama, Y.: Optimum design of cantilevered columns under the combined action of conservative and nonconservative loads Part I: The undamped case. Computers and Structures. 74(4) (2000), 385–398. [15] Langthjem, M. A., Sugiyama, Y.: Optimum design of cantilevered columns under the combined action of conservative and nonconservative loads Part II: The damped case. Computers and Structures. 74(4) (2000), 399–408. [16] Langthjem, M. A., Sugiyama, Y.: Dynamic stability of columns subjected to follower loads: a survey. Journal of Sound and Vibration. 238 (2000), 809–851. [17] Temis, Yu. M., Fedorov, I. M.: Shape optimization of nonconservatively loaded beams with a stability criterion. Problems of Strength and Plasticity. 69 (2007), 24–37. [18] Katsikadelis, J. T., Tsiatas, G. C.: Optimum design of structures subjected to follower forces. International Journal of Mechanical Sciences. 49 (2007), 1 204–1 212. [19] Beck, M.: Die Knicklast des einseitig eingespannten, tangential gedruckten Stabes (The buckling load of the cantilevered, tangentially loaded rod). Zeitschrift fur angewandte Mathematik und Mechanik. 3 (1952), 225–228. [20] Sugiyama, Y., Langthjem, M. A., Ryu, B.-J.: Realistic follower forces. Journal of Sound and Vibration. 225 (1999), 779–782. [21] Sugiyama, Y., Langthjem, M. A., Ryu, B.-J.: Beck’s column as the ugly duckling. Journal of Sound and Vibration. 254 (2002), 407–410. [22] Elishakoff, : Controversy associated with the so-called “follower forces”: critical overview. Applied Mechan- ics Reviews. 58 (2005), 117–142. [23] Gasparini, A. M., Saetta, A. V., Vitaliani, R. V.: On the stability and instability regions of non-conservative continuous system under partially follower forces. Computer Methods in Applied Mechanics and Engineer- ing. 124 (1995), 63–78. [24] Gallina, P.: About the stability of non-conservative undamped systems, Journal of Sound and Vibration. 262(4) (2003), 977–988. [25] Lobas, L. G.: Dynamic behavior of multilink pendulums under follower forces. International Applied Me- chanics. 41(6) (2005), 587–613. [26] Ziegler, H.: Die Stabilitätskriterien der Elastomechanik (Stability criteria in elasticity theory). Ingenieur- Archiv. 20 (1952), 49–56. [27] Kirillov, O. N., Seyranian, A. P.: Overlapping of frequency curves in non-conservative systems. Doklady Physics. 46(3) (2001), 184–189. [28] Kirillov, O. N., Seyranian, A. P.: Metamorphoses of characteristic curves in circulatory systems. Journal of Applied Mathematics and Mechanics. 66(3) (2002), 371–385. [29] Kirillov, O. N., Seyranian, A. P.: A non-smooth optimization problem. Moscow University Mechanics Bulletin. 57(3) (2002), 1–6. [30] Bolotin, V. V.: Non-conservative problems of the theory of elastic stability. Moscow : Fizmatgiz (in Rus- sian), 1961; Oxford : Pergamon, 1963. [31] Butlin, T., Woodhouse, J.: Friction-induced vibration: Should low-order models be believed? Journal of Sound and Vibration. 328 (2009), 92–108. [32] Kirillov, O. N., Verhulst, F.: Paradoxes of dissipation-induced destabilization or who opened Whitney’s umbrella? Zeitschrift fur angewandte Mathematik und Mechanik. 90(6) (2010), 462–488. [33] Wang, G., Lin, Y.: A new extension of Leverrier’s algorithm. Linear Algebra and Applications. 180 (1993), 227–238. 42 Acta Polytechnica Vol. 51 No. 4/2011 [34] Herrmann, G., Jong, I. C.: On the destabilizing effect of damping in nonconservative elastic systems. Transactions of the ASME, Journal of Applied Mechanics. 32(3) (1965), 592–597. [35] Bottema, O.: The Routh-Hurwitz condition for the biquadratic equation. Indagationes Mathematicae. 18 (1956), 403–406. [36] Kirillov, O. N.: Destabilization paradox. Doklady Physics. 49 (2004), 239–245. [37] Bernoulli, J.: Problema novum ad cuius solutionem Mathematici invitantur (A new problem that mathe- maticians are invited to solve). Acta Eruditorum. 15 (1696), 264–269. [38] Lewis, A. S.: The mathematics of eigenvalue optimization. Mathematical Programming, Ser. B 97 (2003), 155–176. [39] Burke, J. V., Henrion, D., Lewis, A. S., Overton, M. L.: Stabilization via nonsmooth, nonconvex optimiza- tion. IEEE Transations on Automatic Control. 51(11) (2006), 1 760–1 769. [40] Ledyaev, Yu. S., Zhu, Q. J.: Nonsmooth analysis on smooth manifolds. Transactions of the American Mathematical Society. 359(8) (2007), 3 687–3 732. Oleg N. Kirillov E-mail: o.kirillov@hzdr.de Magneto-Hydrodynamics Division (FWSH) Helmholtz-Zentrum Dresden-Rossendorf P.O. Box 510119, D-01314 Dresden, Germany 43