# Non-linear vibrations of shallow circular cylindrical panels with complex geometry. Meshless discretization with the R-functions method

Galina Pilgun, Marco Amabili

► To cite this version:

Galina Pilgun, Marco Amabili. Non-linear vibrations of shallow circular cylindrical panels with complex geometry. Meshless discretization with the R-functions method. International Journal of Non-Linear Mechanics, 2012, 47, pp.137 - 152. 10.1016/j.ijnonlinmec.2011.11.011. hal-01403853

HAL Id: hal-01403853

https://hal.science/hal-01403853

Submitted on 28 Nov 2016

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 Attribution 4.0 International License
# Non-linear vibrations of shallow circular cylindrical panels with complex geometry. Meshless discretization with the R-functions method

G. Pilgun, M. Amabili

Department of Mechanical Engineering, McGill University, 817 Sherbrooke Street West, Montreal, Québec, Canada H3A 2K6

Geometrically non-linear forced vibrations of a shallow circular cylindrical panel with a complex shape, clamped at the edges and subjected to a radial harmonic excitation in the spectral neighborhood of the fundamental mode, are investigated. Both Donnell and the Sanders–Koiter non-linear shell theories retaining in-plane inertia are used to calculate the elastic strain energy. The discrete model of the non-linear vibrations is built using the meshfree technique based on classic approximate functions and the R-function theory, which allows for constructing the sequences of admissible functions that satisfy given boundary conditions in domains with complex geometries; Chebyshev orthogonal polynomials are used to expand shell displacements. A two-step approach is implemented in order to solve the problem: first a linear analysis is conducted to identify natural frequencies and corresponding natural modes to be used in the second step as a basis for expanding the non-linear displacements. Lagrange approach is applied to obtain a system of ordinary differential equations on both steps. Different multimodal expansions, having from 15 up to 35 generalized coordinates associated with natural modes, are used to study the convergence of the solution. The pseudo-arclength continuation method and bifurcation analysis are applied to study non-linear equations of motion. Numerical responses are obtained in the spectral neighborhood of the lowest natural frequency; results are compared to those available in the literature. Internal resonances are also detected and discussed.

## 1. Introduction

Circular cylindrical shells (or panels) are largely used in aero-nautics, aerospace and other engineering fields. They are subjected to dynamic loads that can cause vibration amplitude of the order of the shell thickness, giving rise to significant non-linear phenomena. Refs. [1–4] provide the most complete reviews of studies on non-linear vibrations of shells before 2003. First contributions to non-linear dynamics of circular cylindrical shells were made in Refs. [5–7]. A single or two modes expansion of the shell transverse displacement and the Donnell shell theory with tangential and rotary inertia omitted in order to simplify the mathematical formulation were initially used. In many of the past studies, the Donnell non-linear shallow shell theory has been applied to model non-linear vibrations of panels for its simplicity [8–13]. More refined classical non-linear shell theories have been applied to study large-amplitude vibrations of cylindrical shells in Refs. [1, 4, 14–20]. In particular, the Sanders–Koiter (also referred as Sanders), the Flügge–Lur'e–Byrne, the Novozhilov and the Donnell non-linear shell theories retaining in-plane

inertia have been used. A detailed discussion on their strengths and weaknesses has been given by Amabili [4,18].

Another important refinement in studies on non-linear dynamics of shells is achieved increasing the number of modes (i.e. degrees of freedom (dofs)) retained in the approximating functions [10,17–22]. Refs. [19,20] present a comprehensive discussion on the importance of introducing the multimode expansions of shell displacements into analysis.

In general, the majority of the available studies address the non-linear dynamics of circular cylindrical panels or double curved shallow shells with simple geometric shapes; e.g. with a rectangular base [5–13,17–20,22]. The complex geometry of the shell significantly increases the numerical difficulty of the problem. In fact, models with thousands of degrees of freedom obtained for instance by commercial finite-element method (FEM) present large computational problems, with possible divergence of the solution or impossibility to perform a full bifurcation analysis. The discussion on non-linear vibrations of cylindrical shells using FEM can be found, e.g. in Refs. [23–26]. In case of global discretization using admissible functions, in the framework of the Rayleigh–Ritz method or the Lagrangian approach, complex geometry makes their mathematical representation very difficult.

One of the alternatives to the finite-element method is the effective mesh-free discretization technique based on classic approximate functions and the R-functions theory (also referred
as the R-functions method, or shortly RFM). It discretizes the system of non-linear partial differential equations using admissible functions that satisfy given boundary conditions; this is a very difficult task in cases of complex domains. The R-functions theory, developed by the Ukrainian mathematician Rvachev [27], suggests a way for the analytical description of domains with complex boundaries. It allows for the construction of approximate functions that satisfy exactly boundary conditions of different types, including mixed boundary conditions. Therefore, the R-functions method (RFM) is a meshless method that allows all prescribed boundary conditions to be satisfied exactly at all boundary points. When the R-functions method is coupled to the variational Rayleigh-Ritz method, a global discretization of the problem is obtained surpassing the difficulty of expressing admissible functions that satisfy the geometric boundary conditions. A review of mesh-free methods can be found in the book by Liu [28]. The approach based on the R-functions method coupled to the variational Rayleigh-Ritz method has been successfully applied to linear and non-linear vibrations of plates [29-31] and shallow shells [32-35], including circular cylindrical panels with geometry similar or equal to the case of the present study [33,35]. In particular, a single degree of freedom has been used in Ref. [33] to study non-linear free vibrations. This single degree of freedom is associated to the mode shape of the fundamental mode, and it is only a first approximation to the problem solution. Similar limitations are observed in Refs. [34,35] that cannot guarantee the convergence of the solution. For this reason, the present study introduces, for the first time, a multi-mode expansion in the R-functions used to the study of non-linear vibrations of shells.

In the present study, geometrically non-linear forced vibrations of a shallow circular cylindrical panel with a complex shape, clamped at the edges and subjected to a radial harmonic excitation in the spectral neighborhood of the fundamental mode, are investigated. Two different non-linear strain-displacement relationships, according to the Donnell and the Sanders-Koiter shell theories retaining in-plane inertia, are used to calculate the elastic strain energy. It must be observed that, for the panel considered in the section of numerical results, the two theories give practically coincident results, as expected since the panel is shallow and very thin. However, for deep shells the Sanders-Koiter theory has been proven to be more accurate (see e.g. [18]). It must be observed here that the computational effort necessary to obtain the numerical solution is significantly larger for the Sanders-Koiter theory since many additional non-linear terms appear in the discretized equations of motion.

The R-functions theory is applied to construct the sequences of admissible functions that satisfy the given boundary conditions in domains with complex shapes; Chebyshev orthogonal polynomials

[21,36] are used to expand shell displacements. A two-step approach is implemented in order to solve the problem: first a linear analysis is conducted to identify natural frequencies and corresponding natural modes to be used in the second step as a basis for expanding the non-linear displacements. The Lagrange variational approach is applied to obtain a system of ordinary differential equations. This approach introduces multimode expansions of the shell displacements for the first time in studies involving the R-function method to study non-linear vibrations. In the study, different multimodal expansions having from 15 up to 37 generalized coordinates, associated with natural modes of the panel, have been used in order to study the convergence of the solution. The pseudo-arclength continuation method and bifurcation analysis are applied to study the non-linear equations of motion. Numerical responses are obtained in the spectral neighborhood of the lowest natural frequency. Results are compared to those available in the literature. Internal resonances are also identified and discussed.

## 2. Problem formulation

An isotropic shallow circular cylindrical panel of uniform thickness h with complex shape of the base is considered, as shown in Fig. 1. The principal lines of curvature of the middle surface coincide with the coordinates x, y of the Cartesian coordinate system, and z is directed along the normal to the middle surface of the shell, as shown in Fig. 1; the origin is placed at the middle of the panel base in order to take advantage of the symmetry of the structure. The displacements of an arbitrary point of coordinates (x, y) on the middle surface of the panel are denoted by $u(x, y, t)$ in the axial direction, $v(x, y, t)$ in the circumferential direction and $w(x, y, t)$ in the radial direction, respectively; it is assumed that the deflection of the middle surface of the panel is of the same order of the thickness and is taken positive outwards.

### 2.1. Strain energy

