Frictional instabilities in orthotropic hollow cylinders

Frictional instabilities in orthotropic hollow cylinders

Computers and Structures 87 (2009) 1275–1286 Contents lists available at ScienceDirect Computers and Structures journal homepage: www.elsevier.com/l...

2MB Sizes 0 Downloads 30 Views

Computers and Structures 87 (2009) 1275–1286

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Frictional instabilities in orthotropic hollow cylinders A. Pinto da Costa *, M.A. Agwa 1 Departamento de Engenharia Civil e Arquitectura and ICIST, Instituto Superior Técnico, Universidade Técnica de Lisboa, Avenida Rovisco Pais, 1049-001 Lisboa, Portugal

a r t i c l e

i n f o

Article history: Received 17 June 2008 Accepted 3 August 2009 Available online 10 September 2009 Keywords: Coulomb friction Instability Complementarity eigenproblem Cylindrical obstacles Orthotropy Finite elements

a b s t r a c t This study discusses the implications of obstacle curvature, material orthotropy and friction to the occurrence of flutter and divergence instabilities of elastic media. The analytic solution for the divergence instability problem of a hollow cylinder when its principal directions of orthotropy make constant angles with the circumferential and radial directions is closely approximated by finite element solutions of a complementarity eigenproblem in which the coefficient of friction appears as an eigenvalue, quadratically due to the obstacle curvature. Mechanical systems with a finite number of degrees of freedom illustrate some aspects of the onset of instability. In particular, the presented study shows that, under certain conditions, the increase of the obstacle convexity may favour stability. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction The study of divergence and flutter instabilities in frictional contact systems has been receiving some attention. This text represents further progress in these types of Liapunov instabilities of elastic structures in frictional contact with cylindrical obstacles and may be considered an extension of previous works such as [1,22] where flutter and divergence instabilities in infinite layers were studied in some detail and [19,28] where the effects of the obstacle curvature in stability analyses of frictional contact systems were considered. Interfaces between tools or between an insulation device and a container are often curved. The object of this paper is the numerical computation of instabilities, divergence and flutter, in elastic media contacting with curved frictional surfaces. Analytic or finite element models of hollow cylinders and lumped mass systems are used to illustrate these phenomena. Whenever they exist, analytic solutions are used to check the finite element model and to perform parametric studies. The search of unstable solutions is connected with a non-classical eigenproblem which is equivalent to a set of classical generalized eigenproblems whose number grows exponentially with the contact nodes in impending slip and in grazing contact. Some time ago several authors devoted their attention to eigenvalue variational problems motivated by the study of the friction-

* Corresponding author. Tel.: +351 21 841 8408; fax: +351 21 849 7650. E-mail addresses: [email protected] (A. Pinto da Costa), [email protected] (M.A. Agwa). 1 On leave from the Department of Mechanical Engineering, Zagazig University, Egypt. 0045-7949/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2009.08.007

less constrained buckling or just for mathematical purposes [20,21,24,32]. It seems that Ref. [18] is the first to formulate the problem of the stability of a discretized structure in equilibrium in the presence of a frictional obstacle in terms of an algebraic eigenvalue variational (or inclusion) problem, which is equivalent to an eigenvalue complementarity problem with a non-symmetric matrix operator [26,28]. Moirot [22] dedicated attention to frictional instabilities of the divergence and flutter types; special emphasis was given to discretized structures with many degrees of freedom with industrial applications in the automotive industry. More recently Ionescu et al. [11,8] have also done a research effort on frictional instabilities; they deal with discrete nonlinear eigenvalue analyses and propose an algorithm that, in case of convergence, generates a sequence whose limit is a solution of the eigenvalue problem (a divergence instability mode). The complementarity eigenvalue behind these instability analyses problem has also been attracting mathematicians that have developed algorithms for its numerical resolution [12,30,31]. This paper deals with anisotropic (orthotropic) linear elastic hollow cylinders in contact with cylindrical surfaces where the friction law is isotropic. Some authors have addressed anisotropic friction laws, mainly in the three dimensional context [9]. The paper proceeds as follows: Section 2 describes the continuum and the (directionally) linearized problem leading to the instability analysis. The analytic solution and the parametric study for the problem of divergence instability are presented in Section 3. Section 4 is devoted to the finite element based complementarity formulations that are relevant for the numerical computation of divergence and flutter instabilities. Some lumped mass models worked out in Section 5 shed light on the mechanisms of instability and in Section 6 we conclude.

1276

A. Pinto da Costa, M.A. Agwa / Computers and Structures 87 (2009) 1275–1286

nent e33 . The transformation of the constitutive equation into the lo-

2. Statement of the problem

cal reference frame ðer ; eu ; e3 Þ gives r ¼ De where r ¼ TT rp ; D ¼ 2.1. Description of the continuum The present study is mainly concerned with orthotropic elastic materials confined between two cylindrical surfaces of radius Ri and Ro ; Fig. 1a represents schematically a sector of the domain. The particles of the outer surface have kinematic boundary conditions in both directions while the particles of the inner surface satisfy frictional contact conditions. The unit inward normal and tangent vectors to the cylindrical obstacle surface of radius Ri are denoted by n and t, respectively, so that n  t ¼ e3 . The displacement of a particle is denoted by vector u with physical components ður ; uu Þ which are sufficiently smooth functions of the polar coordinates ðr; uÞ. The stress vector ri on the cylindrical surface of radius Ri is given by ri ¼ rrr n þ rru t, where rrr and rru are physical components of the Cauchy stress tensor in the polar coordinate system ðr; uÞ. The hollow cylinder is made of a homogeneous orthotropic transversely isotropic linear elastic material in a state of plane strain. The principal directions of orthotropy are locally defined by p1 ; p2 and e3 (see Fig. 1a). The two arrangements of the principal directions of orthotropy considered here are also portrayed in Fig. 1: (b) pattern A, characterized by a constant angle h between the principal direction p1 and the circumferential direction; (c) pattern B, characterized by two families of mutually orthogonal principal directions of orthotropy with constant orientations with respect to a fixed cartesian reference system. At each point, the plane of isotropy of the material is plane p2  e3 . In the principal orthotropy reference system ðp1 ; p2 ; e3 Þ the constitutive equation may be written as rp ¼ Dp ep , where rp ¼ frp11 ; rp22 ; rp33 ; rp12 ; rp23 ; rp13 gT ; ep ¼ fep11 ; ep22 ; ep33 ; cp12 ; cp23 ; cp13 gT ; Dp ¼ ðCp Þ1 with

2

1 E1

6 m12 6 6 E1 6 m12 6 E 1 p :6 C ¼6 6 0 6 6 6 0 4 0

 mE121

 mE121

0

0

1 E2

 Em2

0

0

 Em2

1 E2

0

0

0

0

1 G12

0

0

0

0

2ð1þmÞ E2

