A new non-linear higher-order shear deformation theory for large-amplitude vibrations of laminated doubly curved shells

Nonlinear vibrations Rotary inertia , A consistent higher-order shear deformation non-linear theory is developed for shells of generic shape, taking geometric imperfections into account. The geometrically non-linear strain–displacement relationships are derived retaining full non-linear terms in the in-plane displacements; they are presented in curvilinear coordinates in a formulation ready to be implemented. Then, large-amplitude forced vibrations of a simply supported, laminated circular cylindrical shell are studied (i) by using the developed theory, and (ii) keeping only non-linear terms of the von Karman type. Results show that inaccurate results are obtained by keeping only non-linear terms of the von Ka´ rma´ n type for vibration amplitudes of about two times the shell thickness for the studied case.


Introduction
Classical shell theories, which neglect shear deformation and rotary inertia, give inaccurate natural frequencies for moderately thick or laminated anisotropic shells and plates (e.g. [1]). In order to overcome this limitation, shear deformation theories have been introduced. These theories can be classified as first-and higherorder shear deformation theories [2]; in the first category, a shear correction factor is required for the equilibrium since a uniform shear strain is assumed through the shell thickness. Higher-order shear deformation theories overcome this limitation since a realistic shear stress distribution through the shell thickness is assumed, which also satisfies the condition of zero shear stresses at both top and bottom shell surfaces.
Several higher-order shear deformation shell theories have been proposed. Librescu [3] developed a non-linear shell theory by expanding the shell displacements with cubic terms in the transverse coordinate. A linear higher-order shear deformation theory of shells has been introduced by Reddy [4] and Reddy and Liu [5]. Arciniega and Reddy [6] have improved the theory developed in [5]. A review of shell theories has been presented by Reddy and Arciniega [7].
Reddy [8] developed the non-linear higher-order shear deformation theory of plates, taking into account von Ká rmá n type non-linear terms. Dennis and Palazotto [9] and Palazotto and Dennis [10] have extended Reddy shell theory [4] to non-linear deformation by introducing the von Ká rmá n type non-linear terms. These theories have been also discussed in the books of Amabili [2] and Reddy [11].
In the existing higher-order shear deformation geometrically non-linear shell theories, the von Ká rmá n type non-linear terms (i.e. those involving the normal displacement only) have been added to the linear equations, so that consistent derivation has not been performed. Moreover, the von Ká rmá n type non-linear terms are known for being accurate for classical shell theories only for small displacements and moderately small rotations, as shown by Amabili [12] in case of vibrations. Therefore, it is important to derive a consistent higher-order shear deformation theory that keeps all the non-linear terms in the normal and inplane displacements.
In the present study, a consistent higher-order shear deformation non-linear theory is developed for shells of generic shape; geometric imperfections are also taken into account. The geometrically non-linear strain-displacement relationships are derived retaining full non-linear terms in the in-plane displacements. They are presented in curvilinear coordinates in a formulation that can be readily implemented. Then, largeamplitude forced vibrations of a simply supported, laminated circular cylindrical shells are studied (i) by using the developed theory, and (ii) keeping only non-linear terms of the von Ká rmá n type. Results show that inaccurate results are obtained by keeping only non-linear terms of the von Ká rmá n type for vibration amplitudes of about two times the shell thickness for the studied case.

