## Abstract

Reduced-order models are essential to study nonlinear vibrations of structures and structural components. The natural mode discretization is based on a two-step analysis. In the first step, the natural modes of the structure are obtained. Because this is a linear analysis, the structure can be discretized with a very large number of degrees of freedom. Then, in the second step, a small number of these natural modes are used to discretize the nonlinear vibration problem with a huge reduction in the number of degrees of freedom. This study finds a recipe to select the natural modes that must be retained to study nonlinear vibrations of an angle-ply laminated circular cylindrical shell that the author has previously studied by using admissible functions defined on the whole structure, so that an accuracy analysis is performed. The higher-order shear deformation theory developed by Amabili and Reddy is used to model the shell.

## 1. Introduction

Reduced-order models are essential to study nonlinear vibrations of structures and structural components. In fact, large models with many thousands or more degrees of freedom, for example obtained by commercial finite-element codes, cannot be used to calculate all the branches of the forced vibration response around a resonance of a structure. This limitation is mainly due to the huge computational problem involved and to the stiff dynamic equations obtained for large dimensional systems, because very different time scales are associated with different degrees of freedom.

A few techniques exist to build reduced-order models for nonlinear vibration problems. The most diffused ones are probably the nonlinear normal modes (under this name are classified the techniques based on the centre manifold theorem, the normal form theory and the inertial manifold) [1,2,3,4,5,6,7,8,9], including the most diffused version with asymptotic approach, the discretization of the equations of motion by using global (i.e. defined on the whole structure) admissible functions [10,11,12,13], the proper orthogonal decomposition method [14,15,16,17,18] and the natural mode discretization [19,20,21,22,23].

The natural mode discretization is based on a two-step analysis. In the first step, the natural modes of the structure are obtained. Because this is a linear analysis, the structure can be discretized with a very large number of degrees of freedom; for example, by using the finite-element method (FEM) or the Rayleigh–Ritz method. Then, in the second step, these natural modes, which are clearly defined on the whole structure (eventually using interpolating functions), are used to discretize the nonlinear vibration problem. The advantage is that many thousands or more degrees of freedom can be used in the linear analysis, but then just a few tens of natural modes can be retained in the nonlinear analysis with a huge reduction in the number of degrees of freedom. In fact, usually, the nonlinear vibration problems do not need many degrees of freedom to be accurately described, but the difficulty is in the choice of them. The use of natural modes, which are quite intuitive, simplifies this choice even if it does not allow the smallest dimensional problem to be obtained, as, for example, is possible by using the nonlinear normal modes. However, the advantage of having a technique that could be implemented by using natural modes obtained by commercial FEM codes is huge, and the dimension of the problem obtained by natural modes is still very small. In addition, it overcomes the limitation of the nonlinear normal modes computed with an asymptotic approach, which give inaccurate responses for large vibration amplitudes [24].

It must be clarified that, at present, there is no technique for the choice of the natural modes to be retained in the nonlinear vibration problem. This study uses the experience of the author in building reduced-order models of shell structures to find a recipe to select the natural modes that must be retained for circular cylindrical shells complete around their circumference.

The recipe works for complete circular cylindrical shells with uniform boundary conditions, isotropic or laminated, even in the presence of moderate geometrical imperfections. However, the case presented is an angle-ply laminated circular cylindrical shell that the author has previously studied by using global admissible functions [25] and here is studied by using natural modes, so that the accuracy of the solution can be studied.

Nonlinear vibrations of shells have been studied by many researchers [26,27,28,29]. Natural modes have never been used before for circular cylindrical shells, and the expansions used to discretize the shell had independent variables for each one of the displacements (and rotations). The difficulty is mostly linked to the identification of the natural modes with prevalent in-plane displacement necessary in the expansions.