0

0

0

0

0

3

7 0 7 7 7 0 7 7 7: 0 7 7 7 0 7 5

ð1Þ

1 G12

TT Dp T and ep ¼ Te. Matrix T is a standard transformation matrix (see Chapter 7 of [6] and Appendix A). An equilibrium state or a steady sliding state of the hollow cylinder is defined by the pair ðu0 ; r0i Þ of a displacement field and a contact stress field on the inner circumference, satisfying: (1) the equilibrium in the interior of the domain

1 r

r0rr;r þ r0ru;u þ

r0rr  r0uu r

2r 1 0 r þ r0ru;r þ r uu;u r

0 ru

þ fr ¼ 0;

þ fu ¼ 0;

ð2Þ ð3Þ

where fr and fu denote the physical components of the volume forces; (2) the generalized Hooke’s law described before; (3) the kinematic boundary conditions on the outer circumference

u0r ðRo Þ ¼ u0u ðRo Þ ¼ 0;

ð4Þ

(4) the frictional contact conditions on the inner circumference

u0r ðRi Þ ¼ 0;

r0ru ¼ lr0rr ;

ð5Þ

corresponding to a counterclockwise state of impending slip of the inner circumference. The coefficient of friction is denoted by l. In order to keep the notation as simple as possible, the dependence of the displacement and stress fields on the coordinates ðr; uÞ was sometimes omitted in this section and in the next section. 2.2. The linearized problem The proof that an equilibrium or steady sliding state is unstable in the sense of Liapunov implies to search for dynamic solutions initiated by an arbitrarily small perturbation, such that, no matter how small the perturbation is, the dynamic solution leaves the neighbourhood of the equilibrium state solution in finite time. Consequently, we study the existence of non-oscillatory (divergence) or oscillatory (flutter) exponentially growing smooth dynamic solutions ðduðtÞ; dri ðtÞÞ starting from initial conditions arbitrarily close to the equilibrium or the steady sliding state ðu0 ; r0i Þ and satisfying the following set of conditions:  The equations of motion in the interior of the hollow cylinder

Matrix Dp has five independent elastic constants; for a state of plane strain there are only four non-dimensional elastic parameters   E1 G12 ; ; m12 ; m while for a state of plane stress the non-dimensional E2 E2   elastic parameters are three EE12 ; GE122 ; m12 . Note however that, for a

drrr ðtÞ  druu ðtÞ 1 €r ðtÞ; drrr;r ðtÞ þ drru;u ðtÞ þ ¼ qdu r r 2drru ðtÞ 1 € u ðtÞ; druu;u ðtÞ þ drru;r ðtÞ þ ¼ qdu r r

state of plane stress, the Poisson’s ratio m relative to plane p2  e3 is needed for the ‘‘a posteriori” computation of the strain compo-

where q stands for mass per unit volume and ðÞ denotes differentiation with respect to time t.

ð6Þ ð7Þ

Fig. 1. (a) Sector of a hollow cylinder with a curved frictional contact. Patterns of the principal directions of orthotropy: (b) each family of principal directions of orthotropy makes the same angle h or h þ p2 with the circumferential direction; (c) constant orientation.

1277

A. Pinto da Costa, M.A. Agwa / Computers and Structures 87 (2009) 1275–1286

 The linear elastic constitutive equations presented in Section 2.1.  The kinematic boundary conditions on the outer circumference in both directions:

in agreement with (8), (9.2) and (10), as well as (12). Upon substituting (12) in the constitutive equation dr ¼ Dde and the result in (6) and (7), the following system of two second order linear ordinary differential equations involving F and G is obtained:

dur ðtÞ ¼ duu ðtÞ  0:



ð8Þ

½dri t ðtÞ þ l

¼ 0;

dur ðtÞ  0;

D24



F 00 ðrÞ



D44 G00 ðrÞ  1 D11 D14 þ 2 r D14 D44   FðrÞ 2 ¼ qk : GðrÞ

D24

 On the inner circumference, the frictional sliding condition in the tangent direction and the kinematic condition in the normal direction:

½dri n ðtÞSignð½r0i t Þ

D22

ð9Þ

½r0i t

where a positive sign of corresponds to a slip in the direct trigonometric sense (counterclockwise). _  Initial conditions duð0Þ  0 and duð0Þ – 0 such that

 D22 1 r D14 þ 2D24   FðrÞ þ

D14



D44

F 0 ðrÞ



G0 ðrÞ

GðrÞ ð15Þ

for all the particles on the circumference of radius Ri . For the study of the stability of a steady sliding state, the initial relative velocity of the contact particles with respect to the moving obstacle must be positive:

Note that the rate of growth of a divergence instability mode is dependent on the (positive) value of eigenvalue k issued from eigenproblem (14) and (15). Alternatively, when the hollow cylinder is submitted to a state of plane stress, conditions r33 ¼ r13 ¼ r23 ¼ 0 must be added. The second row of Table 1 shows the coefficients that ought to be used in system (15) instead of those in the first row when the hollow cylinder is submitted to a state of plane stress.

du_ u ð0Þ  xobst Ri > 0;

3.2. The computation of the stability–divergence instability border

du_ u ð0Þ > 0;

ð10Þ

ð11Þ

where xobst denotes the (possibly) non-vanishing angular velocity of the inner cylinder. In (10) and (11) it is assumed, without loss of generality, that the state of sliding with respect to the cylindrical obstacle is in the counterclockwise sense; negative signs might as well be considered, yielding then symmetric modes.

For a small slide of a volume element on the frictional cylindrical surface r ¼ Ri , the contact stresses variations relevant for the frictional contact conditions depend not only on magnitude changes but also on changes of orientation of the unit vectors n and t:

dri ¼ ðdrrr þ lvduu r0rr Þn þ ðdrru þ vduu r0rr Þt;

The constitutive equation for stress and strain variations written in the cylindrical coordinate system (r; u; x3 ) is dr ¼ TT Dp Tde. See Appendix A for more details on matrix TT Dp T.

where (5.2) was taken into account. All quantities in (16) are evaluated at r ¼ Ri ; v ¼ þ R1i denotes the obstacle curvature (positive for a convex obstacle). The equation yielding the values of the coefficient of friction at the onset of instability, derived from (9.1) and (16), is then quadratic in l,

3. Analytical solution for divergence instability



In this section we consider a pattern of principal directions of orthotropy of type A so that the stress state at equilibrium for constant r coordinate is homogenous.

Assuming that the displacement vectors of each pair of points with the same r are the same, by means of separation of variables we search solutions to (6)–((10) or (11)) in the form

dur ðr; tÞ duu ðr; tÞ



¼ aðtÞ



 FðrÞ ; GðrÞ

aðtÞ ¼