Non-linear higher-order shear deformation theory
A laminated shell of arbitrary shape, made of a finite number of orthotropic layers, oriented arbitrarily with respect to the shell principal curvilinear coordinates (a 1 , a 2 ), is considered, as shown in Fig. 1. The development of the theory remains the same for shells made of isotropic, orthotropic or functionally graded materials. The displacements of an arbitrary point (a 1 , a 2 ) on the middle surface of the shell are denoted by u, v and w, in the a 1 , a 2 and z directions, respectively; w is taken positive outward from the center of the smallest radius of curvature. Initial geometric imperfections of the shell associated with zero initial tension are denoted by displacement w 0 in normal direction, also taken positive outward. The thickness h of the shell is assumed to be small compared to the principal radii of curvature of the shell, so that only moderately thick shells can be considered. The displacements (u 1 , u 2 , u 3 ) of a generic point are related to the middle surface displacements by where f 1 and f 2 are the rotations of the transverse normals at z= 0 about the a 2 and a 1 axes, respectively, and c 1 , c 2 , g 1 , g 2 , y 1 , y 2 and w are functions to be determined in terms of u, v, w, f 1 and f 2 . Thus, the five variables describing the shell deformation are u, v, w, f 1 and f 2 . In Eqs. (1) the in-plane displacements have been expanded up to the 4th order in z while the normal displacement has been assumed to be linear in z; Eqs. (1a,b) give a high-order distribution of shear effects through the thickness. The shear and normal Green's strains for three-dimensional elasticity are [2] Eqs. (2a-c) are non-linear; in Eqs. (2a,b), R 1 and R 2 (functions of the coordinates a 1 and a 2 ) are the principal radii of curvature in a 1 and a 2 directions, respectively, and A 1 and A 2 are the Lamé parameters. The shear deformation, see Eqs. (2a,b), is neglected in classical shell theories, which is a very good approximation for thin and moderately thick isotropic shells and for very thin laminated shells [2,12]. For laminated shells that cannot be considered very thin, shear deformation should be retained in order to obtain accurate results. However, shear deformation plays a smaller role than in-plane strains and bending (since it is even neglected in classical theories), so that it is reasonable to keep only linear displacement terms in Eqs. (2a,b). This assumption leads to the possibility of expressing the displacement field in Eqs. (1a-c) by using only linear expressions in u, v, w, f 1 and f 2 .
Therefore, the following relationships are obtained for the transverse shear strains: It is also assumed that normals to the middle surface of the shell before deformation are not elongated after deformation. This is mathematically expressed by The constraint given by Eq. (4) could be removed and the shell theory will become a six-variable theory, instead of having five variables, w being the additional variable.
The expressions for the transverse shear strains are obtained by substituting Eqs. (1a-c) into (3a,b); the vanishing of the shear strains at the top and the bottom surfaces of the shell requires The following approximation for moderately thick shell is introduced Collecting first-and third-order terms in h from Eq. (5a) into two separate equations, the following expressions are obtained since these expressions change sign if evaluated at z ¼ 7h=2 and must be set equal to zero Eqs. (7a,b) give Then, Eq. (5a) gives Similarly, from Eq. (5b), the following expressions are obtained: By using Eq. (4), the function w can be determined where the approximation in Eq. (10a) is accurate for f 1 ; f 2 51.
Eq. (10a) shows that w is a second-order term with respect to the other five variables, i.e. displacements and rotations, when the constraint given by Eq. (4) is introduced. Therefore, it is possible to approximate wC0: By substituting Eqs. (1a-c) into (3a,b) and using the approximations (6) and Eqs. (8), (9), (10b), the following straindisplacement relationships are obtained for the shear strains keeping terms up to z 3 : 13 þ z 2 k ð2Þ 13 Þ; ð11aÞ Eqs. (11)- (13) show that the strain and stress distribution through the thickness is parabolic and therefore the shear correction factor is no longer required.
The expressions of the other Green's strain tensor e ij in the curvilinear coordinate system, are [2] e 11 ¼ e 1 þ 1 where ; ð15aÞ The strain-displacement equations for the higher-order shear deformation theory to be added to Eqs. (11a,b), keeping terms up to z 3 and using approximations (6) and (10), are obtained by substituting Eqs. (1), (8), (9), (10b) into Eqs. (14) and (15) Eqs. (18a-i) give the changes in curvature and torsion of the middle surface, and they have been obtained retaining only linear terms; in fact, non-linear terms in the changes in curvature and torsion play a very small role, at least for moderate vibration amplitudes of the order of the shell thickness [12]. Eqs. (17a-c), giving the middle surface strains, are coincident with those obtained by using Novozhilov non-linear shell theory [2,13], which neglects shear deformation and rotary inertia. Moreover, it can be observed that k ð1Þ 1 ; k ð1Þ 2 ; k ð1Þ 12 are negligible for shells that are not very thick.