The angle-ply laminated shell presents three complications. The first, which is common to closed shell structures, is that an inward vibration component with twice the excitation frequency must be retained in order to avoid unphysical excessive in-plane stretching. Here, it must be observed that shells bend much more easily than they stretch. In the case of the circular cylindrical shell, this component is axisymmetric and is associated with relatively high natural frequencies. Results will show that a few of these axisymmetric modes must be retained in the model to guarantee accurate results. The second difficulty is common to axisymmetric structures: it is the 1:1 internal resonance that increases the dimension of the nonlinear problem and complicates significantly the nonlinear response to external forcing. The third difficulty is linked to the angle-ply lamination, which introduces a twist in the natural mode shapes (skewed modes) that is not present for isotropic and cross-ply laminates [30]. In the present model, the refined higher-order shear deformation theory developed by Amabili & Reddy [31] is used to model the shell. A harmonic force excitation is applied at mid-length in the radial direction and simply supported boundary conditions are assumed. The equations of motion are obtained by using an energy approach based on Lagrange equations that retains dissipation. Numerical results are obtained by using the pseudo-arclength continuation method and bifurcation analysis.

## 2. First step: natural modes

A circular cylindrical shell made of a finite number of orthotropic layers, oriented arbitrarily with respect to the shell cylindrical coordinates (*x*,*θ*,*z*), is considered, as shown in figure 1. Here, the theory originally developed by Amabili & Reddy [31] for shells of arbitrary shape is used; the theory is given in appendix A. The displacements of an arbitrary point *P* (*x*,*θ*) on the middle surface of the shell are denoted by *u*, *v* and *w*, in the *x*, *θ* and *z* directions, respectively; *w* is taken to be positive outward from the shell. The displacement of a generic point of the shell at distance *z* from the middle surface is denoted by *u*_{x}, *u*_{θ}, *u*_{z}. The kinematic relationships are presented in appendix A. Initial geometrical imperfections of the shell associated with zero initial stress are denoted by the displacement *w*_{0} in the normal (i.e. radial) direction, also taken to be positive outward. The total thickness of the shell is indicated by *h*. Different techniques could be used in order to obtain the natural modes of the shell, the FEM technique being one of them. Here, the Rayleigh–Ritz method is used because the global discretization makes the mode shapes available as continuous functions expressed in closed form, which is not a necessary requirement (see [20,21,23]) but it makes the procedure easier.

Simply supported boundary conditions are assumed at the shell edges, *x*=0, *L*:
2.1a
2.1b
2.1c
2.1d
and
2.1e
where *N*_{x} is the axial stress resultant per unit length, and *M*_{x} is the axial stress moment resultant per unit length, that is,
2.2
Moreover, the three displacements and the two rotations must be continuous in *θ*.

In order to discretize the problem by using the Rayleigh–Ritz method, the middle surface displacements and rotations *u*, *v*, *w*, *ϕ*_{1} and *ϕ*_{2} are expanded by using admissible functions, which satisfy identically only the geometrical boundary conditions in equation (2.1*a*–*c*). The following base of admissible shell displacements, which allows us to describe skewed modes, is used:
2.3a
2.3b
2.3c
and
2.3d
2.3e
for any number *n* of circumferential waves (*n*=0,1,2,…), where *m* is the number of longitudinal half-waves (*m*=1,2,…), *λ*_{m}=*mπ*/*L*, *ω* is the natural frequency in rad per second, *j* is the imaginary unit and *t* is the time; *u*_{m,n}, *v*_{m,n}, *w*_{m,n}, *ϕ*_{1m,n} and *ϕ*_{2m,n} are the unknown coefficients that have to be determined by using the variational approach; the additional subscript ‘c’ or ‘s’ indicates whether the generalized coordinate is associated with the cosine or sine function in *θ* except for *v* and *ϕ*_{2}, for which the notation is reversed. In the case of angle-ply laminated circular shells, natural modes with nodal lines inclined with respect to the *x*-axis (skewed modes) can appear. In order to describe them, terms with odd and even numbers of longitudinal half-waves have been inserted in equations (2.3), and a combination of sin and cos terms in the variable *θ* appears, but all terms will have the same circumferential wavenumber *n*. The presence of cosine and sine modes is the result of the symmetry of the shell and introduces, for any *n*≠0, couples of modes with exactly the same shape but rotated angularly by *π*/(2*n*). These modes, being identical except for the angular position *θ*, have exactly the same natural frequency for a perfect shell, introducing an internal resonance of the type 1:1. Zero initial geometrical imperfections *w*_{0} of the shell are considered in the radial direction in the numerical simulations.