In this section, we investigate how the value of the coefficient of friction at the onset of divergence instability varies with material parameters, the aspect ratio RRoi , a non-dimensional normal stress and the orientation of the principal directions of orthotropy h. Without loss of generality, the computations were done for an inside radius Ri ¼ 250 mm. Table 2 shows the first numerical experiment, with an isotropic material. The (high) numerical values for l

ð12Þ

að0Þ coshðktÞ þ a_ ð0Þ sinhðktÞ; if k > 0; k að0Þ þ a_ ð0Þt; if k ¼ 0

ð13Þ

Table 2 High values of the coefficient of friction that would correspond to the onset of divergence instability yielded by a linear elastic analysis. Data: Isotropic material, r0rr ¼ 0:005; m ¼ 0:4, several aspect ratios RRoi . E2

in which að0Þ P 0 and a_ ð0Þ P 0 are sufficiently small; solutions of this kind diverge from the equilibrium state. Functions FðrÞ and GðrÞ are assumed to be two times continuously differentiable and to satisfy the following boundary conditions

FðRo Þ ¼ GðRo Þ ¼ FðRi Þ ¼ 0;

GðRi Þ ¼ G > 0;

ð17Þ

3.3. Parametric study based on the analytic solution

where the time dependence is governed by the non-decreasing function

