Solution to the problem of Nicolai

Solution to the problem of Nicolai

Journal of Sound and Vibration 333 (2014) 1932–1944 Contents lists available at ScienceDirect Journal of Sound and Vibration journal homepage: www.e...

915KB Sizes 1 Downloads 58 Views

Journal of Sound and Vibration 333 (2014) 1932–1944

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Solution to the problem of Nicolai A.P. Seyranian a, A. Di Egidio b,n, A. Contento b, A. Luongo b a b

Institute of Mechanics, Lomonosov Moscow State University, Russia Dipartimento di Ingegneria Civile, Edile-Architettura e Ambientale, Università dell'Aquila, Italy

a r t i c l e in f o

abstract

Article history: Received 28 August 2013 Accepted 8 November 2013 Handling Editor: W. Lacarbonara Available online 16 December 2013

We consider the problem of Nicolai on dynamic stability of an elastic cantilever rod loaded by an axial compressive force P and a twisting tangential torque L in continuous formulation. The problem is to find the stability region for non-equal principal moments of inertia of the rod in the space of three parameters: P, L and the parameter α for the ratio of principal moments of inertia. New governing equations and boundary conditions, which form the basis for analytical and numerical studies, are derived. An important detail of this formulation is that the pre-twisting of the rod due to the torque L is taken into account. The singular point on the stability boundary at the critical Euler force PE is recognized and investigated in detail. For an elliptic cross-section of a uniform rod the stability region is found numerically with the use of the Galerkin method and the exact numerical approach. The obtained numerical results are compared with the analytical formulas of the asymptotic analysis. & 2013 Elsevier Ltd. All rights reserved.

1. Introduction In 1928 Evgenii L. Nicolai [1] formulated a problem of stability of an elastic rod under action of a tangential twisting torque and an axial force. Assuming that the rod has two equal principal moments of inertia, he found that, for the cantilever boundary conditions, there is no static form of equilibrium except the trivial (straight) one. Then he studied the stability of the trivial form of equilibrium using a dynamic method and came to the conclusion that it is unstable for an arbitrary small magnitude of the twisting torque. This effect is called the paradox of Nicolai. For the dynamic stability study he used a discrete, two degrees of freedom model with a lumped mass attached to the free end of a massless cantilever rod. In the same paper Nicolai introduced a small viscous damping and found that it has a stabilizing effect. In his next paper Nicolai [2] considered the stability problem of a cantilever rod with two different principal moments of inertia loaded only by a tangential twisting torque. He introduced geometrical imperfections related to non-equal principal moments of inertia and used the same discrete dynamic model for the stability study as in [1]. Then Nicolai came to the conclusion that the rod with non-equal moments of inertia possesses a finite non-zero critical twisting torque, i.e., geometrical imperfections have also a stabilizing effect. As an example he considered the rod with elliptic cross-section close to a circle. In the paper [2] Nicolai also promised to write a new paper to take into account the distributed mass of the rod and presence of the axial force but, as we know, this plan did not come true. Later the papers [1,2] were included in the volume of selected works of Nicolai [3]. As it is stated in the book by Bolotin [4], the works of Nicolai [1,2] were the first papers on stability problems with follower (non-conservative) loads. That was

n

Corresponding author. E-mail address: [email protected] (A. Di Egidio).

0022-460X/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jsv.2013.11.017

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

1933

beginning of the era of non-conservative stability problems! Later on the non-conservative stability problems were developed in the books [4–9] and survey papers [10,11]. The problem of Nicolai was reconsidered in the recent paper [12]. For linear vibrational systems of arbitrary degrees of freedom with potential forces having a multiple eigenfrequency it was shown that the addition of arbitrary small nonconservative positional forces generally destabilizes the system. The geometrical interpretation of this effect when multiplicity is equal to two is that the paradox of Nicolai is related to a conical singularity of the stability boundary. The formulas for the stabilizing effect of small internal and external damping and imperfections of the moments of inertia of the rod were derived. However, the pre-twisting effect was not taken into account. In this paper we consider the problem of Nicolai on the dynamic stability of an elastic cantilever rod loaded by a conservative axial compressive force P and twisting tangential torque L, which represents the non-conservative loading, in continuous formulation. The problem is to find the stability region for non-equal principal moments of inertia of the rod in the space of three parameters: P, L and the parameter α for the ratio of principal moments of inertia. It is known that for equal moments of inertia Jx ¼Jy the rod is unstable for any value of twisting torque La0 and all P oPE where PE is the Euler critical force. This effect constitutes the paradox of Nicolai [1]. The new problem is formulated as follows: to find the stability region for non-equal moments of inertia numerically and compare the results with the asymptotic analytical formulas derived for Jx ¼Jy. New governing equations and boundary conditions which form the basis for the analytical and numerical studies are derived. An important detail of this formulation is that the pre-twisting of the rod due to the torque L is taken into account. For an elliptic cross-section of the rod the stability region is found numerically with the use of the Galerkin method and the exact numerical approach. The numerical results are compared with the analytical formulas of the asymptotic analysis. The singular point on the stability region at P¼PE is recognized and investigated in detail. 2. Main equations We consider the problem of dynamic stability of an elastic cantilever rod loaded by a tangential torque L and axial force P. The rod is assumed to be of a uniform cross-section with the constant principal moments of inertia Jx and Jy. In Fig. 1 a local coordinate system x0,y0,z0 of the straight rod is shown in which the axes x0,y0 are in direction of the principal moments of inertia of the rod cross-section and the axis z0 coincides with the rod axis. For a deformed rod we introduce a local coordinate system x,y,z in which the axis z is oriented along normal to cross-section of the rod while x,y coincide with the principal moments of inertia. With the use of D'Alembert principle motion of the rod is governed by Kirchhoff's equations, see Bolotin [4] V ′x r 0 V y  qP þ f x ¼ 0 V ′y þ r 0 V x þpP þ f y ¼ 0 L′x r 0 Ly þ qL V y ¼ 0 L′y þ r 0 Lx  pL þV x ¼ 0

