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

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 build 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: ﬁrst 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.


Introduction
Circular cylindrical shells (or panels) are largely used in aeronautics, 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][2][3][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][6][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][9][10][11][12][13].More refined classical non-linear shell theories have been applied to study large-amplitude vibrations of cylindrical shells in Refs.[1,4,[14][15][16][17][18][19][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][18][19][20][21][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 nonlinear dynamics of circular cylindrical panels or double curved shallow shells with simple geometric shapes; e.g. with a rectangular base [5][6][7][8][9][10][11][12][13][17][18][19][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][24][25][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][30][31] and shallow shells [32][33][34][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 nonlinear 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 nonlinear 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.

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.

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 e x , e y and g xy at an arbitrary point of the panel are related to the middle surface strains e x,0 , e y,0 and g 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: 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] 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] e x;0 ¼ The elastic strain energy U s of a circular cylindrical shell, within the limits of Kirchhoff-Love assumptions, is given by [4,10] where S is the shell middle surface; the stresses s x , s y and t xy are related to the strains for homogeneous and isotropic material (s z ¼0, case of plane stress) by where E is the Young modulus and n is Poisson's ratio.using Eqs.( 1), ( 4) and ( 5), the following expression is obtained: x;0 þe 2 y;0 þ2ne x;0 e y;0 þ 1Àn 2 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.

Kinetic energy, virtual work and damping
The kinetic energy T s of the panel, by neglecting rotary inertia, is given by where r 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 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 f and is given by [4] where o is the excitation frequency, t is the time, d is the Dirac delta function, f is the amplitude of the force applied along the normal to the panel middle surface and taken positive in z direction, x and ỹ 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: The non-conservative damping forces are assumed to be of viscous type and are taken into account using the Rayleigh dissipation function [4] where c is the viscous damping coefficient; it has a different value for each term of the mode expansion.

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 straindisplacement 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.
The studied shell is considered to be clamped, which means 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 t 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] Uðx,yÞ ¼ Oðx,yÞP 1 ðx,yÞ, Vðx,yÞ ¼ Oðx,yÞP 2 ðx,yÞ, where functions P i (i¼1,2,3) are functions expanded in the following way: q j j 2j ðx,yÞ, where q j , j ¼ 1,. ..,N 3 are unknown coefficients to be determined from the corresponding eigenvalue problem and {j ij }, (i¼1,2,3;j ¼1,y,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Þ ¼ Oðx,yÞ Vðx,yÞ ¼ Oðx,yÞ Wðx,yÞ ¼ O 2 ðx,yÞ where T 1i ( Á ) is the ith order Chebyshev polynomial of the first kind, u i,j , v i,j and w i,j are unknown coefficients, which have to be re-arranged into a single indexed sequence q j , j¼ 1,y,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: In Eqs. ( 15) and ( 17) O(x,y)¼0 is the function of the boundary of the panel normalized up to the first order.Normalization of the function O(x,y) is given by [27,29] Oðx,yÞ ¼ 0, 8ðx,yÞ A @S, @Oðx,yÞ @n ¼ À1, 8ðx,yÞ A @S, Oðx,yÞ40, 8ðx,yÞ A S: In the case of the cylindrical panel with the geometry shown in Fig. 1(a), the function of the boundary O(x,y) can be composed in the following form [4,33]: where the functions Z i , (i¼1,y,4), are given by The function O(x,y) for the panel shown in Fig. 1 is plotted in Refs.[4,33] in details.In Eq. ( 19) the R-conjunction (x4 0 y) and the R-disjunction (x3 0 y) have been used.They are defined as [27] x4 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 where o are the natural frequencies and q ¼ ðq 1 ,q 2 ,. ..,qN 3 Þ T are the corresponding eigenvectors giving the mode shapes.The components of the vector q correspond to the reordered sequences of the unknown coefficients in Eq. ( 17).K ¼ fk ij g i,j ¼ 1,...,N 3 and M ¼ fm ij g 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: It is important to point out that coefficients in both matrices K and M are double surface integrals over the domain with complex shape to be numerically computed with accuracy.Eq. ( 21) gives the natural frequencies o m , for m ¼1,y,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.

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 Table 1 Natural frequencies (Hz) and corresponding mode shapes of the panel.

Table 2
Comparison of natural frequencies (Hz) of symmetric modes of the circular cylindrical panel.

Methods
Modes (m,n) 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,y,M 3 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 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] The Lagrange equations of motion are where qT/qp 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 where coefficients z can be obtained only numerically and the linear terms zk,i in Eq. ( 26) are equal to zero for any iak (for 0oirM 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.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.

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 r ¼ 7800 kg=m 3 and Poisson's ratio n¼0.3.The same shell was considered in Ref. [33].In all numerical simulations, a modal damping B 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.

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  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.

The non-linear vibrations
If not differently specified, the non-linear response of the panel to the harmonic force excitation 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 f ¼ 0:02 N has been used when specified.

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: , 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 shallowshell 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.
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 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 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 o/o 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 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 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 o ¼ 0:935o 1,1 , there is a two-to-one internal resonance due to the relationship 2o ¼o 1,5 .This can be 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:94o 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,03o 1,3 )¼o 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 o ¼ 0:85o 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.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 o 1 ¼ 0:937o 1,1 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 o 1 .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][19][20] and for shallow shells with a complex base in [33].In particular, the generalized coordinate w 1,1 has the largest amplitude; the inward displacement of w 1,1 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 1,3 , w 1,5 and w 3,3 are more affected by higher harmonics, while modes w 1,3 and w 1,1 (the last one being the resonant mode, with u 1,1 and v 1,1 ) 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 o 2 ¼ 0:9o 1,1 .The exact location of o 2 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 generalized coordinates in the expansion of w give an important contribution to the shell deflection, as plotted in Fig. 9.This shows that for largeamplitude 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 3,1 , w 1,5 and w 3,3 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 1,1 shows a significant displacement of the mean value (continuous component) in Fig. 9(b), indicating that the displacement inwards is larger than the displacement outward by about 20% for this coordinate.

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 corresponding 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.

