Thin shell structures have an outstanding efficiency in fully utilizing the structural material and have been extensively used in many engineering applications including aircraft structures, pressure vessels, and others. Static and dynamic analysis is essential for the analysis and design of shell structures. Various numerical methods, such as the Finite Difference Method (FDM) and especially the Finite Element Method (FEM) have been used (Lee and Han, 2001) for the dynamic analysis of linear elastic thin shells characterized by complex geometry, loading and boundary conditions. Both methods have been employed successfully for the solution of a variety of static and dynamic shell problems. The Boundary Element Method (BEM) is an efficient alternative to the domain type methods, especially for thin elastic shallow shells (Beskos, 1991), or combined with the AEM for cylindrical shells (Yiotis and Katsikadelis, 2000).

Such methods require the generation of a mesh which can be an incredibly tedious and time-consuming process, while their convergence rate is of 2nd order (Cheng et al., 2003). On the other hand, Meshless Methods (MMs) present an attractive alternative to FEM or BEM, especially for shell structures that are complex regarding both the governing equations and the geometry representation. Comprehensive descriptions of different MMs are presented by Liu (2002); Liu and Gu (2005) and in a review paper by Nguyen et al. (2008).

There are several papers on dynamic analysis of shells using MM. Homogeneous shells are studied using various versions of MM in Liu et al. (2006), Ferreira et al. (2006b), and Dinis et al. (2011). Functionally graded (FG) cylindrical thin shells have also been treated by this method (Ferreira et al., 2006a; Zhao et al., 2009; Roque et al., 2010), as well as thick cylindrical shells (Pilafkan et al., 2013) have been analyzed by this method.

The mesh-free multiquadric radial basis functions (MQ-RBFs) method presented in Kansa (2005) has attracted the interest of researchers, due to its exponential convergence and its easiness of implementation. The significant drawbacks of the method are the ill-conditioning of the coefficient matrix and the inability to accurately approximate the derivatives of the sought solution which renders the method inappropriate for a strong formulation of the problem. These drawbacks of the standard MQ-RBF method, are overcome by a new RBF method presented recently by Katsikadelis (2006, 2008a,b, 2009) and Yiotis and Katsikadelis (2008, 2013, 2015a,b). Another critical issue is the implementation of multiple boundary conditions for equations of order higher than 2nd. In this investigation, the δ-technique is employed (Jang et al., 1989) for the 4th order equation. The problem of multiple boundary conditions is not present when the shell is modeled as a 3D body (Katsikadelis and Platanidi, 2007).

In this paper, the MAEM is extended to the dynamic problem of cylindrical shell panels as described by section MAEM Solution. A first approach to this problem was attempted in a previous work (Yiotis and Katsikadelis, 2015b), where some preliminary results only for the eigenfrequency analysis were presented. In section Problem Statement, the statement of the problem is presented, while several example problems are worked out in section Numerical Examples, which illustrate the applicability of the method and demonstrate its efficiency and section Conclusions contains certain conclusions drawn from this investigation.

We consider a thin cylindrical shell with parametric lines x (s = const.) and s (x = const.) which are assumed to be lines of curvature, as well; x is measured along the x lines of the shell and s along the s lines, while z is measured along the normal to the middle surface of the shell, as shown in Figure 1. R is the radius of curvature and h is the thickness.

In this investigation we use the Flügge equations for the thin shell theory, based on the following assumptions (Love, 1944):

The first assumption defines the meaning of “thin shells,” whereas the second one implies that all calculations refer to the original undeformed configuration and subsequently leads to linear differential equations. Further, the assumption z/R < < 1 is adopted in deriving the stress resultants in integrating the stresses through the thickness of the shell. The 4th assumption is known as Kirchhoff's hypothesis yielding

which implies σ_xz = σ_sz = 0 (Leissa, 1973).

The equations of motion for the case of the thin cylindrical shell can be derived using Hamilton's principle as follows

where Π is the total potential energy given by

in which U₀, is the strain energy

and W₁, W₂ the works produced by the loading and the boundary forces, i.e.,

In Equation (4) σ_x, σ_s are the normal stresses, σ_xs the tangential shear stress, σ_xz, σ_sz the transverse (in the z direction) shear stresses and e_x, e_s, γ_xs, γ_xz, γ_sz the respective strains at an arbitrary point of shell cross-section.

In Equation (5), q_x(t), q_s(t), and q_z(t) are the three components of the loading in the axial, circumferential and normal to the middle surface directions, respectively, while u, ν, and w represent the axial, circumferential and normal displacements at the middle surface of the shell.