Elastic strain and kinetic energies, including rotary inertia, for laminated shells
It is assumed that the transverse normal stress s 3 ¼ 0 is negligible (plane stress assumption); in general, it is verified that s 3 is small compared to t 13 and t 23 , except near the shell edges, so that the hypothesis is a good approximation of the actual behavior of moderately thick shells. The stress-strain relations for the k-th orthotropic lamina of the shell, in the material principal coordinates (x,y,z), are obtained under the hypothesis where G xy , G xz and G yz are the shear moduli in x-y, x-z and y-z directions, respectively, and the coefficients c ij are given in Appendix A and in Ref. [11]; t xz and t yz are the shear stresses and the superscript (k) refers to the k-th layer within a laminate.
where ½Q ðkÞ is the 5 Â 5 matrix of the material properties of the k-th layer expressed in the shell principal coordinates and it is given in Appendix A. Eq. (20) can be rewritten as : In Eq. (21), k ð1Þ 1 ; k ð1Þ 2 ; k ð1Þ 12 have been neglected, supposing a moderately thick shell.
The elastic strain energy U S of the shell is given by ðs ðkÞ 1 e 11 þ s ðkÞ 2 e 22 þ t ðkÞ 12 g 12 þ t ðkÞ 13 g 13 þ t ðkÞ 23 g 23 Þ where K is the total number of layers in the laminated shell and (h (k À 1) , h (k) ) are the z coordinates of the k-th layer. For simplicity a shell of rectangular base of dimension a and b in a 1 and a 2 directions, respectively, has been considered. The kinetic energy T S of the shell, including rotary inertia, is given by The z and z 3 terms vanish after integration on z in the case of a laminate with symmetric density with respect to the z axis.
In particular, for a laminate with the same density for all the layers and uniform thickness, the following simplified expression is obtained:

Boundary conditions and discretization of a circular cylindrical shell
In order to reduce the system to finite dimensions, the middle surface displacements u, v and w are expanded by using approximate functions. Circular cylindrical shells with simply supported boundary conditions are analyzed in the following part of the study, as shown in Fig. 2. In particular, The following boundary conditions are imposed at the shell ends, x=0, L: where N x is the axial stress resultant per unit length and M x is the axial stress moment resultant per unit length, i.e., s ðkÞ Moreover, the three displacements and the two rotations must be continuous in y.
The following base of shell displacements, which satisfy identically the geometric boundary conditions (25a-c), is used to discretize the system: uðx; y; tÞ ¼  The minimal expansion giving accurate results for excitation in the neighborhood of resonance of mode (m=1, n) is uðx; y; tÞ ¼ ½u 1;n;c ðtÞcosðnyÞþu 1;n;s ðtÞsinðnyÞcosðl 1 xÞ This expansion has 23 generalized coordinates (degrees of freedom) and guarantees good accuracy for the calculation performed in the present work. Initial geometric imperfections of the shell are considered only in radial direction. They are assumed to be associated with zero initial stress. The imperfection w 0 is expanded in the same form of w w 0 ðx; yÞ ¼   terms at x=0, L, the following expression is obtained: where the termûðx; y; tÞ has been neglected in second-order terms. In Eq. (31), all the linear terms ð@v=@xÞþð1=RÞð@u=@yÞ, k ð0Þ xy and k ð2Þ xy can be eliminated since they establish a linear relationship which is satisfied by using the minimization of energy in the process of building the Lagrange equations of motion; in fact, this is equivalent to the Rayleigh-Ritz method and therefore it is necessary only to satisfy only geometrical boundary conditions. Therefore Eq. (31a) can be simplified into Â h ðkÞ Àh ðkÀ1Þ þ h ðkÞ2 Àh ðkÀ1Þ2 2R which immediately gives By introducing the notation Eq. (31c) can be expressed in the following form In Eq. (31d) it is necessary to retain only the resonant mode (m,n). The expression ofû can be obtained from Eq. (31d) aŝ Â½w m;n;c ðtÞcosðnyÞþw m;n;s ðtÞsinðnyÞ X N j ¼ 0 where aðtÞ ¼ ðmp=LÞðw 2 m;n;c þ w 2 m;n;s þ v 2 m;n;c þ v 2 m;n;s Þ þðD 12 =D 11 Þ½Ln 2 =ðmpR 2 Þðu 2 m;n;c þ u 2 m;n;s Þ; bðtÞ ¼ ðmp=LÞðw 2 m;n;c Àw 2 m;n;s þ v 2 m;n;s Àv 2 m;n;c Þ þðD 12 =D 11 Þ½Ln 2 =ðmpR 2 Þðu 2 m;n;s Àu 2 m;n;c Þ; cðtÞ ¼ ð2mp=LÞðw m;n;c w m;n;s þv m;n;c v m;n;s Þ À2ðD 12 =D 11 Þ½Ln 2 =ðmpR 2 Þu m;n;c u m;n;s : The boundary condition (25e) is identically satisfied for symmetric laminates if the term z/R is neglected in Eq. (26), i.e. for thin shells. This is due to the expressions of k ð0Þ x and k ð2Þ x given in Appendix B, which are zero at x=0, L for the expansions assumed in Eqs. (27a-e). Additional terms must be added to the expansion of the in-plane displacement u for asymmetric laminates and moderately thick shells. In fact, bending and stretching are coupled for asymmetric laminates.

