Non-linear vibrations of doubly curved shallow shells

International Journal of Non-Linear Mechanics 40 (2005) 683 – 710 www.elsevier.com/locate/nlm

Non-linear vibrations of doubly curved shallow shells M. Amabili∗ Dipartimento di Ingegneria Industriale, Università di Parma, Parco Area delle Scienze 181/A, Parma, 43100, Italy Received 1 March 2004; received in revised form 4 August 2004; accepted 27 August 2004

Abstract Large amplitude (geometrically non-linear) vibrations of doubly curved shallow shells with rectangular base, simply supported at the four edges and subjected to harmonic excitation normal to the surface in the spectral neighbourhood of the fundamental mode are investigated. Two different non-linear strain–displacement relationships, from the Donnell’s and Novozhilov’s shell theories, are used to calculate the elastic strain energy. In-plane inertia and geometric imperfections are taken into account. The solution is obtained by Lagrangian approach. The non-linear equations of motion are studied by using (i) a code based on arclength continuation method that allows bifurcation analysis and (ii) direct time integration. Numerical results are compared to those available in the literature and convergence of the solution is shown. Interaction of modes having integer ratio among their natural frequencies, giving rise to internal resonances, is discussed. Shell stability under static and dynamic load is also investigated by using continuation method, bifurcation diagram from direct time integration and calculation of the Lyapunov exponents and Lyapunov dimension. Interesting phenomena such as (i) snap-through instability, (ii) subharmonic response, (iii) period doubling bifurcations and (iv) chaotic behaviour have been observed. 䉷 2004 Elsevier Ltd. All rights reserved. Keywords: Vibrations; Non-linear; Shells; Curved; Chaos

1. Introduction Doubly curved panels are largely used in aeronautics and aerospace and are subjected to dynamic loads that can cause vibration amplitude of the order of the shell thickness, giving rise to significant non-linear phenomena.

∗ Tel.: +39 0521 905896; fax: +39 0521 905705.

E-mail address: [email protected] (M. Amabili) URL: http://me.unipr.it/mam/amabili/amabili.html. 0020-7462/$ - see front matter 䉷 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijnonlinmec.2004.08.007

An exhaustive literature review of work on the non-linear vibrations of curved panels and shells is given by Amabili and Païdoussis [1], mainly focusing on circular cylindrical shells. Pioneers in the study of large-amplitude vibrations of simply supported, circular cylindrical shallow shells were Reissner [2] and Grigolyuk [3]. Leissa and Kadi [4] studied linear and non-linear free vibrations of doubly curved shallow shells of rectangular boundaries, simply supported at the four edges without in-plane restraints. Donnell’s non-linear shallow-shell theory was used (in case of doubly curved shells, it is a slightly modified version

684

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

of the original Donnell’s theory, originally developed for circular cylindrical shells, to take into account the double curvature). A single mode expansion of the transverse displacement was used. Large amplitude vibrations of shallow shells such as elliptic paraboloids, parabolic cylinders and hyperbolic paraboloids, with zero displacements and rotational springs at the four boundaries were studied by Bhattacharya [5]. Sinharay and Banerjee [6] studied non-linear vibrations of shallow spherical and non-circular cylindrical shells. Furthermore, the same elliptic paraboloid, parabolic cylinder and hyperbolic paraboloid previously studied by Bhattacharya [5] were analysed and results were compared. In the case of cylindrical panels, the edges were immovable and elastically restrained against rotation. The authors used a simple mode expansion and an energy approach. Sinharay and Banerjee [7] also studied shallow spherical and circular cylindrical shells of nonuniform thickness with a single-mode approach. Chia [8] investigated doubly curved panels with rectangular base by using Donnell’s non-linear shallow-shell theory and a single-mode expansion in all the numerical calculations, for both vibration shape and initial imperfection. Only the backbone curves, relative to free vibrations, were obtained. Hui [9] studied free vibrations of rectangular plates and axi-symmetric shallow spherical shells with geometric imperfections. Backbone curves, indicating an initial softening type non-linearity, turning to hardening for larger amplitudes were obtained. A singlemode approach was used. Response to initial condition in presence of viscous damping was also investigated. Librescu and Chang [10] investigated geometrically imperfect, doubly curved, undamped, laminated composite panels. The non-linear theory of shear-deformable shallow panels was used. The nonlinearity was due to finite deformations of the panel due to in-plane loads and imperfections. Only smallamplitude free vibrations superimposed to this finite initial deformation were studied. A single-mode was used to describe the free vibrations and the initial imperfections. Imperfections of the panels, of the same shape as the mode investigated, lowered the vibration frequency significantly. Librescu and Chang [10] also described accurately the post-buckling stability. In fact, curved panels are characterized by an unstable

post-buckling behaviour, in the sense that they are subject to a snap-through-type instability. Kobayashi and Leissa [11] studied free vibrations of doubly curved thick shallow shells; they used the nonlinear first-order shear deformation theory of shells in order to study thick shells. The rectangular boundaries of the panel were assumed to be simply supported at the four edges. A single mode expansion was used for each of the three displacements and two rotations involved in the theory; in-plane and rotational inertia were neglected. The problem was then reduced to one of a single degree of freedom describing the out-of-plane displacement. Numerical results were obtained for circular cylindrical, spherical and paraboloidal shallow shells. Except for hyperbolic paraboloid shells, a softening behaviour was found becoming hardening for vibration amplitudes of the order of the shell thickness. However, increasing the radius of curvature, i.e. approaching a flat plate, the behaviour changed and became hardening. The effect of the shell thickness was also investigated. Free vibrations of doubly curved, laminated, clamped shallow shells with rectangular boundary were studied by Abe et al. [12]. Both first-order shear deformation theory and classical shell theory (analogous to Donnell’s theory) were used. Results obtained neglecting in-plane and rotary inertia are very close to those obtained retaining these effects. Only two modes were considered to interact in the non-linear analysis for higher modes, while a single mode was used for the mode with one circumferential and one longitudinal half-wave. Soliman and Gonçalves [13] studied large amplitude, forced vibrations and stability of axi-symmetric shallow spherical shells (spherical cups, i.e. with circular boundary). Transition to chaos for large harmonic loads, approaching the static stability limit, was investigated. In the present study, large-amplitude (geometrically non-linear) vibrations of doubly curved shallow shells with rectangular base, simply supported at the four edges and subjected to harmonic excitation normal to the surface in the spectral neighbourhood of the fundamental mode are investigated. Two different non-linear strain–displacement relationships, from the Donnell’s and Novozhilov’s shell theories, are used to calculate the elastic strain energy. In-plane inertia and geometric imperfections are taken into account. The

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

solution is obtained by Lagrangian approach. The present theory is based on the studies developed by Amabili [14–16] for circular cylindrical shells, circular cylindrical panels and rectangular plates, respectively, properly adapted to take into account the different geometry. The non-linear equations of motion are studied by using (i) a code based on arclength continuation method that allows bifurcation analysis and (ii) direct time integration. Numerical results are compared to those available in the literature and convergence of the solution is shown. Interaction of modes having integer ratio among their natural frequencies, giving rise to internal resonances, is discussed. Shell stability under static and dynamic load is also investigated by using continuation method, bifurcation diagram from direct time integration and calculation of the Lyapunov exponents and Lyapunov dimension. Interesting phenomena such as (i) snap-through instability, (ii) subharmonic response, (iii) period doubling bifurcations and (iv) chaotic behaviour have been observed; in this case up to four positive Lyapunov coefficients have been found, indicating hyperchaos.

2. Elastic strain energy of the shell A doubly curved shallow (small rise compared to the smallest radius of curvature) shell with rectangular base is considered, as shown in Fig. 1. A curvilinear coordinate system (O; x, y, z), having the origin O at one edge of the panel is assumed; x =
Rx and y = Ry where
and  are the angular coordinates and Rx and Ry are principal radii of curvature (constant); a and b are the curvilinear lengths of the edges and h is the shell thickness. The smallest radius of curvature at every point of the shell is larger than the greatest lengths measured along the middle surface of the shell. The displacements of an arbitrary point of coordinates (x, y) on the middle surface of the shell are denoted by u, v and w, in the x, y and z directions, respectively; w is taken positive outwards the centre of the smallest radius of curvature. Initial imperfections of the shell associated with zero initial tension are denoted by outof-plane displacement w0 , also positive outwards; only out-of-plane initial imperfections are considered. Two different strain–displacement relationships for thin shells are used in the present study, in order to compare results. They are based on Love’s first ap-

685

