# Nonlinear vibrations of deep cylindrical shells by the p-version finite element method

Pedro Ribiero, Bruno Cochelin, Sergio Bellizzi

► To cite this version:

Pedro Ribiero, Bruno Cochelin, Sergio Bellizzi. Nonlinear vibrations of deep cylindrical shells by the p-version finite element method. III European Conference on Computational Mechanics, Solids, Structures and Coupled Problems en Engineering, 2006, Lisbonne, Portugal. 10.3233/SAV-2010-0495. [hal-00097572]

HAL Id: hal-00097572

https://hal.science/hal-00097572v1

Submitted on 23 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Non-linear vibrations of deep cylindrical shells by the p-version finite element method

Pedro Ribeiro<sup>a</sup>, Bruno Cochelin<sup>b</sup> and Sergio Bellizzi<sup>b</sup>

<sup>a</sup>IDMEG/DEMEGI, Faculdade de Engenharia, Universidade do Porto, Porto, Portugal
Laboratoire de Mécanique et d'Acoustique, CNRS, Marseille, France

**Abstract.** A *p*-version shell finite element based on the so-called shallow shell theory is for the first time employed to study vibrations of deep cylindrical shells. The finite element formulation for deep shells is presented and the linear natural frequencies of different shells, with various boundary conditions, are computed. These linear natural frequencies are compared with published results and with results obtained using a commercial software finite element package; good agreement is found. External forces are applied and the displacements in the geometrically non-linear regime computed with the *p*-model are found to be close to the ones computed using a commercial FE package. In all numerical tests the *p*-FE model requires far fewer degrees of freedom than the regular FE models. A numerical study on the dynamic behaviour of deep shells is finally carried out.

**Keywords:** Deep shells, *p*-version finite element method, geometrically non-linear vibrations

## 1. Introduction

Shells vibrating with large amplitudes may present quite rich dynamics, which turns geometrically non-linear shell vibrations into a particularly interesting problem. Furthermore, there are several practical applications – notably, but not only, in aeronautics – where shells experience large amplitude vibrations.

The analysis of the vibration of shells may be carried out with great accuracy using the finite element method. However, the *h*-version of this method – which is the more common version and where refined models are obtained by increasing the number of elements over the domain – is constrained by a high computational cost (references [1, 2], for example).

High order elements and especially elements of the *p*-version, where the accuracy of the approximation is improved by raising the number of shape functions within each element, generally require fewer degrees of freedom and demand a smaller computational effort than low order, simple elements [3–7]. This was verified for example in reference [3], where geometrically non-linear static problems on shells were addressed and a small number of higher-order finite elements were sufficient to achieve accurate results. The authors stated as well that the high order elements provided, over ordinary elements, a considerable reduction of the amount of data preparation and interpretation. In [4] a static analysis indicated that the extension of the *p*-version FEM, with hierarchical basis functions, to geometrically nonlinear analysis of laminated rectangular plates is very successful. Additionally, it was shown in [5,6] that the *p*-version is superior to the *h*-version in linear problems, particularly but not only if the problem is without singularities; in the latter case exponential convergence can be achieved.

<footnote>*Corresponding author. Tel.: +351 22 508 17 16; Fax: +351 22 508 14 45; E-mail: pmleal@fe.up.pt.</footnote>
In what concerns the application of the *p*-version elements to linear vibrations, elements of this type were used in reference [6] to compute the linear natural frequencies of open cylindrical shells and the values obtained were similar to the ones from *h*-version finite element models with many more degrees of freedom. In [7] one *p*-version, finite element was proposed to analyse geometrically non-linear free vibrations of shallow shells and free vibrations were studied. This was a first application of the element and, acknowledging that this approximation may be important in many cases, it was assumed that the oscillations are harmonic; it was verified that the mode shapes of shells change with the vibration amplitude. The same element was used together with the shooting method in [8] to study periodic, forced, non-linear vibrations of shallow shells. The number of degrees of freedom of the *p*-version FEM models was, both in [7,8], quite modest.

A shallow shell is a shell whose smallest radius of curvature is larger than its greatest length measured along the middle surface [9–13], so that the raise is small in comparison with the span. In this case, the displacement field is only slightly more complex than the one of plates. In addition to the already mentioned references [7,8], there are a number of studies on geometrically non-linear vibrations of shallow shells. For example, Abe et al. [14] examined the first and second non-linear modes of clamped laminated shells. Two linear modes were used to analyse the second mode; in the analyses of the first non-linear mode it was assumed that the mode shape does not change and is always equal to the linear one. Shear deformation was considered, whilst membrane and rotary inertias were neglected in most examples. The latter approach was justified by the fact that the membrane and rotary inertias are not very important in the calculation of the first linear natural frequency of thin and thick plates. Amabili [15] analysed the geometrically non-linear vibrations of cylindrical shallow shells, simply supported at the four edges and subjected to harmonic excitations. Two non-linear strain-displacement relationships, from different shell theories, were compared. Internal resonances were found and it was concluded that their study may require models with a quite large number of degrees of freedom.

There are a number of works on large deflections of circular cylindrical shells. For example in reference [16] chaotic vibrations of a cylindrical shell are studied. The governing equations of motion are derived on the basis of single and double mode models. The authors conclude that the single mode method may lead to incorrect results. For other studies on circular cylindrical shells, the reader is referred to [17].

Studies of the geometrically non-linear vibrations of open shells that are not shallow (these will be designated as deep shells) are much more uncommon. There are some applications of multi-degree-of-freedom models based on finite elements to study the dynamics of non-shallow shells, but these are generally restrained to short time spans and to transient vibrations. In [18], a finite element method model for geometrically non-linear free vibrations of thin, closed or open, shells was presented. Before carrying out the numerical calculations, the equations of motion were transformed into modal coordinates and the non-diagonal terms of the transformed non-linear stiffness matrices were discarded. Via this approximation, the non-linear equations of motion were decoupled and the study of each non-linear mode amounted to solve a one degree of freedom problem. A fourth-order Runge-Kutta method was then employed and both hardening and softening non-linearities were found. An implicit time-integration scheme that combines algorithmic dissipation of higher modes, conservation of energy and angular momentum was applied to study snap-through of a shell in [19]. In [20] a facet triangular shell element, designated as TRIC, was proposed for linear and non-linear dynamic problems. Consistent and two lumped mass matrices were compared in test examples, and it was concluded that the consistent mass matrix is the most effective. It was argued that the TRIC element and the Newmark integration scheme can achieve converged and accurate solutions with time steps comparable to other shell elements with more sophisticated integration schemes. It was also concluded that, in the examples investigated, implicit integration is much less expensive than explicit.

It is more or less established that *h*-version finite elements can be used to model a curved shell by using a large number of elements [20]. However, it is not evident and has not yet been demonstrated that the *p*-version FEM will be efficient to follow a similar approach, because most (not all) advantages of the *p*-version rest on the fact that it uses a small number of elements: if many elements are required to approach a curved surface these advantages will be lost. To the best of our knowledge, there is no work that shows that assemblies of *p*-version shallow shell finite elements allow analysing non-linear vibration of deep shells, with a reduced number of degrees of freedom.

In the present paper, a *p*-version, finite element with hierarchic basis functions (also called hierarchical finite element method) is used to study vibrations of open, cylindrical shells with any relation between projected length and curvature. The shells studied have rectangular planform and are constituted by isotropic, linear elastic and
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 1. Shell and shell element.</figcaption></figure>