The kinetic energy *T*_{S} of the shell, including rotary inertia, is given by
2.4
where *ρ*^{(k)}_{S} is the mass density of the *k*th layer of the shell. 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:
2.5
The stress–strain relations for the *k*th orthotropic lamina of the shell, in the material principal coordinates (1,2,3), are obtained under the hypothesis *σ*_{3}=0 [31],
2.6
and
2.7
where *G*_{12}, *G*_{13} and *G*_{23} are the shear moduli in the 1–2, 1–3 and 2–3 directions, respectively, and the superscript (*k*) refers to the *k*th layer within a laminate. Equation (2.6) is obtained (i) under the transverse isotropy assumption with respect to planes orthogonal to axis 1, i.e. assuming fibres in the direction parallel to axis 1, so that *E*_{2}=*E*_{3}, *G*_{12}=*G*_{13} and *ν*_{12}=*ν*_{13}, and (ii) by solving the constitutive equations for *ε*_{33} as a function of *ε*_{11} and *ε*_{22} and then eliminating it.

Equation (2.6) can be transformed to the shell coordinates (*x*,*θ*,*z*) by the following equation [31]:
2.8
where the strain components in the shell coordinates are related to the shell generalized displacement in appendix A by using the Amabili–Reddy nonlinear shell theory; in particular, here all the nonlinear terms in the strain–displacement relationships in appendix A are neglected in the calculation of the natural modes. The matrix [*Q*]^{(k)} is defined in appendix B and is a function of the angle *α* between the shell principal coordinate *x* and the material axis 1.

The elastic strain energy *U*_{S} of the shell is given by
2.9
where *K* is the total number of layers in the laminated shell and (*h*^{(k−1)}, *h*^{(k)}) are the *z* coordinates of the top and bottom surfaces of the *k*th layer.

The vector of the unknown coefficients is introduced for brevity,
2.10
and *n*=0,…,*N*. The dimension of **q** is . The kinetic and potential energies of the shell can be written in vectorial notation
2.11a
and
2.11b
where **M**_{S} and **K**_{S} are the mass and stiffness matrices of the shell, respectively. The natural modes are obtained by solving the following Galerkin equation:
2.12
which gives the natural frequencies *ω* (square root of the eigenvalues) and the natural modes of vibration **q** (eigenvectors).

## 3. Second step: building the reduced-order nonlinear model

A finite number of natural modes, computed by using equation (2.12), are now used to discretize the nonlinear problem. They are defined on the whole structure and allow a reduced-order model to be built with a reasonably small number of degrees of freedom. The small dimension of the reduced-order problem allows continuation technique and bifurcation analysis to be used to obtain the full path of solutions to forced vibration around a resonance. The selection of the natural modes is tricky, and no clear procedure is now available for this operation. In particular, here the selection of the natural modes for the correct simulation of an angle-ply laminated circular cylindrical shell is shown and it is tested versus results obtained by the author with a different discretization technique [25]. This problem is of particular interest because it shows that not only the lowest frequency natural modes have to be retained in the model for accurate simulation, but also axisymmetric terms must be retained for the symmetry of the structure. Moreover, both families of axisymmetric modes with prevalent axial and radial displacements have to be retained. However, for asymmetric modes, only those with prevalent radial displacement have to be retained. Such a non-intuitive conclusion has been reached by developing and integrating many different reduced-order models. They have shown that completely wrong results, not even reported in the present manuscript, are obtained if axisymmetric modes with prevalent axial displacements are neglected.