In the present study, two different non-linear shell theories, namely Donnell and Sanders-Koiter theories, retaining in-plane inertia, are used. Both theories are based on Love's first approximation assumptions. According to them, strain components $\epsilon_x, \epsilon_y$ and $\gamma_{xy}$ at an arbitrary point of the panel are related to the middle surface strains $\epsilon_{x,0}, \epsilon_{y,0}$ and $\gamma_{xy,0}$, and to the changes in the curvature and torsion of the middle surface $k_x, k_y$ and $k_{xy}$ by the following relationships:

$$ \epsilon_x = \epsilon_{x,0} + zk_x, \quad \epsilon_y = \epsilon_{y,0} + zk_y, \quad \gamma_{xy} = \gamma_{xy,0} + zk_{xy} \quad (1) $$

<figure><img src="image_2.png" /><figcaption>Fig. 1. Circular cylindrical panel with complex base. (a) Coordinate system and dimensions; (b) the function of the boundary $\Omega(x,y)$.</figcaption></figure>
where z is the distance of the arbitrary point of the shell from the middle surface. The middle surface strain-displacement relationships and changes in the curvature and torsion have different expressions for the Donnell and the Sanders-Koiter shell theories.

According to the Donnell non-linear shell theory, the middle surface strain-components and the changes in curvature and torsion are given by [4,10]

$$ \epsilon_{x,0} = \frac{\partial u}{\partial x} + \frac{1}{2} \left( \frac{\partial w}{\partial x} \right)^2, \qquad (2a) $$

$$ \epsilon_{y,0} = \frac{\partial v}{\partial y} + \frac{w}{R} + \frac{1}{2} \left( \frac{\partial w}{\partial y} \right)^2, \qquad (2b) $$

$$ \gamma_{xy,0} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial x} \frac{\partial w}{\partial y}, \qquad (2c) $$

$$ k_x = -\frac{\partial^2 w}{\partial x^2}, \qquad (2d) $$

$$ k_y = -\frac{\partial^2 w}{\partial y^2}, \qquad (2e) $$

$$ k_{xy} = -2 \frac{\partial^2 w}{\partial x \partial y}, \qquad (2f) $$

where R is the constant middle surface principal radius of curvature.

The Sanders-Koiter shell theory is a more refined non-linear classical theory developed by Sanders [37] and Koiter [38], separately. It considers finite deformations with small strains and moderately small rotations, which allows for non-linear terms that depend on in-plane displacements to be retained.

According to the Sanders-Koiter non-linear shell theory, the middle surface strain-displacement relationships and changes in the curvature and torsion for a circular cylindrical shell are given by [4,18]

$$ \epsilon_{x,0} = \frac{\partial u}{\partial x} + \frac{1}{2} \left( \frac{\partial w}{\partial x} \right)^2 + \frac{1}{8} \left( \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \right)^2, \qquad (3a) $$

$$ \epsilon_{y,0} = \frac{\partial v}{\partial y} + \frac{w}{R} + \frac{1}{2} \left( \frac{\partial w}{\partial y} - \frac{v}{R} \right)^2 + \frac{1}{8} \left( \frac{\partial u}{\partial y} - \frac{\partial v}{\partial x} \right)^2, \qquad (3b) $$

$$ \gamma_{xy,0} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial x} \left( \frac{\partial w}{\partial y} - \frac{v}{R} \right), \qquad (3c) $$

$$ k_x = -\frac{\partial^2 w}{\partial x^2}, \qquad (3d) $$

$$ k_y = \frac{\partial v}{R \partial y} - \frac{\partial^2 w}{\partial y^2}, \qquad (3e) $$

$$ k_{xy} = -2 \frac{\partial^2 w}{\partial x \partial y} + \frac{1}{2R} \left( 3 \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \right). \qquad (3f) $$

The elastic strain energy $U_s$ of a circular cylindrical shell, within the limits of Kirchhoff-Love assumptions, is given by [4,10]

$$ U_s = \frac{1}{2} \iint_S \int_{-h/2}^{h/2} (\sigma_x \varepsilon_x + \sigma_y \varepsilon_y + \tau_{xy} \gamma_{xy}) dS (1+z/R) dz, \qquad (4) $$

where S is the shell middle surface; the stresses $\sigma_x$, $\sigma_y$ and $\tau_{xy}$ are related to the strains for homogeneous and isotropic material ($\sigma_z=0$, case of plane stress) by

$$ \sigma_x = \frac{E}{1-v^2} (\varepsilon_x + v \varepsilon_y), \quad \sigma_y = \frac{E}{1-v^2} (\varepsilon_y + v \varepsilon_x), \quad \tau_x = \frac{E}{2(1+v)} \gamma_{xy}, \qquad (5) $$

where E is the Young modulus and v is Poisson's ratio. Using Eqs. (1), (4) and (5), the following expression is obtained:

$$ \begin{aligned} U_s = & \frac{1}{2} \frac{Eh}{1-v^2} \iint_S \left( \varepsilon_{x,0}^2 + \varepsilon_{y,0}^2 + 2v\varepsilon_{x,0}\varepsilon_{y,0} + \frac{1-v}{2}\gamma_{xy,0}^2 \right) dS \\ & + \frac{1}{2} \frac{Eh^3}{12(1-v^2)} \iint_S \left( k_x^2 + k_y^2 + 2vk_xk_y + \frac{1-v}{2}k_{xy}^2 \right) dS \\ & + \frac{1}{2R} \frac{Eh^3}{6(1-v^2)} \iint_S \left( \varepsilon_{x,0}k_y + \varepsilon_{y,0}k_x + v\varepsilon_{x,0}k_y + v\varepsilon_{y,0}k_x + \frac{1-v}{2}\gamma_{xy,0}k_{xy} \right) dS + O(h^4), \end{aligned} \qquad (6) $$

where $O(h^4)$ is a higher order term in h according to the Sanders-Koiter theory; the last term in $h^3$ must be omitted, if $z/R$ is neglected with respect to unity in Eq. (4), according to the Donnell non-linear theory [4,10]. The right-hand side of Eq. (6) can be easily interpreted: the first term is the membrane (also referred as stretching) energy, the second one represents the bending energy and, if the last term is retained, membrane and bending energies are coupled.

## 2.2. Kinetic energy, virtual work and damping

The kinetic energy $T_s$ of the panel, by neglecting rotary inertia, is given by

$$ T_s = \frac{1}{2} \rho h \iint_S (\dot{u}^2 + \dot{v}^2 + \dot{w}^2) dS, \qquad (7) $$

where $\rho$ is the mass density of the panel. In Eq. (7) the over-dot denotes a time derivative.

The virtual work W done by the external forces is written as [4]

$$ W = \iint_S (q_x u + q_y v + q_z w) dS, \qquad (8) $$

where $q_x$, $q_y$ and $q_z$ are the distributed forces per unit area acting on the shell in x, y and normal directions, respectively. In the present study, only a single harmonic normal force is considered; therefore $q_x = q_y = 0$. An external, distributed load $q_z$ normal to the shell surface is applied due to the radial concentrated force $\tilde{f}$ and is given by [4]

$$ q_z = \tilde{f} \delta(x-\tilde{x}) \delta(y-\tilde{y}) \cos(\omega t), \qquad (9) $$

where $\omega$ is the excitation frequency, t is the time, $\delta$ is the Dirac delta function, $\tilde{f}$ is the amplitude of the force applied along the normal to the panel middle surface and taken positive in z direction, $\tilde{x}$ and $\tilde{y}$ determine the position of the point of the force application; in present case, the excitation point is located at the origin of the reference system. Eq. (9) is transformed as follows:

$$ W = \tilde{f} \cos(\omega t) w \Big|_{x=0, y=0}, \qquad (10) $$

The non-conservative damping forces are assumed to be of viscous type and are taken into account using the Rayleigh dissipation function [4]

$$ F = \frac{1}{2} c \iint_S (\dot{u}^2 + \dot{v}^2 + \dot{w}^2) dS, \qquad (11) $$

where c is the viscous damping coefficient; it has a different value for each term of the mode expansion.

# 3. The solution method

The problem is solved in two steps: (i) the linear modal analysis is conducted to obtain the natural frequencies and corresponding natural modes; (ii) selected natural modes are used as a basis for expanding the non-linear displacements in the non-linear analysis.
### 3.1. Linear free vibrations: modal analysis using the R-functions method

