Singular divergence instability thresholds of kinematically constrained circulatory systems

Static instability or divergence threshold of both potential and circulatory systems with kinematic constraints depends singularly on the constraints’ coeﬃcients. Particularly, the critical buckling load of the kinematically constrained Ziegler’s pendulum as a function of two coeﬃcients of the constraint is given by the Plücker conoid of degree n = 2. This simple mechanical model exhibits a structural instability similar to that responsible for the Velikhov–Chandrasekhar paradox in the theory of magnetorotational instability.


Introduction
Linear stability analysis of dissipative systems is frequently accompanied by violent behavior of critical parameters at the onset of oscillatory instability (flutter) when dissipation parameters are infinitesimally small. For example, the critical flutter load of the celebrated in structural mechanics Ziegler's pendulum [1,2] exhibits a drop of order O (1) with respect to the undamped case as soon as the dissipation of order o(ε), where 0 < ε ≪ 1, is added [3]. Moreover, the zero dissipation limit generically does not collapse to the undamped result and, if more than one damping parameter is present corresponding to different dissipation mechanisms, the critical flutter load in the limit of vanishing dissipation depends on the order in which the damping coefficients tend to zero. The zero dissipation limits are thus non-commuting-the eff e c tk n o w na l s oi nh y d r o d y n a m i c s [4]. Similar singular phenomena are known for the critical vibration frequency at the flutter threshold [5]. Already in 1956 the structurally unstable behavior of the critical flutter load in the presence of small dissipation was related to the Whitney umbrella singularity on the marginal stability * Corresponding author. Tel.: +49 3512602154; fax: +49 3512602007. boundary by Bottema [6]. Recently, drops in the vibration frequencies of nearly-inextensible post-buckled rods in the limit of vanishing thickness were discovered in a pure conservative problem without dissipation and explained via the Whitney umbrella singularity [7].
In this respect, a question arises, whether a critical buckling (or divergence) load in a conservative or non-conservative system could exhibit the singular behavior described above? Remarkably, a related phenomenon is well-known in magnetohydrodynamics. Namely, in the theory of magnetorotational instability (MRI), which is a static instability of a hydrodynamically stable cylindrical Couette-Taylor flow of an electrically conducting fluid by a magnetic field parallel to the axis of rotation, the divergence threshold in the case of vanishing magnetic field does not collapse to that of Rayleigh's centrifugal instability. This effect, known as the Velikhov-Chandrasekhar paradox, was discovered in 1959-60 in the works [8,9] and subsequently discussed by many authors who employ the assumptions that the fluid is inviscid and perfectly electrically conducting [10][11][12][13][14]. We note that a fluid-mechanical analogue of magnetorotational instability-the so-called elastorotational instability-is known and studied both theoretically and experimentally [15]. determined by the Rossby number at the divergence (MRI) threshold as a function of the Alfvèn frequency and resistive frequency for arbitrary viscosity (including the inviscid case) is given by the surface which is reduced to the Plücker conoid of degree n = 2 [19]. The singular surface with a self-intersection and two Whitney umbrella singularities, one at the Velikhov-Chandrasekhar threshold and another one at the Rayleigh threshold provides a switching between the two static instability criteria when the ratio of the Alfvèn frequency to the resistive one is varying.
Is there a pure mechanical analogue of such a switching between different divergence instability thresholds, for example, in the problems of buckling of simple structures? A positive answer comes from the recent study [20] of mechanical systems with the kinematic constraints motivated by material instability problems of geomechanics [21][22][23][24][25][26]. In this Letter we demonstrate that the kinematic constraints may result in significant changes in the critical buckling load because the limit of vanishing constraints is singular. For the example of the Ziegler pendulum under partial follower load we find that switching between the first and the second buckling loads of the unconstrained pendulum is governed in the presence of a constraint by the Plücker conoid surface.