Assuming harmonic force excitation in the radial direction at the shell mid-length, the model can be built for any resonance of the lowest flexural (i.e. with prevalent radial displacement, which are usually the lowest frequency ones) mode with *n* circumferential waves by using the following natural modes: (i) the first flexural mode **w**_{1,n,c} with *n* circumferential waves (*m*=1, *n*) and its companion **w**_{1,n,s}; (ii) the first three flexural modes with 2*n* circumferential waves and their companions, i.e. **w**_{1,2n,c}, **w**_{2,2n,s}, **w**_{3,2n,c}, **w**_{1,2n,s}, **w**_{2,2n,c}, **w**_{3,2n,s}; and (iii) the axisymmetric modes (only those axisymmetric modes that are symmetric with respect to a plane cutting the shell at half length, denoted, for example, by the subscript 1sym for the first symmetric one) with prevalent radial displacement **w**_{(1sym,0)r}, **w**_{(2sym,0)r}, **w**_{(3sym,0)r}, and with prevalent axial displacement **w**_{(1sym,0)a}, **w**_{(2sym,0)a}, **w**_{(3sym,0)a}. Even if the contribution of the axisymmetric generalized coordinates is small with respect to the global vibration amplitudes, they must be retained in the model in order to have an accurate prediction of the shell nonlinearity. If they are neglected, a wrong strong hardening behaviour is predicted. In fact, these generalized coordinates are associated with the dynamic axisymmetric contraction during large-amplitude vibration and must be retained in the expansion in order to avoid large non-physical in-plane stretching of the shell. In fact, thin shells bend easily but they are very stiff to in-plane stretching. Instead, companion modes have to be retained in the expansion owing to the 1:1 internal resonance. Here, the following notation has been used for natural modes:
3.1
where
3.2a
3.2b
3.2c
3.2d
and
3.2e
where the superscript indicates that the coefficients, introduced in equation (2.3*a*–*e*), have to be taken for the natural mode indicated in the subscript, i.e. for mode (1,*n*,*c*) in this case; the corresponding eigenvector is indicated by **q**^{(1,n,c)}. In equation (3.1), the mode shape has been normalized with respect to the amplitude *w*^{(1,n,c)}_{1,n,c} of the dominant flexural component of mode (1,*n*,*c*). For mode **w**_{(1sym,0)r}, the normalization coefficient is *w*^{(1sym,0)r}_{1sym,0} and so on for other natural modes.

The vector **w** of the middle plane displacements and rotations of the shell is defined as
3.3
and it can be discretized by using the selected natural modes as follows:
3.4
where *g*_{i}(*t*) are the modal coordinates, which are time functions, expressing the degrees of freedom of the discretized reduced-order system. Equation (3.4) gives the recipe to discretize a circular cylindrical shell with uniform boundary conditions and any lamination sequence, even in the case of moderate geometrical imperfections, which have to be introduced in the first step (linear) of the analysis.

Additional modes can be included in the model (3.4), but they will have a very small contribution, whereas those indicated here are important for good accuracy of the nonlinear response. In the case of internal resonance with modes with a different number of circumferential modes, those also have to be included in the model, but this is a special case that can be verified by observing the integer relationship among the natural frequencies and their combinations in the first linear step of the analysis.