For the purpose of carrying out a linear vibration analysis, only linear terms are considered in the shell middle surface strain-displacement relationships for both shell theories, i.e. only quadratic terms are retained in Eqs. (2a)-(2c) and (3a)-(3c).

In the present study, natural frequencies and mode shapes are found using the Lagrange equations and the R-functions method (RFM) in order to expand shell displacements using admissible functions, which satisfy the given boundary conditions.

The main idea of the R-functions method is to build the so called solution structures of the boundary value problem, which are appropriate for domains with arbitrary shapes. Such structures form the basis of the appropriate sequences of admissible functions in a form that satisfies exactly all the boundary conditions and contain functions to be determined in order to satisfy the differential equations governing the problem in an approximate way. Details on the R-function method are given in the appendix.

Vibrations of the shell are assumed to be harmonic in time. The space and time variable are separated, which yields

$$u(x,y,t) = U(x,y)f(t), \quad v(x,y,t) = V(x,y)f(t), \quad w(x,y,t) = W(x,y)f(t), \qquad (12)$$

where $U(x,y), V(x,y), W(x,y)$ represent the modal shape and $f(t)$ is a harmonic function, i.e.

$$f(t) = \cos(\omega t). \qquad (13)$$

The studied shell is considered to be clamped, which means

$$u_n|_{\partial S} = 0, \quad v_n|_{\partial S} = 0, \quad w|_{\partial S} = 0, \quad \frac{\partial w}{\partial n}|_{\partial S} = 0, \qquad (14)$$

where $u_n = ul + vm, v_n = -um + vl$ and $l, m$ are the direction cosines of the edge with respect to the in-plane axes $x$ and $y$ (in case of curved edge, $l$ and $m$ are functions of the position); $n$ and $\tau$ are the normal and the tangent to the boundary, respectively; both of them are contained in the plane tangent to the panel at the boundary, with the normal pointed outward the shell domain. Corner points do not need any special treatment.

According to the R-function theory, the structure of admissible functions for the clamped shallow shell is [29,33]

$$\begin{aligned} U(x,y) &= \Omega(x,y)P_1(x,y), & V(x,y) &= \Omega(x,y)P_2(x,y), \\ W(x,y) &= \Omega^2(x,y)P_3(x,y), & & \end{aligned} \qquad (15)$$

where functions $P_i$ ($i=1,2,3$) are functions expanded in the following way:

$$\begin{aligned} P_1(x,y) &= \sum_{j=1}^{N_1} q_j \varphi_{1j}(x,y), & P_2(x,y) &= \sum_{j=N_1+1}^{N_2} q_j \varphi_{2j}(x,y), \\ P_3(x,y) &= \sum_{j=N_2+1}^{N_3} q_j \varphi_{3j}(x,y), & & \end{aligned} \qquad (16)$$

where $q_j, j=1,...,N_3$ are unknown coefficients to be determined from the corresponding eigenvalue problem and $\{\varphi_{ij}\}$, ($i=1,2,3; j=1,...,N_3$) are any complete sequences of functions, e.g. power polynomials, splines and trigonometric polynomials; in the present study, Chebyshev orthogonal polynomials are used since they have given good results for complete circular shells [39,40]. Thus, mode shapes, expanded in double series in terms of Chebyshev polynomials, are written as

$$U(x,y) = \Omega(x,y) \sum_{i=0}^{M_u} \sum_{j=0}^{M_u-1} u_{ij} T_{1i}(x) T_{1j}(y), \qquad (17a)$$

$$V(x,y) = \Omega(x,y) \sum_{i=0}^{M_u} \sum_{j=0}^{M_u-1} v_{ij} T_{1i}(x) T_{1j}(y), \qquad (17b)$$

$$W(x,y) = \Omega^2(x,y) \sum_{i=0}^{M_w} \sum_{j=0}^{M_w-1} w_{ij} T_{1i}(x) T_{1j}(y), \qquad (17c)$$

where $T_{1i}(\cdot)$ is the ith order Chebyshev polynomial of the first kind, $u_{ij}, v_{ij}$ and $w_{ij}$ are unknown coefficients, which have to be re-arranged into a single indexed sequence $q_j, j=1,...,N_3$, in order to obtain Eq. (16) and the corresponding eigenvalue problem. This can be easily achieved first by expanding the double series over Chebyshev polynomials in Eqs. (17) and then reordering the coefficients in the following way:

$$\mathbf{q} = \{w_{0,0}^{(1)}, w_{0,1}^{(2)}, \dots, w_{M_w, M_w-1}^{(N_1)}, u_{0,0}^{(N_1+1)}, \dots, u_{M_u, M_u-1}^{(N_2)}, v_{0,0}^{(N_2+1)}, \dots, v_{M_u, M_u-1}^{(N_3)}\}. \qquad (18)$$

In Eqs. (15) and (17) $\Omega(x,y)=0$ is the function of the boundary of the panel normalized up to the first order. Normalization of the function $\Omega(x,y)$ is given by [27,29]

$$\begin{aligned} \Omega(x,y) &= 0, & \forall(x,y) \in \partial S, \\ \frac{\partial \Omega(x,y)}{\partial n} &= -1, & \forall(x,y) \in \partial S, \end{aligned}$$

$$\Omega(x,y) > 0, \quad \forall(x,y) \in S.$$

In the case of the cylindrical panel with the geometry shown in Fig. 1(a), the function of the boundary $\Omega(x,y)$ can be composed in the following form [4,33]:

$$\Omega(x,y) = \eta_1 \wedge 0 \eta_2 \wedge 0 (\eta_3 \vee 0 \eta_4), \qquad (19)$$

where the functions $\eta_i$ ($i=1,...,4$), are given by

$$\begin{aligned} \eta_1 &= \frac{1}{2b}(b^2 - x^2) \ge 0, & \eta_2 &= \frac{1}{2a}(a^2 - y^2) \ge 0, \\ \eta_3 &= \frac{1}{2d}(d^2 - y^2) \ge 0, & \eta_4 &= \frac{1}{2c}(x^2 - c^2) \ge 0. \end{aligned}$$

The function $\Omega(x,y)$ for the panel shown in Fig. 1 is plotted in Fig. 1(b).

The complete procedure of building the equation of the domain $\Omega(x,y)=0$ for the panel shown in Fig. 1 is addressed in Refs. [4,33] in details. In Eq. (19) the R-conjunction $(x \wedge 0y)$ and the R-disjunction $(x \vee 0y)$ have been used. They are defined as [27]

$$\begin{aligned} x \wedge 0y &= x + y - \sqrt{x^2 + y^2}, \\ x \vee 0y &= x + y + \sqrt{x^2 + y^2}. \end{aligned}$$

Eqs. (17)-(19) are first inserted into Eq. (12) and then into the linear expression of the strain energy, Eq. (6), and the kinetic energy (7), respectively. Then, for harmonic linear vibration, it is possible to write

$$U_S = \frac{1}{2} \mathbf{q}^T \mathbf{K} \mathbf{q}, \quad T_S = \frac{1}{2} \dot{\mathbf{q}}^T \mathbf{M} \dot{\mathbf{q}}, \quad \dot{\mathbf{q}} = \omega \mathbf{q}, \qquad (20)$$

where $\omega$ are the natural frequencies and $\mathbf{q} = (q_1, q_2, \dots, q_{N_3})^T$ are the corresponding eigenvectors giving the mode shapes. The components of the vector $\mathbf{q}$ correspond to the reordered sequences of the unknown coefficients in Eq. (17). $\mathbf{K} = \{k_{ij}\}_{i,j=1,...,N_3}$ and $\mathbf{M} = \{m_{ij}\}_{i,j=1,...,N_3}$ are the stiffness matrix and the mass matrix, respectively.

The Rayleigh-Ritz method yields the classical non-standard eigenvalue problem:

$$(\mathbf{K} - \omega^2 \mathbf{M}) \mathbf{q} = 0. \qquad (21)$$