Lagrange equations of motion
The virtual work W done by the external forces is written as where q x , q y and q z are the distributed forces per unit area acting in x, y and z directions, respectively, applied at the middle surface. Only a single harmonic force orthogonal to the shell is considered; therefore, q x = q y = 0. The external distributed load q z applied to the shell, due to the normal concentrated forcef , is given by where o is the excitation frequency, t is the time, d is the Dirac delta function,f gives the force magnitude positive in z direction andx andỹ give the position of the point of application of the force. Here, the point excitation is located at the center of shell, that is,x ¼ L=2,ỹ ¼ 0. Eq. (33) can be rewritten in the following form: The non-conservative damping forces are assumed to be of viscous type and are taken into account by using the Rayleigh dissipation function where c has a different value for each term of the mode expansion; in particular In Eq. (37) displacements are non-dimensionalized dividing by h, while rotations are already non-dimensional. The damping coefficient c m,n is related to the modal damping ratio that can be evaluated from experiments by z m;n ¼ c m;n =ð2m m;n o m;n Þ, where o m;n is the natural circular frequency of mode (m, n) and m m,n is the modal mass of this mode.
The following notation is introduced for brevity: q ¼ fu m;n ; v m;n ; w m;n ; f 1m;n ; f 2m;n g T ; m ¼ 1; . . . ; M orM and n The generic element of the time-dependent vector q is referred to as q j ; the dimension of q is N, which is the number of degrees of freedom (dofs) used in the mode expansion. The generalized forces Q j are obtained by differentiation of the Rayleigh dissipation function and of the virtual work done by external forces: The Lagrange equations of motion are d dt where @T P =@q j ¼ 0. The complicated term, derived from the maximum potential energy of the plate, giving quadratic and cubic non-linearities, can be written in the form where the coefficients f have long expressions that include also geometric imperfections. It is interesting to observe that in Eq. (41) there are quadratic and cubic terms.