(



vduu r0rr l2 þ ðdrrr Þl þ ðdrru þ vduu r0rr Þ ¼ 0:

The critical value of l for divergence is given by the smallest positive real root of the previous equation in which one uses the analytic solutions of (12)–(15) for duu ; drrr and drru obtained with the symbolic manipulator program MAPLE [17].

3.1. The search for non-oscillatory exponentially growing solutions



ð16Þ

ð14Þ

Ro Ri

4 3

2

4

l

18.0242669

13.7650351

12.3036967

Table 1 Table of correspondences for the elastic coefficients in Eq. (15) between states of plane strain and stress. Term

11

22

44

12

14

24

Pl. Strain Pl. Stress

D11

D22

D44

D12

D14

D24

D23 D12  D13 D33

D34 D14  D13 D33

D34 D24  D23 D33

D2

D11  D13 33

D2

D22  D23 33

D2

D44  D34 33

1278

A. Pinto da Costa, M.A. Agwa / Computers and Structures 87 (2009) 1275–1286

presented in this table show that, for isotropic materials, one cannot extract conclusions on divergence instability since such high values of l for the onset of instability correspond to unacceptably high values of strains for a linear analysis. For orthotropic materials, the relationship between the critical value of l and h, for several values of EE12 and fixed values of m12 ; m and GE122 is displayed in Fig. 2. It is observed that l ! þ1 when the principal directions are aligned with the directions normal and tangent to the obstacle (for h ? 0, 90° or 180°) as in the infinite layer [1]. Two values of the non-dimensional shear modulus GE122 were considered, the smaller value of GE122 yielding smaller values of l at the onset of instability. The graphs in Fig. 3 show the lines in the friction-curvature space that solve Eq. (17) and separate the stable (below) and the unstable (above) regions for several orientations of the principal directions of orthotropy (type A) of a homogeneous hollow cylinder. Every curve contains a plateau (for sufficiently large curvatures of the obstacle) preceded by an increasing or decreasing part, depending on the inclination of the principal direction of orthotropy: one has an initial increasing (decreasing) section when the stiffer principal direction and the sense of impending slip make an obtuse (acute) angle. 4. The finite element model

of two subvectors: a vanishing part rF relative to the nodes that are not candidates to contact and a potentially non-vanishing part including the reactions of the nodes that may establish contact. The contact candidate nodes have their displacement and reaction expressed in terms of the directions that are normal and tangential to the obstacle. The smooth dynamic evolution of the previously described system beginning at the equilibrium or steady sliding state ðu0 ; r0 Þ is governed by the following set of incremental Lagrange’s equations and conditions involving static and kinematic variables

€ ðtÞ þ KduðtÞ ¼ drðtÞ; Mdu drF ðtÞ ¼ 0;

ð18Þ ð19Þ

p 0 P dupn ðtÞ ? r0p n þ dr n ðtÞ 6 0;

ð20Þ

r0p t

ð21Þ

þ

dr pt ðtÞ

2l

ðr 0p n

þ

drpn ðtÞÞSignðdu_ pt ðtÞÞ;

where M and K stand for the N  N symmetric positive definite mass and stiffness matrices, dð Þ is a small variation of ( ), subscripts F, n and t denote the degrees of freedom that, respectively, are free from any kinematic constraint or are normal or tangent to the obstacle; superscript p is the label identifying the contact candidate particles. In (21), Sign denotes the multi-application defined as follows: SignðxÞ ¼ jxjx if x – 0 or SignðxÞ 2 ½1; 1 if x ¼ 0. The force term does not appear in (18) since forces are to be considered constant for instability studies in the sense of Liapunov.

4.1. Formulation 4.2. Divergence instability We turn, now, to a discretization of the (continuum) solid. The N generalized coordinates are the nodal displacements, grouped in vector u. The vector r of reactions from the obstacle is composed

It may be shown [18] that dynamic solutions to (18)–(21) of the form

4

4 16

0.5 3

μ

3

2

4

μ

2

1 1

8

4

0.5

1

150

0 0

8

16 0 0

30

60

90

θ

120

0.5

2

4

2 2

16

8

1

180

16 30

60

90

θ

120

150

180

Fig. 2. Coefficient of friction at the onset of divergence instability for an orthotropic hollow cylinder of type A in function of h. The values of the ratio EE12 are indicated near the r 0 curves. Data: m12 ¼ 0; m ¼ 0; RRoi ¼ 43 ; E2rr ¼ 0:005.

3

3

−30

80 2

2

μ

μ 1

0 0

1

60

0.02

0.04

χ

0.06

0.08

0.1

0 0

−10

0.02

0.04

χ

0.06

0.08

0.1

Fig. 3. Values of l at the onset of divergence instability for an orthotropic hollow cylinder. The values of h (in degrees) are indicated near each curve. Data: r0 E1 ¼ 4; GE122 ¼ 0:1; m12 ¼ 0:4; m ¼ 0:3; E2rr ¼ 0:005 and R0 ¼ 500 mm. E2

1279

A. Pinto da Costa, M.A. Agwa / Computers and Structures 87 (2009) 1275–1286

duðtÞ ¼ aðtÞv ;

drðtÞ ¼ aðtÞw;

ð22Þ

define escape directions from equilibrium, where vectors v and w are constant right admissible directions of displacement rate and reaction rate, respectively; recall that aðtÞ is the non-decreasing function given by (13). It is now convenient to decompose the set of degrees of freedom at the equilibrium state ðu0 ; r0 Þ for which the instability problem is formulated. In the sequel, the following meanings for subscripts that label kinematic and static quantities hold: f: degrees of freedom without any kinematic constraint in the current equilibrium state; they are the degrees of freedom that have no kinematic constraint in all instants (F) and the degrees of freedom of the contact candidate particles that are currently out of contact (free); d: degrees of freedom of the particles that are currently with the reaction in the interior of the friction cone and will be blocked in the near future (vanishing displacement rate) sn(st): degrees of freedom normal (tangent) to the obstacle of the particles that are currently in a state of impending slip. Assuming that there are no particles in contact with zero reaction (no particles in grazing contact), a necessary and sufficient condition for (18)–(21) to have a solution of the type (22) is the existence of a k P 0 and vectors v ; w 2 RN with v – 0 solving the following complementarity eigenproblem [18,28]

ðk2 M þ KÞv ¼ w;

ð23Þ

wf ¼ 0; v d ¼ 0; v v pst Signðr0p þ lwpsn 6 0; st Þ 6 0; p 0p p ½v pst Signðr 0p st Þ½wst Signðr st Þ þ lwsn  ¼ 0: p sn ¼ 0; p wst Signðr 0p st Þ

ð24Þ ð25Þ ð26Þ

The previous eigenproblem is not a common linear eigenproblem since the right hand side in (23) is not the zero vector and the components of the eigenpair (v, w) have to satisfy inequalities (25) and complementarity conditions (26). The vanishing entities in (24) are the rate of change of the reactions on the degrees of freedom that are not in contact with the obstacle wf , the velocity of the particles that are blocked due to contact reactions strictly inside the friction cone v d , and the velocities normal to the obstacle of impending slip particles (with non-vanishing reactions on the border of the friction cone) v psn . 4.2.1. Two relevant non-classical eigenproblems The elimination of the two degrees of freedom of each blocked particle (d), the elimination of the degree of freedom normal to the obstacle of each particle in impending slip (sn) and the change of variables





xf xst



 ¼

vf Svst

 ;





yf yst



 ¼

 0 ; ðSwst þ lwsn Þ ð27Þ

diagðSignðr0p st ÞÞ,

where S ¼ lead to the following Mixed Complementarity EigenProblem in terms of k2 (MCEIP-k2 ): Compute k P 0 and ðx; yÞ with x – 0, such that

ðk2 M þ K Þx ¼ y; yf ¼ 0;

ð28Þ ð29Þ

0 6 xst ? yst P 0:

ð30Þ

The vector inequalities in (30) are to be interpreted in a componentwise way. The matrix pencils M (linear in l) and K (quadratic in l) are generally non-symmetric. The structures of these mass and stiffness matrix pencils are:

M ¼ M0 þ lM1 ¼



Mf ;f SMst;f





Mf ;st S 0 þl SMst;st S Msn;f

K ¼ K0 þ lK1 þ l2 K2    Kf ;st S 0 Kf ;f þl ¼ SKst;f SKst;st S þ C Ksn;f

0 Ksn;st S



0 Msn;st S

þ l2



0

 ;

0

0 C

ð31Þ  : ð32Þ p 0p rsn Þ,

The curvature terms, incorporated in submatrix C ¼ diagðv affect the diagonal entries of the constant and the quadratic terms of the stiffness matrix pencil connected with the contact nodes in a state of impending slip. After the above mentioned elimination of degrees of freedom, the importance of this geometric stiffness effect, due to the change of orientation of the obstacle reactions when sliding occurs, is quantified by the term vr 0n ð1 þ l2 Þ, for each impending slip node. Note that r 0n 6 0, while v > 0 ðv < 0Þ for a convex (concave) obstacle, so that a convex (concave) obstacle, in general, favours instability (stability) (under particular circumstances this may not be true: recall the curves in Fig. 3a for which the transition to instability occurs for larger values of l as v (convex obstacle) increases). Klarbring and Björkman [4,13,15] had already considered the effect of curvature in the tangent stiffness matrix for the purpose of computing equilibrium states or for analyzing the stability of frictionless equilibrium states. Notice that, due to the change of variables (27), the variables of problem (28)–(30) belong to cones with a much simpler structure than the cones defined by (25, 26). However, the matrices in (28) are no longer symmetric, as opposed to those in (23). A pertinent question related with problem MCEIP-k2 is how to compute the values of l corresponding to the beginning of instability (vanishing k). The answer may be given by the Mixed Complementarity EigenProblem in terms of l2 ðMCEIP-l2 Þ: Compute l P 0 and ðx; yÞ with x – 0, such that

ðK0 þ lK1 þ l2 K2 Þx ¼ y;

ð33Þ

yf ¼ 0;

ð34Þ

0 6 xst ? yst P 0:

ð35Þ

In the previous problem the coefficient of friction plays the role of an eigenvalue. The values of l solving (33)–(35) are homologous, for discrete structures, to those solving (17). It is possible to rewrite the nonlinear eigenproblems MCEIP-k2 and MCEIP-l2 as Mixed nonlinear Complementarity Problems (MCPs), which are then solved with algorithm PATH [7], part of a modelling code for mathematical programming and optimization [2]. This algorithm is based on a non-classic Newton method that builds a piecewise linear path and does a line search along it. The transformation into a MCP is achieved by taking the eigenvalue (k2 or l) as an additional non-negative variable, which is complementary to another non-negative variable, used in a normalizing constraint equation that involves the kinematic components of the eigenvector. In a manner similar to that used in [25] it is possible to show that solutions to problem MCP-l2 also solve problem MCEIP-l2 .

4.2.2. Homogeneous hollow cylinders In Table 3, a comparison is made between the value of the coefficient of friction at the onset of divergence instability computed from the characteristic equation (17) and by the finite element method (33)–(35) for three different values of the aspect ratio RRo . i In all the three cases the total number of nodes is 3024 and the number of contact nodes is 144. It is observed that the value of l increases with the increase of RRoi , since h is between 0° and 90° (recall Fig. 3a).

1280

A. Pinto da Costa, M.A. Agwa / Computers and Structures 87 (2009) 1275–1286

Table 3 The effect of the aspect ratio RRo on the value of the coefficient of friction at the onset of i divergence instability for an orthotropic hollow cylinder (type A). Data: EE12 ¼ 4; GE122 ¼ r0 0:1; m12 ¼ 0:4; m ¼ 0:3; h ¼ 60 ; E2rr ¼ 0:005. Ro Ri

4 3

2

4

Analytic FEM

0.880056 0.880707

1.108571 1.109571

1.510201 1.514582

Fig. 4 illustrates the contact velocity patterns of the instability modes for orthotropic hollow cylinders of type B. As might be expected, the modes are non-uniform in the circumferential direction (the tangent velocities are represented in the radial direction). For the particular set of data used in Fig. 4, an increase of the aspect ratio RRoi corresponded to an increase of the critical value of the l (compare (a) with (b)) while an increase of EE12 corresponded to a decrease of l (compare (b) with (c)). Note that the maximum velocities occur in diametrally opposed nodes that do not make 45° with the horizontal: they are not aligned with the stiffer principal direction of orthotropy. 4.2.3. Searching for non-all-slip instability modes Many combinations of the non-dimensional parameters were used without success in the search of instability modes involving stiction of a part of the contact boundary for homogeneous (orthotropic or not) materials. Such non-all-slip modes occur however for non-homogeneous hollow cylinders as shown in Fig. 5; the hollow cylinder is composed of two homogeneous half cylinders connected at the vertical axis. The nodes that do not slide during the instability are indicated with the label sd while the nodes that slide are labeled with ss. In graphs (a) and (b) of Fig. 5, (total number of nodes, number of contact nodes) is (792, 72) and (3024, 144),

Fig. 5. Solutions of MCEIP-l2 for a non-homogeneous hollow cylinder: the left half is an orthotropic material with h ¼ þ80 (type A), E1 ¼ 8; E2 ¼ 1; m12 ¼ 0:4; m ¼ 0:3; G12 ¼ 0:1; the right half is an orthotropic material with h ¼ þ30 0 (type A), E1 ¼ 6; E2 ¼ 1; m12 ¼ 0:1; m ¼ 0:1; G12 ¼ 0:1; Er2rr ¼ 0:005.

respectively. Fig. 6 depicts the polar representation of the contact variables of problem MCEIP-l2 for the more refined mesh. Non-all-slip divergence instability modes in the context of plane obstacles were obtained in [29,8]. 4.3. Flutter instability of the steady sliding state The study of the occurrence of flutter instability is done for steady sliding equilibrium states with strictly negative normal reactions, i.e., conditions un ¼ rn ¼ 0; r t ¼ lrn ; v st ¼ ðdu_ u ð0Þ xobst Ri Þ < 0 (v st is the relative velocity between each contact node and the obstacle), wst þ lwsn ¼ 0 and v st ðwst þ lwsn Þ ¼ 0 hold. Flutter instability involves the search of non-trivial solutions of the classical generalized eigenproblem ðk2 M þ K Þx ¼ 0 corre-

Fig. 4. Top: divergence instability modes. Bottom: polar representation of the tangent contact velocities. Orthotropic material with u ¼ 45 (type B), 0 m12 ¼ 0:4; m ¼ 0:3; GE122 ¼ 0:1; Er2rr ¼ 0:005.

1281

A. Pinto da Costa, M.A. Agwa / Computers and Structures 87 (2009) 1275–1286

Fig. 6. Polar representation of the absolute value of the tangential velocities ðxst Þ and the complementarity variable measuring the evolution of the reaction rates with respect to the border of the friction cone ðyst Þ, at the onset of divergence instability, for the mode in Fig. 5b.

sponding to k 2 C with ReðkÞ > 0; it occurs for coefficients of friction l such that the matrix pencil k2 M þ K is non-symmetrizable, the system being called circulatory [10,34]. Flutter instability occurs when, for a certain value of the coefficient of friction (possibly zero), an eigenvalue of algebraic multiplicity 2 (a double eigenvalue) has geometric multiplicity equal to 1 (does not correspond to two linearly independent eigenvectors) [33]. For hollow cylinders that are isotropic or have an orthotropy pattern of type A: (i) the modes associated with eigenvalues that have algebraic multiplicity equal to 1 for moderate values of l are uniform in the circumferential direction (graphs (a) of Figs. 7 and 8) and (ii) the modes associated with eigenvalues of algebraic multiplicity 2 and geometric multiplicity 1 at l ¼ 0 are non-uniform in the circumferential direction (graphs (b) of Figs. 7 and 8). The previous two features are shared with flutter modes of infinite layers [1] (substitute ‘‘circumferential direction” by ‘‘the direction

parallel to the obstacle”). Most eigenvalues have algebraic multiplicity 2 and geometric multiplicity 1. Naturally, the flutter modes for orthotropic hollow cylinders of type B are non-uniform in the circumferential direction (Fig. 9); unlike isotropic and orthotropic cylinders of type A, most eigenvalues have algebraic multiplicity 1 and eigenvalues with algebraic multiplicity 2 and geometric multiplicity 1 correspond always to high frequencies. Note that the magnitude of the obstacle’s velocity in (11) does not affect the (exponential) rate of growth of unstable modes; it affects rather the time interval of validity of the unstable solution of (12)–(15). Anyway, once the Liapunov instability of the system is demonstrated, it will abandon the neighbourhood of the steady sliding equilibrium state in finite time and will tend to a limit cycle (if not to chaotic motion): Moirot et al. [23] give examples of waves involving stick, slip (and in some cases separation) that represent the limit cycle towards which hollow cylinders converge. Fig. 10 shows three of such limit cycles for the tangential degree of freedom of a contact candidate node of an isotropic hollow cylinder; they correspond to the obstacle angular velocities xobst ¼ 0:1; 0:5 and 1.0 rad s1. In all the presented limit cycles there is a section of stiction corresponding to the horizontal upper segment of constant velocity equal to 0:25xobst and a longer section of backward slip. It is observed that the horizontal (tangential displacements) and vertical (tangential velocities) dimensions of the limit cycles are exactly proportional to the obstacle’s angular velocity. In linear elasticity, the consequences of divergence instability cannot be as visible as those of flutter instability since flutter is observed when the rigid inner cylinder is rotating (thus permanently injecting a certain amount of power into the system) while the consequences of divergence instabilities depend very much on the amount of elastic potential energy released during the

Fig. 7. Vibration/instability modes of an isotropic hollow cylinder with E ¼ 5;

0

m ¼ 0:3; q ¼ 1:2  103 g mm3 and rE rr ¼ 0:005.

1282

A. Pinto da Costa, M.A. Agwa / Computers and Structures 87 (2009) 1275–1286

Fig. 8. Vibration/instability modes of an orthotropic hollow cylinder with each family of principal directions of orthotropy makes the same angle h ¼ 30 with the 0 circumferential direction (type A). E1 ¼ 4; E2 ¼ 1; m12 ¼ 0:1; m ¼ 0:3; q ¼ 1:2  103 g mm3 and GE122 ¼ 0:25; Er2rr ¼ 0:005.

Fig. 9. First three vibration modes associated with simple ks (for 0 m ¼ 0:3; GE122 ¼ 0:1; q ¼ 1:2  103 g mm3 and Er2rr ¼ 0:005.

l ¼ 0) of an orthotropic hollow cylinder of type B with u ¼ 45 ; E1 ¼ 16; E2 ¼ 1; m12 ¼ 0:4;

instability process (in this case there is no power injected into the system since the inner rigid cylinder is considered to be at rest). Only in the finite elasticity framework and for strongly deformed bodies, which is outside the scope of this study, the consequences of divergence would turn out to be spectacular, because in that case the amount of elastic energy accumulated before the instability process may be capable of deviating considerably the solid from its initial state. 4.4. Locii of the eigenvalues A large number of numerical experiments on the way that the spectrum associated with the pair ðM ðlÞ; K ðlÞÞ in k2 M þ K depends on l showed some patterns that are qualitatively represented in Fig. 11. This figure illustrates the locii of the eigenvalues k in terms of the coefficient of friction. All the eigenvalues corresponding to l ¼ 0 are pure imaginary numbers representing

the natural frequencies of the frictionless unilaterally constrained hollow cylinder. Two kinds of eigenvalues may be identified: simple eigenvalues, denoted by S, with an algebraic multiplicity equal to 1 and double eigenvalues, denoted by D, with an algebraic multiplicity equal to 2. For isotropic hollow cylinders (Fig. 11a), some single eigenvalues depend on l and others do not; on the other hand, for orthotropic hollow cylinders of types A or B all single eigenvalues vary with l (Fig. 11b and c). Divergence instability is in general possible due to the coalescence of two simple pure imaginary eigenvalues at k ¼ 0 for a certain critical value of l, then bifurcating into two mutually symmetric real eigenvalues. Unlike orthotropic hollow cylinders of type B, isotropic and orthotropic hollow cylinders of type A exhibit flutter instability for arbitrarily small values of l due to the bifurcation of double eigenvalues with an algebraic multiplicity equal to 2 into two complex conjugate eigenvalues. Moreover, some simple eigenvalues of

A. Pinto da Costa, M.A. Agwa / Computers and Structures 87 (2009) 1275–1286 0.4

1283

5. Simple models

0.3

Two lumped parameter models will now be worked out with the purpose of illustrating that some aspects of the continuum may be captured by examples with a small number of degrees of freedom.

vt (m s−1)

0.2 0.1 0 −0.1

5.1. Klarbring’s example with a curved obstacle

−0.2 −0.3 −0.4 −0.5

0

0.5

1

ut (m)

1.5

2

2.5

−4

x 10

Fig. 10. Limit cycles (ut nodal coordinate) of a contact candidate node of an isotropic hollow cylinder with E ¼ 5; m ¼ 0:3; q ¼ 1:2  103 g mm3 and r0rr ¼ 0:005, computed with ABAQUS software. The inner (intermediate, outer) E cycle corresponds to an angular velocity of the obstacle of xobst ¼ 0:1 ð0:5; 1:0Þ rad s1 . Coefficient of friction: l ¼ 0:1.

the orthotropic types A and B coalesce in pairs for certain finite values of the coefficient of friction yielding flutter; in isotropic hollow cylinders, flutter instability modes cannot appear at finite values of l due to coalescence of simple eigenvalues, it can only appear at vanishing values of the coefficient of friction. For the orthotropy pattern B, the locii of the double eigenvalues at l ¼ 0 corresponding to high frequency modes do not bifurcate until a certain finite value of the coefficient of friction is reached, a phenomenon that also has no correspondence in orthotropic infinite layers (see Section 4.3 of [1]).

The two degree of freedom system introduced by Klarbring [14] and studied in its original or modified forms in several publications [3,5,19,27] may be used to illustrate the effect of the obstacle curvature on the minimum value of l required for the occurrence of divergence instability. Fig. 12a represents a particle of mass m restrained by two linear elastic springs, one vertical of stiffness K v and another making an angle h with the vertical with stiffness K i . The particle is in impending slip to the right at the top of a cylindrical surface of radius R. By choosing the vector of generalized coordinates u ¼ fu1 u2 gT , the mass and elastic stiffness coefficients are M11 ¼ M22 ¼ m; M12 ¼ M21 ¼ 0 and K 11 ¼ K v þ K i c2 ; K 12 ¼ K 21 ¼ K i sc; K 22 ¼ K i s2 , where s ¼ sinðhÞ and c ¼ cosðhÞ. At the onset of divergence instability, problem CEIP-l2 is formed by the following conditions

½sðs  lcÞ  ð1 þ l2 ÞXx ¼ y;

ð36Þ

0 6 x ? y P 0;

ð37Þ

vr0n

where X ¼  K is the non-dimensional obstacle curvature i (v ¼ 1R > 0 for a convex obstacle). This is not a mixed problem since there are no nodes free from kinematic constraints. For the previous

Fig. 11. Qualitative representation of the eigenvalues locii as a function of the coefficient of friction: (a) isotropic, (b) orthotropy pattern A, (c) orthotropy pattern B.

1284

A. Pinto da Costa, M.A. Agwa / Computers and Structures 87 (2009) 1275–1286

μ

4

4

3

3

μ

2

1

2

1

60 45

−30

30 0 0

0.25

0.5

0.75

0 0

1

−45

0.25

0.5

−60 0.75

1

X

X

Fig. 12. (a) A particle constrained by a system of springs and a curved frictional contact. (b, c) Coefficient of friction at the onset of instability (the angle h is indicated near each curve).

problem to have a non-trivial ðx > 0; y ¼ 0Þ it is required that



solution

involving

sðs  lcÞ : 1 þ l2

sliding

ð38Þ

According to the previous equation, the equilibrium position of the mass on a frictionless obstacle is unstable for a non-dimensional curvature X ¼ s2 . From Fig. 12b and c it may be concluded that,

when the inclined spring is on the same side to which the particle is in a state of impending slip (graph on the right, h < 0), for mod the onset of divergence instability occurs erate values of l l < 1þs c for coefficients of friction that increase with the increase of the nondimensional obstacle curvature. Note that this phenomenon occurs for orientations of the inclined spring for which the absolute value of the normal reaction increases during sliding. The curves in Fig. 12b and c are not similar to those of the continuum case presented in Fig. 3. In some numerical experiments

4

3

μ

2

1

−60 0 0

1

−45 2

−30 3

4

X

6

6 5

0.5

5 1 4

4

μ

μ

3

16 3

2 2

2 4

4

8

1

1 16

0 0

30

60

8 2

90

θ

120

150

180

0 0

1 0.5 30

60

90

120

150

180

θ

Fig. 13. (a) Concentrated masses restrained by springs; the masses are in frictional contact with a cylinder. (b) Values of l at the onset of divergence instability for N ¼ 16; k1 ¼ 4 and k2 ¼ 1. Coefficient of friction at the onset of instability for N ¼ 144 concentrated masses restrained by springs, in frictional contact with a cylinder, for two v r 0 2pr0 values of X ¼ K 3 n ¼ NK 3 rr : (c) 0.05555, (d) 0.55555 (the values of the stiffness ratio kk12 are indicated near the curves).

A. Pinto da Costa, M.A. Agwa / Computers and Structures 87 (2009) 1275–1286

not reported here, curves similar to those in Fig. 12b and c were obtained for the continuous hollow cylinder, but for unacceptably r0 high values of E2rr . We are convinced that dependences of l on the curvature v, similar to those presented in Fig. 12b and c, are obtainable in large deformation analyses of hollow cylinders submitted to high internal pressures r0rr . 5.2. Lumped masses around a cylinder With the purpose of developing a discrete model that is more similar to the continuous hollow cylinder we have developed the axisymmetric model presented next. Fig. 13a illustrates a lumped mass system consisting on an even number N of masses in frictional unilateral contact with a cylindrical surface of radius R. Each mass p is restrained by two mutually orthogonal springs of stiffnesses K p1 and K p2 inclined by an angle hp with respect to the directions tangential and normal to the cylindrical surface. Each two consecutive masses are linked by a spring of stiffness K 3 . This is a 2N degree of freedom system: ðu2p1 ; u2p Þ is the displacement vector of mass p (respectively normal and tangential to the contact surface). All masses are assumed to be in a state of impending slip in the trigonometric direct sense. In the special case where all the stiffnesses and angles hp are equal ðhÞ for all particles, the matrix operators of problems CEIP-k2 and CEIP-l2 become circulant [16] due to the cyclic nature of the structure. In a circulant matrix each row is obtained from the previous row by moving it one place to the right. The non-dimensional version of problem CEIP-k2 is: Find k P 0; x – 0 and y such that

 a circ k2 þ k1 c2 þ k2 s2 þ 2 cos2 þ lðk2  k1 Þsc  ð1 þ l2 ÞX 2 1 a a s s ð39Þ þ l ; 0; 0; . . . ; 0;  cos2  l Ax ¼ y;  cos2 2 |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} 2 2 2 N3 zeros