It is important to point out that coefficients in both matrices $\mathbf{K}$ and $\mathbf{M}$ are double surface integrals over the domain with complex shape to be numerically computed with accuracy. Eq. (21) gives
the natural frequencies $\omega_m$, for $m=1,...,N_3$, and the corresponding eigenvectors $q_m$. The mode shapes are obtained once $q_m$ is inserted back in Eq. (16) and then in (15). The mode shapes $W_m(x,y)$, $U_m(x,y)$ and $V_m(x,y)$ obtained are then normalized by $W_m(x,y)/\max(W_m(x,y))$, $U_m(x,y)/\max(U_m(x,y))$ and $V_m(x,y)/\max(V_m(x,y))$ for any $(x,y)$ in order to have the maximum amplitude in the shell domain $S$ equal to unit.

## 3.2. Non-linear vibrations

In the non-linear analysis, the full expressions of the panel middle surface strain-displacement relationships, Eqs. (2) or (3) and Eq. (6), are used to evaluate the potential energy, which contains terms up to the fourth order. The panel displacements $u, v, w$ are expanded using the normalized natural modes as a

<figure><img src="image_3.png" /><figcaption>**Table 2**

Comparison of natural frequencies (Hz) of symmetric modes of the circular cylindrical panel.</figcaption></figure>

<table><thead><tr><th rowspan="2">Methods</th><th colspan="7">Modes (m,n)</th></tr><tr><th>(1,1)</th><th>(1,3)</th><th>(3,1)</th><th>(1,5)</th><th>(3,3)</th><th>(3,5)</th><th>(1,7)</th></tr></thead><tbody><tr><td>FEM (SOLID186)</td><td>367.39</td><td>451.24</td><td>626.20</td><td>678.22</td><td>923.21</td><td>1057.00</td><td>1120.50</td></tr><tr><td>RFM (Donnell theory)</td><td>368.64</td><td>453.44</td><td>627.84</td><td>689.51</td><td>933.44</td><td>1089.47</td><td>1152.45</td></tr><tr><td>RFM (Sanders-Koiter theory)</td><td>368.64</td><td>453.43</td><td>627.83</td><td>689.50</td><td>933.43</td><td>1089.46</td><td>1152.45</td></tr></tbody></table>
basis

$$w(x,y,t) = \sum_{m=1}^{M_1} W_m(x,y)p_m(t), \qquad (22a)$$

$$u(x,y,t) = \sum_{m=M_1+1}^{M_2} U_m(x,y)p_m(t), \qquad (22b)$$

$$v(x,y,t) = \sum_{m=M_2+1}^{M_3} V_m(x,y)p_m(t), \qquad (22c)$$

where $p_m(t)$ are the generalized coordinates. The total number of degrees of freedom in the non-linear analysis $M_3$ is much smaller than $N_3$ used in the linear analysis. Also, the number of terms necessary for the expansion of the in-plane displacement is larger than for the normal displacement $w$. The advantage of building a reduced-order non-linear model is in the possibility of having a full bifurcation analysis and at the same time to reduce numerical problems like divergence. The generalized coordinates to be retained in the expansion (22) have to be carefully chosen in order to reach accuracy of the solution.

The normalized mode shapes $W_m(x,y)$, $U_m(x,y)$ and $V_m(x,y)$ are known functions expressed in terms of Chebychev polynomials. In Eq. (22a)-(22c) the generalized coordinates $p_m(t), m=1,...,M_3$

<table><caption>**Table 3**

Natural frequencies (Hz) of the circular cylindrical panel.</caption><thead><tr><th>Modes (n,m)</th><th>FEM (SOLID186)</th><th>RFM (Donnell theory)</th></tr></thead><tbody><tr><td>1,1</td><td>367.39</td><td>368.64</td></tr><tr><td>1,2</td><td>384.48</td><td>385.52</td></tr><tr><td>2,1</td><td>393.73</td><td>396.73</td></tr><tr><td>2,2</td><td>408.14</td><td>409.88</td></tr><tr><td>1,3</td><td>451.24</td><td>453.44</td></tr><tr><td>1,4</td><td>544.76</td><td>548.27</td></tr><tr><td>2,3</td><td>570.70</td><td>579.91</td></tr><tr><td>3,1</td><td>626.20</td><td>627.84</td></tr><tr><td>3,2</td><td>627.78</td><td>629.53</td></tr><tr><td>1,5</td><td>678.22</td><td>689.51</td></tr><tr><td>2,4</td><td>680.86</td><td>695.70</td></tr><tr><td>2,5</td><td>811.32</td><td>823.77</td></tr><tr><td>1,6</td><td>843.50</td><td>853.88</td></tr><tr><td>4,2</td><td>897.98</td><td>907.23</td></tr><tr><td>3,3</td><td>923.21</td><td>933.44</td></tr><tr><td>3,5</td><td>1057.00</td><td>1089.47</td></tr><tr><td>1,7</td><td>1120.50</td><td>1152.45</td></tr></tbody></table>

are unknown time functions; due to normalization, the generalized coordinates represent the maximum amplitude of vibration in that specific mode shape since the maximum of $W_m(x,y)$, $U_m(x,y)$ and $V_m(x,y)$ after normalization is equal to unit.

Eqs. (22a)-(22c) are inserted in the expressions of the strain and kinetic energies (6) and (7), virtual work (10) and damping (11). Calculations for damping give

$$F = \frac{1}{2} \left( \sum_{m=1}^{M_1} c_m \dot{p}_m^2 A_m + \sum_{m=M_1+1}^{M_2} c_{m-M_1} \dot{p}_m^2 A_m + \sum_{m=M_2+1}^{M_3} c_{m-M_2} \dot{p}_m^2 A_m \right), \qquad (23a)$$

$$A_m = \begin{cases} \iint_S W_m^2 dS, & \text{if } m = 1, \dots, M_1, \\ \iint_S U_m^2 dS, & \text{if } M_1 + 1 \le m < M_2, \\ \iint_S V_m^2 dS, & \text{if } M_2 + 1 \le m < M_3. \end{cases} \qquad (23b)$$

In Eq. (23a) the damping coefficient $c_m$ is the viscous modal damping coefficient that can be related to the damping ratio.

The generalized forces $Q_j$ are obtained by differentiation of Rayleigh's dissipation function and of the virtual work done by external forces [4]

$$Q_j = -\frac{\partial F}{\partial \dot{p}_j} + \frac{\partial W}{\partial p_j} = -A_j c_j \dot{p}_j + \begin{cases} 0 \\ \tilde{f} \cos(\omega t) \quad \text{if } 1 \le j < M_1, \quad j = 1, \dots, M_3. \end{cases} \qquad (24)$$

The Lagrange equations of motion are

$$\frac{d}{dt} \left( \frac{\partial T}{\partial \dot{p}_i} \right) - \frac{\partial T}{\partial p_i} + \frac{\partial U}{\partial p_i} = Q_i, \quad i = 1, 2, \dots, M_3, \qquad (25)$$

where $\partial T/\partial p_i=0$. These second-order differential equations have very long expressions containing quadratic and cubic non-linear terms. In particular, the complicated term containing the linear terms and quadratic and cubic non-linearities can be written in the form

$$\frac{\partial U}{\partial p_i} = \sum_{k=1}^{M_3} p_k \tilde{z}_{k,i} + \sum_{k,j=1}^{M_3} p_j p_k \tilde{z}_{j,k,i} + \sum_{k,j,l=1}^{M_3} p_j p_k p_l \tilde{z}_{j,k,l,i}, \qquad (26)$$

