Finding the strongest stable weightless column with a follower load and relocatable concentrated masses
aa r X i v : . [ phy s i c s . c l a ss - ph ] A ug FINDING THE STRONGEST STABLE WEIGHTLESSCOLUMN WITH A FOLLOWER LOAD ANDRELOCATABLE CONCENTRATED MASSES by Oleg N. Kirillov and
Michael L. Overton(
Northumbria University, Newcastle upon Tyne, NE1 8ST, UK.Email: [email protected] )( Courant Institute of Mathematical Sciences, New York University,New York, NY, 10012, USA. Email: [email protected] ) [August 2020] Summary
We consider a problem of optimal placement of concentrated masses along a weightlesselastic column (rod) that is clamped at one end and loaded by a nonconservativefollower force at the free end. The goal is to find the largest possible interval suchthat the variation in the loading parameter within this interval preserves stability ofthe structure. When the number of masses is two we find the optimal load analytically.For the case of three or more masses, we obtain conjectures for the globally optimalload that are strongly supported by extensive computational results. To obtain these,we use techniques from nonsmooth optimization, employing the recently developedopen-source software package granso (GRadient-based Algorithm for Non-SmoothOptimization) to maximize the load subject to appropriate nonsmooth stabilityconstraints.
1. Introduction
The Beck column is an elastic Euler-Bernoulli beam clamped at one end and loaded at thetip by a follower force ( , ). The follower force is defined as a force with the line of actionthat always coincides with the tangent line to the neutral axis of the deformed beam at itsfree end, much like a rocket thrust ( ). The follower force does not depend on the velocityof the beam. However, it cannot be derived from a potential: the work done by the followerforce along a closed contour is non-zero ( , ). A straight form of the beam is in a stableequilibrium when the follower force is absent or relatively small. Nevertheless, at somesufficiently large value the follower force excites exponentially growing oscillations of thebeam that are known as the flutter instability ( ). It is well-established that the uniformBeck column is stable against flutter if the follower force, P , is such that 0 ≤ P . . EIl ,where E is Young’s modulus, I is the moment of inertia of a cross-section of the columnand l is the the length of the column ( , , , ).Flutter is an oscillatory instability that is critically important both for safety ofengineering structures interacting with fluid flows and for efficiency of energy harvestingdevices that are based on the fluid-structure interactions. Recent years have seen anincreasing interest in the Beck column in the modelling of biological filaments and their Q. Jl Mech. Appl. Math. (????) ?? (??), 1–27 c (cid:13) Oxford University Press ????
Kirillov and Overton
Stability p P b Fig. 1
The Pfl¨uger column and its stability diagram. The ratio of the end mass to the mass ofthe column, µ = M/ ( ml ), is parameterized by µ = tan β . The Beck column corresponds to thevanishing end mass ( M = 0, so β = 0) and the weightless Pfl¨uger column (or Dzhanelidze’s column( )) to the vanishing mass of the rod ( m = 0, so β = π/ artificial biomimetic analogues, i.e., hair-like slender microscale structures that play animportant part in such biological processes as swimming, pumping, mixing, and cytoplasmicstreaming by performing rhythmic, wave-like motion that usually sets in via flutterinstability ( , , , ).Structures loaded by follower forces have long been questioned for their practicalrealization ( , ), despite an evident example given by flexible missiles ( , , ). Inthe 1970-90s, Sugiyama et al. used solid rocket motors to demonstrate flutter of cantileversunder a follower thrust on relatively short (several seconds) time intervals ( , , ). Amechanism recently invented by Bigoni and Noselli produces a frictional follower force ( )and enables experimental realization of fluttering cantilevered rods under follower loadson virtually infinite time intervals ( , ). These practical realizations differ from theclassical Beck column, however, by the presence of a finite-size loading unit at the tip of thecantilever and therefore are better described by the model of the Pfl¨uger column ( , ),which is the Beck column with a point mass at the loaded end; see the left panel of Fig. 1.To model the Pfl¨uger column we consider an elastic beam of length l , with Young’smodulus E and mass per unit length m , clamped at one end and loaded by a tangentialfollower force P at the other end, where a point mass M is mounted. The moment of inertiaof a cross-section of the column is denoted by I . Small lateral vibrations of the Pfl¨ugercolumn near the undeformed equilibrium are described by the linear partial differentialequation ( , ) EI ∂ y∂s + P ∂ y∂s + m ∂ y∂t = 0 (1.1)where y ( s, t ), is the amplitude of the vibrations and s ∈ [0 , l ] is a coordinate along the he strongest column with a follower load and relocatable masses s = 0) equation (1.1) satisfies the boundary conditions y = 0 , ∂y∂s = 0 , s = 0 , (1.2)while at the loaded end ( s = l ), the boundary conditions are EI ∂ y∂s = 0 , EI ∂ y∂s = M ∂ y∂t , s = l. (1.3)Introducing the dimensionless quantities ξ = sl , τ = tl r EIm , p = P l EI , µ = Mml , (1.4)and separating the time variable through y ( ξ, τ ) = lf ( ξ ) exp( λτ ), we obtain thedimensionless boundary eigenvalue problem ∂ ξ f + p∂ ξ f + λ f = 0 , (1.5) ∂ ξ f (1) = 0 , ∂ ξ f (1) = µλ f (1) ,f (0) = 0 , ∂ ξ f (0) = 0 (1.6)defined on the interval ξ ∈ [0 , , ) f ( ξ ) = A (cosh( g ξ ) − cos( g ξ )) + B ( g sinh( g ξ ) − g sin( g ξ )) (1.7)with g , = s p p − λ ± p , where the subscripts 1 and 2 correspond to the signs + and − , respectively. Imposing theboundary conditions (1.6) on the solution (1.7) yields the characteristic equation ∆( λ ) = 0for the determination of the eigenvalues λ , where∆( λ ) = ∆ − ∆ µλ and ∆ = g g ( g + g + 2 g g cosh g cos g + g g ( g − g ) sinh g sin g )∆ = ( g + g )( g sinh g cos g − g cosh g sin g ) . (1.8)Parameterizing the mass ratio in (1.4) by µ = tan β with β ∈ [0 , π/
2] enables theexploration of all possible ratios µ = M/ ( ml ) of the end mass to the mass of the columnfrom zero ( β = 0) to infinity ( β = π/ ). Kirillov and Overton
It is well known that the Beck column loses its stability at p ≈ .
05 ( ). In contrast,the Dzhanelidze column becomes unstable at p ≈ .
19, which is the smallest positive rootof the equation ( ) tan √ p = √ p. (1.9)These values, representing two extreme situations, are connected by a marginal stabilitycurve in the ( β, p )-plane ( , , , , , ); see the right panel of Fig. 1.For every fixed value β ∈ [0 , π/ p causes the imaginary eigenvalues of two different modes to approach each otherand merge into a double imaginary eigenvalue with one eigenfunction. When p crosses thethreshold, the double eigenvalue splits into two complex eigenvalues, one with positive realpart, which determines a flutter-unstable mode.At β = π/ p & .
19; see ( , , ).Structural optimization of the Beck column against flutter and divergence instabilitiesis usually formulated as a problem on a redistribution of the material of the column of agiven density under an isoperimetric constraint fixing the volume of the column in order tomaximize the range of variation of the follower load corresponding to the stable structure.In the literature many specific numerically optimized shapes of the Beck column have beenreported ( , , , , , ) with the maximal critical load reaching the values of p ≈ .
00 ( ), p ≈ .
30 ( ), p ≈ .
59 ( ), and p ≈ .
62 ( ), which significantlyimprove upon the critical load p ≈ .
05 of the uniform column with the constant cross-section. Nevertheless, none of these designs is proven to be a global or even a local optimizer.Such a proof would be difficult to obtain because the problem of structural optimizationof the critical flutter load for the elastic Beck column is both nonconvex and nonsmooth( , ).Indeed, the elastic Beck column is a time-reversible dynamical system in which thetransition from stability to flutter instability generically happens via the reversible-Hopfbifurcation, i.e., through the formation of a double imaginary eigenvalue with a Jordan blockat the stability boundary and its subsequent splitting when parameters enter the instabilityregion ( ). Codimension-1 parts of the stability boundary are thus smooth hypersurfacescorresponding to double imaginary eigenvalues with a Jordan block (provided that theremaining eigenvalues are simple and imaginary) ( , , ). These hypersurfaces can meeteach other at sets of higher codimension such as intersections, cuspidal edges and points,conical points etc.; see ( ) for a full classification of generic singularities on the stabilityboundary of mechanical systems with non-potential positional forces. The unavoidablesingularities linked to multiple eigenvalues is the main reason for nonsmoothness of themerit functionals in the optimization of such systems, including the Beck column, withrespect to stability criteria ( , ).Many studies report on the phenomenon of overlapping of eigenvalue curves thataccompanies the process of optimization of the Beck column. The eigenfrequencies plottedas functions of the load exhibit sudden crossings during the optimization that lead totransfer of instability between modes and to a discontinuous change in the merit functional( , , , , , , , , , , ). The high sensitivity of the optimized designto variation of parameters is caused by the nonconvexity of the stability domain ( , ). he strongest column with a follower load and relocatable masses , , ).All of the phenomena described above were also observed in simplified settings with theuniform Beck or Pfl¨uger column carrying relocatable lumped masses ( , , , , , , , ). Nevertheless, to the best of our knowledge, no rigorously proven local or globaloptimal solutions or credible numerical guesses exist in the literature even in the problemsof optimal localization of point masses along elastic beams loaded by the follower force.In recent mechanical laboratory experiments with follower forces ( , ), the ratio of theend mass to the mass of the column, µ = M/ ( ml ), was chosen to be very large, approachingthe Dzhanelidze limit µ = ∞ corresponding to a weightless column. The instabilitythresholds obtained in these experiments were in a very good agreement with the theoreticalpredictions based on the Pfl¨uger model. In the Dzhanelidze limit, the mathematical modelis reduced to a system of ordinary differential equations ( , , , ). The works( , , ) considered stability of a weightless Pfl¨uger column with an additional relocatablemass. A recent work ( ) corrected some of the results reported in ( ) and proposedextending the model to incorporate several relocatable masses.The primary purpose of this paper is to study the problem of the weightless Pfl¨uger(Dzhanelidze) column with discrete relocatable masses in depth. We formulate the problemof optimal distribution of the masses along the column with the goal of maximizing theinterval of stability for the follower load parameter. There are two key advantages ofthe discrete mass model over the classical Beck column: first, an analytical treatment ispossible, and second, numerical optimization does not require Galerkin or finite elementsdiscretization, as is typically necessary for continuous nonconservative elastic systems ( ).For the case of two masses we give an analytical derivation of the optimal load. For thecase of three or more masses, we give conjectures for the optimal load that are stronglysupported by extensive computational results. To obtain these, we use techniques fromnonsmooth optimization, employing a recently developed open-source software package, granso (GRadient-based Algorithm for Non-Smooth Optimization) ( , ), to maximizethe load subject to appropriate nonsmooth stability constraints. A second purpose of ourpaper is to introduce these nonsmooth optimization techniques to a broad audience, withthe hopes that they will be useful for other stability optimization problems arising in civiland mechanical engineering.
2. A weightless elastic column with n concentrated masses It is convenient to first consider the simple model of the Pfl¨uger column without relocatablemasses, with zero mass per unit length ( m = 0) and vanishing end mass ( M = 0). Then,the boundary value problem (1.5), (1.6) takes the form ∂ ξ f + κ ∂ ξ f = 0 , (2.1) f (0) = 0 , ∂ ξ f (0) = 0 , ∂ ξ f (1) = 0 , ∂ ξ f (1) = 0 , (2.2)where κ = p. (2.3)Following ( , , ), consider the case when a concentrated constant force F is acting Kirillov and Overton P n i 1s i v i Fig. 2
The weightless ( m = 0) Beck column loaded by the follower force P with n concentratedmasses ( M , M , . . . , M i , . . . , M n ) attached ( , ). in a direction perpendicular to the non-deformed column at the point s = αl . Introducingthe dimensionless version of the force parameter, φ = F l EI , we seek the general solution tothe equation (2.1) in the form ( , , ) f ( ξ ) = u ( ξ ) + , ξ ∈ [0 , α ) v ( ξ ) , ξ ∈ [ α,
1] (2.4)where u ( ξ ) = A sin κξ + B cos κξ + Cξ + D and v ( ξ ) = A sin κξ + B cos κξ + C ξ + D . To determine the coefficients A , B , C , and D , we require that u ( α ) = f ( α ) , ∂ ξ u ( α ) = ∂ ξ f ( α ) ,∂ ξ u ( α ) = ∂ ξ f ( α ) , ∂ ξ f ( α ) − ∂ ξ u ( α ) = φ. (2.5)This yields v ( ξ ) = ( ξ − α ) κ − sin(( ξ − α ) κ ) κ φ. (2.6)Taking (2.6) into account in the general solution (2.4) and then substituting f ( ξ ) into theboundary conditions (2.2), we find the coefficients A , B , C , and D to obtain u ( ξ ) = sin( κα ) − ξκ cos( κα ) + sin(( ξ − α ) κ ) κ φ. (2.7)Let us now assume that the weightless cantilevered column loaded by the follower force atits free end carries n concentrated masses such that the mass M n > M i ≥ i = 1 , . . . , n − he strongest column with a follower load and relocatable masses s i < l from the clamped end of the rod. Let v i be a transversal displacement of the mass M i from the equilibrium configuration, as shown in Fig. 2. Introducing the dimensionlessdisplacements of the masses, w i , the distances, α i , and the mass ratios, µ i , as w i = v i l , α i = s i l , µ i = M i M n , i = 1 , . . . , n, (2.8)we write the equations of motion of the masses ( , , ) w = − γ µ d w dτ − γ µ d w dτ − . . . − γ n µ n d w n dτ ,w = − γ µ d w dτ − γ µ d w dτ − . . . − γ n µ n d w n dτ , ... w i = − γ i µ d w dτ − γ i µ d w dτ − . . . − γ in µ n d w n dτ , ... w n = − γ n µ d w dτ − γ n µ d w dτ − . . . − γ nn µ n d w n dτ , (2.9)where the dimensionless time τ is defined now as τ = t r EIM n l . (2.10)Note that α α . . . α n = 1 and µ n = 1 . The coefficient γ ij is the displacement ofthe mass µ i as a result of application to the column of a unit force φ = 1 at the point α j .According to (2.4) with the functions (2.6) and (2.7) the coefficient γ ij is given by δ ij /κ ,where δ ij = sin( κα j ) − α i κ cos( κα j ) + sin(( α i − α j ) κ )+ , i j ( α i − α j ) κ − sin(( α i − α j ) κ ) , i > j. (2.11)Separating time with the ansatz w i = u i e σκ / τ we arrive at the quadratic eigenvalueproblem ( M σ + K ) u = 0 , (2.12)where u = ( u , u , . . . , u n ), K is the n × n unit matrix, and M = µ δ µ δ · · · µ n δ n µ δ µ δ · · · µ n δ n ... ... . . . ... µ δ n µ δ n · · · µ n δ nn , (2.13) Kirillov and Overton where, as already noted, µ n = 1. The eigenvalues σ k are given by σ k = ± q − λ − k (2.14)where the λ k are the eigenvalues of the matrix M .The trivial equilibrium of the circulatory system (2.9) is stable if and only if theeigenvalues σ k are imaginary and semisimple, or equivalently, the λ k are real, positive andsemisimple. Cases with a multiple imaginary eigenvalue σ k with a Jordan block lie on theboundary between the stability and flutter domains. In the generic case the crossing of thisstability boundary is accompanied by merging of two simple imaginary eigenvalues into adouble imaginary eigenvalue with a Jordan block, indicating the onset of the reversible-Hopfbifurcation or flutter ( , , , ). Non-oscillatory instability or divergence correspondsto one or more positive real eigenvalues σ k and in this model it generically sets in when twoconjugate simple imaginary eigenvalues meet at infinity, split and turn back towards theorigin along the real axis in the complex plane ( , , ).Summarizing, for a given number of masses n , the quadratic eigenvalue problem (2.12)is defined by (2.11) and (2.13), which depend on the dimensionless parameters α i and µ i , i = 1 , . . . , n −
1, defined in (2.8), as α n = µ n = 1, as well as the given load κ . It is convenientto use the parameterization µ i = tan β i , i = 1 , . . . , n − , β i ∈ [0 , π/ . (2.15)There is a slight abuse of notation here, since when β i = π/ µ i is infinite, M n = 0 and M is not defined. However, it is convenient to allow β i to take the value π/
2, and it willalways be understood that in that case the appropriate limit must be taken when referringto the matrix M and the eigenvalues σ k or λ k . Given α i , β i , i = 1 , . . . , n , let us define κ α,β crit as the largest value such that the eigenvalues σ k (which depend on α i , β i and κ ) areimaginary for all κ ∈ [0 , κ α,β crit ]. Our goal is to find the supremum of κ α,β crit over all parameters α i ∈ [0 ,
1] and β i ∈ [0 , π/ i = 1 , . . . , n −
1. We begin with the case n = 2, where we findan analytical solution.
3. Analytical derivation of the optimal load for the weightless column carryingtwo masses
When n = 2, the weightless column carries a relocatable mass M between the clamped endand the free end of the rod with mass M fixed at the free end. There are two parameters, α and β . Expression (2.11) allows us to find the coefficients δ ij in the explicit form, cf.( , ), δ = sin( κα ) − κα cos( κα ) δ = sin( κ ) − κα cos( κ ) − sin( κ (1 − α )) δ = sin( κα ) − κ cos( κα ) + κ (1 − α ) δ = sin( κ ) − κ cos( κ ) . (3.1)As we will see, already in this simplest possible mechanical system, the subdivision of theparameter space into the domains of stability, flutter instability, and divergence instabilityis highly nontrivial. However, we will be able to explore it completely and find an apparent he strongest column with a follower load and relocatable masses , κ α ,β crit ] inthe space of parameters α ∈ [0 , , β ∈ [0 , π/ p ( σ ) = det( M σ + K ) can be obtained with the use of the Gallina criterion ( , , )that is based on the investigation of the discriminant of the polynomial. For n = 2, p ( σ ) isa biquadratic function p ( σ ) = σ tan β { κ ( α − κ − κα cos κ + sin( κα − κ )) − sin( κ ( α − κα ) − κ cos( κα ) − κ ( α − } + σ [tan β (sin( κα ) − κα cos( κα )) − κ cos κ + sin κ ] + 1 . (3.2)Notice that the coefficient at the leading power of σ in the polynomial (3.2) is nothingelse but det M ; see ( ). The system loses stability by divergence as soon as det M = 0,which yields the following equation determining the divergence boundary:sin κ − κα cos κ + sin( κα − κ )sin( κα ) − κ cos( κα ) − κ ( α −
1) = sin( κ ( α − κ ( α − . (3.3)Note that this equation is independent of β . The right panel of Fig. 3 shows the divergenceboundary (3.3) in the ( α , β , κ )-space.The roots of the characteristic polynomial (3.2) are double imaginary if the discriminantof the biquadratic function vanishes:(sin( κα ) − κα cos( κα )) (tan β ) + 2 α κ tan β cos κ [cos( κα ) + 2( α − β sin( κα ) [2 sin( κ ( α − κ ] − κ tan β [7 sin( κ ( α − α −
1) + sin( κ ( α + 1))( α + 1)] − κ tan β [(2 α −
3) sin κ + sin( κ (2 α − κ − κ cos κ ) = 0 . (3.4)For this reason ( , , ) equation (3.4) determines the boundary of the flutter domainthat is shown in the left panel of Fig. 3.For a given ( α , β ), the critical value of the load parameter is given by κ α ,β crit = min { κ : ( κ, α , β ) satisfies either (3 .
3) or (3 . } , as this is the length of the longest vertical line segment rising from the point ( α , β ,
0) thatdoes not enter either the flutter or the divergence domain. Consequently, the quantity κ crit = sup { κ α ,β crit : α ∈ [0 , , β ∈ [0 , π/ } is the critical load for the strongest possible choice of stable column. Note that although thedivergence boundary (3.3) is smooth, the boundary of the flutter domain (3.4) is nonsmooth.Fig. 4 shows cross-sections of the flutter boundary and the divergence boundary in the( α , κ )- and ( β , κ )-planes. In the left panel, for which β is fixed to ˆ β = 1 . α , κ )-plane at α = ˆ α =0 Kirillov and Overton
Fig. 3
The case of n = 2 concentrated masses. (Left) the flutter domain is a finite solid set inthe ( α , β , κ ) space, enclosed within the singular surface defined by (3.4). (Right) The divergencedomain lies above the boundary set defined by (3.3). For a given ( α , β ), the critical value ofthe load parameter, κ α ,β crit , is the minimal value of κ that satisfies either (3.4) or (3.3), as this isthe length of the longest vertical line segment rising from the point ( α , β ,
0) that does not entereither the flutter or divergence domain. Consequently, this is the largest value ˜ κ such that thecolumn is stable for all κ ∈ [0 , ˜ κ ). The strongest stable column is obtained by finding ( α , β ) withthe largest possible critical load κ α ,β crit . . κ = 5 . α is fixed to ˆ α , the flutterboundary has a vertical tangent in the ( β , κ )-plane at β = ˆ β , as is visible in the rightpanel of Fig. 4. Consequently, when α = ˆ α , the maximal stable load κ α ,β crit varies smoothlyfor β ∈ (0 , ˆ β ), but when β reaches ˆ β it jumps up discontinuously from the flutter boundaryto the divergence boundary. For the system under study such jumps were first described inthe work ( ) that corrected the classical result of Bolotin ( ), whose plot in the ( β , κ )-plane did not contain the divergence boundary at all, but provided a correct shape for theflutter boundary. Notice that such overlapping of eigenvalue branches typically accompaniesoptimization of nonconservative systems and was reported in numerous studies ( , , , , , , , , , , , ). The general theory of this effect has been developedin ( , , ).We can obtain a clearer picture of the jump discontinuity by plotting the real andimaginary parts of the eigenvalues σ which describe the flutter boundary, as is done inFig. 5. For α = ˆ α , when β is decreased from the value ˆ β , a bubble of complex eigenvaluescorresponding to flutter appears, but this vanishes for β > ˆ β , resulting in the transition ofthe critical load from the flutter boundary to the divergence boundary.Looking at the discriminant (3.4) we notice that it degenerates into the equation κ cos( κ ) − sin( κ ) = 0 (3.5) he strongest column with a follower load and relocatable masses Fig. 4
Stability diagrams for (left) β = ˆ β =1 . α , κ )-plane and (right) for α =ˆ α =0 . β , κ )-plane. The solid blue curves designate the divergence boundary(3.3) and the solid green curves mark the flutter boundary (3.4). The flutter boundary in theleft panel has a crossing at the saddle point located at α = ˆ α and κ = 5 . β = ˆ β − .
01 and (left and right curves) β = ˆ β +0 .
01. In the right panel, the divergence boundaryis a horizontal blue line with height κ = ˆ κ = 7 . for β = 0 (i.e. when µ = 0) and reduces to the equationsin( κα ) − κα cos( κα ) = 0 (3.6)for β = π/ µ → ∞ ). The sets defined by the equations (3.5) and (3.6) areshown by the solid red lines and curve, respectively, in Fig. 6. The flutter boundary istangent to the planes β = 0 and β = π/ κ = κ ≈ . κ ≈ .
19 which is the optimal load for the Dzhanelidze column in whichthe relocatable mass M is absent: µ = 0. The lines κ = κ at α = 0 and α = 1 formsingularities (edges) of the flutter domain.As soon as β starts deviating from zero, a closed region of flutter instability appearsaround the horizontal red line κ = κ in the ( α , κ )-plane. Further, another region of flutteroriginates above it that touches the divergence boundary. These two regions coalesce when β reaches ˜ β = 0 . β theflutter region in the ( α , κ )-plane is simply connected, as shown in the two upper panelsof Fig. 7 corresponding to β = π/ − . β = π/ − .
15, respectively, until thisparameter passes the value β = 1 . Kirillov and Overton
Fig. 5 (Left) imaginary and (right) real roots of the characteristic polynomial (3.2) for α = ˆ α =0 . β = ˆ β = 1 . β ± .
01. A bubble ofcomplex eigenvalues appears for β = 1 . − .
01 and corresponds to flutter instability. Theblack dotted vertical line at κ = ˆ κ = 7 . β from ˆ β − .
01 to ˆ β + 0 .
01 results in the disappearance of the complex eigenvalues and henceis accompanied by the transition from the overlapping eigenvalue branches to an avoided crossingthat yields a jump in the critical load parameter to the maximal value that is reached at κ = ˆ κ onthe divergence boundary ( ); see also the right panel of Fig. 4. As β approaches π/
2, the upper portion of the flutter region concentrates around thecurve defined by (3.6), as shown in the lower panels of Fig. 7, and coincides with this curveexactly at β = π/
2. At this very limit the critical load κ attains its supremum κ ∗ whichis the intersection point of the red curve (3.6) and blue curve of the divergence boundary(3.3). Solving the equations (3.3) and (3.6) simultaneously, we find κ ∗ ≈ . , α ∗ ≈ . , β ∗ = π . (3.7)Stability diagrams in Fig. 8 presented in the ( β , κ )-plane show the decrease in the criticalload κ when α deviates from the value α ∗ , indicating that the value κ ∗ is a local supremumin the parameter space α ∈ [0 , β ∈ [0 , π/ κ ∗ is actually the global supremum. However, note that M is notdefined at β ∗ = π/
2, so the supremum is not attained. Furthermore, as ( κ, α , β ) → ( κ ∗ , α ∗ , π/ M diverges to ∞ and M converges to 0 (see (3.3)),but M is the product of two quantities, one diverging to ∞ and the other convergingto 0 (see (3.6)). For this reason it is difficult to rigorously state limiting properties of theeigenvalues λ k of M as the supremum is approached, though based on both our symbolic andnumerical calculations, it seems that, under the appropriate assumptions, the eigenvaluesconverge to a double zero eigenvalue with a Jordan block, indicating that the parameters he strongest column with a follower load and relocatable masses Fig. 6
Stability diagrams for (left) β = ˜ β = 0 . α , κ )-plane and (right)for α = ˜ α = 0 . β , κ )-plane. The solid blue curves designate the divergenceboundary (3.3), and the solid green curves mark the flutter boundary (3.4). The flutter boundaryin the left panel has a crossing at the saddle point located at α = ˜ α and κ = 6 . β = ˜ β − .
05 and (left and right curves) β = ˜ β + 0 .
05. In the right panel, the divergenceboundary is the horizontal blue line with height κ = 7 . κ = 4 . β = 0.The other red solid curve in the left panel is the solution to (3.6): the flutter boundary for the case β = π/ are on the boundary of both the flutter and divergence domains, and that the limitingeigenvalues σ k of (2.12) coalesce into a quadruple eigenvalue at ∞ .We conclude this section with an interesting observation. The supremum defined by(3.7) occurs when the divergence boundary meets the flutter boundary, and furthermore β = π/
2. If we substitute (3.6), which is the equation for the flutter boundary in the limit β = π/
2, into the divergence boundary equation (3.3), the latter can be simplified andreduced to κ ( α −
1) sin( κ ( α − κα ) − = 0 . (3.8)Writing sin( κ ( α − κα − κ + πk = 0, k ∈ Z . On the other hand, the relation(3.6) can be written as κα = tan( κα ), yielding κα = κ , where κ ≈ . κ = κ . Combining the results, weobtain κ = κ + πk, with k ∈ Z . For k = 0, we obtain κ = κ , the optimal load when themass M is absent (and the square root of the load for the optimal Dzhanelidze column),while for k = 1, we obtain κ = κ + π ≈ . M and M . Let us therefore set k = n − κ = κ + π ( n − , (3.9)4 Kirillov and Overton
Fig. 7
Stability diagrams in the ( α , κ )-plane for (upper left) β = π/ − .
5, (upper right) β = π/ − .
15, (middle left) β = π/ − .
1, (middle right) β = π/ − .
05, (lower left) β = π/ − .
01, and (lower right) β = π/ − . α ∗ ≈ . κ ∗ ≈ . he strongest column with a follower load and relocatable masses Fig. 8
Stability diagrams for (left) α = α ∗ − . α = α ∗ ≈ . α = α ∗ + 0 .
1. The green and blue curves respectively show the flutter and divergence boundaries.In the left and center panels, the critical load reaches the divergence boundary, but this is higher inthe center panel, and there it is reached only if β = π/
2. In the right panel, the flutter boundaryprevents the critical load from reaching the divergence boundary. and hence, using κα = κ , α = κ κ + π ( n − . (3.10)For n = 1 the expression (3.10) yields α = 1 and for n = 2 we have α = κ ( κ + π ) − ≈ . α ∗ provided by (3.7). Table 1 shows the values of κ and α defined by (3.9) and (3.10) for n = 1 , , . . . ,
6, while Fig. 9 shows these values asdefined by the intersections of equations (3.3) and (3.6), the divergence boundary equationand the flutter boundary equation in the limit β = π/
2, respectively. † Table 1
The approximate values of the parameters κ and α according to (3.9) and (3.10). n κ = κ + ( n − π α = κ [ κ + ( n − π ] − . . . . . . . . . . . Remarkably, the numerical computations reported in the next section for n concentratedmasses, with n = 2 , , ,
5, strongly indicate that the optimal load and the optimal value of α are precisely the values shown in Table 1 and Fig. 9. † Note that, for all n , we have tan( κ + ( n − π ) = tan( κ ) = κ . Kirillov and Overton
Fig. 9
Graphs of (red) equation (3.6) defining the flutter boundary in the limit β = π/ α and (blue) equation (3.3) defining the divergence boundary as a function of α . Theintersection points are given by the expressions (3.9) and (3.10) and listed in Table 1.
4. Numerical derivation of the optimal load for the weightless columncarrying multiple relocatable masses
Recall that, as discussed in Section 2, for a given number of masses n , our stability constraintis defined by the quadratic eigenvalue problem ( M σ + K ) u = 0 (see (2.12)). Here K is theunit matrix while M is defined by (2.11) and (2.13), which depend on the dimensionlessparameters α i and µ i = tan β i , i = 1 , . . . , n −
1, defined in (2.8), as well as a given load κ . Letus write M ( α, β, κ ) for the matrix M defined by α = [ α , . . . , α n − ] T , β = [ β , . . . , β n − ] T and κ . As noted in (2.14), the eigenvalues σ k of ( M ( α, β, κ ) σ + K ) u = 0 are related to λ k ,the eigenvalues of the matrix M ( α, β, κ ), by σ k = ± ( − λ − k ) / .The stability constraint requires that, for given ( α, β, κ ), all eigenvalues σ k should beimaginary, or equivalently, that all eigenvalue reciprocals λ − k are real and nonnegative.Clearly, another equivalent condition is that all eigenvalues λ k are real and nonnegative,interpreting 1 / ∞ . Consequently, we define a stability violation function ˜ v : R n − → R + by (cid:0) α, β, κ (cid:1) max (cid:16) Re p − λ k (cid:17) , (4.1)where the maximum is taken over all eigenvalues of M ( α, β, κ ), using the principal squareroot, hence implying that ˜ v cannot take negative values. Besides avoiding the nonlinearity inthe reciprocal, the stability violation function ˜ v has the virtue that it is continuous, thoughnot Lipschitz continuous, at points in parameter space where a positive eigenvalue λ k passesthrough the origin to the negative real axis, and hence ˜ v changes continuously from the valuezero to a positive value that grows like the square root function at zero. In this case, theparameters cross the divergence boundary, since a conjugate pair of imaginary eigenvalues σ k coalesce at ∞ and split along the real axis. The function ˜ v is also continuous, though notLipschitz continuous, at points in parameter space where two positive real eigenvalues λ k , λ ℓ coalesce and split into a complex conjugate pair, and hence again ˜ v increases from zero to apositive quantity that, generically, increases with the square root of the perturbation. In this he strongest column with a follower load and relocatable masses σ k , σ ℓ (and also their conjugates) coalesce on the imaginary axis and split into a complexpair.We argued in Section 3 that, in the case n = 2, the optimal parameter configurationis simultaneously at both the flutter boundary and the divergence boundary, likely witha double eigenvalue λ at zero (equivalently, a quadruple eigenvalue σ at ∞ ) and, if thisis the case, generically, the stability violation function ˜ v would grow at nearby parameterconfigurations with the fourth root of the perturbation.To compensate for this non-Lipschitz behavior of ˜ v , we define a modified stability violationfunction v : R n − → R by v ( α, β, κ ) = ˜ v ( α, β, κ ) ρ , ˜ v ( α, β, κ ) ∈ [0 , ρ ˜ v ( α, β, κ ) − ( ρ − , ˜ v ( α, β, κ ) ∈ [1 , ∞ ] (4.2)where ρ is a positive integer. In the situations just discussed, the choice ρ = 2 is sufficientto make v generically Lipschitz continuous at points where the parameters cross eitherthe divergence or the flutter boundary separately, and ρ = 4 is sufficient to make v Lipschitz continuous even when the parameters cross the divergence and flutter boundariessimultaneously, at least at the conjectured optimal configuration for n = 2. In ourcomputations, we experimented with choices of ρ from 1 to 5 and we found that ρ = 4 gavesignificantly better results than ρ <
4, but that setting ρ = 5 made no further improvement.Consequently, we chose to use ρ = 4. Note that the specific form of v is chosen so that itdoes not cause blow-up when ˜ v ( α, β, κ ) is large, and so that it is continuously differentiablewhere ˜ v ( α, β, κ ) = 1.However, what makes this problem particularly difficult is that as any β i → π/
2, thecoefficient µ i → ∞ in (2.12). Consider the case n = 2. We already mentioned in Section3 that as α → α ∗ , β → π/ κ → κ ∗ , we have M → ∞ and M →
0, while M is a product of tan( β ) with a second factor that converges to zero. If this second factorconverges to zero more slowly than (tan( β )) − does, so that | M | → ∞ , a change in itssign causes an eigenvalue λ k to discontinuously pass through ∞ from the positive real tothe negative real axis, implying infinitely large growth in the stability violation v as theparameters cross the divergence boundary. This presents a serious difficulty as we shall see.In order to solve our optimization problem, we need to impose the stability constraintnot only at a given point ( α, β, κ ), but also at all points ( α, β, ν ) with ν ∈ [0 , κ ]. Althoughwe could construct an approximation to v ( α, β, · ) on the interval [0 , κ ] using approximationsoftware such as Chebfun ( ), this is computationally expensive. In our optimizationcomputations, we found that a more effective approach is to impose the stability constrainton a coarse grid of ˜ q logarithmically spaced points on (0 , κ ], defining c ( α, β, κ ) = max j ˜ q ( v ( α, β, ν j ) : ν = κ, ν j = (1 − − j ) κ, j = 1 , . . . , ˜ q ) (4.3)and imposing the constraint c ( α, β, κ )
0, or equivalently, c ( α, β, κ ) = 0. Then, after apotential solution is obtained by optimization, we check its stability on a much finer gridof q ≫ ˜ q uniformly spaced points on (0 , κ ), rejecting it if this test is not passed. We foundthat using a coarse grid with ˜ q = 10 points and a fine grid with q = 10 ,
000 points worked8
Kirillov and Overton well, typically with the majority of the solutions obtained by optimization that are feasiblefor the coarse grid also passing the fine grid test.We then pose our optimization problem assup α ∈ R n − ,β ∈ R n − ,κ ∈ R κ (4.4)subject to c ( α, β, κ ) , α . . . α n − , β i π/ , i = 1 , . . . , n − , This is not an easy problem to solve, since the stability constraint is nonconvex andnonsmooth, as well as discontinuous as β i → π/
2. We tackled it using granso (GRadient-based Algorithm for Non-Smooth Optimization), a recently developed open-source softwarepackage for nonsmooth constrained optimization ( , ).As its name suggests, the algorithm implemented in granso is based on employinguser-supplied gradients. This might seem contradictory since it is intended for nonsmoothoptimization problems, but although the constraints are not differentiable everywhere,they are differentiable almost everywhere. Specifically, the stability violation function v is differentiable at ( α, β, κ ) if the following conditions hold:(i) β i < π/ i = 1 , . . . , n − j ∈ (0 , . . . , ˜ q )(iii) the maximum in (4.1) is attained only at one eigenvalue λ k of M ( α, β, ν j )(iv) this eigenvalue λ k is simple and nonzero.Thus, evaluating the gradient of v makes sense almost everywhere in parameter space. Ofcourse, the gradient does not vary continuously, but granso is designed to exploit gradientdifference information, even near points where the gradient varies discontinuously, buildinga model of the constraint function on the parameter space using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton updating method. For more details, see ( ), andfor application of BFGS in other stability optimization problems, see ( ) and the paperscited there.To derive the gradient of v , we need to differentiate an eigenvalue λ k with respect tochanges in the matrix M . Let us write M ( t ) = M + t ( ∆M ) and let λ ( t ) denote theeigenvalues of M ( t ). It is well known ( ) that, if λ k = λ (0) is a simple eigenvalue of M = M (0) satisfying the right and left eigenvector equations M u = λu and w ∗ M ∗ = λw ∗ ,where the asterisk denotes complex conjugate transpose, then ddt λ ( t ) (cid:12)(cid:12)(cid:12)(cid:12) t =0 = w ∗ ( ∆M ) uw ∗ u . With this in mind, deriving the gradient of v with respect to the 2 n − α, β, κ ) is straightforward, employing the chain rule to incorporate the variation in thepower function in (4.2), the square root in (4.1), and the formulas (2.13), (2.11) and (2.15).We now describe our experiments using granso (version 1.6.4), running in matlabhe strongest column with a follower load and relocatable masses
50 100 150 200 250 300 350 400 450 50002468 n = 2, conjectured optimal loadcomputed loadstermination codes
50 100 150 200 250 300 350 400 450 50000.51 distances: 0 conjectured
50 100 150 200 250 300 350 400 450 50000.511.5 mass ratio arctangents: 0 /2 /2
50 100 150 200 250 300 350 400 450 50010 -20 eigenvalues of M Fig. 10
Solving (4.4) for n = 2 (release R2020a) on a MacBook Air laptop, to solve (4.4). We used the default choice ofparameters with the following exceptions: we set maxit , the limit on the iteration count, to500, and we set the tolerances opt_tol and feas_tol to zero, to obtain the highest possibleaccuracy. We added bound constraints on the load variable formulated as 0 κ κ max with κ max = 1 . × ( κ + ( n − π ), that is, with a lower bound of zero and an upper boundset to 10% higher than the conjectured optimal value of κ given in Table 1. Since granso may generate iterates violating these bounds or the other bound constraints in (4.4), wedefined v to be zero if κ β i in (2.15) by pi/2 , the 16 digit rounded valueof π/
2, if β i exceeds pi/2 , to avoid the discontinuity in the tangent function at π/ tan(pi/2) ≈ . × has the desired positive sign). Because of the difficulty ofthe problem, we ran the code from many randomly generated starting points, with theinitial values for κ , α i and β i generated from the uniform distribution on [0 , κ max ], [0 , , π/
2] respectively, with the α i then sorted into increasing order.0 Kirillov and Overton
Results for n = 2Our analytical discussion of the case n = 2 was given in Section 3; the results here stronglysupport our claim that the optimal configuration is given by (3.7). Fig. 10 shows the resultsobtained by running granso from 1000 randomly generated starting points. Of the 1000candidate solutions generated by granso , 734 satisfied the bound and coarse grid stabilityconstraints imposed by granso , and of these, 691 also passed the fine grid stability testdescribed above. The top panel in the figure shows the computed optimal loads κ for thebest 500 of these feasible solutions, sorted into decreasing order, while the second and thirdpanels show the associated final values of α and β computed by these same 100 runs. Thefourth panel shows the eigenvalues of the final associated matrix M ( α , β , κ ). The highesttwo final values of κ agree with each other, and with the conjectured optimal load κ + π given in Table 1 and in (3.7), to 10 digits. The final values for α and β for these sametwo best results agree with the conjectured optimal values κ [ κ + π ] − (see Table 1) and π/
2, to 10 and 12 digits, respectively. It’s also worth noting that the top 100 final valuesfor the computed optimal load agree with κ + π to 4 digits.Looking at all four panels of Fig. 10, we see that the top 500 results come in severalclearly distinct flavours. The first flavour is exhibited by the best 180 or so runs which allgive good approximations to the conjectured optimal value κ + π . However, starting withthe 284th result, we find a very different flavour: many runs find that the computed optimalload is about κ = 4 . κ , the optimal load for the Dzhanelidze column,to four digits. Clearly, this is a locally maximal value for (4.4); otherwise, it would not befound so frequently. If we look at the associated computed α and β values, usually α isclose to zero, but if not, then β is close to zero. It is easily checked that, regardless of thevalue of β , if α = 0 then M ( α , β, κ ) is the zero matrix, with a double semisimple zeroeigenvalue, so this locally optimal parameter configuration, like the conjectured globallyoptimal configuration (3.7), is on both the flutter and divergence boundaries. Physically,this corresponds to the mass M being fixed at the mounted end of the column. On theother hand, regardless of the value of α , if β = 0, then M ( α , β , κ ) has all zero entriesexcept for M , and hence has a double zero eigenvalue with a Jordan block. Again, thisparameter configuration is on both the flutter and divergence boundaries, and physically,it corresponds to the mass M having zero weight. Note that the computed eigenvalues forthis second flavour of solutions are relatively small.The third flavour of results is exhibited by the results numbered approximately 180 to280. In these cases, granso terminated prematurely, without approximating a globally orlocally maximal value, and we can see also that, on average, the larger κ is, the closer α isto its conjectured optimal value. Investigation of these cases shows that termination occursbecause of the discontinuity in the stability constraint that we described above. This isalso supported by the enormous associated eigenvalues of M shown in the fourth panel.Note also that as κ increases towards its optimal value, these eigenvalues decrease, but theyneither converge to specific values, nor do they become very small. In fact, the matrix M associated with the best computed optimal κ is " . × − . × − . × − . × − Although its eigenvalues are not close to each other or to zero, they are small relative to the he strongest column with a follower load and relocatable masses
50 100 150 200 250 300 350 400 450 50002468 n = 2, conjectured optimal loadcomputed loadstermination codes
50 100 150 200 250 300 350 400 450 50000.51 distances: 0 conjectured
50 100 150 200 250 300 350 400 450 50000.511.5 mass ratio arctangents: 0 arctan(100) /2
50 100 150 200 250 300 350 400 450 50010 -20 eigenvalues of M
50 100 150 200 250 300 350 400 450 50002468 n = 2, conjectured optimal loadcomputed loadstermination codes
50 100 150 200 250 300 350 400 450 50000.51 distances: 0 conjectured
50 100 150 200 250 300 350 400 450 50000.511.5 mass ratio arctangents: 0 arctan(10) /2
50 100 150 200 250 300 350 400 450 50010 -20 eigenvalues of M Fig. 11
Solving (4.4) for n = 2 with (left) the constraint β arctan(100) and (right) theconstraint β arctan(10) norm of the matrix, and their associated right eigenvectors are almost identical, indicatingthe nearby presence of a double eigenvalue. Furthermore, the diagonal and upper triangularelements are very small compared to the norm of the matrix, implying that a relatively smallperturbation removing them yields a Jordan block with a double zero eigenvalue.Finally there is a fourth flavour of results: those that did not even reach a goodapproximation to the locally optimal value κ .A final comment on Fig. 10: the granso termination codes are plotted at the bottomof the first panel. The value 1 means that granso terminated because the limit of 500iterations was reached, while the value 2 means that it terminated because it could not finda higher feasible value for the load. Observe that the latter termination always occurred forthe runs which approximated the conjectured globally optimal load κ + π well (the firstflavour) and the runs that obtained loads higher than κ but terminated without reachinga good approximation to κ + π (the third flavour). Thus, increasing the iteration limitwould not have improved any of these values. On the other hand, the runs that provideda good approximation to the locally maximal value κ (the second flavour) or terminatedbefore reaching that value (the fourth flavour) sometimes, but not always, terminated byexceeding the maximum iteration limit.The physical interpretation of the conjectured supremum (3.7) is that, since the massratio µ is infinite, the mass M mounted on the free end of the rod is zero. It’s of interestto consider what happens if we disallow this case, putting an upper limit on µ . Fig. 112 Kirillov and Overton
50 100 150 200 250 300 350 4000510 n = 3, conjectured optimal loadcomputed loadstermination codes
50 100 150 200 250 300 350 40000.51 distances: 0 conjectured
50 100 150 200 250 300 350 40000.511.5 mass ratio arctangents: 0 i /2, i=1,2 /2
50 100 150 200 250 300 350 40010 -20 eigenvalues of M
50 100 150 200 250 300 350 400 450 5000510 n = 3, conjectured optimal loadcomputed loadstermination codes
50 100 150 200 250 300 350 400 450 50000.51 distances: = =
50 100 150 200 250 300 350 400 450 50000.511.5 mass ratio arctangents: 0 = /2 = /2
50 100 150 200 250 300 350 400 450 50010 -20 eigenvalues of M Fig. 12
Solving (4.4) for n = 3 with (left) no additional constraints and (right) with the constraints α = α ∗ = κ [ κ + 2 π ] − , β = π/ shows the results when we introduce the constraint µ
100 (left) or µ
10 (right) bylimiting β to arctan(100) or arctan(10) respectively. For µ . µ
10, we find the optimal load 7 . κ + π . The biggest difference we observefrom comparing Fig. 11 with Fig. 10 is that the eigenvalues of the final computed M arenow dramatically reduced, from more than 10 to less than 100 and 15 respectively. Thus,we obtain an only slightly reduced optimal load while introducing a much more physicallyreasonable model with much better numerical properties.4.2 Results for n = 3The left panel in Fig. 12 shows the results for solving (4.4) for n = 3. There are five variables: α , α , β , β and κ . Of the 1000 candidate solutions generated by granso , 517 satisfied thebound and coarse grid stability constraints imposed by granso , and of these, 416 passed thefine grid test. We see immediately that the problem for n = 3 is significantly harder than for n = 2, with not many runs approximating the conjectured optimal value well. Nonetheless,the two best runs generate κ ≈ . κ + 2 π to five digits. These two runs also generate α ≈ . β ≈ . α ∗ = κ [ κ + 2 π ] − and π/ α = α ∗ and β = π/ α , β and κ . Then the best two he strongest column with a follower load and relocatable masses
50 100 150 200 250 3000510 n = 3, conjectured optimal loadcomputed loadstermination codes
50 100 150 200 250 30000.51 distances: 0 conjectured
50 100 150 200 250 30000.511.5 mass ratio arctangents: 0 i arctan(100), i=1,2 /2
50 100 150 200 250 30010 -20 eigenvalues of M
50 100 150 200 250 3000510 n = 3, conjectured optimal loadcomputed loadstermination codes
50 100 150 200 250 30000.51 distances: 0 conjectured
50 100 150 200 250 30000.511.5 mass ratio arctangents: 0 i arctan(10), i=1,2 /2
50 100 150 200 250 30010 -20 eigenvalues of M Fig. 13
Solving (4.4) for n = 3 with the constraints (left) β i arctan(100), i = 1 , β i arctan(10), i = 1 , runs generate κ agreeing with κ + 2 π to 12 digits, and the best 100 runs agree with thisto 10 digits. Together, the results reported in the left and right panels of Fig. 12 make aconvincing argument that the values shown in Table 1 are indeed the optimal values for κ and α when n = 3, and that the optimal β is again π/ µ i
100 (left) or µ i β i arctan(100), i = 1 , β i arctan(10), i = 1 , µ i κ + 2 π .However, when we constrain µ i
10, the best optimal load found is only 7.59, which is a30% reduction from κ + 2 π . When we repeat these runs with 10,000 starting points insteadof 1000, these numbers are only slightly improved.4.3 Results for n = 4 and n = 5The optimization problem is so much harder for n = 4 and n = 5 that we needed 10,000starting points to get good results, even when we set α to its conjectured optimal valueand β to π/
2, optimizing over the remaining 5 and 7 variables, respectively. The resultsare shown in the left and right panels of Fig. 14. For n = 4, the best 5 results agree with ourconjectured optimal load κ + 3 π to 8 digits, while for n = 5, the best 25 results agree with κ + 4 π to 7 digits. These results strongly support our conjecture regarding the optimalload κ for n masses given in Table 1.4 Kirillov and Overton
50 100 150 200 250 300 350 400 450 5000510 n = 4, conjectured optimal loadcomputed loadstermination codes
50 100 150 200 250 300 350 400 450 50000.51 distances: = =
50 100 150 200 250 300 350 400 450 50000.511.5 mass ratio arctangents: 0 i = /2 (i=2,3) = /2
50 100 150 200 250 300 350 400 450 50010 -20 eigenvalues of M
50 100 150 200 250 300 350 400 450 500051015 n = 5, conjectured optimal loadcomputed loadstermination codes
50 100 150 200 250 300 350 400 450 50000.51 distances: = =
50 100 150 200 250 300 350 400 450 50000.511.5 mass ratio arctangents: 0 i = /2 (i=2,3,4) = /2
50 100 150 200 250 300 350 400 450 50010 -20 eigenvalues of M Fig. 14 (Left) solving (4.4) for n = 4 with the constraints α = α ∗ = κ [ κ + 3 π ] − , β = π/ n = 5 with the constraints α = α ∗ = κ [ κ + 4 π ] − , β = π/
5. Concluding Remarks
We believe we have made a convincing case that the supremal load for the strongest stableweightless column with a follower load and n relocatable masses is, in the dimensionlessmodel defined in Section 2, κ +( n − π , where κ is the smallest positive root of tan( κ ) = κ .This conjecture has not previously appeared in the literature as far as we know, except inthe case n = 1 where it has been known to be true since 1963 ( ). We have given a detailedanalytical derivation of this result for n = 2, and presented extensive computational resultsthat support it for n = 2 , , ,
5, using numerical nonsmooth optimization.With this model problem effectively solved, we believe it would be interesting to applyour nonsmooth optimization techniques to more realistic columns with follower loads ( ),such as the Beck, Pfl¨uger and Leipholz columns with a single free end as well as to free-freebeams both with distributed and concentrated masses to get new insights about the natureof the optimal solution to these long-standing optimization problems. We believe it is alsoimportant to consider extending traditional stability constraints to more robust stabilityconstraints based on pseudospectra ( ), a topic that is beyond the scope of this paper. Acknowledgements.
The authors thank Tim Mitchell, the author of granso , for manyhelpful discussions and suggestions regarding the formulation of the stability constraint.They also thank the London Mathematical Society for supporting the second author’s visitto Northumbria through the Scheme 4 Research in Pairs grant No 41820. The second authorwas supported in part by the U.S. National Science Foundation Grant DMS-2012250. he strongest column with a follower load and relocatable masses References1.
Beck M (1952) Die Knicklast des einseitig eingespannten, tangential gedr¨uckten St¨abes.Z angew Math Mech 3:225–228. Carr J and Malhardeen M Z M (1979) Beck’s problem. SIAM J Appl Math 37:261–262. Sugiyama Y, Langthjem M, Katayama K (2019) Dynamic stability of columns undernonconservative forces: theory and experiment. Solid mechanics and its applications.Vol 255 Springer, Berlin. Ziegler H (1953) Linear elastic stability. A critical analysis of methods, First part. ZAMPZ angew Math Phys 4:89–121. Ziegler H (1953) Linear elastic stability. A critical analysis of methods, Second part.ZAMP Z angew Math Phys 4:167–185. Bolotin V V (1963) Nonconservative problems of the theory of elastic stability. PergamonPress, Oxford. Kirillov O N (2013) Nonconservative stability problems of modern physics. De Gruyterseries in mathematical physics. Vol 14 De Gruyter, Berlin, Boston. Bayly P V and Dutcher S K (2016) Steady dynein forces induce flutter instability andpropagating waves in mathematical models of flagella. J R Soc Interface 13:20160523. De Canio G, Lauga E, Goldstein R E (2017) Spontaneous oscillations of elastic filamentsinduced by molecular motors. J R Soc Interface 14:20170491.
Zhu L, Stone H A (2019) Propulsion driven by self-oscillation via an electrohydrody-namic instability. Phys Rev Fluids 4:061701.
Zhu L, Stone H A (2020) Harnessing elasticity to generate self-oscillation via anelectrohydrodynamic instability. J Fluid Mech 888:A31
Sugiyama Y, Langthjem M, Ryu B-J (1999) Realistic follower forces. J Sound Vibr225:779–782.
Sugiyama Y, Langthjem M, Ryu B-J (2002) Beck’s column as the ugly duckling. JSound Vibr 254:407–410.
Sundararajan C (1975) Optimization of a nonconservative elastic system with stabilityconstraint. J Opt Theory Appl 16(3/4):355–378.
Park Y P, Mote C D (1985) The maximum controlled follower force on a free-free beamcarrying a concentrated mass. J Sound Vibr 98(2):247–256.
Kirillov O N, Seyranian A P (1998) Optimization of stability of a flexible missile underfollower thrust. AIAA Paper 98-4969:2063–2073.
Sugiyama Y, Matsuike J, Ryu B-T, Katayama K, Kinoi S, Enomoto N (1995) Effect ofconcentrated mass on stability of cantilevers under rocket thrust. AIAA J. 33(3):499–503.
Sugiyama Y, Langthjem M A, Iwama T, Kobayashi M, Katayama K, Yutani H (2012)Shape optimization of cantilevered columns subjected to a rocket-based follower forceand its experimental verification. Struct. Multidisc. Opt. 46:829–838.
Bigoni D, Noselli G (2011) Experimental evidence of flutter and divergence instabilitiesinduced by dry friction. J Mech Phys Sol 59:2208–2226.
Bigoni D, Kirillov O N, Misseroni D, Noselli G, Tommasini M (2018) Flutter anddivergence instability in the Pfl¨uger column: Experimental evidence of the Zieglerdestabilization paradox. J Mech Phys Sol 116:99–116.
Bigoni D, Misseroni D, Tommasini M, Kirillov O N, Noselli G (2018) Detecting singular6
Kirillov and Overton weak-dissipation limit for flutter onset in reversible systems. Phys Rev E 97(2):023003.
Pfl¨uger A (1955) Zur Stabilit¨at des tangential gedr¨uckten St¨abes. Z Angew Math Mech35(5):191.
Deineko K S, Leonov M Ia (1955) A dynamic method for the investigation of the stabilityof a compressed bar. Prikl Mat Mekh 19:738–744.
Tommasini M, Kirillov O N, Misseroni D, Bigoni D (2016) The destabilizing effect ofexternal damping: Singular flutter boundary for the Pfl¨uger column with vanishingexternal dissipation. J Mech Phys Sol 91:204–215.
Oran C (1972) On the significance of a type of divergence. J Appl Mech 39:263–265.
Sugiyama Y, Kashima K, Kawagoe H (1976) On an unduly simplified model in thenon-conservative problems of elastic stability. J Sound Vib 45(2):237–247.
Chen L W, Ku D W (1992) Eigenvalue sensitivity in the stability analysis of Beck’scolumn with a concentrated mass at the free end. J Sound Vib 153(3):403–411.
Gajewski A, Zyczkowski M (1988) Optimal structural design under stability constraints.Kluwer, Dordrecht.
Claudon J L (1975) Characteristic curves and optimum design of two structuressubjected to circulatory loads. J. de Mecanique 14(3):531–543.
Hanaoka M, Washizu K (1980) Optimum design of Beck’s column. Comp Struct11(6):473–480.
Kounadis A N, Katsikadelis J T (1980) On the discontinuity of the flutter load forvarious types of cantilevers. Intern J Solids Struct 16:375–383.
Bogacz R, Frischmuth K (2018) On optimality of column geometry. Arch Appl Mech88:317–327.
Mahrenholtz O, Bogacz R (1981) On the shape of characteristic curves for optimalstructures under non-conservative loads. Arch Appl Mech 50:141–148.
Langthjem M A, Sugiyama Y (2000) Optimum design of cantilevered columns underthe combined action of conservative and nonconservative loads Part I: The undampedcase. Comp Struct 74(4):385–398.
Ringertz U T (1994) On the design of Beck’s column. Struct Opt 8(2):120–124.
Temis Yu M, Fedorov I M (2007) Shape optimization of nonconservatively loaded beamswith a stability criterion. Probl Strength Plast 69:24–37.
Kirillov O N, Seyranian A P (2002) Metamorphoses of characteristic curves incirculatory systems. J Appl Math Mech 66(3):371–385.
Kirillov O N, Seyranian A P (2002) A non-smooth optimization problem. Moscow UnivMech Bulletin 57(3):1–6.
O’Reilly O M, Malhotra N K, Namachchivaya N S (1996) Some aspects of destabilizationin reversible dynamical systems with application to follower forces. Nonlin Dynamics10:63–87.
Seyranian A P, Kirillov O N (2001) Bifurcation diagrams and stability boundaries ofcirculatory systems. Theor Appl Mech 26: 135–168.
Kirillov O N, Seyranian A P (2004) Collapse of the Keldysh chains and stability ofcontinuous non-conservative systems. SIAM J Appl Math 64(4):1383–1407.
Kirillov O N, Overton M L (2013) Robust stability at the swallowtail singularity.Frontiers in Physics 1: 24.
Kirillov O N (2011) Singularities in structural optimization of the Ziegler pendulum.Acta Polytechn 51(4):32–43. he strongest column with a follower load and relocatable masses Langthjem M A, Sugiyama Y (1999) Optimum shape design against flutter of acantilevered column with an end-mass of finite size subjected to a non-conservativeload. J Sound Vibr 226(1):1–23.
Katsikadelis J T, Tsiatas G C (2007) Optimum design of structures subjected to followerforces. Intern J Mech Sci 49:1204–1212.
Kordas Z, Zyczkowski M (1963) On the loss of stability of a rod under a super-tangentialforce. Arch Mech Stos 1(15):7–31.
Lee H P (1997) Flutter of a cantilever rod with a relocatable lumped mass. ComputMethods Appl Mech Engrg 144:23–31.
Hauger W (1967) Bemerkungen zu dem Einfluss der Massenverteilung bei nichtkonser-vativen Stabit¨atsproblemen elastischer St¨abe. Ing.-Arch. 35:283–291.
Leipholz H, Lindner G (1970) ¨Uber den Einfl¨uss der Massenverteilung auf dasnichtkonservative Knicken von St¨aben. Ing.-Arch. 39:187–194.
Kapoor R N, Leipholz H (1974) Stability analysis of a damped polygenic systems withrelocatable mass along its length. Ing-Arch 43:233–239.
Kapoor R N, Leipholz H (1974) On mass distribution, rotary inertia and externaldamping of a viscoelastic polygenic system. Z Angew Math Mech 54(3):205–208.
Kounadis A N (1977) Stability of elastically restrained Timoshenko cantilevers withattached masses subjected to a follower force. ASME J Appl Mech 44(4):731–736.
Bazant Z P, Cedolin L (2010) Stability of structures: elastic, inelastic, fracture anddamage theories. World Scientific, Singapore.
Kagan-Rosenzweig L M (2001) Quasi-static approach to non-conservative problems ofthe elastic stability theory. Int J Solids Struct 38:1341–1353.
Ingerle K (2013) Stability of massless non-conservative elastic systems. J Sound Vibr332:4529–4540.
Kagan-Rosenzweig L M (2014) Topics in nonconservative stability theory. (St.-Petersburg Univesity Press, St. Petersburg. In Russian.
Curtis F E, Mitchell T, Overton M L (2017) A BFGS-SQP method for nonsmooth,nonconvex, constrained optimization and its evaluation using relative minimizationprofiles. Optim Methods Softw 32(1):148–181.
Gallina P (2003) About the stability of non-conservative undamped systems. J. SoundVibr. 262:977–988.
Bulatovic R M (2011) A sufficient condition for instability of equilibrium of non-conservative undamped systems. Phys Lett A 375:3826–3828.
Driscoll T A, Hale N, Trefethen L N (2014), Chebfun Guide, Pafnuty Publications,Oxford.
Overton M L (2014) Stability optimization for polynomials and matrices. In: NonlinearPhysical Systems: Spectral Analysis, Stability and Bifurcations (O. Kirillov, D.Pelinovsky, eds.), 351–375. John Wiley & Sons Inc., New York.
Greenbaum A, Li R-C, Overton M L (2020). First-order perturbation theory foreigenvalues and eigenvectors. SIAM Rev 62 (1) 463–482.64.