Inertial coupling in the equations of motion
For shells with rotary inertia, inertial coupling arises in the equations of motion (see Eqs. (23) and (24)) so that they cannot be immediately transformed in the form required for numerical integration. In particular, the equations of motion take the following form: where M is the non-diagonal mass matrix of dimension N Â N (N being the number of degrees of freedom), C is the damping matrix, K is the linear stiffness matrix, which does not present terms involving q, N 2 is a matrix that involves linear terms in q, therefore giving the quadratic non-linear stiffness terms, N 3 is a matrix that involves quadratic terms in q, therefore giving the cubic non-linear stiffness terms, f 0 is the vector of excitation amplitudes and q is the vector of the N generalized coordinates, defined in Eq. (38). In particular, by using Eq. (41), the generic elements k j,i , n 2 j;i and n 3 j;i , of the matrices K, N 2 and N 3 , respectively, are given by Eq. (42) is pre-multiplied by M À1 in order to diagonalize the mass matrix, as a consequence that the matrix M is always invertible; the result is ð44Þ which can be rewritten in the following form I € q þC _ q þ½K þM À1 N 2 ðqÞþM À1 N 3 ðq; qÞq ¼f 0 cosðotÞ; ð45Þ Eq. (45) is in the form suitable for numerical integration.

Numerical results
The dimensions and material properties of a simply supported, imperfection-free, laminated circular cylindrical shell made of graphite/epoxy plain weave fabric are: L= 520 mm, R= 150 mm, h=1 mm, E 1 = E 2 =55 Â 10 9 Pa, G 12 =4.48 Â 10 9 Pa, G 13 = G 23 =4.14 Â 10 9 Pa, n 12 = n 21 = 0.04 and r=1500 kg/m 3 ; the fundamental mode (m =1, n= 4) is considered, with natural frequency of 309.7 Hz; it has n =4 circumferential waves. Numerical results have been obtained by using the software AUTO [14] for continuation and bifurcation analysis by using the pseudo-arc length continuation method. Fig. 3 presents the frequency-response curve for forced vibrations of the shell around the resonance of the fundamental mode computed by using the 3rd order shear deformation non-linear theory introduced in Section 2. A modal damping coefficient z 1;n ¼ 0:001 and harmonic point force excitation of magnitude 3.04 N and frequency o applied at shell mid-length are assumed. The main branch 1 in Fig. 3 corresponds to vibration with zero amplitude of the companion mode w 1;n;s ðtÞ. This branch has pitchfork bifurcations (BP) at o=o 1;n ¼ 0:9905 and at 0.9997, where branch 2 appears. This new branch corresponds to participation of both w 1;n;c ðtÞ and w 1;n;s ðtÞ, giving a traveling-wave response. The companion mode presents a node at the location of the excitation force and therefore it is not directly excited; its amplitude is different than zero only for largeamplitude vibrations, due to non-linear coupling. In the narrow frequency region where both w 1;n;c ðtÞ and w 1;n;s ðtÞ are different from zero, they give rise to a traveling wave around the shell; phase shift between the two coordinates is almost equal to p/2. The appearance of branch 2 is related to the 1:1 internal resonance of w 1;n;c ðtÞ and w 1;n;s ðtÞ, which is due to the axial symmetry of the circular cylindrical shell. This branch appears for sufficiently large excitation, and it can be observed for vibration amplitude of the order of 1/10 of the shell thickness. Branch 2 undergoes two Neimark-Sacker (torus) bifurcations (TR), at o=o 1;n ¼ 0:9909 and 0.9965. Amplitudemodulated (quasi-periodic) response is indicated in Fig. 3 on branch 2 for 0:9909 o o=o 1;n o0:9965, that is, bracketed by the two Neimark-Sacker bifurcations.
A comparison between the results obtained by using the present theory and those obtained keeping only von Ká rmá n nonlinear terms in Eqs. (B.1)-(B.3), or in Eqs. (17a-c), is shown in Fig. 4. Results indicate that, for the studied case, it is important to retain non-linearities in the in-plane displacements in order to obtain accurate results. In fact, the difference between the two calculations is significant from the quantitative point of view, even if the qualitative behavior is the same.

Conclusion
The proposed theory is the first 3rd order shear deformation theory retaining non-linearities in the in-plane displacements. Moreover, it has been derived with consistency without introducing nonlinearities of the von Kármán at the end of the derivation of the theory. Numerical results show that, for the studied case, these nonlinear terms are important to predict with accuracy the non-linear forced response of a laminated graphite/epoxy plain weave fabric circular cylindrical shell.