In Equation (6) the quantities N̄x, N̄xs, Q̄x, M̄x, M̄xs, and N̄sx, N̄s, Q̄s, M̄sx, M̄s denote prescribed boundary forces along an edge (x = const.) and an edge (s = const.) respectively; u, ν and w represent the axial, circumferential and normal displacements at the boundary and θ_x, θ_s are the rotations of the normal to the middle surface about the s and x axes respectively.

Furthermore, in Equation (2) the quantity K is the kinetic energy of the body and is given regarding the shell variables as

where ρ is the mass density of the material of the shell.

The quantity δW_nc represents the work of the damping forces, non-conservative forces, due to the virtual displacements and is given by the relation

where η is the damping coefficient.

Neglecting the contribution from the rotatory inertia terms ρh3θx/12 and ρh3θs/12 in Equation (7), inserting Equation (3) and taking the variation (Katsikadelis, 2016), we obtain the Flügge type differential equations (Flügge, 1962; Kraus, 1967), in terms of the displacements as well as the associated boundary and initial conditions

where g_i(x), h_i(x) (i = 1, 2, 3) are specified functions and (x) = (x, s).

The stress resultants N_x, N_s, N_xs, N_sx, M_x, M_s, M_xs, M_sx, Q_x, Q_s are expressed in terms of the displacements as

where D = E/12(1−ν²) and, T_xs V_x the effective tangential membrane and transverse shear forces at the edges x = 0, l given as

Similarly, T_sx and V_s represent the effective tangential membrane, and transverse shear force at the edges s = 0, a and are given as

MAEM (Katsikadelis, 2002) is used for the solution of the initial boundary problem (9), (10), (11), as shown in the following. Let u, ν and w be the solution to the problem. Since Equations (9) are of the 2nd order with respect to u, ν and of the 4th order with respect to w, the analog equations which are convenient to use are

where b_i = b_i(x, t) (i = 1, 2, 3) are unknown fictitious sources depending on time, which, however, is treated as a parameter, i.e., Equations (15a,b,c) are quasi-static, treated as static at each instant. The fictitious sources can be established as follows.

The fictitious sources are approximated with MQ-RBFs series. Thus we have

where c is the shape parameter; K, L represent the number of collocating points inside Ω and on Γ, respectively; fj(r)=r2+c2, r = |x−x_j| (see Figure 2), and aj(1), aj(2), aj(3) time-dependent coefficients to be determined. Note that, while the derivatives of the membrane displacements u, ν are collocated at K domain and L boundary points, the derivatives of the normal displacement w according to the δ-technique (Ferreira et al., 2005) are collocated in K domain and 2L boundary nodal points placed on the auxiliary boundary Γ~ at a small distance δ from the actual one.

Equations (16) can be directly integrated to yield

where ûj(r),ν^j(r),ŵj(r) are solutions of the equations

Since the functions f_j(r) depend only on the radial distance r, the solution of Equations (18) can be obtained after writing them in polar coordinates. For the 2nd order equations, we have

which after integration gives

Similarly, we have

The regularity condition at r = 0 demands G₁ = G₂ = 0. The remaining constants H₁, H₂ together with the shape parameter c, if not arbitrarily specified, can be determined with an optimization procedure, such as to ensure the regularity of coefficients matrix (control of the condition number) and the error minimization. It has been shown that the coefficient matrix resulting from the new RBFs is always invertible (Sarra, 2006), and as a result, we take in this analysis H₁ = H₂ = 0 for convenience. Thus only c, the shape parameter is involved in the error minimization procedure (Katsikadelis, 2008a).

For the 4th order equation, one can write

Integrating Equation (22), and removing the singular terms and the terms including the arbitrary constants (Yao et al., 2010) yield

Direct differentiation of Equations (17) obtains the derivatives of the displacements involved in the governing equations (9a,b,c).

where i, k, l stand for x, s.

Furthermore, the derivatives of the displacements with respect to time can also be obtained by direct differentiation of Equations (17). Thus we have

Collocating Equations (9) at the K nodal points inside Ω and the four boundary conditions, Equations (10), at the L boundary nodal points (Figure 2) using the well-known δ-technique for multiple boundaries (Yiotis and Katsikadelis, 2015a), and inserting Equations (17) and (24) to (26) in the resulting expressions, a system of ordinary differential equations is obtained, namely