Fig. 1. Geometry of the doubly curved shallow shell and coordinate system.

proximation assumptions. These theories are: (i) Donnell’s and (ii) Novozhilov’s non-linear shell theories. According to these two theories, the strain components εx , εy and xy at an arbitrary point of the panel are related to the middle surface strains εx,0 , εy,0 and xy,0 , and to the changes in the curvature and torsion of the middle surface kx , ky and kxy by the following three relationships: εx = εx,0 + zk x , εy = εy,0 + zk y , xy = xy,0 + zk xy ,

(1)

where z is the distance of the arbitrary point of the panel from the middle surface. The middle surface strain–displacement relationships and changes in the curvature and torsion have different expressions for the Donnell’s and Novozhilov’s theories. According to Donnell’s non-linear shell theory, the middle surface strain–displacement relationships and changes in the curvature and torsion for a doubly curved shallow shell are εx,0 =


 ju w 1 jw 2 + + Rx j
Rx 2 Rx j
jw j w 0 + , Rx j
Rx j


(2a)

686

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

εy,0 =


 jv w 1 jw 2 + + Ry j Ry 2 Ry j jw j w 0 + , Ry j Ry j

ju jv jw jw + + Ry j Rx j
Rx j
Ry j jw j w 0 jw0 jw + + , Rx j
Ry j Rx j
Ry j

xy,0 =

kx = −

ky = −

j2 w Rx2 j
2

j2 w Ry2 j2

kxy = −2

∼ = (2b)

(2c)

,

(2d)

,

(2e)

j2 w . Rx Ry j
j

(2f)

According to Novozhilov’s non-linear theory, displacements u, ˜ v˜ and w˜ of points at distance z from the middle surface are introduced; they are related to displacements on the middle surface by the relationships [17] u˜ = u + z,

(3a)

v˜ = v + z,

(3b)

w˜ = w + w0 + z,

(3c)


 j(w + w0 ) jv u w 1+ = − − + Rx j
Rx Ry j Ry
 j(w + w0 ) ju v + − Ry j R Ry j  y
jw u w0 − − , (3d) Rx j
Rx Ry



 j(w + w0 ) ju v w 1+ − + Ry j Ry Rx j
Rx
 j(w + w0 ) jv u + − Rx j
R Rx j
 x
jw v w0 − , (3e) − Ry j Ry Rx


(3f)

The strain–displacement relationships for a generic point of the shell are [17] 
 2
1 1 ju˜ ju˜ εx = +w˜ + +w˜ Rx + z j
j
2(Rx +z)2
2
2  jv˜ jw˜ , (4a) + + − u˜ j
j

 
jv˜ ju˜ 2 1 1 εy = + w˜ + Ry + z j j 2(Ry + z)2
2
2  jv˜ jw˜ , (4b) + + w˜ + − v˜ j j

xy =

jv˜ ju˜ 1 1 1 1 + + Rx + z j
Ry + z j Rx + z Ry + z 

  ju˜ ju˜ jv˜ jv˜ × + w˜ + + w˜ j j
j
j

 jw˜ jw˜ + − u˜ − v˜ . (4c) j
j

The following thinness approximations are now introduced 1/(R + z) = (1/R){1 − z/R + O[(z/R)2 ]},

where

= −


 1 1 ju jv + + (w + w0 ) + Rx j
Ry j Rx Ry

 ju jv w + w0 w + + + R j
Rx Ry j Ry
x  jv ju ju w w0 + + − . Rx j
Rx Ry Ry j Rx j


(5a)

1/(R + z)2 = (1/R 2 ){1 − 2z/R + O[(z/R)2 ]}, (5b) where O((z/R)2 ) is a small quantity of order (z/R)2 . By introducing Eqs. (3) into Eqs. (4) and using approximations (5), the middle surface strain–displacement relationships, changes in the curvature and torsion are obtained for the Novozhilov theory of shells 
 2
2 ju ju jv w 1 εx,0 = + + + w + 2 Rx j
Rx 2Rx j
j


2   1 jw 0 jw jw + 2 + −u −u j
Rx j
j

 ju +w0 +w , (6a) j


M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710




jv w 1 εy,0 = + + Ry j Ry 2Ry2
2  1 jw + 2 + −v j Ry
 jv , +w0 w + j xy,0 =


2  ju 2 jv + +w j j 
 jw0 jw −v j j (6b)

ju jv + Ry j Rx j

 ju ju 1 + + w + w0 Rx Ry j j

 jv jv + + w + w0 j
j

 jw jw + −u −v j
j

  jw 0 jw jw 0 jw + −v + −u . j
j j j


case of plane stress) by [18] E (εx + εy ), 1 − 2 E =  , 2(1 + ) xy

x = xy

ky = −

j2 w Rx2 j
2

US =

j2 w Ry2 j2

+

jv ju w + + , Rx Ry Rx Ry j Rx2 j


(6d)

+

ju jv w + + , Rx Ry Rx Ry j
Ry2 j

(6e)

j2 w




1 jv 1 − + Rx Ry j
j Rx j
Rx Ry
 1 ju 1 + . − Ry j Rx Ry

kxy = − 2



+

(6f)