where coefficients $\tilde{z}$ can be obtained only numerically and the linear terms $\tilde{z}_{k,i}$ in Eq. (26) are equal to zero for any $i \neq k$ (for $0 < i \le M_1$) since normal modes are used to expand the unknown non-linear displacements. Linear coupling is present only among the components $W_m(x,y)$, $U_m(x,y)$ and $V_m(x,y)$ of the same mode shape, i.e. with the same $m$.

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Amplitude of the vibration response of the panel versus the excitation frequency; excitation frequency around the fundamental mode (1,1), nondimensional excitation $f=0.02$, where $f = \tilde{f}/(h\omega_{1,1}\mu_{1,1})$, $\mu_{1,1}$ is the modal mass of the mode (1,1) and $\zeta_{1,1}=0.004$; model with 27 dofs; Donnell theory. ——, stable solution; –, unstable solution; --, backbone curve from Kurpa et al. [33]. (a) Maximum of the generalized coordinate $w_{1,1}$ (in a vibration period); (b) minimum of the generalized coordinate $w_{1,1}$ changed of sign.</figcaption></figure>
Gaussian quadrature in the Mathematica software [41] has been used in order to compute double surface integrals in Eqs. (6) and (7). The system of Lagrange equations is pre-multiplied by the inverse of the mass matrix **M** and then each equation is transformed in two equations of the first order. For computational convenience, variables are non-dimensionalized in the following way: the frequencies are non-dimensionalized dividing by the natural frequency of the resonant mode (i.e. the one excited near its resonance frequency) and the vibration amplitudes are divided by the shell thickness *h*. The resulting $M_3$ equations are studied using the software AUTO [42] for numerical continuation and bifurcation analysis of non-linear ordinary differential equations. The software AUTO is capable of continuation of the solution, bifurcation analysis and branch switching using pseudo-arclength continuation method. In particular, the shell response under harmonic excitation has been studied using an analysis in two steps: (i) first the excitation frequency has been fixed far enough from resonance and the magnitude of the excitation has been used as bifurcation parameter; the solution has been started at zero force where the solution was the trivial undisturbed configuration of the shell and has been continued up to reach the desired force magnitude and (ii) when the desired magnitude of excitation has been reached, the solution has been continued using the excitation frequency as bifurcation parameter.

## 4. Numerical results and discussion

Numerical calculations are performed for the clamped cylindrical panel with the complex base shown in Fig. 1(a), with the following dimensions and material properties: overall length $2b = 0.199$ m, curvilinear width $2a = 0.132$ m, length of the cut $2c = 0.041$ m, curvilinear width at the cut $2d = 0.092$ m, radius of curvature $R = 2$ m, thickness $h = 0.00028$ m, Young modulus $E = 195$ GPa, mass density $\rho = 7800$ kg/m$^3$ and Poisson's ratio $\nu = 0.3$. The same shell was considered in Ref. [33]. In all numerical simulations, a modal damping $\zeta_{1,1} = 0.004$ is assumed. If not differently specified, the calculations have been performed using the Donnell non-linear shell theory retaining in-plane inertia. Where specified, the results obtained using the Sanders-Koiter non-linear shell theory.

### 4.1. Linear analysis

The frequency range around the fundamental mode, i.e. the lowest natural frequency mode (1,1), is investigated. Here (1,1) indicates that the mode shape has one half-wave in both x and y directions, as shown in Table 1. The fundamental frequency $\omega_{1,1}$ obtained with the Donnell shell theory is 368.64 Hz, while $\omega_{1,3} = 453.44$ Hz, $\omega_{3,1} = 627.84$ Hz, $\omega_{1,5} = 689.51$ Hz and $\omega_{3,3} = 933.44$ Hz. Table 1 gives the natural frequencies of interest and corresponding mode shapes obtained with the Sanders-Koiter shell theory. The relative difference between natural frequencies obtained using the Donnell and the Sanders-Koiter shell theories is less than 0.01%, i.e. practically the same results for both theories.

In order to verify that our linear solution is practically converged, a comparison of the present results to a commercial FEM code was performed. The FEM results have been obtained after mesh refinement for convergence study using the ANSYS finite element code [42,43]. The panel was discretized using the standard three-dimensional non-layered structural solid element SOLID186 [44]. The element is defined by 20 nodes having three degrees of freedom per node: translations in the nodal x, y and z directions. The coordinate system for this element and its full description can be found in ANSYS Element Manual [44].

Natural frequencies obtained by the finite element method (FEM) are shown in Table 2 (for symmetric modes) and in Table 3. The convergence of the finite element model was studied by increasing the total number of elements in discretization until the relative difference reached 0.01%. Convergence was reached with a total number of 24,576 elements. Table 2 compares natural frequencies obtained using RFM models with the Donnell and the Sanders-Koiter shell theories to FEM results. The difference between natural frequencies obtained using the RFM and FEM is the range from 0.33% to 1.6% and becomes noticeable for higher frequencies.

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

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Amplitude of the vibration response of the panel versus the excitation frequency; mode (1,1), $f = 0.034$ and $\zeta_{1,1} = 0.004$; (a) Donnell theory. –, back-bone curve from Kurpa et al. [33]. (a) Maximum of the generalized coordinate $w_{1,1}$; (b) minimum of the generalized coordinate $w_{1,1}$ changed of sign.</figcaption></figure>
## 4.2. The non-linear vibrations

If not differently specified, the non-linear response of the panel to the harmonic force excitation $\tilde{f} = 0.034$ N applied at the center of the panel in z direction is investigated; this gives the advantage of having a symmetric excitation. A different amplitude $\tilde{f} = 0.02$ N has been used when specified.

### 4.2.1. Calculations with the Donnell shell theory

Fig. 2(a) and (b) presents the maximum and the minimum of the amplitude-frequency response of the present study model with 27 dofs, i.e. retaining only one term, mode (1,1), in the panel deflection w and 13 terms in expansions for each one of the in-plane displacements u and v in Eqs. (22); specifically the generalized coordinates used in the model are: $w_{1,1}, u_{1,1}, v_{1,1}, u_{1,3}, v_{1,3}, u_{3,1}, v_{3,1}, u_{1,5}, v_{1,5}, u_{3,3}, v_{3,3}, u_{3,5}, v_{3,5}, u_{1,7}, v_{1,7}, u_{5,1}, v_{5,1}, u_{3,7}, v_{3,7}, u_{5,3}, v_{5,3}, u_{5,5}, v_{5,5}, u_{1,9}, v_{1,9}, u_{7,1}, v_{7,1}$. Due to the curvature of the panel in the x direction, in a vibration period, the amplitude of the displacement inward the curvature (minimum displacement since it is in the negative z direction) is larger than the amplitude outward. Therefore, the panel in a vibration period is for a longer time deflected inward than outward.

In the previous study [33], a single-mode expansion was used for the transversal displacement w and the Donnell non-linear shallow-shell theory, which neglects in-plane inertia, was applied to condensate the in-plane displacements u and v into the Airy potential function. The potential function was directly obtained by solving a partial differential equation in w with Galerkin method. The present model with 27 degrees of freedom (dofs) includes one mode in the panel deflection w and thirteen modes for each in-plane displacement, so that it can be directly compared to Ref. [33]. However, it must be clarified that the present 27 dofs model gives slightly different results since in-plane inertia is retained. In Fig. 2, the forced vibration response obtained with the present model is very close to the backbone curve from Ref. [33]. In fact, the backbone curve of [33], indicating the free vibration response for different vibration amplitudes, approximately passes through the middle of the forced response curve (and its peak) computed in the present study. The softening type non-linearity, turning to hardening type for amplitudes of vibration much larger than the shell thickness, can be clearly observed in Fig. 2.

The convergence of the solution has been checked by calculating the amplitude-frequency response of the panel with larger models, including 33 dofs, 35 dofs and 37 dofs. The model with 37 dofs retains seven generalized coordinates in the panel deflection

<figure><img src="image_7.png" /><figcaption>Fig. 4. Amplitude-frequency response of the panel with respect to the different numbers of in-plane modes around the mode (1,1); Donnell theory; $f = 0.034$ and $\zeta_{1,1} = 0.004$. In all models 3 transversal modes are considered: (a) maximum of the generalized coordinate $w_{1,1}$; (b) minimum of the generalized coordinate $w_{1,1}$ changed of sign.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Amplitude of the response of the panel versus the excitation frequency; mode (1,1), $f = 0.034$ and $\zeta_{1,1} = 0.004$; Donnell theory. All models keep five transversal modes. (a) Maximum of the generalized coordinate $w_{1,1}$; (b) minimum of the generalized coordinate $w_{1,1}$ changed of sign.</figcaption></figure>
$w_{1,1}, w_{1,3}, w_{3,1}, w_{1,5}, w_{3,3}, w_{3,5}, w_{1,7}$ and 15 terms in the expansions for both the in-plane displacements $u$ and $v$; specifically $u_{1,1}, u_{1,3}, u_{3,1}, u_{1,5}, u_{3,3}, u_{3,5}, u_{1,7}, u_{5,1}, u_{3,7}, u_{5,3}, u_{5,5}, u_{1,9}, u_{7,1}, u_{9,1}, u_{3,9}; v_{1,1}, v_{1,3}, v_{3,1}, v_{1,5}, v_{3,3}, v_{3,5}, v_{1,7}, v_{5,1}, v_{3,7}, v_{5,3}, v_{5,5}, v_{1,9}, v_{7,1}, v_{9,1}, v_{3,9}$. The model with 35 dofs has five generalized coordinates in the panel deflection $w_{1,1}, w_{1,3}, w_{3,1}, w_{1,5}, w_{3,3}$ and the model with 33 dofs keeps three flexural modes in the expansion of deflection, $w_{1,1}, w_{1,3}, w_{3,1}$, respectively, in addition to the same in-plane coordinates. All these models exclude non-symmetric modes, as e.g. $w_{1,2}$ and $w_{2,1}$, from the analysis since the problem is symmetric. This is correct excluding the case of internal resonances with the non-symmetric modes. The frequency range of internal resonances can be reasonably predicted by linear