(1)

€ f y ¼  mv€ are the inertial forces, m is where Vx, Vy and Lx, Ly are the forces and bending moments, respectively, f x ¼  mu; the mass per unit length, u and v are the displacements in x and y directions, respectively, p and q are the curvatures, and r0 is the pre-twisting. Primes and dots denote differentiation with respect to coordinate s along the rod axis and time t, respectively. It is assumed that the pre-twisting due to the torque L can affect deformation of the rod and we take it into

Fig. 1. Elastic cantilever rod loaded by a tangential torque L and axial force P.

1934

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

account putting r0 ¼L/GJd with GJd being twisting stiffness. It means that the local coordinate system x0,y0,z0 is associated with the twisted rod in its straight position. We differentiate the last two equations (1) and substitute there the derivatives V ′x and V ′y from the first two equations. Then we obtain L″x r 0 L′y þ q′L þr 0 V x þ pP þf y ¼ 0 L″y þ r 0 L′x  p′L þr 0 V y þ qP  f x ¼ 0

(2)

now we substitute in these equations the forces Vx and Vy from the last two equations (1) and find L″x  r 0 L′y þ q′L r 0 ðL′y þr 0 Lx  pLÞ þ pP þf y ¼ 0 L″y þ r 0 L′x  p′L þ r 0 ðL′x r 0 Ly þ qLÞ þ qP f x ¼ 0

(3)

using relations for the bending moments Lx ¼ EJxp and Ly ¼EJyq, where E is Young's modulus, and Jx and Jy are the moments of inertia of cross-section of the rod, we obtain equations for the curvatures as EJ x p″ þ ðL  2r 0 EJ y Þq′þ ðP  r 20 EJ x þ r 0 LÞp þf y ¼ 0 EJ y q″ ðL  2r 0 EJ x Þp′ þ ðP r 20 EJ y þr 0 LÞq  f x ¼ 0

(4)

The kinematic relations yield, see Bolotin [4]  φ ¼ v′ þ r 0 u;

ψ ¼ u′ r 0 v

p ¼ φ′ r 0 ψ;

q ¼ ψ′þ r 0 φ

(5)

q ¼ u″ 2r 0 v′  r 20 u

(6)

From these equations we find the curvatures p ¼ v″  2r 0 u′ þr 20 v;

we assume that the twisting torque L and the pre-twisting r0 are small quantities of the first order and in the following calculations will neglect higher order terms. Substituting relations (6) into (4) and keeping only first order terms we get EJ y u′′′′ þðL  2r 0 EðJ x þ J y ÞÞv‴ þ Pu″ 2Pr 0 v0 þmu€ ¼ 0 EJ x v′′′′ ðL 2r 0 EðJ y þ J x ÞÞu‴ þ Pv″ þ 2Pr 0 u0 þmv€ ¼ 0

(7)

The boundary conditions at the clamped end of the rod mean that the displacements u, v and angles φ, ψ vanish. Using (5) we obtain s¼0 :

u ¼ v ¼ 0;

u′ ¼ v′ ¼ 0

(8)

the boundary conditions at the free end of the rod can be written as Lx ¼ Ly ¼ 0 and V y ¼  Pφ;

V x ¼ Pψ

(9)

which mean absence of bending moments and shear forces. With the use of Lx ¼EJxp, Ly ¼EJyq and (6) the first two conditions (9) yield s¼l :

u″  2r 0 v′ ¼ 0; v″ þ2r 0 u′ ¼ 0

(10)

to derive third and fourth boundary conditions (9) at the free end we express forces Vx, Vy from the last two equations (1) and obtain L′x  r 0 Ly þ qL þPφ ¼ 0;

L′y þ r 0 Lx  pL þPψ ¼ 0

(11)

then using (5), (6), and (10) and neglecting higher order terms we finally get s¼l :

EJ y u‴ þ Pðu′ r 0 vÞ ¼ 0;

EJ x v‴ þPðv′ þr 0 uÞ ¼ 0

(12)