Appendix. R-functions
Let us consider a boundary value problem of continuum mechanics in the form: where A is a differential operator, L i are boundary operators and u is the solution field function in the domain S. In Eq. (A.1b) qS i are the boundary sections (or the whole boundary); f is the field exciter [48,49].
In the 1950s Kantorovich [51] introduced a numerical technique based on the observation that the solution of a differential equations Eq. (A.1a) with homogeneous boundary conditions, known as the Dirichlet problem u9 @S ¼ 0 ðA:2aÞ can be represented in the form u ¼ OP: ðA:2bÞ In the context of structural analysis, the solution field function u is the displacement, and Eq.(A.2a) represents a boundary that has zero displacements.In Eqs.(A.2a) and (A.2b) qS is the boundary of the domain S; the function O vanishes at the boundary (O9 qS ¼0) and is positive in the interior of S, and P is an unknown function that allows us to satisfy the differential equations governing the problem.Because O is identically zero at the boundary qS, then no matter what the indefinite component P (differentiable up to the required order) is, the function u given by Eq. (A.2b) will satisfy the boundary condition (A.2a) exactly.The function P is determined in order to satisfy the differential equations; therefore, in general, it can be expressed in approximate way using a series expansion of basis functions c i where C i are scalar coefficients.Thus, the solution takes the form The undetermined coefficients C i can be found numerically, for example, using Rayleigh-Ritz or Galerkin method.The function P originally was represented in a global polynomial basis, but can be also constructed using many other basis functions, such as B-splines, radial basis functions, or trigonometric polynomials, among others.In spite of the intrinsic advantage of this method, which was able to give a clean and modular separation of the geometric information represented by the function O from the differential equation and numerical procedure used to determine the analytic component P, Kantorovich's idea did not have particular success.In fact, (i) no available technique existed to construct such real functions O for domains of complex shape, and (ii) the same solution was not applicable to other types of boundary value problems.Rvachev found the way to overcome both of the obstacles by creating the R-functions theory.He recognized Kantorovich's representation of the field as a special form of the Taylor series expansion of the function O and showed that the solution structure can be practically generalized to any type of engineering problem [27,50].He also developed a theory of R-functions specifically as a method for constructing the functions O for arbitrary shapes [27,29,49].R-functions are defined as functions whose sign is completely defined by the sign of their arguments.Functions that satisfy these properties are, for example, xyz, x þ yþ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi xy þ x 2 þ y 2 p and xyþz þ9z À xy9.The best known system of R-functions is given by the R-conjunction (x4 a y) and the R-disjunction (x3 a y), which are defined as x þyÀ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi x 2 þ y 2 À2axy q , ðA:3Þ where a is a continuous function satisfying the condition À 1oa(x,y)r1; the denial function, simply given by a minus, must be added to complete this system of R-functions.R-functions behave as continuous analogs of logical Boolean functions.Every Boolean function has infinity analog R-functions.The R-conjunction (A.3) is an R-function whose companion Boolean function is logical ''and''(4), whereas Eq. (A.4) has companion Boolean function logical ''or'' (3).Note that the precise value of a is not important in many applications, and it is often set to a constant.If a¼1 is taken, Eqs.The domain and its characteristic function are usually denoted with the same symbol.
In general, the characteristic function of domain S is obtained by applying Boolean operations of intersection \, union [ and absolute complement -, respectively, to subdomains S i , where a characteristic function S i for each subdomain S i can be defined by Eq. (A.5), when S is substituted by S i .Because S i AB 2 {0,1}, that is, to the space of Boolean functions, then S i may be used as an argument of the Boolean function F S ¼ FðS 1 ,S 2 ,. ..,S n Þ, ðA:6Þ which is the characteristic function of domain S. It is obvious that domain S is determined not only by the shape of subdomains S i , but also by the type of Boolean functions involved in F (union, intersection and absolute complement).
As next step in the determination of O, the functions Z i are introduced as the continuous analog of the characteristic functions S i .They may be lines (surfaces) described by known equations and are defined as Z i ðx,yÞ 40 8ðx,yÞ A S i , Z i ðx,yÞ o0 8ðx,yÞ= 2S i , Z i ðx,yÞ ¼ 0 8ðx,yÞ A @S i : ðA:7Þ Let us assume that domain S is defined by the characteristic function (two-valued predicate) represented in Eq. (A.6), where F(S 1 ,y,S n ) is a known Boolean function.Then, the inequality f(Z 1 ,y,Z n ) Z0 describes domain S, where O¼f(Z 1 ,Z 2 ,y,Z n ) is an R-function that corresponds to the Boolean function F(S 1 ,y,S n ).
To construct the function O(x,y), it is sufficient to perform a formal substitution of the characteristic function S i with the continuous function Z i and the Boolean operations S i \ S j , S i [ S j , S i (intersection, union, absolute complement) with the corresponding symbols of R-operations, Z i 4 a Z j , Z i 3 a Z j , ÀZ i .Following this substitution, the continuous function O of the domain S is given by O ¼ f ðZ 1 ,Z 2 ,. ..,Z n Þ: ðA:8Þ

Fig. 1 .
Fig. 1.Circular cylindrical panel with complex base.(a) Coordinate system and dimensions; (b) the function of the boundary O(x,y).

Fig. 1 (
Fig. 1(b).The complete procedure of building the equation of the domain O(x,y)¼0 for the panel shown in Fig.1is addressed in

Fig. 4 .
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 z 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.
to 3:1 internal resonances of mode (1,1) with mode (3,3) for o/o 1,1 ¼ 0.85 and with mode (1,7) for o/o 1,1 ¼1.04.The natural frequencies in Table 2 divided by three times o 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.

Fig. 6 (
b) indicates that this internal resonance is achieved around 0:9o 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 terms, respectively.Internal resonances could be further complicated by introducing larger models retaining non-symmetric modes, e.g.modes(2,4) and (4,2).

4. 2 . 2 .
Calculations with the Sanders-Koiter shell theory Fig.11shows 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 1,1 , 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][19][20].
(A.3) and (A.4) become the functions Min(x, y) and Max(x, y), respectively.Setting a¼0 in (A.3) and (A.4) gives the simpler equations introduced in Section 3.1.R-functions are closed under composition.Therefore, the function O(x,y) can be obtained for domains of complicate shape whose boundary consists of sections of lines (surfaces) described by known equations of the kind Z i (x,y)¼0 or Z i (x,y,z)¼0 in the three-dimensional case, where Z i are elementary functions (or other functions when necessary).Using R-operations, such as R-disjunctions (x3 a y), which has analog Boolean union [, R-conjunction (x4 a y), with analog Boolean intersection \ and R-denial ( Àx), with analog Boolean absolute complement -(indicated with an overline), it is possible to construct the analytical expression of O(x,y) for almost any domain S with no concern of its connectivity.The first step in order to build O(x,y) for the domain S, it is necessary to construct the so-called characteristic function (also named two-valued predicate) of domain S, which is defined as S ¼ Sðx,yÞ ¼ 0, 8ðx,yÞ= 2S, 1, 8ðx,yÞ A S: ( ðA:5Þ

Table 3
Natural frequencies (Hz) of the circular cylindrical panel.