0 6 x ? y P 0;

ð40Þ k1 ¼ KK 13

k2 ¼ KK 23

vr0n

where s ¼ sinðhÞ; c ¼ cosðhÞ; a ¼ 2Np ; ; and X ¼  K 3 > 0 ðv ¼ þ 1R > 0; r 0n < 0 is the static normal reaction from the obstacle on each particle at equilibrium). For a sufficiently high number N of particles the relationship between the static equilibrium normal stress of the continuous struc ture ðr0rr Þ and the normal reaction of the lumped model r0n is 0 0 N rrr ¼ vrn 2p. As observed in the modified example of Klarbring, in this example similar curves of l versus the non-dimensional curvature X were obtained containing a turning point below which the transition between stability and instability is a line corresponding to decreasing l for decreasing non-dimensional curvature (see Fig. 13b); this happens for a sufficiently high degree of orthotropy and when the orientation of the stiffer principal direction of orthotropy makes an acute angle with the direction towards which the masses are in a state of impending slip. Note that in this lumped mass model the prescription of a normal stress r0rr (as we did for obtaining Fig. 3) corresponds to prescribing a single value for the non-dimensional curvature X in (39) and consequently to obtain a single value for l, not a curve of l versus v as obtained for the continuum case. For a set of N ¼ 144 masses and k2 ¼ 1, Fig. 13c and d plots the coefficient of friction at the onset of instability versus the inclination angle h that simulates the orientation of the principal directions of orthotropy; two values of the normal stress r0rr and several values of the stiffness ratio kk12 were considered. According to the figure, an increase of the absolute value of r0rr implies a decrease of l. There are two orientations of the principal directions of orthotropy for which l at the onset of divergence instability is