conditions (10) and (12) mean absence of bending moments and shear forces at the free end, respectively. Combining boundary conditions at both ends of the rod (8), (10) and (12) we obtain s¼0 : s¼l :

u ¼ v ¼ 0;

u′ ¼ v′ ¼ 0;

u″  2r 0 v′ ¼ 0; v″ þ 2r 0 u′ ¼ 0;

EJ y u‴ þ Pðu′ r 0 vÞ ¼ 0;

EJ x v‴þ Pðv′ þr 0 uÞ ¼ 0

(13)

in derivation of linear equations (7) and (13) it was assumed that the displacements u and v, angles between new and initial coordinate systems φ and ψ, curvatures p and q, twisting torque L and pre-twisting r0 are small quantities.

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

1935

3. Asymptotic analysis With the substitution u ¼ UexpðiωtÞ; v ¼ VexpðiωtÞ, where ω is an eigenfrequency, we can write equations and boundary conditions (7) and (13) in a vector form as EJU″″ þξLGU″′ þ PU′′ þ 2Pr 0 GU′ μmU ¼ 0

(14)

Uð0Þ ¼ U′ð0Þ ¼ 0

(15)

U′′ðlÞ þ2r 0 GU′ðlÞ ¼ 0;

EJU‴ðlÞ þ PU′ðlÞ þ Pr 0 GUðlÞ ¼ 0

(16)

2

where μ¼ω and  U¼

U V

 ; J¼

Jy

0

0

Jx

!



; G¼

0

1

1

0

 ; ξ¼

2EðJ x þJ y Þ L  1; r 0 ¼ GJ d GJ d (17)

the system governed by Eqs. (14)–(16) is stable if all eigenvalues μ are positive and simple or semi-simple (with the number of eigenvectors equal to the algebraic multiplicity of μ), otherwise it is unstable, see Seyranian and Mailybaev [9]. We assume that in the initial (undisturbed) state the rod possesses two equal principal moments of inertia Jx ¼Jy, axial force P is less than the critical Euler value PoPE ¼π2EJx/4l2, and no tangential torque L is applied. Then we introduce the disturbance: small tangential torque L, variation of axial force δP, imperfections of the principal moments of inertia δJx, δJy and mass distribution δm. Our goal is to find the influence of the disturbance to the stability region. First, we consider the undisturbed rod with Jx ¼ Jy ¼J0, L¼0 and fixed PoPE. Then the governing equations and boundary conditions (14)–(16) take the form ″ EJU″″ 0 þ PU0  μmU0 ¼ 0

(18)

U0 ð0Þ ¼ U′0 ð0Þ ¼ 0

(19)

U″0 ðlÞ ¼ 0;

′ EJU″′ 0 ðlÞ þ PU0 ðlÞ ¼ 0

this system has double eigenvalues μ with two linearly independent eigenvectors     w 0 U1 ¼ ; U2 ¼ 0 w

(20)

(21)

where the scalar function w(s) satisfies the equations and boundary conditions EJ 0 w′′′′ þ Pw′′  μmw ¼ 0

(22)

wð0Þ ¼ w′ð0Þ ¼ 0

(23)

w″ðlÞ ¼ 0;

EJ 0 w‴ðlÞ þ Pw′ðlÞ ¼ 0

(24)

with the normality condition Z

l

m 0

w2 ðsÞds ¼ 1

with this equality the eigenvectors U1 and U2 satisfy the orthonormality condition Z l UTi Uj ds ¼ δij ; i; j ¼ 1; 2 m

(25)

(26)

0

where δij is Kronecker's delta. Now we disturb the initial system by small imperfections such as δm, δJx, δJy, variation δP and twisting torque L. Note that the pre-twisting r0 ¼ L/GJd was also assumed to be small. Due to these perturbations the eigenvalues and eigenvectors take variations δμ, δU, and we obtain the variational equation and boundary conditions as ″′ ″ ′ EδJU″″ 0 þξLGU0 þ δPU0 þ 2Pr 0 GU0 þ ½EJδU′′′′ þ PδU′′  μmδU  δmμU0  mδμU0 ¼ 0

δUð0Þ ¼ δU′ð0Þ ¼ 0;

(27) (28)

δU′′ðlÞ þ 2r 0 GU′0 ðlÞ ¼ 0; ′ EδJU″′ 0 ðlÞ þδPU0 ðlÞ þ Pr 0 GU0 ðlÞ þ EJδU‴ðlÞ þ PδU′ðlÞ ¼ 0

(29)

For the eigenvector U0 of the undisturbed system we have U0 ¼ γ 1 U1 þ γ 2 U2

(30)

1936

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

where U1 and U2 are the eigenvectors from (21) with unknown coefficients γ1 and γ2, see Seyranian and Mailybaev [9,12]. We substitute (30) into (27), multiply it by UT1 and UT2 from left side, and integrate over [0, l]. With the use of boundary conditions (28) and (29) we carry out integration by parts, and the term related to δU drops out. Then we obtain the system of two linear equations with respect to the coefficients γ1 and γ2 as ða11  δμÞγ 1 þ a12 γ 2 ¼ 0 a21 γ 1 þ ða22  δμÞγ 2 ¼ 0