The middle surface of the shell is loaded by distributed forces per unit area *q*_{x}, *q*_{θ} and *q*_{z} acting in the *x*, *θ* and *z* directions, respectively. In this study, only a single harmonic force normal to the shell surface is considered; therefore, *q*_{x}=*q*_{θ}=0. The external-distributed load *q*_{z} applied to the shell, owing to the normal concentrated force , is given by
3.5
where *ω* is the excitation frequency, *t* is the time, *δ* is the Dirac delta function, gives the force amplitude positive in the *z*-direction, and and give the position of the point of application of the force. Here, the point excitation is located at the centre of shell, that is, , . Equation (3.5) is correct assuming that the angle 0≤*θ*≤2*π*; actually, equation (3.5) should work for any angle, because the shell is periodic in *θ* with period 2*π*, then equation (3.5) can be replaced with the more complicated expression
This specific dynamic force excites directly the mode **w**_{1,n,c}, which takes the name of the driven mode. On the other hand, the mode **w**_{1,n,s} is not directly excited, because the force is applied at a node for this mode. Therefore, **w**_{1,n,s} is named the companion mode, because it is excited only by nonlinear interaction through a 1:1 internal resonance between **w**_{1,n,c} and **w**_{1,n,s}. The virtual work *W* carried out by the external forces is given by
3.6
The non-conservative damping forces are assumed to be of viscous type and are taken into account by using the Rayleigh dissipation function
3.7
where *c* has a different value for each natural mode; in particular
3.8
In equation (3.8), displacements are non-dimensionalized dividing by *h*, whereas rotations are already non-dimensional. The damping coefficient *c*_{i} is related to the modal damping ratio *ζ*_{i} that can be experimentally evaluated by *ζ*_{i}=*c*_{i}/(2*μ*_{i}*ω*_{i}), where *ω*_{i} is the natural circular frequency of mode *i* and *μ*_{i} is the modal mass of the same mode.

The time-dependent vector of the modal coordinates is defined by
3.9
The generalized forces *Q*_{i} are obtained by differentiation of the Rayleigh dissipation function and of the virtual work carried out by external forces,
3.10
The Lagrange equations of motion are
3.11
where ∂*T*_{S}/∂*g*_{i}=0. The complicated term, derived from the potential energy of the shell, giving quadratic and cubic nonlinearities, can be written in the form
3.12
where the coefficients *γ* have long expressions that also include geometrical imperfections. In equation (3.12), there are quadratic and cubic terms.

For shells with rotary inertia, inertial coupling arises in the equations of motion (see equation (2.5)), so that they cannot be immediately transformed in the first-order form required for numerical integration. In particular, the equations of motion take the following form:
3.13
where **M** is the non-diagonal mass matrix of dimension 14×14 (14 being the number of degrees of freedom), **C** is the damping matrix, **f**_{0} is the vector of excitation amplitudes, **g** is the vector of the modal coordinates, defined in equation (2.6), and
3.14
is the vector containing the stiffness terms (including linear and nonlinear terms) appearing in each one of the 14 equations of motion. Equation (3.13) 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
3.15
Equation (3.15) is in the form suitable for numerical integration.

## 4. Numerical results