homogeneous materials. Comparisons with established finite element codes are made, showing that the p-version, hierarchical, approach here proposed requires fewer degrees of freedom than the models based on h-version finite elements. To solve the equations of motion in geometrically non-linear vibrations Newmark implicit integration scheme is employed. The influence of the span to height relation and of the thickness of the shells on their dynamics is discussed.

## 2. Finite element model

The p-version element here employed combines a shallow shell theory which may be found, for example, in reference [9], with the first order shear deformation theory, which is usually attributed to Reissner [21], being Mindlin [22] credited with its first application to vibrations. A recent book where this theory is very clearly presented, albeit for laminates, is reference [10]. If the shell is shallow, that is if the its raise is small in comparison with its spans, then the curvilinear coordinates employed in shell theories can be replaced by the Cartesian coordinates x and y, and the Lamé parameters are $A = B = 1$ [9]. Hence, within each element, the displacement components u and v - along the local x and y directions, respectively - are functions of the element's middle surface membrane displacements $u^0$ and $v^0$, and of the rotations of the normal to the middle surface about the local x and y axis. The latter are here denoted by $\theta_x^0$ and $\theta_y^0$ and follow the right-hand rule (Fig. 1). $w^0$ represents the displacement of a particle along the z direction and with respect to the initial curved configuration. The initial curved surface is represented by $w^i$ and for a cylindrical element with curvature radius equal to R it is:

$$w^i(x, y) = -\frac{1}{2} \left(\frac{x^2}{R}\right) \qquad (1)$$

It is assumed that the middle surface displacements $u^0$, $v^0$, $w^0$ and the rotations $\theta_x^0$ and $\theta_y^0$ do not depend on z.

Deep shells will be analysed assembling several elements, but the commonly accepted limit of shallow shells theory, a/R $\le$ 0.5 [13], is respected in each element. This theory of shallow shells is based upon the assumption that the squares and products of $\partial w^i(x, y)/\partial x$ and $\partial w^i(x, y)/\partial y$ are small [11]. In the particular case of cylindrical shells defined by Eq. (1) $\partial w^i(x, y)/\partial y$ is, evidently, zero, so only the restriction $(\partial w^i(x, y)/\partial x)^2 \ll 1$ remains. The strain field is of the following form [10]

$$\varepsilon_x = \varepsilon_x^0 - z\kappa_x \quad \varepsilon_y = \varepsilon_y^0 - z\kappa_y \quad \gamma_{xy} = \gamma_{xy}^0 - z\kappa_{xy} \quad \gamma_{zx} = \gamma_{zx}^0 \quad \gamma_{yz} = \gamma_{yz}^0 \qquad (2)$$

where $\varepsilon_x^0$, $\varepsilon_y^0$ and $\gamma_{xy}^0$ are the membrane strain components at $z = 0$ defined by (Von Kármán relations):

$$\varepsilon_x^0 = u_{,x}^0 + \frac{w^0}{R} + \frac{1}{2} (w_{,x}^0)^2, \quad \varepsilon_y^0 = v_{,y}^0 + \frac{1}{2} (w_{,y}^0)^2, \quad \gamma_{xy}^0 = u_{,y}^0 + v_{,y}^0 + w_{,x}^0 w_{,y}^0 \qquad (3)$$

The comma represents partial derivation and $\kappa_x$, $\kappa_y$ and $\kappa_{xy}$ are the curvature changes, which are given by:
$$ \kappa_x = -\frac{\partial \theta_y^0}{\partial x}, \quad \kappa_y = \frac{\partial \theta_x^0}{\partial y}, \quad \kappa_{xy} = -\frac{\partial \theta_y^0}{\partial y} + \frac{\partial \theta_x^0}{\partial x} \qquad (4) $$

The transverse shear strains are considered to be constant along the element thickness and to be given by

$$ \gamma_{zx}^0 = w_{,x}^0 + \theta_y^0, \quad \gamma_{yz}^0 = w_{,y}^0 - \theta_x^0 \qquad (5) $$

The term $u^0/R$ which could appear in $\gamma_{zx}^0$ will not be here considered, to be consistent with the assumption that $R$ is large [10] and because $u^0$ will be small in the applications.

In each element, the middle surface displacements are expressed in the form:

$$ \left\{ \begin{matrix} u^0(x, y, t) \\ v^0(x, y, t) \\ w^0(x, y, t) \\ \theta_x^0(x, y, t) \\ \theta_y^0(x, y, t) \\ \theta_z^0(x, y, t) \end{matrix} \right\} = \left[ \begin{matrix} \mathbf{N}^u(x,y)^T & 0 & 0 & 0 & 0 & 0 \\ 0 & \mathbf{N}^v(x,y)^T & 0 & 0 & 0 & 0 \\ 0 & 0 & \mathbf{N}^w(x,y)^T & 0 & 0 & 0 \\ 0 & 0 & 0 & \mathbf{N}^{\theta_x}(x,y)^T & 0 & 0 \\ 0 & 0 & 0 & 0 & \mathbf{N}^{\theta_y}(x,y)^T & 0 \\ 0 & 0 & 0 & 0 & 0 & \mathbf{N}^{\theta_z}(x,y)^T \end{matrix} \right] \left\{ \begin{matrix} \mathbf{q}_u(t) \\ \mathbf{q}_v(t) \\ \mathbf{q}_w(t) \\ \mathbf{q}_{\theta_x}(t) \\ \mathbf{q}_{\theta_y}(t) \\ \mathbf{q}_{\theta_z}(t) \end{matrix} \right\} \quad (6) $$

where $\mathbf{q}_u(t)$, $\mathbf{q}_v(t)$ are the vectors of generalised membrane displacements, $\mathbf{q}_w(t)$ are the vectors of transverse displacements, and $\mathbf{q}_{\theta_x}(t)$, $\mathbf{q}_{\theta_y}(t)$ and $\mathbf{q}_{\theta_z}(t)$ are the vectors of generalised rotations.

The rotation about $z$, $\theta_z$, was introduced to facilitate the transformation from local to global coordinates [23], however, it is not present in the strain displacement relations or in the acceleration terms and its consideration would result in the appearance of zero blocks in the stiffness and mass matrices. In order to avoid singularities in these matrices, artificial stiffness and mass sub-matrices connected with $\theta_z$ are introduced. These artificial matrices are computed by multiplying the sub-matrices associated with the rotation about $\theta_y$ by parameters, the influence of which is discussed in the applications. Either than that, the rotations about $z$ are ignored in the formulation.

The matrix of shape functions

$$ [N] = \begin{bmatrix} \mathbf{N}^{u^T} & 0 & 0 & 0 & 0 \\ 0 & \mathbf{N}^{v^T} & 0 & 0 & 0 \\ 0 & 0 & \mathbf{N}^{w^T} & 0 & 0 \\ 0 & 0 & 0 & \mathbf{N}^{\theta_x^T} & 0 \\ 0 & 0 & 0 & 0 & \mathbf{N}^{\theta_y^T} \end{bmatrix} \qquad (7) $$

is constituted by the row vectors of bi-dimensional surface or membrane ($\mathbf{N}^{u^T}$ and $\mathbf{N}^{v^T}$), transverse ($\mathbf{N}^{w^T}$) and rotational ($\mathbf{N}^{\theta_x^T}$, $\mathbf{N}^{\theta_y^T}$) shape functions. To facilitate the elaboration of the code and the assembling process the same set and the same number of shape functions will be employed for generalised displacements related with $u$, $v$ and $w$, i.e., $\mathbf{N}^{u^T} = \mathbf{N}^{v^T} = \mathbf{N}^{w^T}$. The sets and number of shape functions related with the rotations $\theta_x$ and $\theta_y$ are also equal, $\mathbf{N}^{\theta_x^T} = \mathbf{N}^{\theta_y^T}$, but the shape functions used for the rotations may be different from the ones employed for the displacements.