(31)

where R l ′2 w″2 ds  μδm m  δP 0 w ds; Rl ′2 a12 ¼  a21 ¼ L 0 w‴wds ¼  L w 2ðlÞ

a11 ¼ EδJ y

Rl

0

a22 ¼ EδJ x

Rl 0

w″2 ds  μδm m  δP

Rl 0

w′2 ds

(32)

here we used integration by parts with the boundary and normality conditions (23)–(26). For the nontrivial solution of system (31) its determinant must vanish. Thus, we obtain the quadratic equation ðδμÞ2 ða11 þ a22 Þδμ þa11 a22 a12 a21 ¼ 0 The two solutions of this equation are δμ ¼

a11 þ a22 7

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ða11  a22 Þ2 þ4a12 a21 2

(33)

(34)

this means that, due to perturbation, the double positive eigenvalue μ splits into two different simple eigenvalues. If solutions of (34) are real the disturbed system remains stable for rather small perturbations, and if they are complex the system becomes unstable. Thus, the instability condition implies that the discriminant of Eq. (34) must be negative. Since a12 ¼ a21 we obtain the instability condition as ða11  a22 Þ2 o4a212

(35)

Using expressions (32) for the coefficients aij we obtain Z EðδJ y  δJ x Þ

l

!2 w″2 ds

0

oL2 w′4 ðlÞ

(36)

thus, the instability condition yields Rl EjδJ y δJ x j 0 w″2 ds jLj 4 w′2 ðlÞ

(37)

we see that the asymptotic formula for the instability region depends on difference between principal moments of inertia and deflection function w(s) but it does not depend on pre-twisting r0 ¼L/GJd, i.e., on twisting stiffness of the rod. If there are no changes of the moments of inertia δJx ¼δJy ¼ 0 then according to (37) the system is unstable for any La0. This is the paradox of Nicolai [1–4] according to which the rod with equal moments of inertia Jx ¼Jy is unstable for any infinitely small twisting torque L. Thus, we have shown that taking into account pre-twisting does not eliminate the paradox of Nicolai. 4. Singularity at the critical force P¼ PE At the critical force P¼ PE the undisturbed system (18)–(20) possesses the double zero eigenvalue μ¼0 of the undisturbed system (18)–(20). Due to perturbation it splits into two different simple eigenvalues δμ. In this case the stability conditions take the form ða11  a22 Þ2 4a212 4 0

(38)

a11 a22 þa212 4 0

(39)

a11 þ a22 40

(40)

condition (38) means that the two eigenvalues δμ of characteristic equation (33) are real while conditions (39) and (40), according to Vieta's theorem, imply that both δμ are positive. Inequalities (39) and (40) can be also represented as ða11  a22 Þ2 =4 4a212 4 a11 a22

(41)

a11 þ a22 40

(42)

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

1937

thus, according to equalities (32) we obtain the stability conditions in the form !2 Z l

EðδJ y  δJ x Þ Z 4 4 EδJ y

l

0

Z

w″2 ds δP

w″2 ds

0

l

4 L2 w′4 ðlÞ

!

Z

w′2 ds

EδJ x

0

l 0

w″2 ds  δP

Z

l

! w′2 ds

(43)

0

and Z EðδJ y þ δJ x Þ

l 0

w″2 ds  2δP

Z

l

w′2 ds4 0

(44)

0

relations (43) and (44) depend on the function w(s) which at P ¼PE is the first buckling mode of the cantilever rod loaded by the axial force, see Eqs. (22)–(25). Up to a scaling factor the first buckling mode is wðsÞ ¼ 1  cos ðπs=2lÞ

(45)

with the constants Z I1 ¼

0

l

w″2 ds ¼

π4 32l

Z ; 3

I2 ¼

l

0

w′2 ds ¼

π2 ; 8l

w′ðlÞ ¼

π 2l

(46)

Thus, the stability conditions (43) and (44) with relations (46) represent restrictions on the parameters: the tangential torque L, variation of axial force δP and imperfections of the principal moments of inertia δJx, δJy. 5. Elliptic cross-section of the rod Let us consider an elliptic cross-section of the rod with semi-axes a and b. For such cross-section we have 3

Jx ¼

πab ; 4

Jy ¼

3

πa3 b ; 4

Jd ¼

πa3 b

2

a2 þ b

;

m ¼ πabρ

(47)

where ρ is the mass density per unit volume. We introduce a new parameter α ¼b/a and assume that the semi-axis a remains constant while the semi-axis b is changing around b¼a (α¼1) where we have J0 ¼πa4/4 and m0 ¼πa2ρ. Note that according to (47) α2 ¼Jx/Jy. We will study the influence of the parameter α on the stability region of the rod along with the parameters L and P. For the sake of convenience we introduce non-dimensional quantities s sn ¼ ; l

Un ¼

U ; l

2

Pn ¼

Pl ; EJ 0

Ln ¼

Ll ; EJ 0

4

μn ¼

ω2 m 0 l α EJ 0