A simply supported, imperfection-free, laminated circular cylindrical shell made of graphite/epoxy layers has been investigated: *R*=0.15 m, *L*=0.52 m, *h*=0.003 m, *E*_{1}=50×10^{9} Pa, *E*_{2}=2×10^{9} Pa, *G*_{12}=*G*_{13}=1×10^{9} Pa, *G*_{23}=0.4×10^{9} Pa, *ν*_{12}=*ν*_{23}=0.25 and *ρ*=1500 kg m^{−3}. This is a thin shell, being *R*/*h*=50. The shell is made of two layers 0/30^{°} of the same thickness (*α*=0 for the internal layer); this choice has been made to have significant skewness in the mode shapes. The linear natural frequencies of the shell are given for modes with *m*=1 longitudinal half-waves as a function of the number *n* of circumferential waves in figure 2. The natural frequencies calculated with models with a different number of longitudinal terms *M* in equations (2.3) are given in table 1; the first mode (*m*=1) for a different number of circumferential waves *n* is given. These natural frequencies are for flexural modes (i.e. with prevalent radial displacement *w*) except for the torsional mode for *n*=1. The fundamental mode is obtained for *n*=4 and its natural frequency is 241.4 Hz (with model *M*=4).

The most important mode shapes used in equation (3.4) to discretize the nonlinear vibration problem for the radial mode with *n*=3 circumferential waves are presented in figure 3. It is possible to observe in figure 3*a*, which is for mode (*m*=1, *n*=3) with prevalent radial displacement, that some twisting of the nodal lines parallel to the longitudinal axis is present; this is due to the angle-ply lamination sequence, as discussed in Amabili [25]. Modes with *n*=6 circumferential waves are shown in figure 3*b*,*c* because they appear in equation (3.4). Also axisymmetric modes (*n*=0) are shown in figure 3*d*–*f*, but only those symmetric with respect to the mid-length because they appear in equation (3.4). In figure 3, only modes with prevalent radial displacement are shown.

The nonlinear forced vibration response under harmonic excitation in the spectral neighbourhood of the resonance of mode (*m*=1, *n*=3) with prevalent radial displacement is presented in figure 4. For convenience, a non-dimensionalization of the harmonic force excitation acting in the radial direction is introduced. The non-dimensional force *f* is defined as where *μ* is the modal mass of mode (*m*=1, *n*=3) for the radial displacement *w*, given by *μ*=*ρhRπL*/2. A non-dimensional harmonic point force excitation of amplitude *f*=0.00437 and frequency *ω* applied at shell mid-length is assumed. A modal damping coefficient *ς*_{1,n}=0.001 has been used for all the generalized coordinates having *n*≠0. For the axisymmetric coordinates (*n*=0), a damping proportionally increasing with frequency has been assumed. The nonlinear equations of motion have been integrated by using the software auto [32] for continuation and bifurcation analysis by using the pseudo-arclength continuation method, starting at zero force from the trivial solution. The solution has been initially continued with the excitation amplitude as the parameter and fixed excitation frequency. Once the desired excitation level has been reached, the solution has been continued by using the excitation frequency as the continuation parameter. Pitchfork bifurcations (BPs) and Neimark–Sacker bifurcations are detected in figure 4. Results show the response computed by using the reduced-order model based on natural modes with the 14 degrees of freedom indicated in equation (3.4), compared with results obtained in Amabili [25] with 37 degrees of freedom with an ad hoc constructed reduced-order model. As will be shown later, the present results are more accurate because they take into account a larger number of axial terms (*M*=4), but at the same time they are obtained with a smaller model saving computational time.

The branch 1 of the nonlinear response in figure 4 corresponds to vibration with zero amplitude of the companion mode g_{2}(*t*). Branch 1 has two BPs at *ω*/*ω*_{1,n}=0.964 and at 1.0001, where branch 2 appears. This new branch corresponds to participation of both g_{1}(*t*) and g_{2}(*t*), giving a travelling-wave response, i.e. nodal lines travel around the shell during vibrations. The companion mode presents a node at the location of the excitation force, and therefore it is not directly excited; its amplitude is different from zero only for large-amplitude vibrations, owing to nonlinear coupling through 1:1 internal resonance. In the frequency region where both g_{1}(*t*) and g_{2}(*t*) are different from zero, they give rise to a travelling wave around the shell; phase shift between the two coordinates is practically equal to *π*/2 when the two generalized coordinates have almost the same amplitude, as can be observed comparing figure 4*a*,*b*. This branch appears for sufficiently large excitation. Branch 2 undergoes two Neimark–Sacker (torus) bifurcations, at *ω*/*ω*_{1,n}=0.9641 and 0.9941. The amplitude-modulated (quasi-periodic) response is indicated in figure 4 on branch 2 for 0.9641<*ω*/*ω*_{1,n}<0.9941, that is, bracketed by the two Neimark–Sacker bifurcations.

Another comparison of the present reduced-order model versus two larger models presented in Amabili [25] is shown in figure 5 for the generalized coordinate *g*_{1}(*t*) only. As anticipated, figure 5 shows that the present results are more accurate than those in Amabili [25] for *M*=2 and they seem to be better than those with *M*=4 in Amabili [25] because more terms with 2*n* circumferential waves have been considered in the present reduced-order model. An important reduction in the number of degrees of freedom is obtained, even with respect to the ad hoc reduced-order model presented in Amabili [25].