The vectors of shape functions $\mathbf{N}^{u^T} = \mathbf{N}^{v^T} = \mathbf{N}^{w^T}$ are defined using products of functions from a set that will be designated as the $f$ set of functions. This is constituted by the four Hermite cubics and by the polynomials given by the following formula

$$ f_{r-4}(\xi) = \sum_{n=0}^{INT(r/2)} \frac{(-1)^n (2r - 2n - 7)!!}{2^n n! (r - 2n - 1)!} \xi^{r-2n-1}, \quad r > 4 \qquad (8) $$

where $\xi$ is an adimensional coordinate. The vectors of shape functions $\mathbf{N}^{\theta_x^T} = \mathbf{N}^{\theta_y^T}$ are either defined using the $f$ set of functions or, alternatively, another "g" set of functions constituted by two linear functions and by

$$ g_{r-2}(\xi) = \sum_{n=0}^{INT(r/2)} \frac{(-1)^n (2r - 2n - 5)!!}{2^n n! (r - 2n - 1)!} \xi^{r-2n-1}, \quad r > 2 \qquad (9) $$

More details concerning the shape functions are given in reference [7].
The constitutive equations for linear elastic isotropic materials are used, with E representing the Young modulus and ν the Poisson ratio. The value λ = 5/6 [24] is employed for the shear correction factor.

The direct and membrane shear strain components are

$$
\begin{equation}
\begin{bmatrix} \epsilon_x \\ \epsilon_y \\ \gamma_{xy} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & -z & 0 & 0 \\ 0 & 1 & 0 & 0 & -z & 0 \\ 0 & 0 & 1 & 0 & 0 & -z \end{bmatrix} \boldsymbol{\epsilon} = \begin{bmatrix} \mathbf{I} & -z\mathbf{I} \end{bmatrix} \boldsymbol{\epsilon} \tag{10}
\end{equation}
$$

with

$$
\epsilon = \begin{Bmatrix} \epsilon_o^p \\ \epsilon_o^b \end{Bmatrix} + \begin{Bmatrix} \epsilon_L^p \\ 0 \end{Bmatrix} \qquad (11)
$$

The linear membrane and bending strains, $\epsilon_0^p$ and $\epsilon_0^b$, and the geometrically non-linear membrane strain, $\epsilon_L^p$, may be
written in vector form as

$$
\epsilon_0^p = \begin{bmatrix} \mathbf{N}_{,x}^{u^T} & 0 \\ 0 & \mathbf{N}_{,y}^{u^T} \\ \mathbf{N}_{,y}^{u^T} & \mathbf{N}_{,x}^{u^T} \end{bmatrix} \left\{ \begin{matrix} \mathbf{q}_u \\ \mathbf{q}_v \end{matrix} \right\} + \left\{ \begin{matrix} \frac{1}{R} \mathbf{N}^{w^T} \mathbf{q}_w \\ 0 \\ 0 \end{matrix} \right\}, \quad \epsilon_o^b = \begin{bmatrix} -\mathbf{N}_{,x}^{\theta^T} & 0 \\ 0 & \mathbf{N}_{,y}^{\theta^T} \\ -\mathbf{N}_{,y}^{\theta^T} & \mathbf{N}_{,x}^{\theta^T} \end{bmatrix} \left\{ \begin{matrix} \mathbf{q}_{\theta_y} \\ \mathbf{q}_{\theta_x} \end{matrix} \right\} \tag{12}
$$

$$
\epsilon_L^p = \begin{cases} \frac{1}{2} \mathbf{q}_w^T \mathbf{N}_{,x}^w \mathbf{N}_{,x}^{w^T} \mathbf{q}_w \\ \frac{1}{2} \mathbf{q}_w^T \mathbf{N}_{,y}^w \mathbf{N}_{,y}^{w^T} \mathbf{q}_w \\ \mathbf{q}_w^T \mathbf{N}_{,x}^w \mathbf{N}_{,y}^{w^T} \mathbf{q}_w \end{cases}
$$

The transverse shear strains are given by

$$
\left\{ \begin{matrix} \gamma_{zx} \\ \gamma_{yz} \end{matrix} \right\} = \begin{bmatrix} \mathbf{N}_{,x}^{w^T} & 0 & \mathbf{N}_{,y}^{\theta^T} \\ 0 & -\mathbf{N}_{,x}^{\theta^T} & 0 \end{bmatrix} \left\{ \begin{matrix} \mathbf{q}_w \\ \mathbf{q}_{\theta_x} \\ \mathbf{q}_{\theta_y} \end{matrix} \right\} \qquad (13)
$$

For an element of uniform thickness h, the membrane stress resultants {$T_x, T_y, T_{xy}$} and couples {$M_x, M_y, M_{xy}$},
all per unit length, are defined by

$$
\{T_x, T_y, T_{xy}\} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \{\sigma_x, \sigma_y, \tau_{xy}\} dz \qquad (14)
$$

$$
\{M_x, M_y, M_{xy}\} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \{\sigma_x, \sigma_y, \tau_{xy}\} z dz \qquad (15)
$$

The shear stress resultants per unit length are

$$
\{Q_x, Q_y\} = \int_{-\frac{h}{2}}^{\frac{h}{2}} \{\tau_{xz}, \tau_{yz}\} dz \qquad (16)
$$

Using Hooke's generalised law and Eqs (14), (15) and (16), one arrives at the constitutive relations of the shallow
shell finite element:

$$
\begin{Bmatrix} \mathbf{T} \\ \mathbf{M} \end{Bmatrix} = \begin{bmatrix} \mathbf{A} & 0 \\ 0 & \mathbf{D} \end{bmatrix} \boldsymbol{\epsilon} \qquad (17)
$$

$$
\left\{ \begin{matrix} \mathbf{Q}_x \\ \mathbf{Q}_y \end{matrix} \right\} = [\mathbf{C}] \gamma \qquad (18)
$$

Matrices **A**, **D** and **C** may be found in [7].

The element matrices are derived by the principle of virtual work that may be written as: δW<sub>in</sub> = δW<sub>ex</sub> + δW<sub>j</sub>,
i.e., the virtual work of the elastic restoring forces, or internal forces (δW<sub>in</sub>), should be equal to the sum of the virtual
works of the external ($\delta W_{ex}$) and inertia forces ($\delta W_j$). Inertia forces are defined using the well known principle of d'Alembert.

The virtual work of the internal forces is

$$ \delta W_{in} = \int_{\Omega} \left( \left\{ \delta\varepsilon_{0}^{p} \right\}^{T} + \left\{ \delta\varepsilon_{L}^{p} \right\}^{T} \right) \begin{bmatrix} \mathbf{A} & 0 \\ 0 & \mathbf{D} \end{bmatrix} \left( \left\{ \varepsilon_{0}^{p} \right\} + \left\{ \varepsilon_{L}^{p} \right\} \right) d\Omega + \int_{\Omega} \delta\gamma^{T} C \gamma d\Omega \quad (19) $$