analysis. Natural frequencies of symmetric and non-symmetric modes are given in Table 3.

Comparison of the computed responses for both maximum and minimum of the generalized coordinate $w_{1,1}$ is given in Fig. 3(a) and (b), where the backbone curve from Ref. [33] is also shown. The response computed with 35 dofs model (with 5 deflection coordinates) is slightly moved to the right with respect to the one of the 33 dofs (3 deflection coordinates) model. In order to test if the convergence is reached, simulations with the model of 37 dofs (with 7 deflection coordinates) have been performed. The response of the 37 dofs models is very close to the one computed for the model with 5 deflection coordinates and is also presented in Fig. 3, showing solution convergence. On

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Amplitude-frequency responses of the panel around the mode (1,1); $f = 0.034$ and $\zeta_{1,1}=0.004$; Donnell theory. —, stable solution; -, unstable solution. (a) Maximum of the generalized coordinate $w_{1,3}$; (b) maximum of the generalized coordinate $w_{3,1}$; (c) maximum of the generalized coordinate $w_{1,5}$; (d) maximum of the generalized coordinate $w_{3,3}$; (e) maximum of the generalized coordinate $u_{1,1}$; (f) maximum of the generalized coordinate $v_{1,1}$.</figcaption></figure>
the other hand, Fig. 3 shows that the 27 dofs model, with a single deflection mode, is not accurate for vibration amplitude larger than about half the shell thickness.

Fig. 3(a) and (b) shows an initial softening behavior of the panel, turning to hardening for vibration amplitude around 0.8h.

While Fig. 3(a) shows the displacement outwards (positive), Fig. 3(b) shows the displacement inwards (negative). It is interesting to observe that the displacement outwards is slightly smaller than inwards, as observed in general for curved panels. However, Fig. 3(b) shows that the inward displacement presents

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Time-domain response of the panel at the excitation frequency $\omega=0.937\omega_{1,1}$; $f = 0.034$ and $\zeta_{1,1}=0.004$; Donnell theory; model with five transversal modes (35 dofs); (a) Force excitation; (b) generalized coordinate $w_{1,1}$; (c) generalized coordinate $w_{1,3}$; (d) generalized coordinate $w_{3,1}$; (e) generalized coordinate $w_{1,5}$; (f) generalized coordinate $w_{3,3}$; (g) generalized coordinate $u_{1,1}$; (h) generalized coordinate $v_{1,1}$.</figcaption></figure>
almost a discontinuity in the slope of the response at the transition from softening to hardening behavior, which is observed for vibration amplitude equal to 0.8h. This transition is instead smooth in Fig. 3(a) where the outward displacement is presented.

Fig. 3 shows also two small peaks for the 37 dofs model for $\omega/\omega_{1,1}=0.85$ and $1.04$. These peaks do not appear for the 33 dofs model, while for the 35 dofs model only the peak at 0.85 is observed. These two peaks can be easily attributed to 3:1 internal resonances of mode (1,1) with mode (3,3) for $\omega/\omega_{1,1}=0.85$ and with mode (1,7) for $\omega/\omega_{1,1}=1.04$. The natural frequencies in Table 2 divided by three times $\omega_{1,1}$ give exactly the peak values. The presence or absence of these two internal resonances in a model is also perfectly justified by the corresponding expansion: mode (3,3) is used in both the 35 and 37 dofs models while mode (1,7) is used only in the 37 dofs model.

Fig. 4(a) and (b) shows the effect of the number of dofs taken into account in the in-plane displacements on the panel response.

Here three bending modes are used in the panel deflection $w$ in all models; the model with 15 dofs has 6 terms for both the in-plane displacements $u$ and $v$; the model with 21 dofs keeps 9 terms in expansions of both $u$ and $v$; the model with 27 dofs has 12 terms of both $u$ and $v$; finally, the model with 33 dofs retains 15 terms in each one of the in-plane displacements. Results of the 33 dofs model are close to the ones of the 27 dofs model, while the two smaller models are significantly moved to the right. Therefore, Fig. 4 shows that 9 in-plane terms for both $u$ and $v$ are not enough to guarantee accurate results, while 12 terms give good results.

A similar analysis has been performed for the more refined models with 25 dofs and 35 dofs. In particular, the 35 dofs model has been previously introduced and retains 5 deflection terms on $w$ and 15 terms for each one of $u$ and $v$. The 25 dofs model retains the same 5 terms for $w$, but the number of in-plane terms is reduced to 10 for both $u$ and $v$. The maximum and minimum of the generalized coordinate $w_{1,1}$ in a vibration period are given in Fig. 5(a) and (b). Results show that 10 terms for the in-plane

<figure><img src="image_11.png" /><figcaption>Fig. 8. Phase plane diagram of the panel response; excitation frequency $\omega=0.937\omega_{1,1}$; $f=0.034$ and $\zeta_{1,1}=0.004$; Donnell theory; 35 dofs model. (a) Generalized coordinate $w_{1,1}$; (b) generalized coordinate $w_{1,3}$; (c) generalized coordinate $w_{3,1}$; (d) generalized coordinate $w_{3,3}$.</figcaption></figure>
displacements are not sufficient to give very accurate results, as shown by comparison with the model with 15 in-plane terms. This is in agreement with the result in Fig. 4.

The more complicated amplitude-frequency response of the models retaining 5 bending modes in Fig. 5 with respect to Fig. 4,

obtained with only 3 bending modes, is due to the presence of internal resonances that involves the additional bending modes $w_{1,5}$ and $w_{3,3}$, which are excluded in the models in Fig. 4. In fact, at the excitation frequency $\omega = 0.935\omega_{1,1}$, there is a two-to-one internal resonance due to the relationship $2\omega = \omega_{1,5}$. This can be

<figure><img src="image_12.png" /><figcaption>Fig. 9. Time-domain response of the panel; the excitation frequency $\omega=0.9\omega_{1,1}$; $f = 0.034$ and $\zeta_{1,1}=0.004$; model with 35 dofs; Donnell theory. (a) Force excitation; (b) generalized coordinate $w_{1,1}$; (c) generalized coordinate $w_{1,3}$; (d) generalized coordinate $w_{3,1}$; (e) generalized coordinate $w_{1,5}$; (f) generalized coordinate $w_{3,3}$; (g) generalized coordinate $u_{1,1}$; (h) generalized coordinate $v_{1,1}$.</figcaption></figure>
observed in Fig. 5(a) and (b) for $w_{1,1}$ and Fig. 6(a)–(f) for the other most significant generalized coordinates, with indication of response stability. In fact, the complication in the non-linear response of the model with 5 bending modes is due to a loop appearing around excitation frequency $0.94\omega_{1,1}$. This loop is more evident in Fig. 6(c) relative to the coordinate $w_{1,5}$. The difference between the 35 dofs model in Fig. 5 and the 33 dofs model in Fig. 4 is mostly limited to a small frequency range around the internal resonance. Another internal resonance is $2(1,03\omega_{1,3})=\omega_{3,3}$ in the frequency range shown in Fig. 6, but this does not affect the resonant generalized coordinate $w_{1,1}$.