The additional generalized coordinates used in the present model and showing a significant amplitude are given in figure 6*a*–*d*. Out of those, *g*_{8}(*t*) is the one with the largest displacement and is associated with the third axisymmetric mode (symmetric with respect to the shell mid-length) with prevalent axial displacement. For axisymmetric modes with prevalent radial displacement, *g*_{3}(*t*) is the coordinate with the largest amplitude. In comparison, the contribution of modes with six circumferential waves is small, as indicated by the coordinate *g*_{9}(*t*).

## 5. Conclusion

A reduced-order model for nonlinear vibrations of angle-ply laminated circular cylindrical shells has been developed by using natural modes of vibration. A recipe for the selection of the modes to be retained in the discretization process is given, which is valid for circular shells with uniform boundary conditions and generic lamination sequence even in the case of moderate geometrical imperfections. The model presents very good results compared with an ad hoc reduced model developed to study the same problem, but it gives a further reduction in the number of degrees of freedom. Reduced-order models based on natural modes of vibration do not have the limitation of losing accuracy increasing the vibration amplitude, as is the case in asymptotic nonlinear normal modes. The method gives models of small size that can be studied by using bifurcation analysis and continuation methods. It seems to have a huge potential to be applied in conjunction with commercial finite-element codes for linear analysis.

## Acknowledgements

M.A. acknowledges the financial support of the NSERC Discovery Grant, Canada Research Chair and Canada Foundation for Innovation (Leader Opportunity Fund) programmes of Canada and the PSR-SIIRI programme of Québec.

## Appendix A. The Amabili–Reddy nonlinear higher-order shear deformation theory

The displacements (*u*_{x},*u*_{θ},*u*_{z}) of a generic point of the shell at distance *z* from the middle surface are related to the middle surface displacements by [31]
A1a
A1b
and
A1c
where *ϕ*_{1} and *ϕ*_{2} are the rotations of the transverse normals at *z*=0 about the *y* (*y* is the curvilinear coordinate (*R*+*z*)*θ*, which is orthogonal to *x* and *z*) and *x* axes, respectively, and *ψ*_{1}, *ψ*_{2}, *γ*_{1}, *γ*_{2}, *θ*_{1}, *θ*_{2} and *χ* are functions to be determined in terms of *u*, *v*, *w*, *ϕ*_{1} and *ϕ*_{2}, the five variables describing the shell deformation.

Calculations give
A2a
A2b
A2c
A3a
A3b
A3c
A4
Equations (A.1*a*–*c*) are linear in the five variables describing the shell deformation *u*, *v*, *w*, *ϕ*_{1} and *ϕ*_{2}.

The strain–displacement equations for the higher-order shear deformation theory, keeping terms up to *z*^{3}, are
A5a
A5b
A5c
A5d
and
A5e
where
A6a
A6b
A6c
A6d
A6e
and
A7a
A7b
A7c
A7d
A7e
A7f
A7g
A7h
A7i
and
A8a
A8b
A8c
A8d
A8e
and
A8f
Equations (A.5*d*,*e*) show that the strain and stress distribution through the thickness is parabolic and therefore the shear correction factor is no longer required. Equations (A.6*a*–*c*), giving the middle surface strains, are coincident with those obtained by using Novozhilov nonlinear shell theory [26], which neglects shear deformation and rotary inertia.

## Appendix B. From material to shell reference system

Usually, the lamina material axes (1,2) do not coincide with the shell reference axes (*x*,*θ*), whereas the 3 axis is coincident with *z*. Then, the strains and stresses on material axes can be related to the reference axes by using the following invertible expressions:
B1a
B1b
where
B2
and
B3
*α* being the angle between the shell principal coordinate *x* and the material axis 1.

It can be shown that
B4
Therefore, it is convenient to introduce the matrix [*Q*]^{(k)} given by
B5
where **C** is the 5×5 matrix of *c*_{ij} and *G*_{ij} coefficients in equation (2.6). As a consequence of the discontinuous variation of the stiffness matrix [*Q*]^{(k)} from layer to layer, the stresses may be discontinuous from layer to layer.

## Footnotes

One contribution of 17 to a Theme Issue ‘A celebration of mechanics: from nano to macro’.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.