From the products involving linear strains and their variation, the constant stiffness matrices are defined. Because these constant matrices originate linear terms in the equations of motion, they are commonly designated as linear stiffness matrices. The membrane and the bending linear stiffness matrices are represented by $K_1^p$ and $K_1^b$; the matrices due to coupling between the transverse response and the membrane displacements are $K_1^{sp}$ and $K_1^{ps}$; the shear stiffness matrix is $K_1^\gamma$. These matrices are given in [7]. As explained before, an artificial matrix connected with $\theta_z$ was introduced, by means of a parameter $k_{\theta z}$. Thus, the complete element linear stiffness matrix has the following form

$$ K_{el} = \begin{bmatrix} K_1^{p11} & K_1^{p12} & K_1^{pS13} & 0 & 0 & 0 \\ K_1^{p21} & K_1^{p22} & K_1^{pS23} & 0 & 0 & 0 \\ K_1^{SP31} & K_1^{SP32} & K_1^{ss33} + K_1^{\gamma33} & K_1^{\gamma34} & K_1^{\gamma35} & 0 \\ 0 & 0 & K_1^{\gamma43} & K_1^{b44} + K_1^{\gamma44} & K_1^{b45} & 0 \\ 0 & 0 & K_1^{\gamma53} & K_1^{b54} & K_1^{b55} + K_1^{\gamma55} & 0 \\ 0 & 0 & 0 & 0 & 0 & k_{\theta z} (K_1^{b55} + K_1^{\gamma55}) \end{bmatrix} \quad (20) $$

The dimension of this square elemental matrix is $p_{el} = 3p_o^2 + 3p_b^2$; the superscript numbers indicate the position of each sub-matrix in Eq. (20).

The matrices that give rise to quadratic non-linear terms in the equations of motion, matrices $K_2$ and $K_3$, result from terms containing either $\delta\varepsilon_L^p$ or $\varepsilon_L^p$, and the matrix that originates cubic non-linear terms, matrix $K_4$, results from terms that contain products of $\delta\varepsilon_L^p$ by $\varepsilon_L^p$ [7].

The virtual work of the inertia forces is:

$$ \delta W_j = - \int_{-\frac{h}{2}}^{\frac{h}{2}} \int_{\Omega} \rho (\delta u \ddot{u} + \delta v \ddot{v} + \delta w \ddot{w}) d\Omega dz \quad (21) $$

where $\rho$ represents the mass per unit volume and $\Omega$ represents the finite element area. It is noted that rotary inertia is included in expression Eq. (21) because $u = u^0 + z\theta_y^0$ and $v = v^0 - z\theta_x^0$. Hence, the consistent mass matrix of each element would be given by

$$ M = \rho h \begin{bmatrix} \int_{\Omega} u u^{T} d\Omega & 0 & 0 & 0 & 0 & 0 \\ 0 & \int_{\Omega} u u^{T} d\Omega & 0 & 0 & 0 & 0 \\ 0 & 0 & \int_{\Omega} w w^{T} d\Omega & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{h^2}{12} \int_{\Omega} \theta_y \theta_y^{T} d\Omega & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{h^2}{12} \int_{\Omega} \theta_x \theta_x^{T} d\Omega \end{bmatrix} \quad (22) $$

However, an artificial mass matrix for the rotations about $z$ was introduced and the elemental mass matrix has the following form