Mechanical systems with kinematic constraints
In [21] it has been proposed to interpret failure in geomaterials (its relation to landslides is discussed, e.g., in [22][23][24][25])a sal o s so f stability of a stress-strain equilibrium state. A usual tool to investigate instability occurrence from a mechanical viewpoint is the Hill stability condition which is based on the sign of the second-order work [26]. This criterion corresponds to the loss of definite positiveness of the elasto-plastic stiffness matrix (the constitutive 6 × 6 matrix is defined in the usual 6-dimensional stress rate and strain rate spaces). It is thus related to a spectral analysis of the symmetric part of the elasto-plastic matrix [23], distinguishing essentially the second-order work criterion applied to plastic non-associative materials (like geomaterials) from the classical plastic failure theory related to the singularity of the stiffness matrix.
In geotechnical tests, one often controls linear combinations of stresses or strains. The active loading to the boundary of a specimen is usually applied by means of static or kinematic control variables that are linear combinations of the stress or strain components. For example, the classical undrained axisymmetric triaxial test is an isochoric loading with the sum of the incremental strains maintained equal to zero during the loading.
Since the critical load depends on the scenario of the loading, ap r o b l e ma r i s e sw i t ht h ed e t e r m i n a t i o no ft h i sc r i t i c a ll o a dw i t h respect to the constraints imposed during the loading process such as the preservation of the volume. Such linear combinations correspond to what is called in structural mechanics a "kinematic constraint" [27,28]. According to the experiments and in agreement with the second-order work criterion, for these mixed loading paths failure can appear in the bifurcation domain (thus strictly inside the plastic limit surface) when the second-order work is vanishing at least for a given loading direction.
Therefore, the second-order work criterion is related to a critical kinematic constraint, which can be interpreted as an instability direction when the material stability analysis is considered or as a direction inside the isotropic cone of the elasto-plastic matrix. This indicates existence of close parallels between material and structure stability analyses [20,21].
In structural mechanics, the stiffness matrix K encompasses both the constitutive behavior of the materials and the balance equations where M and K are the n × n mass and stiffness matrices and x is the n-vector of perturbation around the equilibrium. Note that the kinematic constraints prescribed to a structure can be regarded as a notion similar to that of the control (or loading) parameter in geomechanics [20].
Let us introduce the kinematic constraint by means of the n × q-matrix α,a sα T x = 0. Then, the constrained system is described via the q-vector of Lagrange multipliers, Λ,a s The two equations form an extended (n + q) × (n + q) circulatory system [29,30] Assuming x ∼ exp(µt) in Eq. (3), we obtain an eigenvalue problem for a matrix polynomial Mµ 2 + K. The eigenvalues, µ,o ft h i s problem determine stability of the extended system (3).
A circulatory system is marginally stable when all its eigenvalues are pure imaginary and simple. When there exists a real positive eigenvalue, the system is statically unstable (divergence or non-oscillatory instability). A complex eigenvalue with the positive real part yields dynamical instability (flutter or oscillatory instability).
Transition between stability and flutter in circulatory systems with the change of parameters generically happens through the merging of two pure imaginary eigenvalues into one double pure imaginary eigenvalue with the Jordan block and its subsequent splitting into a complex conjugate pair (a non-semisimple 1 : 1r e sonance or a reversible Hopf bifurcation) [3].
Transition between stability and divergence is accompanied by the merging of a pair of pure imaginary eigenvalues into a double zero eigenvalue with the Jordan block that then splits into two real eigenvalues of different sign. Thus, the boundary between the divergence and stability is characterized by the double zero eigenvalues µ.
Since the characteristic equation det( Mµ 2 + K) = 0i sg i v e nb y a real polynomial with respect to µ 2 , the divergence boundary is determined by the vanishing free term of the characteristic polynomial, which yields µ 2 = 0. According to the modified Leverrier's algorithm [32], the free term of the characteristic polynomial of a matrix pencil Mµ 2 + K is exactly det K.F o rd e tK = 0S c h u r ' sd eterminant identity [31] gives Therefore, the divergence instability boundary of the extended system (3) is If the matrix K depends on the loading parameter, say p,b yt h e implicit function theorem it can be expressed as a function of the components of the matrix of constraints, α: p = p(α).  Note that the kinematically constrained stability problems were previously discussed in other areas of mechanics, e.g. in [27,28]. In [20] the effect of additional constraints on the divergence instability of the Ziegler pendulum with partial follower load (a nonconservative problem with coexistence of divergence instability with flutter instability) has been treated. When the divergence instability controls the problem, it was found that additional constraints may both increase and decrease the divergence load of the free non-conservative system [20]. The constraint related to the lowest critical buckling load is given exactly by the secondorder work criterion. Excluding flutter instabilities, the secondorder work criterion thus forms an envelope of the divergence boundaries calculated for every particular kinematic constraint and therefore provides a sufficient stability condition for a constrained pendulum [20]. As a consequence, for all kinematic constraints, there necessarily exists a constraint which destabilizes the free Ziegler pendulum [20].