Eqs. (6a)–(6f) are an improved version of Novozhilov’s non-linear shell theory because approximations (5) have been used instead of neglecting terms in z and z2 with respect to R as done in [17]. Non-linearities in changes of curvature and torsion have been neglected in this case for simplicity. The elastic strain energy US of the shell, neglecting z as stated by Love’s first approximation assumptions, is given by    1 a b h/2 US = (x εx + y εy + xy xy ) 2 0 0 −h/2 × (1 + z/Rx )(1 + z/Ry ) dx dy dz, (7) where the stresses x , y and xy are related to the strains for homogeneous and isotropic material (z =0,

y =

E (εy + εx ), 1 − 2 (8)

where E is the Young’s modulus and is the Poisson’s ratio. By using Eqs. (1), (7) and (8), the following expression is obtained:

(6c) kx = −

687

1 Eh 2 1 − 2



a



b



2 2 + εy,0 εx,0  1− 2 +2 εx,0 εy,0 + xy,0 dx dy 2   a b Eh3 1 kx2 + ky2 + 2 kx ky + 2 12(1 − 2 ) 0 0  
1 1 1− 2 1 kxy dx dy + + + 2 2 Rx Ry  a b 3 Eh εx,0 kx +εy,0 ky + εx,0 ky × 6(1− 2 ) 0 0  1− + εy,0 kx + xy,0 kxy dx dy + O(h4 ), 2 (9) 0

0

where O(h4 ) is a higher-order term in h and the last term in h3 disappears if z/R is neglected with respect to unity in Eq. (7), as it must be done for Donnell’s theory. If this term is neglected, the right-hand side of Eq. (9) can be easily interpreted: the first term is the membrane (also referred to as stretching) energy and the second one is the bending energy. If the last term is retained, membrane and bending energies are coupled.

3. Mode expansion, kinetic energy and external loads The kinetic energy TS of the shell, by neglecting rotary inertia, is given by TS =

1

h 2 S



a 0



b

(u˙ 2 + v˙ 2 + w˙ 2 ) dx dy,

(10)

0

where S is the mass density of the shell. In Eq. (10) the overdot denotes a time derivative.

688

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

The virtual work W done by the external forces is written as  a b (qx u + qy v + qz w) dx dy, (11) W= 0

M

N

v(x, y, t) =

× cos(n y/b),

0

where qx , qy and qz are the distributed forces per unit area acting in x, y and normal directions, respectively. Initially, only a single harmonic normal force is considered; therefore qx = qy = 0. The external, normal, distributed load qz applied to the shell, due to the radial concentrated force f˜, is given by qz = f˜(x − x) ˜ (y − y) ˜ cos( t),

(12)

where is the excitation frequency, t is the time,  is the Dirac delta function, f˜ gives the normal force amplitude positive in w direction, x˜ and y˜ give the position of the point of application of the force; here, the point excitation is located at x˜ = a/2, y˜ = b/2. Eq. (11) can be rewritten in the following form: W = f˜ cos( t)(w)x=a/2,y=b/2 .

(13)

Eq. (13) specialized for the expression of w used in the present study is given in Section 4, together with the virtual work of loads in the x, y plane and uniform normal pressure. In order to reduce the system to finite dimensions, the middle surface displacements u, v and w are expanded by using approximate functions. The classical simply supported boundary conditions at the four edges of the shell are assumed in this study:

(14a–f)

u = w = w0 = Ny = My = j w0 /jy = 0, at y = 0, b,

(15a–f)

2

2

2

where N is the normal force and M is the bending moment per unit length. A base of panel displacements is used to discretize the system. The displacements u, v and w can be expanded by using the following expressions, which satisfy identically the geometric boundary conditions (14a,b, 15a,b): u(x, y, t) =

N M



um,n (t) cos(m x/a)

m=1 n=1

× sin(n y/b),

w(x, y, t) =

(16a)

N M



(16b)

wm,n (t) sin(m x/a)

m=1 n=1

× sin(n y/b),

(16c)

where m and n are the number of half-waves in x and y direction, respectively; um,n (t), vm,n (t) and wm,n (t) are the generalized coordinates that are unknown functions of t. By using a different number of terms in Eqs. (16), it is possible to study the convergence and accuracy of the solution. A sufficiently accurate model for the mode (m = 1, n = 1) in case of no internal resonances has been shown to have 9 degrees of freedom (dofs) (for the case studied in Section 6). In particular, the following terms in Eq. (16) have been used: m = 1, 3 and n = 1, 3 in Eqs. (16a, b) and m = 1 and n = 1 in Eq. (16c). 3.1. Geometric imperfections Initial geometric imperfections of the shell are considered only in the normal direction, z. They are associated with zero initial stress. The normal imperfection w0 is expanded in the same form of w, i.e. in a double Fourier sine series satisfying the boundary conditions (14c,f, 15c,f) at the shell edges w0 (x, y) =

v = w = w0 = Nx = Mx = j w0 /jx = 0, at x = 0, a, 2

vm,n (t) sin(m x/a)

m=1 n=1

N˜ M˜



Am,n sin(m x/a)

m=1 n=1

× sin(n y/b),

(17)

where Am,n are the modal amplitudes of imperfections; N˜ and M˜ are integers indicating the number of terms in the expansion. 3.2. Boundary conditions Eqs. (14, 15) give the boundary conditions for a simply supported shell. Eqs. (14a–c,f, 15a–c,f) are identically satisfied by the expansions of u, v, w and w0 . On the other hand, Eqs. (14d,e, 15d,e) can be rewritten in the following form [18]: Mx =

Eh3 (kx + ky ) = 0, 12(1 − 2 )

(18)

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

My =

Eh3 (ky + kx ) = 0, 12(1 − 2 )

(19) (20)

juˆ jw jw0 + + Ry j Ry j Rx j


×

k=1 s=1

N

M

k=1 s=1

1 wm,n (t) sin(m x/a) 2

k ws,k (t) sin(s x/a) n+k

N˜

M˜

j Ai,j sin(i x/a) n+j j =1 i=1 × sin[sin(n + j ) y/b] .

×

1 u(t) ˆ = − [am,n (t) + bm,n (t) cos(2n y/b)] 8 × sin(2m x/a) − (m /a) × wm,n (t) sin(n y/b)

M˜ N˜

j =1 i=1

 = 0,

(23)

=0,b/Ry

s ws,k (t) sin(k y/b) m+s

M˜ N˜



i Ai,j sin(j y/b) m+i j =1 i=1 × sin[(m + i) x/a] ,

(24a)

(25a)

A simplified version of Eqs. (24a, 25a) that can be used in case of small interaction among mode (m, n) and other modes is

i Ai,j m+i

× sin(j y/b) sin[(m + i) x/a],

× sin[(m + s) x/a] + wm,n (t) sin(n y/b) ×

(n /b)

× sin[(n + k) y/b] + wm,n (t) sin(m x/a)


=0,a/Rx

where uˆ and vˆ are terms that must be added to the expansion of u and v, given in Eqs. (16a, b), in order to satisfy exactly the in-plane boundary conditions Nx =0 and Ny = 0. As a consequence that uˆ and vˆ are second order terms in the panel displacement, they have not been inserted in the second order terms that involve u and v in Eqs. (22, 23). Non-trivial calculations give (see Appendix A)  M N

1 wm,n (t) sin(n y/b) u(t) ˆ = − (m /a) 2 n=1 m=1 M N



×

(21)

Eqs. (18, 19) are identically satisfied for the expressions of kx and ky given in Eqs. (2d, e) and (6d, e) for Donnell’s and Novozhilov’s non-linear shell theories, respectively. Eqs. (20, 21) are not identically satisfied for both the non-linear shell theories. According to Donnell’s theory, see Eqs. (2a, b), and eliminating null terms at the panel edges, Eqs. (20), (21) can be rewritten as  
juˆ 1 jw 2 + Rx j
2 Rx j
 jvˆ jw jw0 + + = 0, (22) Ry j Rx j
Rx j

jvˆ 1 jw 2 + Ry j 2 Ry j



n=1 m=1

Eh Nx = (εx,0 + εy,0 ) = 0, 1 − 2 Eh Ny = (εy,0 + εx,0 ) = 0. 1 − 2



v(t) ˆ = −

M N



689

(24b)

1 v(t) ˆ = − [cm,n (t) + dm,n (t) cos(2m x/a)] 8 × sin(2n y/b) − (n /b) × wm,n (t) sin(m x/a)

N˜

M˜

j =1 i=1

j Ai,j n+j

× sin(i x/a) sin[(n + j ) y/b],

(25b)

where 2 am,n (t) = (m /a)wm,n (t),

(26)

2 (t), bm,n (t) = −(m /a)wm,n

(27)

2 (t), cm,n (t) = (n /b)wm,n

(28)

2 (t). dm,n (t) = −(n /b)wm,n

(29)

The simplified expressions (24b) and (25b) are obtained by assuming that all the generalized coordinates, except the three ones associated to the resonant mode (m, n), which are um,n (t), vm,n (t) and wm,n (t), are negligible because they are an infinitesimal of higher order. For Novozhilov’s non-linear shell theory, the approximate equations (24b, 25b) are still correct, but

690

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

Eqs. (26–29) are substituted by (boundary conditions are only approximately satisfied in this case) am,n (t) =

m 2 a u2m,n 2 )+ (wm,n + vm,n a m Rx2 a  n 2 2 2wm,n um,n − − um,n , 8 m b Rx m 2 a u2m,n 2 )− (wm,n − vm,n a m Rx2  a n 2 2 2wm,n um,n + + um,n Rx 8 m b 1 m 2 + v , 4 a m,n

(30)

The virtual work done by the concentrated normal force f˜, given by Eq. (13) specialized for the expression of w given in Eq. (16c), is W = f˜ cos( t)(w)x=a/2,y=b/2 = f˜ cos( t)   Nˆ Mˆ

(−1)i+1 (−1)j +1 w2i−1,2j −1 (t) , × i=1 j =1

bm,n (t) = −

cm,n (t) =

2 n 2 b vm,n (wm,n + u2m,n ) + b n Ry2 b  m 2 2 2wm,n vm,n − − vm,n , Ry 8 n a 2 n 2 b vm,n (wm,n − u2m,n ) − b n Ry2  b m 2 2 2wm,n vm,n + + vm,n Ry 8 n a 1 n 2 u . + 4 b m,n

(36)

(31)

(32)

where Mˆ = Integer Part(M/2 + 1) and Nˆ = Integer Part(N/2 + 1). In the presence of in-plane loads and pressure acting on the panel, additional virtual work is done by the external forces. Let us consider a time-dependent load N (t) per unit length, applied at both the shell edges x=0, a; N (t) is positive in the x direction. In particular −N (t) is applied at x = 0 and N (t) is applied at x = a. The distributed force qx has the following expression: qx = N (t)[−(x) + (x − L)],

dm,n (t) = −

(37)

where  is the Dirac delta function. The virtual work done by the load is  (33)

4. Lagrange equations of motion

a

W=



0

b

N (t)[−(x) + (x − a)]u dx dy

0

Mˆ

Nˆ

4b N (t) u2i−1,2j −1 (t)/(2j − 1). = −

i=1 j =1

The nonconservative damping forces are assumed to be of viscous type and are taken into account by using Rayleigh’s dissipation function  a b 1 F= c (u˙ 2 + v˙ 2 + w˙ 2 ) dx dy, (34) 2 0 0 where c has a different value for each term of the mode expansion. Simple calculations give F=

1 ab 2 4

M N

n=1 m=1

2 2 cm,n (u˙ 2m,n + v˙m,n + w˙ m,n ). (35)

The damping coefficient cm,n is related to modal damping ratio, that can be evaluated from experiments, by m,n = cm,n /(2m,n m,n ), where m,n is the natural circular frequency of mode (m, n) and m,n is the modal mass of this mode, given by m,n = S h(ab/4).

(38) In case of uniform internal time-varying pressure pz (t), the normal distributed force qz is obviously qz = pz (t).

(39)

The virtual work done by pressure is  W= 0

a



b 0

pz (t)w dx dy

Nˆ Mˆ

w2i−1,2j −1 (t) 4ab . = 2 pz (t)

(2i − 1)(2j − 1)

(40)

i=1 j =1

Eqs. (38) and (40) show that only modes with odd number of half-waves in x and y direction are directly excited by uniform in-plane loads and pressure.

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

The following notation is introduced for brevity: q = {um,n , vm,n , wm,n } , and n = 1, . . . N. T

m = 1, . . . M (41)

The generic element of the time-dependent vector q is referred to as qj and is referred as generalized coordinate; the dimension of q is dofs, which is the number of degrees of freedom used in the mode expansion. The generalized forces Qj are obtained by differentiation of Rayleigh’s dissipation function and of the virtual work done by external forces

jF jW + = −(ab/4)cj q˙j jq˙j jqj  0 if qj = um,n , vm,n ; or wm,n   with m or n even, + (42) ˜ cos( t) if qj = wm,n f   with both m, n odd. The Lagrange equations of motion are
 d jT S jT S jUS − + = Qj , dt jq˙j jqj jq j j = 1, . . . dofs, (43)

Qj = −

where jTS /jqj = 0. These second-order equations have very long expressions containing quadratic and cubic non-linear terms. In particular, 
d jT S = S h(ab/4)q¨j , (44) dt jq˙j which shows that no inertial coupling among the Lagrange equations exists for the shell with simply supported edges with the mode expansion used. The very complicated term giving quadratic and cubic non-linearities can be written in the form dofs dofs



jUS = qk z˜ j,k + qi qk z˜ j,i,k jqj k=1

+

i,k=1

dofs

qi qk ql z˜ j,i,k,l ,

(45)

i,k,l=1

where coefficients z˜ have long expressions that include also geometric imperfections. 5. Numerical techniques The equations of motion have been obtained by using the Mathematica 4 computer software [19] in order

691

to perform analytical integrals of trigonometric functions. The generic Lagrange equation j is divided by the modal mass associated with q¨j and is then transformed into the following two first-order equations by using the dummy variables yj : q˙j = yj , y˙j = − 2j j yj −

dofs

i=1

−

dofs

dofs dofs

i=1 k=1 l=1

+ f cos( t)

zj,i qi −

dofs

dofs

zj,i,k qi qk

i=1 k=1

zj,i,k,l qi qk ql for j = 1, . . . dofs,

(46)

where  0 if qj = um,n , vm,n ; or   wm,n with m or n even, f= ˜

hab/4] if qj = wm,n f /[  S  with both m, n odd, and coefficients z are obtained by those in Eq. (45), z = z˜ /[ S hab/4]. In Eq. (46) each generalized coordinate qj (and therefore j and j ) has to be referred to mode (m, n), as stated in Eq. (41). For computational convenience a non-dimensionalization of variables is also performed: the time is divided by the period of the resonant mode and the vibration amplitudes are divided by the shell thickness h. The resulting 2 × dofs first-order non-linear differential equations are studied by using (i) the software AUTO 97 [20] for continuation and bifurcation analysis of non-linear ordinary differential equations, and (ii) direct integration of the equations of motion by using the DIVPAG routine of the Fortran library IMSL. Continuation methods allow following the solution path, with the advantage that unstable solutions can also be obtained; these are not ordinarily attainable by using direct numerical integration. The software AUTO 97 is capable of continuation of the solution, bifurcation analysis and branch switching by using arclength continuation and collocation methods. In particular, the shell response under harmonic excitation has been studied by using an analysis in two steps: (i) first the excitation frequency has been fixed far enough from resonance and the magnitude of the excitation has been used as bifurcation parameter; the solution has been started at zero force where the solution is the trivial undisturbed configuration of the shell and has been continued up to reach

692

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

the desired force magnitude; and (ii) when the desired magnitude of excitation has been reached, the solution has been continued by using the excitation frequency as bifurcation parameter. Direct integration of the equations of motion by using Gear’s BDF method (routine DIVPAG of the Fortran library IMSL) has also been performed to check the results and obtain the time behaviour. The Adams Gear algorithm has been used due to the relatively high dimension of the dynamical system. Indeed, when a high-dimensional phase space is analysed, the problem can present stiff characteristics, due to the presence of different time scales in the response. In simulations with adaptive step-size Runge–Kutta methods, spurious non-stationary and divergent motions can be obtained. Therefore, the Adams Gear method, designed for stiff equations, was used. The bifurcation diagram of the Poincaré maps was also used in case of non-stationary response. It has been constructed by using the time integration scheme and by varying the force amplitude.

In order to evaluate the maximum Lyapunov exponent, which is useful to characterize regular or chaotic motion of the system, it is necessary to assume a reference trajectory xr (t) in the phase plane (q, q˙ plane) and observe a neighbouring trajectory originated at infinitesimal initial perturbation x(t0 ) from the reference trajectory (see Argyris et al. [21]). The evolution of the perturbation during time, x(t), is governed by the following variational equations directly obtained from Eq. (46) d qj = yj , dt dofs

zj,i,k qn (k,n qi + i,n qk )

n=1 i=1 k=1

−

zj,i,k,l qn (i,n qk ql

(48)

1 = lim sup t→∞

1 ln |x(t)|. t

(49)

Then, by restoring at each integration time step k, the amplitude of x(t) to its original unitary measure by the following re-normalization

x¯ (t)k =

x(t)k , dk

(50)

where |x(t)|k = dk , the following formula for the maximum Lyapunov exponent, evaluated at step k, is obtained k 1

ln di . k
t

(51)

In the numerical calculation of the maximum Lyapunov exponent, the non-dimensional time previously introduced has been used. A perturbation vector f1 is determined with the maximum Lyapunov exponent 1 . 5.2. All Lyapunov exponents

n=1 i=1 k=1 l=1

+ k,n qi ql + l,n qk qi ) for j = 1, . . . dofs,

1 |x(t)| ln . t |x(t0 )|

i=1

i=1

dofs

dofs

dofs dofs



t→∞

1,k =

d yj = − 2j j yj − zj,i qi dt −

1 = lim sup

Assuming the initial perturbation of unitary amplitude, Eq. (48) is simplified into

5.1. Maximum Lyapunov exponent

dofs

dofs dofs



where k,n is the Kronecker delta. Assuming qj and yj as new variables, the simultaneous integration of the 4 × dofs first-order differential equations (Eqs. (46) are non-linear and are integrated by using DIVPAG IMSL routine and Eqs. (47) are linear, but with time-varying coefficients, and are integrated by using the adaptive step-size fourth–fifth order Runge–Kutta method) has been performed. The excitation period has been divided in 10 000 integration steps
t in order to have accurate evaluation of the time-varying coefficients in Eqs. (47) that are obtained at each step by integration of Eqs. (46). To find a reference trajectory, 6×106 steps are skipped in order to eliminate the transient and 1×106 steps are skipped to eliminate the transitory on the variational Eqs. (47). Then 1 × 106 steps are used for evaluation of the maximum Lyapunov exponent 1 for the reference trajectory xr (t), which is given by

(47)

The 2 × dofs numbers designating the spectrum of the Lyapunov exponents [21] can be ordered according

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

to their magnitude

Then,

1
2
· · ·
2dofs .

(52)

It is clear that, for any perturbation, the largest Lyapunov exponent 1 will emerge and the perturbation vector will orient as f1 . In order to calculate all the Lyapunov exponents, it is necessary to introduce the Lyapunov exponent of pth order, indicated as (p) , which is the mean exponential growth rate of the volume V (p) of a p-dimensional parallelepiped (p = 2, . . . , 2dofs) included in the 2dofs space

(p) = lim sup t→∞

1 V (p) (t) ln (p) . t V (t0 )

(53)

This p-dimensional parallelepiped is spanned by p orthogonal perturbation vectors f1 , . . . , fp . As a consequence that the volume V (p) is obtained as a product of the length of these orthogonal perturbation vectors, which are directly related to the Lyapunov exponents, the Lyapunov exponent of pth order is equal to the sum of the corresponding Lyapunov exponents, i.e.

(p) = 1 + · · · + p .

(54)

Therefore once (p) is obtained,

p

(55)

− 1 − · · · − p−1 .

Therefore in order to obtain all the Lyapunov exponents, the Lyapunov exponent of pth order must be computed up to p = 2dofs. The main computational problem is to keep all the perturbation vectors orthogonal to each other. The set of p orthogonal perturbations is indicated with {f1 , . . . , fp }. After choosing casually an initial set, the following operation is performed at each integration time step k in order to restore the amplitude of each vector to its original unitary measure and keeping all of them orthogonal to each other [21] (f¯1 )k = (f1 )k /|f1 |k , (56a)     i−1 i−1  



  ¯ ci,m fm ci,m fm  (fi )k = fi +  fi +   m=1

k

(p)

Vk

= |f1 |k × · · · × |fp |k .

(57)

By using this computational procedure, the Lyapunov exponent of pth order is given by (p) k

k 1

(p) = ln Vi . k
t

(58)

i=1

By using Eqs. (57) and (55), all the Lyapunov exponents can be calculated. The Fortran computer program developed to calculate all the Lyapunov exponents has been initially validated by computing the two Lyapunov exponents of Duffing’s equation for different parameters and results have been very satisfactorily compared to those reported in Table 6-1 at page 313 of Moon’s book [22]. In particular, 1 × 107 steps have been used to evaluate the Lyapunov exponents (10 times larger than for calculation of the maximum Lyapunov exponent, that has been evaluated for all the bifurcation diagram, with high computational cost). 5.3. Lyapunov dimension

2 = (2) − 1 , .. . = (p)

693

m=1

for i = 2, . . . , p, with (ci,m )k = (−fi · f¯m )k .

k

(56b)

The long-term behaviour of dissipative systems is characterized by attractors of the most varied types if the trajectories are not drawn towards infinity. After a transient state, in which some modes of motion finally vanish due to damping, the state of the system approaches an attractor where the number of independent variables, which determine the dimension of the phase space, is generally reduced considerably [21]. The fractal dimension is a measure of the strangeness of an attractor [22] and indicates the number of effective independent variables determining the long-term behaviour of a motion. There exist several measures of the fractal dimension, including the well-known Lyapunov dimension, which is defined as [21–23] dL = s +

s

r /|s+1 |,

(59)

r=1

where the Lyapunov exponents are ordered by their magnitude, and s is obtained by satisfying the

694

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

following conditions: s

r=1

r > 0 and

s+1

r < 0.

(60)

r=1

6. Numerical results for Harmonic excitation Numerical calculations have been initially performed for doubly curved shallow shells, simply supported at the four edges, having the following dimensions and material properties: curvilinear dimensions a = b = 0.1 m, radius of curvature Rx = 1 m, thickness h = 0.001 m, Young’s modulus E = 206 × 109 Pa, mass density = 7800 kg/m3 and Poisson ratio = 0.3. Shallow shells with the same dimension ratios (Rx /a = 10, h/a = 0.01, a/b = 1, = 0.3) were previously studied by Kobayashi and Leissa [11]. In all the numerical simulations a modal damping 1,1 = 0.004 and harmonic force excitation at the centre of the shell in z direction are assumed. If not explicitly specified, all calculations have been performed by using Donnell’s shell theory. 6.1. Case with Rx /Ry = 1, spherical shell A spherical shallow shell (Ry = 1 m) is considered. The frequency range around the fundamental frequency (mode (m=1, n=1) in this case, where m and n are the numbers of half-waves in x and y direction, respectively) is investigated. The fundamental frequency 1,1 is 952.31 Hz according to Donnell’s shell theory and 952.26 Hz according to Novozhilov’s shell theory, i.e. practically the same results for both theories. Other natural frequencies useful in the present study are (according with Donnell’s theory): 1,3 = 3,1 = 2575.9 Hz, 3,3 = 4472.3 Hz. The amplitude of the harmonic force is f˜ = 31.2 N. Fig. 2 shows the maximum (in the time period; this is positive, i.e. outwards) and the minimum (negative, i.e. inwards) of the shell response in z direction in the spectral neighbourhood of the fundamental mode (m= 1, n = 1) versus the excitation frequency. Calculations have been performed with the 9 dofs model. This model includes the following terms in Eqs. (16a–c): w1,1 , u1,1 , v1,1 , u3,1 , v3,1 , u1,3 , v1,3 , u3,3 and v3,3 . Results are compared to those obtained by Kobayashi and Leissa [11], where only the backbone

Fig. 2. Amplitude of the response of the shell (generalized coordinate w1,1 ) versus the excitation frequency; Rx /Ry = 1 (spherical shallow shell); fundamental mode (m = 1, n = 1), f˜ = 31.2 N and 1,1 = 0.004; model with 9 dofs; Donnell’s theory. ——, present stable results; - -, present unstable results; − · −, backbone curve from Kobayashi and Leissa [11]. (a) Maximum of the generalized coordinate w1,1 (in a vibration period); and (b) minimum of the generalized coordinate w1,1 .

curve is given. The present results with 9 dofs are quite close to those in Ref. [11] (considering that the backbone curve, indicating the maximum of the

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

Fig. 3. Amplitude of the response of the shell versus the excitation frequency; Rx /Ry = 1 (spherical shell); fundamental mode (m = 1, n = 1), f˜ = 31.2 N and 1,1 = 0.004; model with 9 dofs. ——, Donnell’s theory; - -, Novozhilov’s theory.

response for different force excitations, approximately passes through the middle of the forced response curve computed in the present calculation) and shows a softening type non-linearity, turning to hardening for vibration amplitudes about two times larger than the shell thickness. Comparing Figs. 2(a) and (b), it is evident that displacement inwards is about two times larger than displacement outwards during a vibration period. Comparison of results obtained by using Novozhilov’s and Donnell’s non-linear shell theories is given in Fig. 3. Results are practically coincident in this case (m = 1, n = 1) and have been computed by using 9 dofs. A similar result has been recently obtained by Amabili [15] for a circular cylindrical panel also in case of multi-mode response with internal resonances. Results indicate that, at least for the thin shallow shell studied here, there is no advantage in using the more refined Novozhilov’s shell theory, which is much more time-consuming in computations. In order to check the convergence of the solution, the response has been calculated with a larger model, including 22 dofs. This model has the following additional generalized coordinates with respect to the 9

695

Fig. 4. Amplitude of the response of the shell versus the excitation frequency; Rx /Ry = 1 (spherical shell); fundamental mode (m = 1, n = 1), f˜ = 31.2 N and 1,1 = 0.004; Donnell’s theory. ——, 22 dofs model; - -, 9 dofs model; − · −, backbone curve from Kobayashi and Leissa [11].

dofs model: w1,3 , w3,1 , w3,3 , u1,5 , u3,5 , u5,1 , u5,3 , u5,5 , v1,5 , v3,5 , v5,1 , v5,3 and v5,5 . Comparison of the response computed with the 22 dofs and the 9 dofs models is given in Fig. 4, where the backbone curve of Kobayashi and Leissa [11] is also shown. The results of the 22 dofs model are moved slightly to the left with respect to the smaller 9 dofs model, and present a more complicate curve, specially in the frequency region around 0.9 1,1 . In fact, for excitation frequency = 0.9 1,1 , there is a 3:1 internal resonance with modes (m = 3, n = 1) and (m = 1, n = 3), giving 3 = 3,1 = 1,3 . A second relationship between natural frequencies that leads to internal resonances is for = 0.77 1,1 where 6 = 3,3 . This more complicate response, due to internal resonances, can be observed in Fig. 5(a)–(c) where the generalized coordinates involved are shown (for the symmetry of the spherical shell, w3,1 = w1,3 ). In synthesis, the main difference between the 9 dofs model and the 22 dofs model is related to the internal resonances that cannot be studied with the 9 dofs model because w1,3 , w3,1 and w3,3 are not included. Otherwise the trend of non-linearity does not

696

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

change much if investigated with the 9 or 22 dofs model. It is of great importance to analyse the results in the time domain by using direct integration of the equations of motion. Fig. 6 shows the most significant (i.e. those with larger amplitude) generalized coordinates for excitation frequency = 0.8 1,1 (with the 22 dofs model and initial conditions to catch the higher amplitude stable solution, which is in the softening type region). The generalized coordinates related to normal displacement w clearly indicate that there is a strong asymmetry between vibration inwards and outwards. In fact, the vibration is much larger inwards (negative) than outwards (positive), as previously observed, e.g. in Refs. [4,11]; this is related to the curvature of the shell. In particular, for the generalized coordinate w1,1 , which has the largest amplitude, the inward displacement is about 1.8 times bigger than the outward displacement. Moreover, there is a large contribution of higher harmonics, as clearly indicated by the frequency spectra in Fig. 7. The generalized coordinate w1,1 (associated to the resonant mode, with u1,1 and v1,1 ) is less affected by these higher harmonics. Phase plane diagrams in Fig. 8 are also provided to better understand the asymmetry of inward and outward vibration; loops are related to higher harmonics. 6.2. Case with Rx /Ry = 0.5

Fig. 5. Amplitude of the response of the shell versus the excitation frequency; Rx /Ry = 1 (spherical shell); fundamental mode (m = 1, n = 1), f˜ = 31.2 N and 1,1 = 0.004; 22 dofs model; Donnell’s theory. ——, stable solutions; - -, unstable solutions. (a) Maximum of w1,1 ; (b) maximum of w3,1 ; and (c) maximum of w3,3 .

A shallow shell with Rx /Ry =0.5 (Ry =2 m) is considered. The fundamental frequency 1,1 is 784.0 Hz according with Donnell’s shell theory. The magnitude of the harmonic force is f˜ = 22.85 N. Fig. 9 shows the maximum (outwards) and the minimum (inwards) of the shell response in z direction in the spectral neighbourhood of the fundamental mode (m = 1, n = 1) versus the excitation frequency. Calculations have been performed with the 9 dofs model. Results are compared to those obtained by Kobayashi and Leissa [11], where only the backbone curve is given. The present results with 9 dofs are quite close to those in Ref. [11] and shows a softening type nonlinearity, turning to hardening for vibration amplitudes about 1.5 times larger than the shell thickness. Comparing Figs. 9(a) and (b), it is evident that displacement inwards is larger than displacement outwards during a vibration period.

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

697

Fig. 6. Time-domain response of the spherical shallow shell; excitation frequency = 0.8 1,1 , f˜ = 31.2 N and 1,1 = 0.004; 22 dofs model; Donnell’s theory. (a) Excitation; (b) generalized coordinate w1,1 ; (c) generalized coordinate w3,1 ; (d) generalized coordinate w3,3 ; and (e) generalized coordinate u1,1 .

6.3. Case with Rx /Ry = 0, circular cylindrical shell A circular cylindrical shell is considered. The fundamental frequency 1,1 is 636.9 Hz according with Donnell’s shell theory. The magnitude of the harmonic force is f˜ = 10.68 N. This case has been extensively studied by Amabili [15] with the same approach; in particular, in reference [15] many additional results are given for this geometry.

Fig. 10 shows the maximum (outwards) and the minimum (inwards) of the shell response in z direction in the spectral neighbourhood of the fundamental mode (m = 1, n = 1) versus the excitation frequency. Calculations have been performed with the 9 dofs model. Results are compared to those obtained by Kobayashi and Leissa [11], where only the backbone curve is given. The present results with 9 dofs are very close to those in references [11] and shows a softening

698

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

Fig. 7. Spectrum of the response of the spherical shallow shell; excitation frequency = 0.8 1,1 , f˜ = 31.2 N and 1,1 = 0.004; 22 dofs. (a) Generalized coordinate w1,1 ; (b) generalized coordinate w3,1 ; (c) generalized coordinate w3,3 ; and (d) generalized coordinate u1,1 .

Fig. 8. Phase plane diagram of the response of the spherical shallow shell; excitation frequency = 0.8 1,1 , f˜ = 31.2 N and 1,1 = 0.004; 22 dofs. (a) Generalized coordinate w1,1 ; (b) generalized coordinate w3,1 ; and (c) generalized coordinate w3,3 .

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

699

Fig. 9. Amplitude of the response of the shell versus the excitation frequency; Rx /Ry = 0.5; fundamental mode (m = 1, n = 1), f˜ = 22.85 N and 1,1 = 0.004; model with 9 dofs; Donnell’s theory. ——, present stable results; - -, present unstable results; −·−, backbone curve from Kobayashi and Leissa [11]. (a) Maximum of w1,1 ; and (b) minimum of w1,1 .

Fig. 10. Amplitude of the response of the shell versus the excitation frequency; Rx /Ry = 0 (circular cylindrical shallow shell); fundamental mode (m = 1, n = 1), f˜ = 10.68 N and 1,1 = 0.004; model with 9 dofs; Donnell’s theory. ——, present stable results; -, present unstable results; − · −, backbone curve from Kobayashi and Leissa [11]. (a) Maximum of w1,1 ; and (b) minimum of w1,1 .

type non-linearity, turning to hardening for vibration amplitudes larger than the shell thickness. A comparison of Figs. 10(a) and (b) evidences that displacement inwards is larger than displacement outwards during a vibration period. Also, in this case results obtained by using Novozhilov’s and Donnell’s non-linear shell theories are practically coincident for shallow and deeper shells [15]. It can be also observed that for a circular cylindrical shell complete around the circum-

ference (closed shell), Amabili [14] has shown that significant differences arise between Donnell’s and Novozhilov’s theories only when in-plane inertia is neglected in Donnell’s theory. However, differences between Donnell’s and Novozhilov’s theories retaining in-plane inertia for complete and deep shells are slightly larger (but still negligible) than in the case the shallow shells, due to the fact that in these cases larger in-plane displacements arise.

700

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

tion in the spectral neighbourhood of the fundamental mode (m = 1, n = 1) versus the excitation frequency. Calculations have been performed with the 9 dofs model. Results are compared to those obtained by Kobayashi and Leissa [11], where only the backbone curve is given. The present results with 9 dofs are quite close to those in Ref. [11] and shows a hardening type non-linearity. Comparing Figs. 11(a) and (b) it can be observed that the displacement inwards (with respect to the centre of the smallest radius of curvature) is larger than displacement outwards during a vibration period. 6.5. Case with Rx /Ry = −1, hyperbolic paraboloid shell

Fig. 11. Amplitude of the response of the shell versus the excitation frequency; Rx /Ry = −0.5; fundamental mode (m = 1, n = 1), f˜ = 5.44 N and 1,1 = 0.004; model with 9 dofs; Donnell’s theory. ——, present stable results; - -, present unstable results; − · −, backbone curve from Kobayashi and Leissa [11]. (a) Maximum of w1,1 ; and (b) minimum of w1,1 .

6.4. Case with Rx /Ry = −0.5 A shallow shell with Rx /Ry = −0.5 (Ry = −2 m) is considered. In this case, the centres of principal curvatures are placed on opposite sides of the shell surface. The fundamental frequency 1,1 is 529.3 Hz according to Donnell’s shell theory. The magnitude of the harmonic force is f˜ = 5.44 N. Fig. 11 shows the maximum (outwards) and the minimum (inwards) of the shell response in z direc-

A shallow hyperbolic paraboloid shell with Rx /Ry = −1 (Ry = −1 m) is considered. The fundamental frequency 1,1 is 488.1 Hz according to Donnell’s shell theory. Other natural frequencies useful in the present study are (according with Donnell’s theory): 1,3 = 3,1 = 2528.7 Hz and 3,3 = 4396.6 Hz. The magnitude of the harmonic force is f˜ = 4.37 N. Fig. 12 shows the maximum and the minimum of the shell response in z direction in the spectral neighbourhood of the fundamental mode (m=1, n=1) versus the excitation frequency. Calculations have been performed with the 9 dofs model. Results are compared to those obtained by Kobayashi and Leissa [11], where only the backbone curve is given. The present results with 9 dofs are very close to those in Ref. [11] and shows a hardening type non-linearity. Comparing Figs. 12(a) and (b) it can be observed that the displacement in z and −z direction is perfectly symmetric during a vibration period. In order to check the convergence of the solution, the response has been calculated with the 22 dofs model. Comparison of the response computed with the 22 dofs and the 9 dofs models is presented in Fig. 13, where the backbone curve of Kobayashi and Leissa [11] is also given. The results of the 22 dofs model are moved slightly to the left with respect to the smaller 9 dofs model, showing good accuracy of the 9 dofs model in this case where no internal resonances appear. The response of the more significant generalized coordinates is given in Figs. 14(a)–(g) with stability indication.

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

701

Fig. 13. Amplitude of the response of the shell versus the excitation frequency; Rx /Ry = −1 (hyperbolic paraboloid shallow shell); fundamental mode (m = 1, n = 1), f˜ = 4.37 N and 1,1 = 0.004; Donnell’s theory. ——, 22 dofs model; - -, 9 dofs model; − · −, backbone curve from Kobayashi and Leissa [11].

Fig. 12. Amplitude of the response of the shell versus the excitation frequency; Rx /Ry = −1 (hyperbolic paraboloid shallow shell); fundamental mode (m = 1, n = 1), f˜ = 4.37 N and 1,1 = 0.004; model with 9 dofs; Donnell’s theory. ——, present stable results; -, present unstable results; − · −, backbone curve from Kobayashi and Leissa [11]. (a) Maximum of w1,1 ; and (b) minimum of w1,1 .

6.6. Effect of different curvature Fig. 15 synthesizes all the maximum responses presented in Sections from 6.1 to 6.4 for the 9 dofs model for different shell curvature aspect ratios Rx /Ry . It is clearly shown that for Rx /Ry = 1 (spherical), 0.5 and 0 (circular cylindrical) the shallow shell considered exhibits a softening type behaviour turning to hardening type for vibration amplitude of the order of magnitude of the shell thickness. The softening behaviour

becomes weaker with the decrement of the curvature aspect ratio Rx /Ry . In particular, the softening behaviour of the spherical shallow shell for vibration amplitude around 1.3 times the shell thickness is very strong. On the other hand, for Rx /Ry = −0.5, −1 (hyperbolic paraboloid) the shell has a strong hardening type behaviour without the presence of the softening type region. The effect of the curvature on the maximum response of spherical shallow shells is investigated in Fig. 16, obtained with the 9 dofs model and Donnell’s shell theory. Here the radius of curvature Rx = Ry is varied with respect to the value 1 m previously used, while all the other geometric and material characteristics are kept constant. In particular, three responses for Rx /a = 10, 20 and 100 are presented (solid line) and compared to the response of a flat square plate (dashed line) of the same dimension and material, also computed with the same theory [16] and the 9 dofs model. It must be observed that the flat plate is the limiting case for a spherical shallow shell when Rx /a → ∞. Fig. 16 shows that increasing the curvature ratio Rx /a the shell response changes gradually from (i) strongly

702

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

Fig. 14. Amplitude of the response of the shell versus the excitation frequency; Rx /Ry = −1 (hyperbolic paraboloid shallow shell); fundamental mode (m = 1, n = 1), f˜ = 4.37 N and 1,1 = 0.004; 22 dofs model; Donnell’s theory. ——, stable solutions; - -, unstable solutions. (a) Maximum of w1,1 ; (b) minimum of w1,1 ; (c) maximum of w3,1 ; (d) maximum of w1,3 ; (e) maximum of w3,3 ; (f) maximum of u1,1 ; and (g) maximum of v1,1 .

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

Fig. 15. Effect of the curvature aspect ratio Rx /Ry on the shell response (maximum of the generalized coordinate w1,1 ) versus the excitation frequency; fundamental mode (m = 1, n = 1); 1,1 = 0.004; 9 dofs model; Donnell’s theory.

softening turning to hardening for vibration amplitude about two times the shell thickness, to (ii) fully hardening, which is the typical behaviour of flat plates. It is curious that a flat plate (Rx /a → ∞) presents a slightly weaker hardening type behaviour than the spherical shallow shell with Rx /a =100. However, the difference is practically negligible.

7. Numerical results for buckling under static load The same shallow spherical shell studied in Section 6.1 is investigated here with the 22 dofs model. A static load is applied at the shell centre in −z direction and the shell deformations are reported in Figs. 17(a)–(c). A snap-through-type buckling instability is clearly identified around 350 N, where the load–displacement curve in Fig. 17(a) presents an unstable region and the shell deformation jumps from a relatively small value to a significantly larger one without any load increment. The shell deformations just before and after the snap-through are shown in Figs. 17(d) and (e).

703

Fig. 16. Effect of the curvature ratio Rx /a on the response (maximum of the generalized coordinate w1,1 ) of spherical shallow shells versus the excitation frequency; fundamental mode (m = 1, n = 1); 1,1 = 0.004; 9 dofs model; Donnell’s theory. For Rx /a = 10, 1,1 = 952.3 Hz and f˜ = 31.2 N; for Rx /a = 20, 1,1 =637.1 Hz and f˜ =11.16 N; for Rx /a =100, 1,1 =495.4 Hz and f˜ = 5.56 N; for Rx /a = ∞ (square plate), 1,1 = 488.6 Hz and f˜ = 5.4 N.

8. Bifurcation diagrams varying the amplitude of the harmonic excitation The same shallow spherical shell studied in Section 6.1 is considered here and the 22 dofs model is used. Poincarè maps have been computed by direct integration of the equations of motion. The excitation frequency has been kept constant, = 0.8 1,1 , and the excitation amplitude has been varied between 0 and 1400 N. The force range has been divided into 800 steps, so that the force is varied of 1.75 N at each step. 600 periods have been waited each time the force is changed of a step in order to eliminate the transient motion. The initial condition at the first step is zero displacement and velocity for all the variables; at the following steps the solution at the previous step, with addition of a small perturbation in order to find stable solution, is used as initial condition. The bifurcation diagrams obtained by all these Poincarè maps are shown in Figs. 18 and 19 . In particular, in Fig. 18 the load is increased from 0 to 1400 N; in Fig. 19 the

704

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

Fig. 17. Deformation of the spherical shallow shell under static concentrated load at the shell centre in −z direction; 22 dofs model. (a) Generalized coordinate w1,1 ; (b) generalized coordinate w3,1 ; (c) generalized coordinate w3,3 ; (d) deformed shell shape for load of 348.9 N (just before the snap-through instability; and (e) deformed shell shape for load of 350 N (just after the snap-through instability).

load is decreased from 1400 N to 0. Simple periodic motion, period doubling bifurcation, subharmonic response, amplitude modulations and chaotic response have been detected, as indicated in Figs. 18(a) and 19(a) and (c). This indicates a very rich and complex non-linear dynamics of the spherical shallow shell subject to large harmonic excitation. Comparing Figs. 18 and 19, where only stable solutions are reported, it is clearly visible that only in some regions the shell presents the same behaviour when the load (bifurcation parameter) is increased or decreased. This indicates that different stable solutions coexist for the same set of system parameters, so that the solution is largely affected by initial conditions. In particular, Figs. 19(b) and (d) give the maximum Lyapunov exponent 1 associated to the bifurcation diagram. It can be easily

observed that (i) for periodic forced vibrations 1 < 0, (ii) for amplitude modulated response 1 = 0 and (iii) for chaotic response 1 > 0. Therefore 1 can be conveniently used for identification of the system dynamics. A three dimensional representation of the bifurcation diagram in the displacement, velocity and load space is shown in Fig. 19(e), while in the other figures the bifurcation diagrams have been projected on a plane orthogonal to the velocity axis. All the Lyapunov exponents have been evaluated for the case with excitation f˜ = 1396 N in Fig. 19, corresponding to chaotic response, and are given in Table 1 and Fig. 20. In this case, four positive Lyapunov exponents have been identified, allowing to classify this response as hyperchaos. The Lyapunov dimension in this case is dL = 24.59.

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

705

Fig. 18. Bifurcation diagram of Poincaré maps for the spherical shallow shell under increasing harmonic load with frequency = 0.8 1,1 ; 1,1 = 0.004; 22 dofs model; Donnell’s theory. (a) Generalized coordinate w1,1 ; T = response period equal to excitation period; PD = period doubling bifurcation; 2T = periodic response with two times the excitation period; M = amplitude modulations; C = chaos; (b) generalized coordinate w3,1 ; and (c) generalized coordinate w3,3 .

The case of excitation of 910 N gives chaotic behaviour increasing or decreasing the load, as shown by Figs. 18 and 19; the maximum Lyapunov exponent for this case is 22.7. The time response for this case is re-

ported in Fig. 21. The frequency spectra and Poincarè maps are given in Figs. 22 and 23, respectively. Figs. 21–23 show typical chaotic behaviour of the spherical shallow shell for this specific harmonic excitation.

706

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

Fig. 19. Bifurcation diagram of Poincaré maps and maximum Lyapunov exponent for the spherical shallow shell under decreasing harmonic load with frequency = 0.8 1,1 ; 1,1 = 0.004; 22 dofs model; Donnell’s theory. (a) Bifurcation diagram: generalized coordinate w1,1 ; T = response period equal to excitation period; PD = period doubling bifurcation; M = amplitude modulations; C = chaos; (b) maximum Lyapunov exponent; (c) bifurcation diagram: generalized coordinate w1,1 , enlarged scale; T = response period equal to excitation period; PD = period doubling bifurcation; 2T = periodic response with two times the excitation period; 5T = periodic response with five times the excitation period; M = amplitude modulations; and C = chaos; (d) maximum Lyapunov exponent, enlarged scale; (e) 3-D representation of the bifurcation diagram: generalized coordinate w1,1 , enlarged scale; (f) bifurcation diagram: generalized coordinate w3,1 ; and (g) bifurcation diagram: generalized coordinate w3,3 .

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

707

Fig. 19. (continued). Table 1 All the Lyapunov exponents for the spherical shallow shell; excitation frequency = 0.8 1,1 , f˜ = 1396 N and 1,1 = 0.004; 22 dofs Positive Lyapunov exponents

Negative Lyapunov exponents

34.2, 23.9, 14.0, 5.29

−3.64, −3.67, −3.68, −3.69, −3.73, −3.74, −3.74, −3.76, −3.76, −3.76, −3.76, −3.77, −3.77, −3.77, −3.80, −3.81, −3.81, −3.82, −3.82, −3.82, −3.83, −3.83, −3.84, −3.84, −3.85, −3.85, −3.86, −3.86, −3.86, −3.87, −3.87, −3.89, −3.90, −3.91, −3.95, −4.02, −13.1, −21.5, −31.5, −41.8

Fig. 20. All the 44 Lyapunov exponents for the spherical shallow shell; excitation frequency = 0.8 1,1 , f˜ = 1396 N and 1,1 = 0.004; 22 dofs.

708

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

Fig. 22. Spectrum of the response of the spherical shallow shell; excitation frequency = 0.8 1,1 , f˜ = 910 N and 1,1 = 0.004; 22 dofs; Donnell’s theory. (a) Generalized coordinate w1,1 ; and (b) generalized coordinate w1,3 .

9. Conclusions

Fig. 21. Time-domain response of the spherical shallow shell; excitation frequency = 0.8 1,1 , f˜ = 910 N and 1,1 = 0.004; 22 dofs model; Donnell’s theory. (a) Excitation; (b) generalized coordinate w1,1 ; and (c) generalized coordinate w1,3 .

The present study introduces a multi-mode expansion and uses accurate shell theories retaining in-plane inertia to study large-amplitude, forced vibrations of doubly curved shallow shells. This overcomes two frequent limitations in previous studies: (i) the use of mode expansions with one or two degrees of freedom; and (ii) the use of less accurate, but simpler, Donnell’s shallow-shell theory, which neglects in-plane inertia. The occurrence of internal resonances in the problem studied is a clear indication that this important non-linear phenomenon has fundamental importance in the study of curved shells. Internal resonances can be studied only with multi-mode expansions, in some cases with a quite large number of degrees of freedom. Bifurcation diagrams constructed by Poincarè maps and Lyapunov exponents show period doubling bifurcations and highly complex non-linear behaviour of

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

709

Acknowledgements This work was partially supported by the FIRB 2001 and COFIN 2003 grants of the Italian Ministry for University and Research (MIUR).

Appendix A. Calculation of additional terms to satisfy in-plane boundary conditions Integration of Eq. (22) gives three contributions, related to the 2nd, 3rd and 4th term, respectively, in Eq. (22) u(t) ˆ = uˆ 1 (t) + uˆ 2 (t) + uˆ 3 (t),

(A.1)

where 1 uˆ 1 (t) = − 2 = − ×

 


jw jx

2  dx x=0,a

M N  1

(m /a) wm,n (t) sin(n y/b) 2

n=1 m=1 N M

k=1 s=1

s ws,k (t) sin(k y/b) m+s

× sin[(m + s) x/a]} ,    jw jw0 dx uˆ 2 (t) = − jx jx x=0,a = − ×

M N

n=1 m=1 M˜ N˜

j =1 i=1

Fig. 23. Poincaré maps for the spherical shallow shell; excitation frequency = 0.8 1,1 , f˜ = 910 N and 1,1 = 0.004; 22 dofs. (a) Generalized coordinate w1,1 ; and (b) generalized coordinate w1,3 .

a spherical shallow shell under very large harmonic excitation. Up to four positive Lyapunov coefficients have been found for the studied case, indicating hyperchaos.

(A.2)

 (m /a) wm,n (t) sin(n y/b) i Ai,j sin(j y/b) m+i

× sin[(m + i) x/a]} ,    jv dx. uˆ 3 (t) = − jy x=0,a