where M, C, and K are known square matrices having dimension 3K+4L; g is a vector including the 3K values of the external load g(x, t) and a is the vector of the 3K+4L values of the unknown time-dependent coefficients aj(1)(t), aj(2)(t), aj(3)(t).

Equation (27) is the semi-discretized equation of motion of the cylindrical shell with M, C, and K representing the generalized mass, damping and stiffness matrices, respectively. It can be solved numerically, using any time step integration technique to establish the time-dependent unknown coefficients. Here the method presented in Katsikadelis (2014a,b) is employed. The initial conditions of Equation (27) result from Equations (17) and (25) on the basis of Equations (11) as follows:

Once the coefficients aj(1)(t), aj(2)(t), aj(3)(t) have been computed, the field functions u, ν, w, their derivatives, and the stress resultants can be evaluated from Equations (17), (24) to (26) and (12) to (14).

For free vibrations it is C = g(x, t) = 0 and the equation of motion, Equation (27), takes the form

while the essential boundary conditions, Equations (10), become homogeneous.

By setting

Equation (29) results in the eigenvalue problem

which gives the eigenfrequencies ω_i and the eigenvectors α = [α⁽¹⁾, α⁽²⁾, α⁽³⁾]T, where α(1)=[α1(1),α2(1),…,αK+L(1)]T, α(2)=[α1(2),α2(2),…,αK+L(2)]T, α(3)=[α1(3),α2(3),…,αK+2L(3)]T. The elements of these vectors are the three sets of coefficients corresponding to the functions u, ν, and w, respectively. Subsequently, the mode shapes are obtained by substituting α = [α⁽¹⁾, α⁽²⁾, α⁽³⁾]T in Equations (17).

The accuracy of the approximation (9) depends heavily on c (see Equations 20-21-23). This was also verified in the problem at hand. Thus we come across to the problem of selecting a “good” value for c, that is, a value of the shape parameters that produces results of acceptable accuracy. Several methods have been suggested (Hardy, 1971; Franke, 1982; Foley, 1994; Rippa, 1999; Katsikadelis, 2009) for selecting a good value for c in 2D problems. Katsikadelis (2006, 2008b) proposed the minimization of the functional (total potential) for obtaining an optimal value for c. For the present problem, the optimal value is obtained by the search method as the value of c which yields the minimum value of the eigenfrequencies ω_i. It was observed that the optimum c giving the minimum 1st eigenfrequency differs negligibly from that yielding the higher minimum eigenfrequencies. Therefore the same optimum value of c can be used to avoid the search method for higher eigenfrequencies.

On the basis of the above analysis, a Fortran program has been written. The expressions of the derivatives involved in Equations (9) to (11) and Equations (12) to (14) have been obtained using the symbolic language MAPLE. Though the method applies to cylindrical shell of variable radius of curvature, for reasons of simplicity, the efficiency and accuracy of the developed method are demonstrated by studying the forced and free vibrations of circular cylindrical panels, (Figure 3), under different sets of boundary conditions. The NASTRAN FEM code and a model with 400 rectangular elements (Figure 4B) are used to compare the results. In all examples the employed material constants are: E = 21 × 10⁶kN/m², ν = 0.30. The results have been obtained running the relevant programs on an Intel Core 1.6 GHz with RAM 4 GB computer.

The Meshless Analog Equation Method, a truly meshless method, has been applied to the dynamic analysis of thin cylindrical shell panels in the present study. MAEM is based on the principle of the analog equation, converting the original equations into three substitute equations, two Poisson's and one biharmonic, which are solved using a meshless method. The use of integrated MQ-RBFs to approximate the fictitious sources allows the approximation of the sought solutions by new RBFs, which approximate both the solution and its derivatives accurately. This way the strong formulation of the problem avoids the drawbacks inherent in the conventional MQ-RBFs, while maintaining all the advantages of a truly mess-free method. A method is presented to obtain optimum values for the shape parameter, eliminating the uncertainty in its choice. It was observed that the optimum value of the shape parameter for the 1st mode differs negligibly from those of higher modes and therefore the same value can be used to obtain the eigenfrequencies of higher modes. The solution algorithm is straightforward and quite reasonably easy to program. The numerical examples presented demonstrate the efficiency and accuracy of the proposed method and show that MAEM can be used as an efficient solver for challenging problems in engineering analysis.

AY has implemented the MAEM (Meshless Analog Equation Method) for the solution of problems governing by fourth order differential equations (plates-shells etc.). JK has invented the AEM (Analog Equation Method), which is a general method for the solution of various engineering problems.