1285

insensitive to the stiffness ratio kk12 . The case kk12 ¼ 1 corresponds to isotropy: it yields a value of l independent of h.

6. Conclusions An analytic and finite element study showing how obstacle curvature, normal stress, elastic coefficients and orientation of the principal directions of orthotropy affect the onset of instability of hollow cylinders and discrete models was presented. Special attention was paid to hollow cylinders with two patterns of principal directions of orthotropy: type A (constant orientation with respect to radial and circumferential directions) and type B (constant orientation with respect to a fixed cartesian reference system). In this paper we have presented a complementarity eigenvalue problem, whose coefficient matrix is a pencil of matrices, quadratic in the coefficient of friction; this special eigenproblem is able to consider all the potential tangent behaviours of a finite dimensional system in frictional contact with curved obstacles. For certain orientations of the principal directions of orthotropy, values of the obstacle curvature and some common relations between the elastic constants, the value of the coefficient of friction at the onset of divergence instability may easily get below 0.5. Modes with a non-uniform velocity profile at contact, sometimes involving stiction of nodes that were in impending slip at equilibrium, were obtained for non-homogeneous hollow cylinders. It was also concluded that, for some combinations of material properties, a non-intuitive phenomenon may occur: the onset of divergence instability for coefficients of friction that increase for increasing obstacle convexity. This study further showed that (1) hollow cylinders of type A are unstable by flutter for arbitrarily small values of the coefficient of friction and (2) hollow cylinders of type B are not unstable by flutter for arbitrarily small values of the coefficient of friction but may be unstable by flutter at a finite value of l by bifurcation of double pure imaginary eigenvalues into complex conjugate eigenvalues. Two lumped mass models provide additional physical insight into the role of the obstacle curvature, namely that the increase of the curvature of a convex obstacle may, in some circumstances, delay the onset of divergence instability. Both models may exhibit divergence instabilities at vanishing coefficients of friction: in this case instability occurs due to the combined effect of curvature and normal stress, not friction. A similar phenomenon is expected to occur in finite elasticity analyses of solids strongly compressed against frictionless convex obstacles.