Fig. 6(b) shows that the amplitude of the generalized coordinate $w_{3,1}$ is as big as the main coordinate $w_{1,1}$ around the resonance peak. The response of the coordinate $w_{3,1}$ grows very quickly after the response turns from softening to hardening behavior for vibration amplitude around 0.25h (amplitude similar to the one of $w_{1,3}$, that anyway is not growing so fast after the turning point), with a transfer of energy from $w_{1,1}$, that in fact

presents the sudden change of slope already observed in Fig. 3(b), to $w_{3,1}$. The sudden increase of the response $w_{3,1}$ is attributed to a 1:2 internal resonance happening at excitation frequency $\omega = 0.85\omega_{1,1}$ according to linear frequencies. Anyway this internal resonance relationship is moved to different frequency by the system non-linearity. In fact, internal resonances must be checked on the non-linear frequency-amplitude relationship so they do not occur at the exact linear frequency relationship [45,46]. An accurate way of identifying internal resonances is the construction of frequency energy plots, as discussed by Lee et al. [47] and Vakakis et al. [48]; this is beyond the scope of the present study. Fig. 6(b) indicates that this internal resonance is achieved around $0.9\omega_{1,1}$ in the present case.

Responses given in Fig. 6(a)–(f) have been computed using two large models with 33 dofs and 35 dofs, with 3 and 5 bending modes, respectively. Internal resonances could be further complicated by introducing larger models retaining non-symmetric modes, e.g. modes (2,4) and (4,2).

<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Phase plane diagram of the panel response; excitation frequency $\omega=0.9\omega_{1,1}$; $f = 0.034$ and $\zeta_{1,1}=0.004$; Donnell theory; 35 dofs model. (a) Generalized coordinate $w_{1,1}$; (b) generalized coordinate $w_{1,3}$; (c) generalized coordinate $w_{3,1}$; (d) generalized coordinate $w_{1,5}$; (e) generalized coordinate $w_{3,3}$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Amplitude of the response of the panel versus excitation frequency; fundamental mode (1,1), *f* = 0.034 and ζ<sub>1,1</sub> = 0.004; model with 35 dofs. —, Donnell theory; –, Sanders–Koiter theory. (a) Maximum of the generalized coordinate w<sub>1,1</sub>; (b) minimum of the generalized coordinate w<sub>1,1</sub>.</figcaption></figure>

It is important to analyze the results in the time domain. Fig. 7 (a)–(h) shows the time responses of the most significant (i.e. those with larger amplitude) generalized coordinates for excitation frequency ω₁ = 0.937ω₁, with initial conditions to catch the stable solution after the peak of the response, next to the loop in Fig. 5. Fig. 3 (zoomed part) depicts the exact location of ω₁. Results of the generalized coordinates related to radial displacement *w* clearly indicate that there is a strong asymmetry between vibration inwards and outwards. In fact, the vibration is much larger inwards (negative) than outwards (positive). Such effect, which is due to the curvature of the panel, has been previously observed for shallow cylindrical shells with simple shapes in Refs. [13,18–20] and for shallow shells with a complex base in [33]. In particular, the generalized coordinate w<sub>1,1</sub> has the largest amplitude; the inward displacement of w<sub>1,1</sub> is about 30% larger than the outward displacement in this case. A significant contribution of higher harmonics can be also observed in Fig. 7. In particular, the generalized coordinates w<sub>1,3</sub>, w<sub>1,5</sub> and w<sub>3,3</sub> are more affected by higher harmonics, while modes w<sub>1,3</sub> and w<sub>1,1</sub> (the last one being the resonant mode, with u<sub>1,1</sub> and v<sub>1,1</sub>) are almost perfectly sinusoidal. In Fig. 8 phase plane diagrams are provided for better understanding the asymmetry of inward and outward vibration; the effect of higher harmonics is indicated by loops.

The same analysis has been performed at the excitation frequency
ω₂ = 0.9ω₁,¹. The exact location of ω₂ is shown in the zoomed part
of Fig. 3. Appropriate initial conditions have been chosen in order to
catch the highest amplitude solution. Results are given in Figs. 9
and 10 for time responses and phase-plane plots that show an even
more significant contribution of higher harmonics. All the 5 general-
ized coordinates in the expansion of w give an important contribution
to the shell deflection, as plotted in Fig. 9. This shows that for large-
amplitude vibrations, an interaction among the different generalized
coordinates involving w is important and this is the reason for
inaccurate results obtained using a single mode in the expansion of
w. Fig. 9 shows also that the generalized coordinates w₃,₁, w₁,₅ and
w₃,₃ present an important contribution of super-harmonics in the
shell response. Higher harmonics give loops in the phase-plane plots
in Figs. 8 and 10. Instead, the main coordinate w₁,₁ shows a significant
displacement of the mean value (continuous component) in Fig. 9(b),
indicating that the displacement inwards is larger than the displace-
ment outward by about 20% for this coordinate.

4.2.2. Calculations with the Sanders–Koiter shell theory

Fig. 11 shows the comparison of the panel response in the
frequency neighborhood of the fundamental mode obtained using
the Sanders–Koiter and the Donnell non-linear shell theories.

Results have been computed for the refined models with 35 dofs.
Fig. 11(a) and (b) gives the maximum and the minimum of the
generalized coordinate w₁,₁, respectively, and show practically
coincident results. This indicates that, for the thin shallow shell
considered in the present study, there is no advantage in using a
more refined as the non-linear Sanders–Koiter shell theory; the
use of the Donnell non-linear shell theory retaining in-plane
inertia saves computational time in the present case. Similar
results comparing the Sanders–Koiter, the Donnell and several
other non-linear shell theories for shallow shells with simple base
has been observed and discussed by Amabili in Refs. [4,18–20].

**5. Conclusions**

An effective approach to build meshfree discrete models to
study forced large-amplitude vibrations of circular cylindrical
panels with arbitrary shapes is considered in the present study.
It consists in expanding the displacements using approximating
functions and the R-functions theory in order to satisfy the
boundary conditions on a complex domain. A two-step approach
is implemented in order to solve the problem: first a linear
analysis is conducted to identify natural frequencies and corre-
sponding natural modes to be used in the second step as a basis
for expanding of the non-linear displacements and building a
reduced-order model. Lagrange variational approach is applied to
obtain a system of ordinary differential equations on both steps.

The present approach has the advantage of being suitable for
constructing the system of admissible basis functions in an
analytical way in domains of arbitrary shape. The use of the
R-functions theory makes it also very flexible to be applied to a
variety of different boundary conditions.

Once the reduced-order model has been built, this has been
numerically studied using the pseudo-arclength continuation
method and bifurcation analysis. The convergence of the present
solution by increasing the number of degrees of freedom has been
studied and results from two different non-linear shell theories
have been compared.

Acknowledgments

The authors acknowledge the financial support of the NSERC Discovery Grant, Canada Research Chair and Canada Foundation for Innovation programs of Canada and the PSR-SIIRI Program of Québec.
$f(\eta_1, \dots, \eta_n) \ge 0$ describes domain $S$, where $\Omega=f(\eta_1, \eta_2, \dots, \eta_n)$ is an R-function that corresponds to the Boolean function $F(S_1, \dots, S_n)$. To construct the function $\Omega(x,y)$, it is sufficient to perform a formal substitution of the characteristic function $S_i$ with the continuous function $\eta_i$ and the Boolean operations $S_i \cap S_j$, $S_i \cup S_j$, $\bar{S}_i$ (intersection, union, absolute complement) with the corresponding symbols of R-operations, $\eta_i \wedge \alpha \eta_j$, $\eta_i \vee \alpha \eta_j$, $-\eta_i$. Following this substitution, the continuous function $\Omega$ of the domain $S$ is given by

$$ \Omega = f(\eta_1, \eta_2, \dots, \eta_n). \qquad (A.8) $$

## References

[1] Ya.M. Grigorenko, V.I. Guljaev, Nonlinear problems of the shells theory and methods of their solution (review), Prikladnaya Mekhanika 27 (1991) 3-23 (in Russian).

[2] V.D. Kubenko, P.S. Koval'chuk, Nonlinear problems of the vibration of thin shells (review), International Applied Mechanics 34 (1998) 703-728.

[3] M. Amabili, M.P. Païdoussis, Review of studies on geometrically nonlinear vibrations and dynamics of circular cylindrical shells and panels, with and without fluid-structure interaction, Applied Mechanics Reviews 56 (2003) 349-381.