(A.3) (A.4)

It can be observed that solution (A.2) is not unique, but the one proposed satisfies boundary conditions and deformation energy considerations. Then, in order to satisfy boundary conditions (15a), Eq. (A.4) gives    jv uˆ 3 (t) = − dx = 0. (A.5) jy x=0,a

710

M. Amabili / International Journal of Non-Linear Mechanics 40 (2005) 683 – 710

Similar expressions are obtained by integration of Eq. (23). Finally Eqs. (24a) and (25a) are obtained, which satisfy condition (A.5). References [1] M. Amabili, M.P. Païdoussis, Review of studies on geometrically nonlinear vibrations and dynamics of circular cylindrical shells and panels, with and without fluid-structure interaction, Appl. Mech. Rev. 56 (2003) 349–381. [2] E. Reissner, Nonlinear effects in vibrations of cylindrical shells, Ramo-Wooldridge Corporation Report AM5-6, 1955. [3] E.I. Grigolyuk, Vibrations of circular cylindrical panels subjected to finite deflections, Prikl. Mat. Mekh. 19 (1955) 376–382 (in Russian). [4] A.W. Leissa, A.S. Kadi, Curvature effects on shallow shell vibrations, J. Sound Vibr., 16 (1971) 173–187. [5] A.P. Bhattacharya, Nonlinear flexural vibrations of thin shallow translational shells, ASME J. Appl. Mech. 43 (1976) 180–181. [6] G.C. Sinharay, B. Banerjee, Large-amplitude free vibrations of shallow spherical shell and cylindrical shell—a new approach, Int. J. Non-Linear Mech. 20 (1985) 69–78. [7] G.C. Sinharay, B. Banerjee, Large-amplitude free vibrations of shells of variable thickness—a new approach, AIAA J. 24 (1986) 998–1004. [8] C.Y. Chia, Nonlinear analysis of doubly curved symmetrically laminated shallow shells with rectangular platform, Ing.-Arch. 58 (1988) 252–264. [9] D. Hui, Accurate backbone curves for a modified-Duffing equation for vibrations of imperfect structures with viscous damping, ASME J. Vibr. Acoust. 112 (1990) 304–311. [10] L. Librescu, M.-Y. Chang, Effects of geometric imperfections on vibration of compressed shear deformable laminated composite curved panels, Acta Mech. 96 (1993) 203–224.