Acknowledgements The beginning of this work has been developed with the support of a scholarship from project number 1117, ‘‘Programa de Financiamento Plurianual Programático” (through ICIST) and its final part was possible due to scholarship number SFRH/BD/31173/ 2006, both from Fundação Portuguesa para a Ciência e a Tecnologia (FCT), awarded to M.A. Agwa. The authors wish to express their thanks to the editor and to the referees for the pertinent comments that led to the improvement of this work.

Appendix A. Constitutive matrix in the coordinate system ðr; u; x3 Þ The inverse Dp of the compliance matrix has the following components in the local coordinate system parallel to the principal directions of orthotropy

1286

A. Pinto da Costa, M.A. Agwa / Computers and Structures 87 (2009) 1275–1286

2

Dp11 6 Dp 6 12 6 p 6 D12 : Dp ¼6 6 0 6 6 4 0 0

3

Dp12

Dp12

0

0

0

Dp22 Dp23

Dp23 Dp22

0 0

0 0

0

0

G12

0

0

0

0

1 E2 2 1þm

7 7 7 7 7; 0 7 7 7 0 5

0

0

0

0

0 0

and D ¼ 1  m 

2m212 E2 . E1

[13]

G12

Dp11 ¼ ð1DmÞE1 ; Dp12 ¼ m12DE2 ; Dp22 ¼

where

[12]

ð41Þ

[14] E2 ðE1 m212 E2 Þ ; E1 ð1þmÞD

Dp23 ¼