(48)

for the elliptic cross-section we have ξ ¼ ð1 þνÞð2 þ α2 þ α  2 Þ 1;

r0 ¼

2Lð1 þ νÞð1 þα2 Þ Eπa4 α3

(49)

here we have used the idendity G ¼E/2(1 þν). In terms of non-dimensional quantities and parameter α the equations and boundary conditions (14)–(16) with (48) and (49) take the form ″′ ″ ′ Jn U″″ n þ ξLn GUn þ P n Un þ 2ηP n Ln GU′Un  μn Un ¼ 0

(50)

Un ð0Þ ¼ U′n ð0Þ ¼ 0

(51)

U″n ð1Þ þ 2ηLn GU′n ð1Þ ¼ 0;

′ Jn U″′ n ð1Þ þ P n Un ð1Þ þηP n Ln GUn ð1Þ ¼ 0

where primes denote differentiation with respect to sn , and   α 0 ð1 þνÞðα  1 þ α  3 Þ ; η¼ Jn ¼ 3 2 0 α

(52)

(53)

thus, the system of differential equations (50)–(52) depends on three non-dimensional parameters, P n , Ln and α (Poisson's ratio ν is assumed to be fixed). Let us consider a small imperfection of the circular cross-section (a¼ b)which becomes elliptic due to the small imperfection δb of the semi-axis b. The principal moments of inertia take increments δJ x ¼

π3a3 δb ; 4

δJ y ¼

πa3 δb 4

(54)

1938

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

the asymptotic formula (37) in non-dimensional form reads Ln ¼

R1 2jδbj 0 w″2 dsn 4jLjl 4 πEa4 aw′2 ð1Þ

(55)

the right hand side of (55) implicitly depends on the parameter P through the eigenmode w(s). Note that the critical torque Ln does not depend on normalization condition of eigenmodes and the ratio in (55) can be replaced by δα¼δb/a. For P ¼0 the vibration eigenmodes for the cantilever rod (up to the normalizing constants) are wk ðsn Þ ¼ ðsinh βk þ sin βk Þðcosh βk sn  cos βk sn Þ ðcosh βk þ cos βk Þðsinh βk sn  sin βk sn Þ

(56)

where the constants βk satisfy the transcendental equation cosh βk cos βk ¼  1

(57)

with the solutions β1 ¼1.875, β2 ¼4.694, β3 ¼ 7.854, … Calculating the integrals in (55) for the first vibration mode we get jLn j 4 3:26jδαj

(58)

4jLjl 3:26jδbj 4 a πEa4

(59)

This estimate in dimensional form is

This critical value of the twisting torque can be compared with the result of Nicolai [2] who found 4jLjl 8jδbj 4 3a πEa4

(60)

Note that this estimate was obtained for the dynamic model with two degrees of freedom developed for P¼0, and it does not depend on Poisson's ratio ν. Thus, we see that the critical value of instability of the twisting torque obtained by Nicolai [2] is by 18.4 percent lower than our result (59) taken for the continuous dynamic model of the rod. Now let us study the stability region near the critical point P ¼PE. Then, according to inequalities (43) and (44) and constants (46) with Eqs. (48) and (54) we find the stability conditions as 

2   2  π 2 δα 3π δα π 2 δα δP n 4 L2n 4 δP n  4 4 4 π 2 δα  δP n 4 0 2

(61)

(62)

The two inequalities (61) represent conical singularities in two α, Ln and three-dimensional parameter space P n , Ln , α, respectively, while (62) is the semi-space. In Fig. 2 the stability region around the singular point W(P ¼PE, L ¼0, α¼ 1) analytically described by the three inequalities (61) and (62), is shown. The curved surfaces are a portion of the boundary of the conical stability region cut by the planes.

Fig. 2. First approximation of the stability region for the rod with the elliptic cross-section taken at the singular point W(P¼ PE, L ¼0, α ¼1).

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

1939

6. Numerical results for the elliptic cross-section Two types of numerical analyses are performed in this section to study the stability of the uniform rod with elliptic crosssection. The first one, named Approximate Numerical Approach, is based on the Galerkin discretization; the second one, called Exact Numerical Approach, attacks the main equations and the relevant boundary conditions of the problem with no approximation, see [13,14]. 6.1. Approximate numerical approach For the Galerkin discretization with (N þN) shape functions, we use the eigenmodes wk of the cantilever rod, with no compression load, P ¼0, taken from (56) and (57) N

U n ¼ ∑ ck wk ðsn Þ; k¼1

N

V n ¼ ∑ dk wk ðsn Þ

(63)

k¼1

where ck and dk are the unknown constants. These eigenmodes have been chosen since they do not depend on the compression load P and allow an easier study of stability of the system in the parameter space (Pn, Ln, α). Note that the eigenmodes wk satisfy the boundary conditions (51) and (52) for P* ¼0. It is convenient to write expansions (63) in a vector form ! ! N ck Un Un ¼ ∑ Ck wk ðsn Þ; Un ¼ ; Ck ¼ (64) dk Vn k¼1 we substitute (64) into Eq. (50) with the use of (53). Multiplying by the eigenmode wj and integrating over [0,1] we obtain ! Z 1 Z 1 Z 1 Z 1 Z 1 N ″″ ″′ ″ ′ wk wj dsn þξLn G wk wj dsn þ P n I wk wj dsn þ 2ηP n Ln G wk wj dsn  μn I wk wj dsn Ck ¼ 0; j ¼ 1; 2; 3; …; N; ∑ Jn k¼1

0

0

0

0

0

(65) where I is the unit matrix of dimension 2. This is a system of homogeneous linear algebraic equations on unknown coefficients Ck. To have a non-trivial solution the determinant of the system matrix must vanish R1 R1 R1 detðJn 0 w″k w″j dsn  ξLn G 0 w″j w′k dsn þ P n I 0 w″k wj dsn (66) R1 ′ R1 þ2ηP n Ln G 0 wk wj dsn  μn I 0 wk wj dsn Þ ¼ 0 where, with the use of the boundary conditions, an integration by parts has been performed. Due to the orthogonality of the eigenfunctions, we have Z 1 Z 1 wk wj dsn ¼ δjk w2k dsn ; (67) 0

0

where δjk is Kronecker's delta. Eq. (66) is the characteristic equation for eigenvalues μ. As we pointed out above, the system is stable if all eigenvalues μ are positive and simple or semi-simple, otherwise the system is unstable. In Fig. 3 different viewpoints of the stability region for the discretized system, obtained by using 4 eigenmodes for each displacement components (4 þ4), are shown. By performing a zoom in the parameter space, in Fig. 4(a) it is possible to observe the great similarity around the critical point W, with the first approximation of the region reported in Fig. 2. In calculations Poisson's ratio was always taken as ν¼0.3. To distinguish the different types of loss of stability some labels are added in Fig. 4(b), where a larger windows of the parameter space is adopted to show the critical boundaries. In particular across the cone, a static bifurcation (SB) occurs, across the other walls of the solid domain a dynamic bifurcation (DB) takes place, along the ‘edges’ a quadruple-zero (4Z) eigenvalue exists. To better illustrate the role of the three parameters Pn, Ln, α, several sections of the stability region have been plotted. In particular, in Fig. 5, sections at different values of Pn are shown. The regions where the system exhibits a stable behavior are the gray ones. In Fig. 5(a) it is possible to observe the Paradox of Nicolai. In fact, when a small torque Ln (very close to 0) acts on a rod with the circular cross-section (α¼1), even if the compression load is below the critical one, the system is unstable. The introduction of a small perturbation that makes the cross-section elliptic (αa1) destroys the Paradox, since it makes the system stable (see [2,4,6]). Over the Euler critical load, Pnπ2/4E2.46, the increasing of the compression load causes the decreasing of the stability region (Fig. 5(b and c)). However, it is always possible to find a couple of values Ln, α that makes the system stable. In Fig. 6 sections of the stability region at different values of the torque Ln are shown. By the comparison of the three graphs in Fig. 6, it is possible to observe that the presence of a small torque Ln (Fig. 6(b and c)) introduces a discontinuity in the stability region around the circular cross-section (α¼1) that it is not observable in the case with Ln ¼0 (Fig. 6(a)). The accuracy of the used numerical approach, based on the Galerkin discretization of the equations, is now discussed. In Fig. 7(a) the convergence to the exact solution described in Section 6.2, versus the number of shape functions used in the discretization approach, is reported. A zoom of the graph in Fig. 7(a) around the point (Ln,α)¼(0,1) is taken as a reference to compare the results obtained at various levels of approximation. It is possible to observe that the stability region decreases

1940

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

Fig. 3. Stability region for the elliptic cross-section of the rod: (a)–(c) different viewpoints (singular point W (Pn ¼ 2.46, Ln ¼ 0, α¼ 1)).

Fig. 4. Zoom of the stability region for the elliptic cross-section of the rod: (a) and (b) different zooms in the parameters space.

when the number of the shape functions increases. Moreover, it virtually does not change when more than (3 þ3) shape functions are taken, which is therefore expected to be the exact solution. Dashed lines refer to the first approximation of the boundaries. Results obtained in Section 4 for Pn ¼ 0 are shown in Fig. 7(b). The curve labeled with ‘First-order approx.’ refers to the case where pre-twist is taken into account and described by Eq. (59), while the thick curve in Fig. 7(b) represents the stability boundary of the discretized system accounting for (4 þ4) shape functions. When Pn is the Euler critical load, the ‘First-order approx.’ curve is perfectly tangent to the critical boundaries. In Fig. 7(a), for a Pn very close to the Euler load (PE ¼ 2.46), the critical boundaries and the ‘First-order approx.’ curve appear to be almost tangent to each other, while, when Pn becomes very different to the Euler load, the tangency of the two curves is lost (see Fig. 7(b)). 6.2. Exact numerical approach With Exact Numerical Approach we mean to find the exact solution of the eigenproblem of Eqs. (50)–(52). An exact solution can be formally written as ! g βsn Un ¼ De ; D ¼ ; (68) f

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

1941

1.10 1.05 1.00 0.95 0.90

-0.2

-0.1

0.0

0.1

0.2

1.10 1.05 1.00

W

0.95 0.90

-0.2

-0.1

0.0

0.1

0.2

-0.2

-0.1

0.0

0.1

0.2

1.10 1.05 1.00 0.95 0.90

Fig. 5. Sections of the stability region at different values Pn: (a) below the Euler load; (b) at the Euler load Pn ¼ 2.46; (c) over the Euler load.

with D as a vector of constants. By substituting Eq. (68) in Eq. (50) the following first eigenproblem is obtained: ðβ4 Jn þβ 3 ξLn G þβ2 P n I þ2βηP n Ln G  μn IÞ D ¼ 0 ;

(69)

where I is the unit matrix of dimension 2. By vanishing the determinant of the matrix detðβ4 Jn þ β3 ξLn G þβ2 P n I þ 2βηP n Ln G μn IÞ ¼ 0

(70)

it is possible to obtain the eight values of the spatial frequency βi ¼β(μ*,α,L*,P*) and the eight corresponding eigenvectors Di ¼Di(μn,α,Ln,Pn) (i¼1,…,8) as functions of the eigenvalue μn and of the other parameters. The eigenproblem (69) furnishes the solution Un ¼ k1 D1 eβ1 sn þ ⋯ þk8 D8 eβ8 sn

(71)

where ki (i¼1,…,8) are the unknown coefficients of the linear combination. By introducing Eq. (71) in the boundary conditions (51) and (52) the following second eigenproblem is then obtained: Fðμn ; α; Ln ; P n Þ k ¼ 0

(72)

where k is the vector of unknown coefficients ki, i ¼1,…,8 and F is a square matrix of dimension 8. The vanishing of the determinant of the matrix Mðμn ; α; Ln ; P n Þ : ¼ det½Fðμn ; α; Ln ; P n Þ ¼ 0

(73)

1942

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

4.0 3.0

P*

W

2.0 1.0 0.0

0.0

0.5

1.0

1.5

2.0

0.0

0.5

1.0

1.5

2.0

0.0

0.5

1.0

1.5

2.0

4.0 3.0

P*

2.0 1.0 0.0

4.0 3.0

P*

2.0 1.0 0.0

Fig. 6. Sections of the stability region at different values Ln: (a) Ln ¼ 0; (b) Ln ¼ 0.025; (c) Ln ¼0.05.

1.02

1.02

4+4 SF

3+3 = 4+4 SF

1.0

0.98

2+2 SF

First order approx.

P*=2.0 -0.04

-0.02

0

0.02

0.04

1.0

0.98

L*

First order approx.

P*=0 -0.04

-0.02

0

0.02

0.04

L*

Fig. 7. Accuracy of the solution: convergence of the solution: (a) Pn ¼2.0; (b) Pn ¼0.

furnishes the infinite number of the eigenvalues μ* which correspond to the infinite number of the eigenvectors k. As we stated above in Section 2, the stability boundaries are the locus where the double eigenvalues μ manifest themselves. Therefore, also the following condition for the double root μn must be satisfied: ∂Mðμn ; α; Ln ; P n Þ ¼0 δμn

(74)

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

1943

The procedure, due to the high order of the polynomial algebraic equation (70), can not be solved in closed form, then a numerical iterative procedure has been developed in Mathematicas environment. The objective was to describe the stability boundary, shown in Fig. 4, by one or more sections in the parameter plane (Ln,α). The parameter Pn was always taken as a fixed value. The stability boundaries in Fig. 7 are described by bullets representing the numerical exact solution obtained for a finite number of values of the parameter Ln. When a specific value of Ln is fixed, the system depends only on the other two parameters: the ratio α and the eigenvalue μn. Two nested loops, representing the body of the algorithm, able to obtain the single bullet, are briefly described in the following.

 In the first external loop the parameter α is changed in a specific range, close to the value furnished by the discretized approach.

 In the second internal loop the parameter μn is changed in a specific range close to the value furnished by the discrete approach; in correspondence of the couple (α,μn) it is possible to obtain the eight frequencies βi (from Eq. (70)) and the eight corresponding eigenvectors Di (i¼ 1,…,8) (by the solution of the eigenproblem given by Eq. (69)); the exact solution Un (Eq. (71)) can be now evaluated unless the unknown coefficient ki. Finally, the determinant of the matrix F of the boundary conditions (73) and its derivative (74) are numerically evaluated. When both are smaller than a pre-fixed tolerance, the procedure stops and the couple of parameters (α,μn) are registered as values that reach the goal.

To obtain a numerical exact solution with a sufficient precision the incremental steps of the parameters α and μn must be very small. To give an idea of the computational time to obtain the single bullet, by using a PC with a multicore CPU (six cores, 12 threads, clock 3.2 GHz, 32 GB RAM), the average time required to reach the goal is about 12 hours. At small values of the parameter Ln, the bullets shown in Fig. 7 are almost perfectly superimposed to the stability boundaries obtained by the discretized system using (3 þ3) or more shape functions. In particular when Pn a0 (see Fig. 7(a)) the ‘exact’ solution is very close to the discretized curves obtained using (4þ 4) shape functions. Also when Pn ¼0 (see Fig. 7 (b)), the exact solutions are almost perfectly superimposed to the discretized curves. The exact numerical results confirm the correctness of the results obtained by the approximate method and discussed in the previous sections.

7. Conclusion We considered the problem of Nicolai [2] on dynamic stability of an elastic cantilever rod with different principal moments of inertia loaded by an axial compressive force P and a twisting tangential torque L in the continuous formulation. In simplified form this problem was considered and solved by Nicolai [2] assuming that the dynamics of the rod can be described by a two degrees of freedom system with a massless rod and a lumped mass attached at the free end of the rod, and no axial compressive force is applied. In [2] Nicolai also promised to write a new paper to take into account the distributed mass of the rod and presence of the axial force but, as we know, this plan did not come true so far. We emphasize that the problem of Nicolai [1,2] was the first non-conservative stability problem with follower loads. We derived the governing equations and boundary conditions, which form the basis for analytical and numerical studies, and formulated the problem: to find the stability region for non-equal principal moments of inertia of the rod in the space of three parameters: P, L and the parameter α for the ratio of principal moments of inertia. An important detail of this formulation is that the pre-twisting of the rod due to the torque L is taken into account. The singular point on the stability boundary at the critical Euler force PE is recognized and investigated in detail. For an elliptic cross-section of a uniform rod the stability region is found numerically with the use of the Galerkin method and the exact numerical approach. The numerical results agree with the analytical formulas of the asymptotic analysis demonstrating the accuracy and efficiency of the obtained results. We compared the critical value of instability of the twisting torque obtained by Nicolai [2] for the dynamic model with two degrees of freedom taken at P¼ 0 and found that it is by 18.4 percent lower than our result (59) obtained for the continuous dynamic model of the rod. Geometrical interpretation of the obtained results is that the boundary of the stability region represents a conical surface in the three-dimensional parameter space cut by the planes. We confirmed and extended the results of Nicolai showing that the uniform cantilever column with equal principal moments of inertia loaded by an axial force loses stability under the action of an arbitrary small tangential torque, but it is stabilized for non-equal moments of inertia. The analytical and numerical results for the stability region obtained in this paper hold also for the case when the tangential twisting torque is substituted by the axial vertical twisting torque [4], since the corresponding eigenvalue problems, as it is shown in [12], are adjoint and hence, the stability problems are equivalent.

Conflict of interest statement None.

1944

A.P. Seyranian et al. / Journal of Sound and Vibration 333 (2014) 1932–1944

Acknowledgment The first author expresses his gratitude for the scientific visit to University of L'Aquila in June–July 2012 supported by Regional Funds PO FSE Abruzzo 2007–2013. References [1] E.L. Nicolai, On stability of the straight form of equilibrium of a column under axial force and torque, Izvestia Leningradskogo Politekhnicheskogo Instituta 31 (1928) 201–227. (in Russian). [2] E.L. Nicolai, On stability of a column under torsion, Vestnik Mekhaniki i Prikladnoi Matematiki 1 (1929) 41–58. (in Russian). [3] E.L. Nicolai, Works on Mechanics, Gos. Izd. Tekniko-Teoreticheskoi Literatuti, Moscow, 1955 (in Russian). [4] V.V. Bolotin, Nonconservative Problems of the Theory of Elastic Stability, Pergamon, New York, 1963. [5] H. Ziegler, Principles of Structural Stability, Blaisdell, Waltham, MA, 1968. [6] K. Huseyin, Vibrations and Stability of Multiple Parameter Systems, Sijthoff & Noordhoff, Alphen aan den Rijn, 1978. [7] H. Leipholz, Stability of Elastic Systems, Sijthoff & Noordhoff, Alphen aan den Rijn, 1980. [8] Ya.G. Panovko, I.I. Gubanova, Stability and Vibrations of Elastic Systems, Nauka, Moscow, 1987. [9] A.P. Seyranian, A.A. Mailybaev, Multiparameter Stability Theory With Mechanical Applications, World Scientific, New Jersey, 2003. [10] G. Herrmann, Stability of equilibrium of elastic systems subjected to noncorservative forces, Applied Mechanics Reviews 20 (1967) 103–108. [11] M.A. Langthjem, Y. Sugiyama, Dynamic stability of columns subjected to follower loads: a survey, Journal of Sound and Vibration 238 (2001) 809–851. [12] A.P. Seyranian, A.A. Mailybaev, Paradox of Nicolai and related effects, Zeitschrift für Angewandte Mathematik und Physik 62 (2011) 539–548. [13] R.L. Burden, J.D. Faires, Numerical Analysis, Thomson Higher Education, Belmont, USA, 2005. [14] E. Isaacson, H.B. Keller, Analysis of Numerical Methods, Dover Publication, New York, USA, 1994.