[11] Y. Kobayashi, A.W. Leissa, Large amplitude free vibration of thick shallow shells supported by shear diaphragms, Int. J. Non-Linear Mech. 30 (1995) 57–66. [12] A. Abe, Y. Kobayashi, G. Yamada, Non-linear vibration characteristics of clamped laminated shallow shells, J. Sound Vibr. 234 (2000) 405–426. [13] M.S. Soliman, P.B. Gonçalves, Chaotic behaviour resulting in transient and steady state instabilities of pressure-loaded shallow spherical shells, J. Sound Vibr. 259 (2003) 497–512. [14] M. Amabili, Comparison of shell theories for large-amplitude vibrations of circular cylindrical shells: Lagrangian approach, J. Sound Vibr. 264 (2003) 1091–1125. [15] M. Amabili, Nonlinear vibrations of circular cylindrical panels, J. Sound Vibr., in press. [16] M. Amabili, Nonlinear vibrations of rectangular plates with different boundary conditions: theory and experiments, Comput. Struct., in press. [17] V.V. Novozhilov, Foundations of the Nonlinear Theory of Elasticity, Graylock Press, Rochester, USA, 1953 Now available from Dover, N.Y., USA. [18] A.W. Leissa, Vibration of Shells, NASA SP-288, Government Printing Office, Washington, DC, 1973. Now available from The Acoustical Society of America, 1993. [19] S. Wolfram, The Mathematica Book, 4th ed., Cambridge University Press, Cambridge, UK, 1999. [20] E.J. Doedel, A.R. Champneys, T.F. Fairgrieve, Y.A. Kuznetsov, B. Sandstede, X. Wang, AUTO 97: Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont). Concordia University, Montreal, Canada, 1998. [21] J. Argyris, G. Faust, M. Haase, An Exploration of Chaos, North-Holland, Amsterdam, 1994. [22] F.C. Moon, Chaotic and Fractal Dynamics, Wiley, New York, 1992. [23] T. Yamaguchi, K. Nagai, Chaotic vibrations of a cylindrical shell-panel with an in-plane elastic-support at boundary, Nonlinear Dynamics 13 (1997) 259–277.