E2 ðmE1 þm212 E2 Þ E1 ð1þmÞD

The transformation matrix T, taken from

Ref. [6], involved in the transformation of Dp into the polar coordinate system ðr; u; x3 Þ, has the form

2

c2

6 s2 6 6 0 :6 T¼6 6 2sc 6 6 4 0 0

s2

0

sc

c2

0

sc

0

1

0

0

0

3

0

0

0

0 7 7 7 0 0 7 7; 0 0 7 7 7 c s 5

0

0

0

s

2sc

0 c 2  s2

[15]

[16] [17] [18]

[19]

0

ð42Þ

c

with s ¼ sinðhÞ and c ¼ cosðhÞ. The constitutive matrix in the coordinate system ðr; u; x3 Þ is D ¼ TT Dp T which is not shown explicitly because its elements are trigonometric polynomials up to power 4. References [1] Agwa MA, Pinto da Costa A. Instability of frictional contact states in infinite layers. Eur J Mech, A/Solids 2008;27:487–503. [2] Brooke A, Kendrick D, Meeraus A, Raman R, Rosenthal RE. GAMS: a user’s guide. Optimization methods and software. GAMS Development Corporation, 1998. . [3] Basseville S, Leger A, Pratt E. Investigation of the equilibrium states and of their stability for a simple model with unilateral contact and Coulomb friction. Arch Appl Mech 2003;73(5–6):409–20. [4] Björkman G. The solution of large displacement frictionless contact problems using a sequence of Linear Complementarity Problems. Int J Numer Methods Eng 1991;31:1553–66. [5] Cho H, Barber JR. Dynamic behavior and stability of simple frictional systems. Math Comput Model 1998;28(4–8):37–53. [6] Cook RD, Malkus DS, Plesha ME. Concepts and applications of finite element analysis. John Wiley & Sons; 1989. [7] Dirkse SP, Ferris MC. The PATH solver: a nonmonotone stabilization scheme for mixed complementarity problems. Optimiz Methods Software 1995;5:123–56. [8] Hassani R, Ionescu IR, Sakki ND. Unstable perturbation of the equilibrium under Coulomb friction: nonlinear eigenvalue analysis. Comput Methods Appl Mech Eng 2007;196:2377–89. [9] Hjiaj M, Feng Z-Q, de Saxcé G, Mróz Z. Three dimensional finite element computations for frictional contact problems with non-associated sliding rule. Int J Numer Methods Eng 2004;60:2045–76. [10] Huseyin K. Vibrations and stability of multiparameter systems. Sijthoff & Noordhoff International Publishers; 1978. [11] Ionescu IR, Hassani R. Equilibrium configurations with Coulomb friction. Existence, multiplicity and stability. Solids, structures and coupled problems

[20]

[21] [22] [23] [24]

[25]

[26]

[27]

[28]

[29]

[30]

[31] [32] [33] [34]

in engineering. In: Mota Soares CA, et al. Proceedings of the III European conference on computational mechanics, Laboratrio Nacional de Engenharia Civil, Lisboa, Portugal, 5–8 June, 2006. Júdice JJ, Ribeiro IM, Sherali HD. The eigenvalue complementarity problem. Comput Optimiz Appl 2007;37:139–56. Klarbring A. On discrete and discretized non-linear elastic structures in unilateral contact (stability, uniqueness and variational principles). Int J Solids Struct 1988;24(5):459–79. Klarbring A. Examples of non-uniqueness and non-existence of solutions to quasistatic contact problems with friction. Ing Arch 1990;56:529–41. Klarbring A, Björkman G. Solution of large displacement contact problems with friction using Newton’s method for generalized equations. Int J Numer Methods Eng 1992;34:249–69. Lütkepohl H. Handbook of matrices. John Wiley & Sons; 1996. Maple Waterloo Maple Inc., 2001. . Martins JAC, Barbarin S, Raous M, Pinto da Costa A. Dynamic stability of finite dimensional linearly elastic systems with unilateral contact and Coulomb friction. Comput Methods Appl Mech Eng 1999;177:289–328. Martins JAC, Pinto da Costa A. Stability of finite-dimensional nonlinear elastic systems with unilateral contact and friction. Int J Solids Struct 2000;37(18):2519–64. Miersemann E, Mittelmann HD. On the continuation for variational inequalities depending on an eigenvalue parameter. Math Methods Appl Sci 1989;11:95–104. Mittelmann HD. An efficient algorithm for bifurcation problems of variational inequalities. Math Comput 1983;41(164):473–85. Moirot F. Étude de la Stabilité d’un Équilibre en Présence de Frottement de Coulomb. PhD thesis. École Polytechnique, Paris, France, 1998. Moirot F, Nguyen Q–S, Oueslati A. An example of stick–slip and stick–slipseparation waves. Eur J Mech, A/Solids 2002;22:107–18. Naumann J, Wenk HU. On eigenvalue problems for variational inequalities. An application to nonlinear plate buckling. Rend di Matematica, Univ Roma 1976;6:439–63. Pinto da Costa A, Figueiredo IN, Júdice JJ, Martins JAC. Contact avec frottement en terme de complementarité: instabilité et non-unicité. Colloque Instabilité du Frottement, Laboratoire de Mathématiques de l’Université de Savoie, Le Bourget-du-Lac, France, 27–29 September; 1999. Pinto da Costa A, Figueiredo IN, Júdice JJ, Martins JAC. A complementarity eigenproblem in the stability analysis of finite dimensional elastic systems with frictional contact. In: Ferris MC, Mangasarian OL, Pang JS, editors. Complementarity: applications, algorithms and extensions, vol. 50. Dordrecht, Boston, London: Kluwer Academic Publishers; 2001. p. 67–83. Pinto da Costa A, Martins JAC. The evolution and rate problems and the computation of all possible evolutions in quasi-static frictional contact. Comput Methods Appl Mech Eng 2003;192:2791–821. Pinto da Costa A, Martins JAC, Figueiredo IN, Júdice JJ. The directional instability problem in systems with frictional contacts. Comput Methods Appl Mech Eng 2004;193:357–84. Pinto da Costa A, Martins JAC. A numerical study on multiple rate solutions and onset of directional instability in quasi-static frictional contact problems. Comput Struct 2004;82:1485–94. Pinto da Costa A, Seeger A. Cone-constrained eigenvalue problems: theory and algorithms. Comput Optimiz Appl, in press, doi:10.1007/s10589-0089167-8. Queiroz M, Júdice JJ, Humes Jr C. The symmetric eigenvalue complementarity problem. Math Comput 2004;73:1849–63. Quittner P. Spectral analysis of variational inequalities. Comm Math Univ Carolinae 1986;27:605–29. Strang G. Linear algebra and its applications. Hartcourt Brace Jovanovich College Publishers; 1988. Ziegler H. Principles of structural stability. Blaisdell Publishing Company; 1968.