[4] M. Amabili, Nonlinear Vibrations and Stability of Shells and Panels, Cambridge University Press, New York, USA, 2008.

[5] E. Reissner, Nonlinear Effects in Vibrations of Cylindrical Shells, Ramo-Woolridge Corporation Report AM5-6, 1955.

[6] E.I. Grigolyuk, Vibrations of circular cylindrical panels subjected to finite deflections, Prikladnaya Matematika i Mekhanika 19 (1955) 376-382 in Russian.

[7] B.E. Cummings, Large-amplitude vibration and response of curved panels, AIAA Journal 2 (1964) 709-716.

[8] A.W. Leissa, A.S. Kadi, Curvature effects on shallow shell vibrations, Journal of Sound and Vibration 16 (1971) 173-187.

[9] A.S. Volmir, Nonlinear Dynamics of Plates and Shells, Nauka, Moscow, Russia, 1972 (in Russian).

[10] A.S. Vol'mir, A.A. Logvinskaya, V.V. Rogalevich, Nonlinear natural vibrations of rectangular plates and cylindrical panels, Soviet Physics—Doklady 17 (1973) 720-721.

[11] C.Y. Chia, Nonlinear analysis of doubly curved symmetrically laminated shallow shells with rectangular platform, Ingenieur-Archiv 58 (1988) 252-264.

[12] A.A. Popov, J.M.T. Thompson, J.G.A. Croll, Bifurcation analyses in the parametrically excited vibrations of cylindrical panels, Nonlinear Dynamics 17 (1998) 205-225.

[13] Y. Kobayashi, A.W. Leissa, Large amplitude free vibration of thick shallow shells supported by shear diaphragms, International Journal of Non-Linear Mechanics 30 (1995) 57-66.

[14] J.H. Ginsberg, Large amplitude forced vibrations of simply supported thin cylindrical shells, Journal of Applied Mechanics 40 (1973) 471-477.

[15] M. Ganapathi, T.K. Varadan, Large amplitude vibrations of circular cylindrical shells, Journal of Sound and Vibration 192 (1996) 1-14.

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

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

[18] M. Amabili, Comparison of shell theories for large-amplitude vibrations of circular cylindrical shells: Lagrangian approach, Journal of Sound and Vibration 264 (2003) 1091-1125.

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

[20] M. Amabili, Non-linear vibrations of doubly curved shallow shells, International Journal of Non-Linear Mechanics 40 (2005) 683-710.

[21] Ye. Kurylov, M. Amabili, Polynomial versus trigonometric expansions for nonlinear vibrations of circular cylindrical shells with different boundary conditions, Journal of Sound and Vibration 329 (2010) 1435-1449.

[22] T. Yamaguchi, K. Nagai, Chaotic vibration of a cylindrical shell-panel with an in-plane elastic-support at the boundary, Nonlinear Dynamics 13 (1997) 259-277.

[23] C. Sansour, P. Wriggers, J. Sansour, Nonlinear dynamics of shells: theory, finite element formulation, and integration schemes, Nonlinear Dynamics 13 (1997) 279-305.

[24] K.Y. Lam, Wu Qian, Free vibration of symmetric angle-ply thick laminated composite cylindrical shells, Composites Part B: Engineering 31 (2000) 345-354.

[25] C. Sansour, W. Wagner, P. Wriggers, J. Sansour, An energy-momentum integration scheme and enhanced strain finite elements for the non-linear dynamics of shells, International Journal of Non-Linear Mechanics 37 (2002) 951-966.

[26] F. Boumediene, A. Miloudi, J.M. Cadou, L. Duigou, E.H. Boutyour, Nonlinear forced vibration of damped plates by an asymptotic numerical method, Computers and Structures 87 (2009) 1508-1515.

[27] V.L. Rvachev, Theory of R-Functions and Some of Its Applications, Nauka, Dumka, Kiev, Ukraine, 1982 (in Russian).

[28] G.R. Liu, Mesh Free Methods, second ed., CRC Press, USA, 2009.

[29] V.L. Rvachev, L.V. Kurpa, The R-Functions in Problems of Plate Theory, Nauka, Dumka, Kiev, Ukraine, 1987 (in Russian).

[30] L.V. Kurpa, V.L. Rvachev, E. Ventsel, The R-function method for the free vibration analysis of thin orthotropic plates of arbitrary shape, Journal of Sound and Vibration 261 (2003) 109-122.

[31] L.V. Kurpa, G. Pilgun, E. Ventsel, Application of the R-function method to nonlinear vibrations of thin plates of arbitrary shape, Journal of Sound and Vibration 284 (2005) 379-392.

[32] L.V. Kurpa, K.I. Lyubitska, A.V. Shmatko, Solution of vibration problems for shallow shells of arbitrary form by the R-function method, Journal of Sound and Vibration 279 (2005) 1071-1084.

[33] L. Kurpa, G. Pilgun, M. Amabili, Nonlinear vibrations of shallow shells with complex boundary: R-functions method and experiments, Journal of Sound and Vibration 306 (2007) 580-600.

[34] L. Kurpa, T. Shmatko, G. Timchenko, Free vibration analysis of laminated shallow shells with complex shape using the R-functions method, Composite Structures 93 (2010) 225-233.

[35] I.D. Breslavsky, K.V. Avramov, Nonlinear modes of cylindrical panels with complex boundaries. R-function method, Meccanica 46 (2011) 817-832.

[36] D. Zhou, Y.K. Cheung, S.H. Lo, F.T.K. Au, 3D vibration analysis of solid and hollow circular cylinders via Chebyshev-Ritz method, Computer Methods in Applied Mechanics and Engineering 192 (2003) 1575-1589.

[37] J.L. Sanders Jr., Nonlinear theories for thin shells, Quarterly of Applied Mathematics 21 (1963) 21-36.

[38] W.T. Koiter, On the nonlinear theory of thin elastic shells. I, II, III, Proceedings Koninklijke Nederlandse Akademie van Wetenschappen B 69 (1966) 1-54.

[39] F. Pellicano, Vibrations of circular cylindrical shells: theory and experiments, Journal of Sound and Vibration 303 (2007) 154-170.

[40] Y. Kurylov, M. Amabili, Nonlinear vibrations of clamped-free circular cylindrical shells, Journal of Sound and Vibration 330 (2011) 5363-5381.

[41] S. Wolfram, The Mathematica Book, fifth edition, Wolfram Media, USA, 2003.

[42] E.J. Doedel, A.R. Champneys, T.F. Fairgrieve, Y.A. Kuznetsov, B. Sandstede, X. Wang, AUTO 97:Continuation and Bifurcation Software for Ordinary Differential Equations (with HomCont), Concordia University, Montreal, Canada, 1998.

[43] G. Nygus, Numerical Analysis Using Finite Element Method, Ph.D. Thesis, NTU Mech. Eng. Dept., Lagos, 1983.

[44] ANSYS, Engineering Analysis System Theoretical Manual, Ansys Version 9.0. <http://www.Ansys.com>, 2004.

[45] R. Lewandowski, Computational formulation for periodic vibration of geometrically nonlinear structures—part 2: numerical strategy and examples, International Journal of Solids and Structures 34 (1997) 1949-1964.

[46] G. Kerschen, M. Peeters, J.C. Golinval, A.F. Vakakis, Nonlinear normal modes, part I: a useful framework for the structural dynamicist, Mechanical Systems and Signal Processing 23 (2009) 170-194.

[47] Y.S. Lee, G. Kerschen, A.F. Vakakis, P. Panagopoulos, L.A. Bergman, D.M. McFarland, Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment, Physica D 204 (2005) 41-69.

[48] A.F. Vakakis, O.V. Gendelman, L.A. Bergman, D.M. McFarland, G. Kerschen, Y.S. Lee, Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems I, Springer, New York, 2008.

[49] V.L. Rvachev, Boolean Algebra Methods in Mathematical Physics, Nauka, Dumka, Kiev, Ukraine, 1974 (in Russian).

[50] V.L. Rvachev, T.I. Sheiko, R-functions in boundary value problems in mechanics, Applied Mechanics Reviews 48 (4) (1995) 151-187.

[51] L.V. Kantorovich, V.I. Krylov, Approximate Methods of Higher Analysis, Interscience Publisher, New York, 1958.