$$ M_{el} = \begin{bmatrix} M_p^{11} & 0 & 0 & 0 & 0 & 0 \\ 0 & M_p^{22} & 0 & 0 & 0 & 0 \\ 0 & 0 & M_b^{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & M_{Rx}^{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & M_{Ry}^{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & m_{\theta z} M_{Ry}^{55} \end{bmatrix} \quad (23) $$

$M_p$ and $M_b$ are the membrane and transverse inertia matrices, and $M_{Ry}$ and $M_{Rx}$ are due to the rotary inertia. Common simplifications of the mass matrix, like neglecting the rotary or the membrane inertia, or resorting to mass lumping, will not be adopted here.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Sketch of angles between tangents to the element and the global axis, represented by $x_G$ and $z_G$.</figcaption></figure>

All the above matrices are derived in the local adimensional coordinates $\xi$ and $\eta$. Before assembling them, the matrices are transformed into global coordinates by introducing the physical dimensions of each element and “rotating” the matrices. For example the element global stiffness matrix is computed as:

$$ K_{elGlobal} = T_{S_1S_0} K_{el} T_{S_0S_1} \qquad (24) $$

where $T_{S_1S_0}$ transforms local to global coordinates. In the case represented in Fig. 2, $T_{S_1S_0}$ is defined as:

$$ T_{S_1S_0} = \begin{bmatrix} \text{COS} & 0 & -\text{SIN} \\ 0 & \text{I} & 0 \\ \text{SIN} & 0 & \text{COS} \end{bmatrix} \qquad (25) $$

where matrix COS and SIN have the following form

$$ \text{COS}_{p_o^2 \times p_o^2} = \begin{bmatrix} \cos \phi_{i_{p_o \times p_o}} & 0 & 0 \\ 0 & \cos \theta_{i_{(p_o^2 - 2 \times p_o) \times (p_o^2 - 2 \times p_o)}} & 0 \\ 0 & 0 & \cos \phi_{i+1_{p_o \times p_o}} \end{bmatrix} \qquad (26) $$

$$ \text{SIN}_{p_o^2 \times p_o^2} = \begin{bmatrix} \sin \phi_{i_{p_o \times p_o}} & 0 & 0 \\ 0 & \sin \theta_{i_{(p_o^2 - 2 \times p_o) \times (p_o^2 - 2 \times p_o)}} & 0 \\ 0 & 0 & \sin \phi_{i+1_{p_o \times p_o}} \end{bmatrix} \qquad (27) $$

Each sub-matrix in Eqs (26) and (27) is a diagonal matrix. The symbol $\phi_i$ indicates the angle of the tangent to the i side (Fig. 2), and affects shape functions that contain only Hermite cubics or products of Hermite cubics by Legendre polynomials. The symbol $\theta_i$ represents the angle of the shell element planform and affects functions that contain only Legendre polynomials. In this way shallow shell elements can be used to approximate a deep shell.

The generalised coordinates internal to each element, which are due to the p enrichment, do not couple with any other coordinates, since these internal coordinates multiply by displacement shape functions that disappear along the element boundaries and nodes. Continuity between elements is achieved by matching the generalised coordinates that affect not only common nodes, i.e., the coordinates that are linked with products of Hermite cubics, but also common edges, which are associated with products of Hermite cubics by Legendre polynomials (example shown in Fig. 3). Thus, the larger the number of shape functions employed, the higher the degree of deformation that may occur not only inside each element but also at the elements boundaries. Continuity, as understood in the finite element method, is of type $C_0$, that is, it is only enforced on the generalised displacements, which include the independent rotations.

The elemental vector of external forces in local coordinates is also derived applying the principle of virtual work [8], this vector is then transformed into global coordinates and assembled to the vectors of other elements.
<figure><img src="image_5.png" /><figcaption>Fig. 3. Example of shape function defined by a product of a Hermite cubic by a higher order polynomial.</figcaption></figure>

After assemblage, the time domain equations of motion in generalised global coordinates are obtained; they are of the following form:

$$ M\ddot{q}_G + C\dot{q}_G + Klq_G + Knl(q_G)q_G = P \quad (28) $$

In the former equation Kl and Knl(q_G) stand for the constant stiffness matrix and for the stiffness matrix that originates non-linear terms, respectively; M and C represent the assembled mass and damping matrices. The damping matrix here used is proportional to the linear stiffness [24]. The vectors P and q_G are the vectors of external forces and generalised displacements in global coordinates.

The non-linear equations of motion are solved in the time domain by an implicit Newmark method [24,25], which allows one to correct the non-linear stiffness matrix until the equation

$$ (a_0 M + a_1 C + Kl + Knl(q_{G_{t_i+\Delta t}})) q_{G_{t_i+\Delta t}} = P_{t_i+\Delta t} + M (a_0 q_{G_{t_i}} + a_2 \dot{q}_{G_{t_i}} + a_3 \ddot{q}_{G_{t_i}}) + C (a_1 q_{G_{t_i}} + a_4 \dot{q}_{G_{t_i}} + a_5 \ddot{q}_{G_{t_i}}) \quad (29) $$

is satisfied below a desired error condition. The Newmark parameters δ = 0.5, γ = 0.25 were used and the constants a₀–a₅ are given in [25].

In each time step, an inner cycle is necessary to update the non-linear terms until convergence is achieved. This requires transforming back from global to local coordinates, generating the non-linear matrices of each element in local coordinates, returning to global coordinates and re-assembling the matrix. Naturally, this process is time consuming, but it benefits from the fact that the p-element model is not too large. Moreover, the time spent in the assembly stage is reduced because the number of elements is much smaller in the p-version than in the h-version and several coordinates are internal to the elements in the p-version. Finally, the fact that the p-version models generally have less degrees of freedom than h-version models, should mean that less computer memory and less computational time is required to solve the system of equation that arises in each step of the implicit iterative procedure.

# 3. Numerical applications

## 3.1. Linear natural frequencies

In order to carry out a convergence study and to verify the accuracy of the approach presented, the linear natural frequencies of vibration of different shells will be now computed and compared with published values. The dimension and material properties of the shells considered in most of this section and in the study of non-linear vibrations that follows are the following: a = b = 0.3 m, E = 7 10¹⁰ N/m², ρ = 2778 kg/m³, ν = 0.3. In order to compare with reference [26], other material parameters are employed in one particular occasion. As represented in Fig. 1, a and b are the width and length of the projection of the shell on a plane. E represents the Young modulus, ρ the mass per
<table><caption>Table 1
First four frequency parameters of a completely free cylindrical shell ($a/h = 100, a/b = 1, a/R = 0.5$)</caption><thead><tr><th rowspan="2">Number of<br>elements</th><th rowspan="2">Number of shape functions (p-version)<br>or reference</th><th colspan="4">Mode number</th><th rowspan="2">DOF</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>One</td><td>p<sub>o</sub> = 8 f functions; p<sub>θ</sub> = 11 g functions</td><td>13.487</td><td>22.065</td><td>34.804</td><td>48.643</td><td>434</td></tr><tr><td>One</td><td>p<sub>o</sub> = 13 (f); p<sub>θ</sub> = 16(g)</td><td>13.466</td><td>22.063</td><td>34.767</td><td>48.596</td><td>1019</td></tr><tr><td>Four</td><td>p<sub>o</sub> = 5 (f); p<sub>θ</sub> = 5 (g)</td><td>13.410</td><td>21.650</td><td>34.636</td><td>49.425</td><td>510</td></tr><tr><td>Four</td><td>p<sub>o</sub> = 6 (f); p<sub>θ</sub> = 6 (g)</td><td>13.358</td><td>21.645</td><td>34.342</td><td>49.301</td><td>756</td></tr><tr><td>Four</td><td>p<sub>o</sub> = 8 (f); p<sub>θ</sub> = 10 (g)</td><td>13.339</td><td>21.628</td><td>34.284</td><td>48.763</td><td>1806</td></tr><tr><td>Four</td><td>p<sub>o</sub> = 8, p<sub>θ</sub> = 10 (f)</td><td>13.339</td><td>21.628</td><td>34.284</td><td>48.763</td><td>1806</td></tr><tr><td>Four</td><td>p<sub>o</sub> = 16, p<sub>θ</sub> = 18 (f)</td><td>13.321</td><td>21.626</td><td>34.253</td><td>48.720</td><td>6654</td></tr><tr><td>100 × 100</td><td>Ref. [6] Ansys</td><td>13.403</td><td>21.473</td><td>34.148</td><td>48.913</td><td>61 206</td></tr><tr><td>1</td><td>Ref. [6] HFEM, p<sub>o</sub> = p<sub>i</sub> = 10</td><td>13.403</td><td>21.473</td><td>34.147</td><td>48.908</td><td>300</td></tr><tr><td>—</td><td>Ref. [13]</td><td>13.508</td><td>22.073</td><td>34.868</td><td>48.703</td><td>75</td></tr></tbody></table>

<table><caption>Table 2
First five frequency parameters of a FCFC cylindrical shell ($a/h = 100, a/b = 1, a/R = 0.5$)</caption><thead><tr><th rowspan="2">Number of shape functions<br>or reference</th><th colspan="5">Mode number</th><th rowspan="2">DOF</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>p<sub>o</sub> = p<sub>θ</sub> = 10 (f functions)</td><td>36.901</td><td>39.679</td><td>82.203</td><td>84.095</td><td>84.793</td><td>2100</td></tr><tr><td>p<sub>o</sub> = p<sub>θ</sub> = 14 (f functions)</td><td>36.879</td><td>39.645</td><td>82.184</td><td>84.030</td><td>84.734</td><td>4284</td></tr><tr><td>Reference [27]</td><td>36.952</td><td>39.745</td><td>82.244</td><td>84.975</td><td>85.026</td><td>—</td></tr></tbody></table>

<table><caption>Table 3
First five frequencies (in Hz) of a CFCF cylindrical shell $a = 0.25$ m, $b = 0.4$ m, $h = 0.001$ m</caption><thead><tr><th rowspan="2">Number of<br>elements</th><th rowspan="2">number of shape<br>functions or ref.</th><th colspan="5">Mode number</th><th rowspan="2">DOF</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Four</td><td>p<sub>o</sub> = p<sub>θ</sub> = 6</td><td>216.55</td><td>397.88</td><td>407.19</td><td>562.32</td><td>715.91</td><td>684</td></tr><tr><td>Six</td><td>p<sub>o</sub> = p<sub>θ</sub> = 6</td><td>215.84</td><td>396.59</td><td>405.89</td><td>561.09</td><td>713.28</td><td>1044</td></tr><tr><td>h-version FEM, Abaqus</td><td>Ref. [26]</td><td>215.09</td><td>398.75</td><td>406.43</td><td>562.66</td><td>721.73</td><td>16 144</td></tr></tbody></table>

unit volume and $v$ is the Poisson ratio. The letter $h$ represents the thickness, which is often $h = 0.01$ a, but different values will be assumed occasionally.

Table 1 gives the natural frequency parameter – defined as $\Omega = \omega b^2 \sqrt{\rho h/D}$, where $D$ is the flexural rigidity equal to $E h^3 / (12(1 - \nu^2))$ – of a shell with free boundaries and where $a/R = 0.5$. The values are compared with results from Bardell et al. [6], who used both a commercial finite element software and a thin shell p-version element, and from Leissa and Narita (Table 1 of reference [13]), who employed the Rayleigh-Ritz method and thin shell theory. The lower six eigenvalues computed with the present approach are not shown in the table, but they are very close to zero, value they should have in the case of free boundaries. Our results agree with the ones published in the former references. It is worth noting that exactly the same values are calculated using $f$ and $g$ functions for the rotations. In the computation of linear natural frequencies, the parameter $k_{\theta z}$ was generally made equal to $10^{-7}$ and $m_{\theta z}$ equal to $10^{-2}$. However, it was verified that a change in these values results in minor alterations in the natural frequencies. The shell is shallow and one element provides rather accurate results with a small number degrees of freedom. Nevertheless, the good results computed with four elements provide a first validation of the accuracy of the assembling procedure implemented.<fnref n="1" />

Table 2 gives the first natural frequencies of a cylindrical shell, also with $a/R = 0.5$, but with the two straight edges ($x = -a/2$ and $x = a/2$) free (F) and the two curved edges ($y = -b/2$ and $y = b/2$) clamped (C). Four elements were used in the two p-version models employed, which differ in the number of functions used per element. The results are quite close to the ones of reference [27].

Table 3 shows values calculated using the p-version model and the FEM package Abaqus [26]; in this case the straight edges ($x = -a/2$ and $x = a/2$) are clamped and the curved edges ($y = -b/2$ and $y = b/2$) are free. The

<footnote n="1">¹The occasional variation in the number of p-elements would allow us to designate the finite element method here employed as an hp-finite element method.</footnote>
<table><caption>Table 4
First five frequency parameters of a CFCF cylindrical shell a/b = 1, a/R = 1.25, a/h = 100</caption><thead><tr><th rowspan="2">Number of<br>elements</th><th rowspan="2">Number of shape functions<br>for ref.</th><th colspan="5">Mode number</th><th rowspan="2">DOF</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Four</td><td>p<sub>o</sub> = p<sub>θ</sub> = 5 (f functions)</td><td>50.477</td><td>93.667</td><td>99.228</td><td>126.77</td><td>169.43</td><td>450</td></tr><tr><td>Four</td><td>p<sub>o</sub> = 6, (f functions), p<sub>θ</sub> = 6 (g)</td><td>50.357</td><td>93.283</td><td>98.477</td><td>125.84</td><td>168.24</td><td>684</td></tr><tr><td>Four</td><td>p<sub>o</sub> = p<sub>θ</sub> = 6 (f functions)</td><td>50.449</td><td>93.452</td><td>98.696</td><td>126.07</td><td>168.61</td><td>684</td></tr><tr><td>Eight</td><td>p<sub>o</sub> = p<sub>θ</sub> = 7 (f)</td><td>49.984</td><td>92.627</td><td>98.012</td><td>125.65</td><td>167.05</td><td>1974</td></tr><tr><td>Eight</td><td>p<sub>o</sub> = 8 (f), p<sub>θ</sub> = 8 (g)</td><td>49.982</td><td>92.623</td><td>97.820</td><td>125.47</td><td>167.04</td><td>2640</td></tr><tr><td>Nine</td><td>p<sub>o</sub> = 7 (f), p<sub>θ</sub> = 7 (g)</td><td>49.958</td><td>92.581</td><td>97.966</td><td>125.62</td><td>166.96</td><td>2226</td></tr><tr><td>Ten</td><td>p<sub>o</sub> = 7 (f), p<sub>θ</sub> = 7 (g)</td><td>49.939</td><td>92.548</td><td>97.932</td><td>125.60</td><td>166.90</td><td>2478</td></tr><tr><td>Mesh 30 × 30</td><td>h-version FEM</td><td>48.759</td><td>91.162</td><td>96.175</td><td>124.24</td><td>165.33</td><td>—</td></tr></tbody></table>

<table><caption>Table 5
First five frequency parameters of a CFCF cylindrical shell a/b = 1, a/R = 1.5, a/h = 100</caption><thead><tr><th rowspan="2">Number of<br>elements</th><th rowspan="2">number of shape functions<br>or ref.</th><th colspan="5">Mode number</th><th rowspan="2">DOF</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Four</td><td>p<sub>o</sub> = p<sub>θ</sub> = 6 (f functions)</td><td>45.067</td><td>85.313</td><td>89.978</td><td>126.32</td><td>153.42</td><td>684</td></tr><tr><td>Four</td><td>p<sub>o</sub> = p<sub>θ</sub> = 12 (f functions)</td><td>44.958</td><td>85.112</td><td>89.700</td><td>125.97</td><td>153.00</td><td>3096</td></tr><tr><td>Four</td><td>p<sub>o</sub> = 12 (f) p<sub>θ</sub> = 12 (g)</td><td>44.958</td><td>85.112</td><td>89.700</td><td>125.97</td><td>153.00</td><td>3096</td></tr><tr><td>Six</td><td>p<sub>o</sub> = 8 (f), p<sub>θ</sub> = 8 (g)</td><td>44.578</td><td>84.432</td><td>89.037</td><td>125.34</td><td>151.79</td><td>1968</td></tr><tr><td>Eight</td><td>p<sub>o</sub> = 8 (g), p<sub>θ</sub> = 8 (f)</td><td>44.438</td><td>84.182</td><td>88.783</td><td>125.10</td><td>151.34</td><td>2640</td></tr><tr><td>Ten</td><td>p<sub>o</sub> = 7 (f), p<sub>θ</sub> = 7 (g)</td><td>44.377</td><td>84.072</td><td>88.861</td><td>125.22</td><td>151.14</td><td>2478</td></tr><tr><td>Twelve</td><td>p<sub>o</sub> = 7 (f), p<sub>θ</sub> = 7 (g)</td><td>44.342</td><td>84.009</td><td>88.797</td><td>125.16</td><td>151.02</td><td>2982</td></tr><tr><td>Twenty</td><td>p<sub>o</sub> = 7 (f), p<sub>θ</sub> = 7 (g)</td><td>44.291</td><td>83.919</td><td>88.704</td><td>125.07</td><td>150.86</td><td>5040</td></tr><tr><td>Mesh 30 × 30</td><td>h-version FEM</td><td>42.737</td><td>82.086</td><td>86.436</td><td>122.98</td><td>148.88</td><td>—</td></tr></tbody></table>

geometric properties are the following: a = 0.25 m, b = 0.4 m, h = 0.001 m. The curvature radius is equal to the projected length, R = a, hence this shell is not shallow. The material properties are still typical of aluminium, but they are $E = 7 \times 10^{10}$ N/m$^2$, $\rho = 2657$ kg/m$^3$ and $v = 0.3$ in order to make a direct comparison with results from [26]. As seen in Table 3, the values of the natural frequencies computed using the different models are very similar and the dimension of the p-version models is far smaller than the one of the h-version model.

The Abaqus shell element employed in reference [26] is element S8R5, which was chosen after a convergence study where it required less degrees of freedom than other Abaqus elements for shells. Element S8R5 is a thin shell element [28] and therefore should not be used when the transverse shear deformation is important. This element has eight nodes at the boundary and an internally generated midbody node, all nodes with five degrees of freedom, three displacements and two rotations. According to the Abaqus manual [28], reduced integration is employed, whilst in the p-version finite element here used complete integration is carried out. Although damping is only introduced in the next section, we also report here another difference between Abaqus and the p-model. We will define the stiffness proportional damping matrix as $C = \alpha KI$; Abaqus uses the same definition in linear models but not exactly the same when the system is non-linear (reference [28]). These differences between the p-element and element S8R5 will lead to close but not equal results.

Returning to shells with straight edges clamped and curved edges free, but now with $a/R = 1.25$ and $a/R = 1.5$, Tables 4 and 5 show some of the results of convergence studies that were carried out and comparisons with results calculated using the commercial FEM package Abaqus. Note that the length a of the line segment with the direction of the x axis, that results from the projection of each shell on the plane Oxy, is the same for all shells. The Abaqus h-version finite element used was again the S8R5 and the mesh is constituted by 30 per 30 elements. Concerning the agreement between different models, the case $a/R = 1.5$ is inferior to the other numerical tests: nevertheless the difference between the p-version and the Abaqus frequencies attains only the maximum value of 3.6%.

The frequency parameters computed by the p-version model change more with the number of elements as the shells become deeper. This occurs because the larger the number of elements employed, the better the geometry is approached. In any case, four elements provide rather reasonable approximations for the shells here studied. No difference resulted from using g instead of f functions for the rotations.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Section of shell, $y = 0$, and points (1, 2 and 3) where displacements and velocities are computed.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Displacements of point 2 in the z direction computed with the *p*-version model (—) and with Abaqus (- -).</figcaption></figure>

## 3.2. Non-linear forced vibrations

To cause vibrations with large displacement, a localised vertical force sinusoidal in time is applied transversely to the shells at their apex, that is at point (0, 0, 0), as represented in Fig. 4. Four equal elements are employed, with shape functions of the *f* type and $p_o = p_θ = 7$; this model was chosen after convergence studies in linear and non-linear vibrations. The equations of motion are solved by Newmark's method. In all the following examples the dimensions *a* and *b* are equal to 0.3 m and the material properties are $E = 7 \ 10^{10}$ N/m$^2$, $ρ = 2778$ kg/m$^3$ and $ν = 0.3$, which are typical properties of aluminium. The shells are clamped at their straight sides and free at their curved sides. For several curvature radii, the central points of the shells belong to nodal lines of their first modes and are points where the second modes experience large displacements; therefore, the excitation frequency $ω_e$ will be equal to the second linear natural frequency of each shell.

Next, the results of the *p*-version model are compared with the ones computed using a Abaqus model with a 30 per 30 mesh. The geometric properties of the shell respect the relations *a*/*R* = 1.25 and *b*/*h* = 100. In the *p*-version model, viscous damping proportional to the stiffness of the form *C* = *α* *KI*, with a proportionality factor *α* = 10⁻⁵, was used and the same factor *α* was introduced in the Abaqus code. The excitation amplitude is 5000 N.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. Time plots of oscillations of Shells ♦a/R = 0.5; — a/R = 1; Δa/R = 1.5, at the following points and directions: (a) Point 2, direction z; (b) Point 1, direction x; (c) Point 1, direction y.</figcaption></figure>

The curves shown in Fig. 5 follow the same trend but do not exactly coincide. Because in numerical integration the solution at time *t* depends on the solution at time *t* − Δ*t*, the difference between the Abaqus and the *p*-version model grows as time increases. Nevertheless, taking into account the differences in the models that were briefly explained when element S8R5 of Abaqus was introduced in this text, a rather satisfactory agreement is achieved.

In the following numerical examples, the *p*-model is used to study the variation of the response with the curvature of the shell. The shells are excited by a vertical force with amplitude *A* = 2500 N. A stiffness proportional damping model is employed with a damping factor equal to 10<sup>-4</sup>.

The displacements and velocities are computed at Points 1, 2 and 3 represented in Fig. 4. The coordinates of Point 2 are the same for the three shells: (x, y, z) = (0, 0, 0); the local, adimensional coordinates of Points 1 and 3, are
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 7. Phase plots of oscillations of Shells ♦a/R = 0.5; — a/R = 1; Δa/R = 1.5, at the following points and directions: (a) Point 2, direction z; (b) Point 1, direction x; (c) Point 1, direction y.</figcaption></figure>

also the same for the three shells, but their global, physical coordinates naturally depend on the curvature radius of each shell.

The displacements and velocities are significantly influenced by the curvature of the shell, as shown in Figs 6 and in Fig. 7. The central point, Point 2, of the deepest shell, with the lowest curvature radius *a/R* = 1.5, experiences the lowest displacement and the smallest velocity. However, at Point 1 the inverse occurs and the shallow shell – *a/R* = 0.5 – experiences lower displacement and velocity than the other shells. All the oscillations are periodic and visibly the shallow shell is the one where the response further deviates from a harmonic motion.

In order to gain some insight on how the response of deep shells changes with the thickness, different thicknesses
<table><caption>Table 6
Maximum displacements and velocities of shells with different thicknesses</caption><thead><tr><th>Shell</th><th></th><th>h/a = 0.05</th><th>h/a = 0.01</th><th>h/a = 0.005</th><th>h/a = 0.001</th></tr></thead><tbody><tr><td rowspan="2">Point 2</td><td>Maximum displacement w(m)</td><td>8.38 × 10<sup>-5</sup></td><td>1.13 × 10<sup>-3</sup></td><td>2.05 × 10<sup>-3</sup></td><td>1.05 × 10<sup>-3</sup></td></tr><tr><td>Maximum velocity ṉ (m/s)</td><td>1.12</td><td>6.83</td><td>9.07</td><td>6.29 × 10<sup>-1</sup></td></tr><tr><td rowspan="2">Point 1</td><td>Maximum displacement w (m)</td><td>1.78 × 10<sup>-5</sup></td><td>4.07 × 10<sup>-04</sup></td><td>7.81 × 10<sup>-04</sup></td><td>4.24 × 10<sup>-4</sup></td></tr><tr><td>Maximum velocity ṉ (m/s)</td><td>2.56 × 10<sup>-1</sup></td><td>2.17</td><td>3.43</td><td>3.62 × 10<sup>-1</sup></td></tr><tr><td rowspan="2"></td><td>Maximum displacement u (m)</td><td>1.08 × 10<sup>-6</sup></td><td>2.12 × 10<sup>-04</sup></td><td>4.68 × 10<sup>-04</sup></td><td>2.32 × 10<sup>-4</sup></td></tr><tr><td>Maximum velocity ṉ (m/s)</td><td>1.56 × 10<sup>-2</sup></td><td>1.07</td><td>1.46</td><td>1.63 × 10<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 8. Time plots of oscillations of Shells ♦h/a = 0.05, —h/a = 0.01, Δ/a = 0.005, o h/a = 0.001, at the following points and directions: (a) Point 2, direction z; (b) Point 1, direction x; (c) Point 1, direction z.</figcaption></figure>

were considered. In these numerical tests, the radius of all shells are equal to the projected length, i.e. R = a, and the thicknesses are: h/a = 0.05, h/a = 0.01, h/a = 0.005 and h/a = 0.001. The excitation force has similar characteristics to the one of the previous examples, except in what the amplitude is concerned. In the case of shells
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 9. Phase plots of oscillations of Shells ♦h/a = 0.05, —h/a = 0.01, Δ/a = 0.005, o h/a = 0.001, at the following points and directions: (a) Point 2, direction z; (b) Point 1, direction x; (c) Point 1, direction z.</figcaption></figure>

h/a = 0.05, h/a = 0.01 and h/a = 0.005, the amplitude of excitation is 2500 h/0.03 N, i. e., the amplitude of the excitation is directly proportional to the thickness. In the very thin shell, h/a = 0.001, which is analysed as a limit case, the former rule was not followed and the amplitude of the excitation is 5 N. The reason for this choice is that an amplitude equal to 2500 h/0.03 would cause too large displacements in this shell and the present linear elastic model might fail to provide a reasonable approximation to reality.

The maximum amplitudes of positive displacement and velocity achieved in each shell are given in Table 6. As the thickness is decreased from h/a = 0.05, to h/a = 0.01 and h/a = 0.005 the displacement and velocity amplitudes increase, in spite of the fact that the amplitude of excitation is decreasing proportionally to the thickness. In the case
of the very thin shell $h/a = 0.001$ this trend is altered only because a very small force was employed.

Figures 8 and 9 give the displacement time histories and the phase plots in two points. As we just saw, the maximum displacements and velocities of the diverse shells are of very different magnitudes, and, in order to facilitate the comparison, the displacement and velocities in these figures were divided by their maximum values.

The thicker shell, where $h/a = 0.05$, experiences a harmonic motion with a very small amplitude of vibration; for example, in the middle point the transverse displacement oscillates between $+1.21 \ 10^{-3}h$ and $-1.21 \ 10^{-3}h$. Accordingly, the phase plot of this oscillation is an ellipse. As the thickness of the shells decreases, the amplitude of displacement inwards becomes greater than the one outwards and higher harmonics appear in the oscillations, that is non-linear effects become more important. This is particularly true in shells $h/a = 0.005$ and $h/a = 0.001$, where the displacements are larger in comparison with the shell thickness and, therefore, where the oscillations are more pronouncedly non-linear.

## 4. Conclusions

An approach that employs a p-version shallow shell finite element was suggested to model shells of any depth and to study their geometrically non-linear vibrations. Linear natural frequencies and displacements in geometrically non-linear regime originated by external forces were computed using the p-version model and compared with published results or obtained using the finite element software Abaqus and good agreement was found. The main advantage of the method here proposed in comparison with other finite element approaches is that, as demonstrated, it requires far fewer degrees of freedom.

Sinusoidal forces were employed to investigate the variation of the response of the shells with their curvature. In the conditions of this study, the response of the different shells was always periodic. It was found that shallower shells achieve larger displacements at the central point and are the shells where the response further deviates from a harmonic motion. Therefore, linear models will find a greater domain of applicability in deep shells. The variation of the response of deep shells with the thickness was also analysed. In this case, it is highly noticeable that even if the force is proportional to the thickness, the displacements increase very much as the thickness decreases. Therefore, thinner deep shells are more prone to experience non-linear oscillations then thicker ones. In any case, all the responses computed were periodic, whatever the shell thickness

## Acknowledgments

The partial support of this work by GRICES, Portugal, and PAI PESSOA, France, in the framework of project "Large Amplitude Vibrations of Structures" is gratefully acknowledged.

## References

[1] A.A. Popov, J.M.T. Thompson and F.A. McRobie, Low dimensional models of shell vibrations. Parametrically excited vibrations of cylindrical shells, *Journal of Sound and Vibration* **209** (1998), 163-186.

[2] A. Przekop, M. Salim Azzouz, X. Guo and C. Mei, Finite Element Multiple-Mode Approach to Nonlinear Free Vibrations of Shallow Shells, *AIAA Journal* **42** (2004), 2373-2381.

[3] H.C. Chan and W.C. Chung, Geometrically nonlinear analysis of shallow shells using higher order finite elements, *Computers and Structures* **31** (1989), 329-338.

[4] W. Han, M. Petyt and K.-M. Hsiao, An investigation into geometrically nonlinear analysis of rectangular laminated plates using the hierarchical finite element method, *Finite Elements in Analysis and Design* **18** (1994), 273-288.

[5] B.A. Szabo and I. Babuska, *Finite Element Analysis*, Wiley, New York, 1991.

[6] N.S. Bardell, J.M. Dunsdon and R.S. Langley, On the free vibration of completely free, open, cylindrically curved, isotropic shell panels, *Journal of Sound and Vibration* **207** (1997), 647-669.

[7] P. Ribeiro, A hierarchical finite element for geometrically non-linear vibration of doubly curved, moderately thick isotropic shallow shells, *International Journal for Numerical Methods in Engineering* **56** (2003), 715-738.

[8] P. Ribeiro, Forced Large Amplitude periodic vibrations of cylindrical shallow shells, *Finite Elements in Analysis and Design* **44** (2008), 657-674.
[9] E. Ventsel and T. Krauthammer, *Thin Plates and Shells. Theory Analysis and Applications*, Marcel Dekker, New York, 2001.

[10] M.S. Qatu, *Vibration of Laminated Shells and Plates*, Elsevier, Amsterdam, 2004.

[11] W. Flügge, *Stresses in Shells*, Springer-Verlag, Berlin, 1973.

[12] K.M. Liew, M.K. Lim and S. Kitipornchai, Vibration of shallow shells: a review with bibliography, *Applied Mechanics Review* **50** (1997), 431–444.

[13] A.W. Leissa and Y. Narita, Vibrations of completely free shallow shells of rectangular planform, *Journal of Sound and Vibration* **96** (1984), 207–218.

[14] A. Abe, Y. Kobayashi and G. Yamada, Non-linear vibration characteristics of clamped laminated shallow shells, *Journal of Sound and Vibration* **234** (2000), 405–426.

[15] M. Amabili, Nonlinear vibrations of circular cylindrical panels, *Journal of Sound and Vibration* **281** (2005), 509–535.

[16] L.M. Dai, Q. Han and M.Z. Dong, A single and double mode approach to chaotic vibrations of a cylindrical shell with large deflection, *Shock and Vibration* **11** (2004), 533–546.

[17] M. Amabili, *Nonlinear Vibrations and Stability of Shells and Plates*, Cambridge University Press, Cambridge, 2008.

[18] A. Selmane and A.A. Lakis, Influence of geometric non-linearities on the free vibrations of orthotropic open cylindrical shells, *International Journal for Numerical Methods in Engineering* **40** (1997), 1115–1137.

[19] D. Kuhl and E. Ramm, Constraint Energy Momentum Algorithm and its application to non-linear dynamics of shells, *Computer Methods in Applied Mechanics and Engineering* **136** (1996), 293–315.

[20] J. Argyris, M. Papadrakakis and Z.S. Mouroutis, Nonlinear dynamic analysis of shells with the triangular element TRIC, *Computer Methods in Applied Mechanics and Engineering* **192** (2003), 3005–3038.

[21] E. Reissner, The effect of transverse shear deformation on the bending of elastic plates, *Journal of Applied Mechanics* **12** (1945), 69–77.

[22] R.D. Mindlin, Influence of rotatory inertia and shear on flexural vibrations of isotropic, elastic plates, *Journal of Applied Mechanics* **18** (1951), 31–38.

[23] O.C. Zienkiewicz and R.L. Taylor, *The Finite Element Method*, 4th Edition, McGraw-Hill, London, 1988.

[24] M. Petyt, *Introduction to Finite Element Vibration Analysis*, Cambridge University Press, Cambridge, 1990.

[25] K.-J. Bathe, *Finite Element Procedures*, Prentice Hall, Upper Saddle River, 1996.

[26] E. Blanchard, M. Jouet and S. Bellizzi, Vibrations Forcées de Coques Minces Sous Fortes Sollicitations, Internal Report, Ecole Généraliste d'Ingénieurs de Marseille, 2005.

[27] C.W. Lim and K.M. Liew, A pb-2 Ritz Formulation for Flexural Vibration of Shallow Cylindrical Shells of Rectangular Planform, *Journal of Sound and Vibration* **173** (1994), 343–375.

[28] ABAQUS/Standard, *User's Manual*, Volume II, version 5.6. Hibbit, Karlson & Sorensen, Inc., 1996.