Kinematically constrained Ziegler's pendulum
The free pendulum consists of two weightless bars of equal length l which carry concentrated masses m 1 = 2m and m 2 = m, one bar attached via a planar joint to the basement at the point O , see Fig. 1(a). Two rotation angles at the joints constitute the state vector x = (θ 1 ,θ 2 ) T . The stiffness coefficients at the joints are assumed to be equal k.T h ef o r c eF is applied at the free end of the upper bar and its direction constitutes an angle γ θ 2 with the vertical, Fig. 1(a).
We , where p = Fl/k is the loading parameter and γ is the parameter characterizing the proportion between the follower and the dead components of the load: γ = 0 corresponds to the conservative case and γ = 1-to the pure follower loading.
Assuming solutions to Eq. (6) have the form x = u exp(λt),w e arrive at the eigenvalue problem where u is an eigenvector at the eigenvalue λ. When two pure imaginary eigenvalues λ of different signs of the non-constrained pendulum merge at zero, Fig. 2(a), we have det K = 0, which yields the critical divergence (buckling) loads The pendulum is stable against buckling for the loads in the interval p < p 1 . In the interval p 1 p p 2 the pendulum exhibits buckling by the first mode. For p > p 2 and small γ 0t h es e cond mode becomes unstable by divergence, see Fig. 2.W i t ht h e increase of γ , p 2 becomes a point on the boundary between divergence and stability domains [20].A tγ = 5/9 the upper and lower divergence loads coincide: p 1 = p 2 = 3/2.
Consider a point P on the upper bar of the pendulum and assume that it can move only along the vertical axis, which can be realized by means of the sliding support as is shown in Fig. 1(b). This implies the following constraint on the lateral displacement of the point [20] l sin θ 1 + a sin θ 2 = 0.
Taking into account the holonomic constraint (10) Calculating the determinant of the extended stiffness matrix we find the critical divergence load of the constrained pendulum [20] p d (α 1 , α 2 , γ ) = α 2 1 + 2α 1 α 2 + 2α 2 Kinematic constraints reduce the number of degrees of freedom. While the free Ziegler's pendulum is a two-degrees-offreedom system, the constrained system is effectively a singledegree-of-freedom one. Assuming that we have θ 1 ∼ exp(µt), θ 2 ∼ exp(µt), and Λ ∼ exp(µt), we find that the eigenvalues of the con- Curiously enough, in a recent study of optimal mass distribution in the classical Ziegler's pendulum [33], the optimal solutions correspond to vanishing of one of the point masses which yields equations similar to (3). By this reason the optimal solutions obtained in [33] have in fact only one degree of freedom and thus can lose stability only by divergence despite the initial design of the pendulum is susceptible only to the flutter instability.
Therefore, when γ is within the interval 0 γ 1, the constrained pendulum is stable for p < p d and buckles when p p d .
Naturally, there exists only one buckling load p d for a given γ and a given pair (α 1 > 0, α 2 > 0). However, varying the coefficients α 1 and α 2 , we can switch between the critical loads of the nonconstrained system. As we will see in the following, the switching surface (11) is not smooth.

Singularities at the threshold of divergence of the constrained pendulum
We first consider the Ziegler pendulum under a pure conservative load (γ = 0). According to Eq. (8), the first and second divergence loads of the free pendulum are corresponds to the second-order work criterion [20].
so that the non-constrained conservative system is stable for 0 < p < p 1 , otherwise being statically unstable (divergence).
In this notation, the expression for the critical load (16) has the form Therefore, in the (α 1 , α 2 , p) space the threshold surface separating the domains of stability and divergence possesses the following parametric representation (ρ,φ) → ρ cos(φ), ρ sin(φ), which is a canonical form for the ruled surface known as the Plücker conoid of degree 2 [19].
Resolving Eq. (14) with respect to α 2 ,w efi n dt h ee x p r e s s i o n s for the generators The generators (20) are two straight lines that constitute the boundary between the domains of divergence and stability in the (α 1 , α 2 ) plane for a given p from the interval p 1 p p 2 ,s e e Fig. 3. Therefore, the divergence domain in the (α 1 , α 2 ) plane is an angle-shaped region. The angle between the generators is not zero when the load is within the interval p 1 < p < p 2 .H o w e v e r ,a tt h e end points of the stability interval the generators merge into one and the angle between them vanishes, see Fig. 3

(b).
In the (α 1 , α 2 , p) space the generators that are always perpendicular to the p axis sweep out the Plücker conoid surface that has a self-intersection along the interval p 1 < p < p 2 with the Whitney umbrella singularities [6,3] at the buckling loads p 1 and p 2 of the free pendulum, as shown in Fig. 3(a). We see that introducing kinematic constraints we can switch between the first and second buckling loads of the free system through the continuum of intermediate critical buckling loads parameterized by the angle φ or, equivalently, by the ratio α 2 /α 1 .
Remarkably, the minimal and maximal buckling loads of the constrained pendulum can be attained for non-trivial choices of the α 2 /α 1 -ratios that correspond to the two orthogonal directions in the (α 1 , α 2 ) plane. Therefore, we demonstrate that the problem of the optimal choice of the kinematic constraint inevitably leads to the non-smooth and non-convex merit functional, p d (α 1 , α 2 ), represented by the ruled surface (19). In general case such optimization problems can be solved numerically by the algorithms of optimization on Stiefel manifolds [34].
We note that the same ruled surface appears as the divergence threshold in the problem of destabilization of a hydrodynamically stable Couette-Taylor (CT) flow by an external magnetic field that is directed along the axis of rotation of the cylinders of the CTcell [16]. This indicates that the singular shape of the divergence threshold that we observe also in the case of the constrained Ziegler pendulum under the dead load, may be characteristic for a rather broad class of divergence instability problems both in solid and fluid mechanics.

Singular divergence instability threshold in the non-conservative case
Similarly, we can treat the case of a partial follower force in the situation when 0 < γ < 1. Denoting A = 10γ 2 − 14γ + 5, we find that the critical divergence load (11) varies between the two extrema p b p p u , where the minimum (21) corresponds in the (α 1 , α 2 ) plane to the direction while the maximum is attained when The critical divergence load (11) as a function of α 1 and α 2 is a ruled surface with the generators When γ = 0, the generators (25) coincide with those given by Eq. (20).
The divergence domain in the (α 1 , α 2 ) plane is an angle-shaped region. The angle between the generators is not zero when the plane where the critical divergence load of the constrained nonconservative system either exceeds the upper divergence load p 2 of the non-constrained one (stabilization by constraints) or is smaller than the lower divergence load p 1 of the non-constrained pendulum (destabilization by constraints), see Fig. 4(a). It is worth noting that Liu et al. [35], following the results of Rozvany and Mroz [37] and Olhoff and Akesson [36] had found that the optimal locations of internal supports for maximizing the buckling load of a column are at the nodal points of an appropriate higher-order buckling mode. Our study confirms this observation, because the maximal buckling load of the Ziegler pendulum with respect to the kinematic constraints corresponds to the higherorder (second) buckling mode of the non-constrained system while the minimal load is attained at the lower-order (first) one. This effect holds also in case of the non-conservative systems, which were not considered, however, in [35][36][37]. The possible destabilizing effect of additional constraints was not anticipated in these works, neither the fact that the optimal constraints in non-conservative systems are chosen by the second-order work criterion.

Conclusions
In various multiparameter stability problems of solid and fluid mechanics one frequently encounters with the effect of noncommuting limits of the instability threshold for a selected parameter when other parameters (that can correspond to both conservative and non-conservative forces) tend to zero. One of the complicated problems possessing the non-commuting limits of the divergence instability threshold is the Velikhov-Chandrasekhar paradox in the theory of magnetorotational instability. In order to understand this effect better it would be desirable to have a simple mechanical model that would demonstrate such a discontinuous switching between the critical buckling loads. In this work we have presented such a system, which is the Ziegler pendulum with kinematic constraints. The very idea to impose the constraints to this non-conservative mechanical system came, however, from the ma-terial instability problems of geomechanics. We have established that the divergence critical load as a function of the two parameters of a constraint is represented by a singular surface with the interval of self-intersection and two Whitney umbrella singular points at its ends. There exist directions in the plane of the constraint coefficients such that approaching the origin along them results in higher or lower values of the divergence threshold than in the case of the free pendulum. The extrema of the singular surface correspond to the values of the divergence load given by the second-order work criterion. The singular surface clearly illustrates how a switching between the first and second buckling loads of the free system happens when the constraints are varying. The simple but paradigmatic toy model studied in the Letter provides a link between the structural instability effects arising in seemingly disconnected applications.