DYNAMIC STIFFNESS METHOD IN THE VIBRATION ANALYSIS OF CIRCULAR CYLINDRICAL SHELL

Ljuske se često koriste kao elementi inženjerskih konstrukcija, naročito u mašinstvu, građevinarstvu, avio i brodskom inženjeringu. Razlog je u tome što ljuske mogu da se koriste za konstrukcije velikih raspona, bez potrebe za međuosloncima, i što ovakve konstrukcije imaju povoljan odnos krutosti i težine. Ljuske se koriste za konstrukcije rezervoara, delova trupova aviona, delova konstrukcija brodova i slično. Tokom životnog veka, odnosno tokom njihove primene, ove konstrukcije izložene su složenim uslovima okruženja koji podrazumevaju najrazličitije granične uslove i opterećenja, kao što su npr. dinamički uticaji velikog intenziteta koji mogu da dovedu do kolapsa konstrukcije. Potpuno poznavanje dinamičkih karakteristika ovih konstruktivnih elemenata od velike je važnosti da bi se obezbedilo sigurno, uspešno i ekonomski isplativo projektovanje. Metod dinamičke krutosti (MDK), koji u poslednje vreme privlači pažnju sve većeg broja istraživača, predstavlja alternativu metodu konačnih elemenata (MKE) u analizi vibracija. Kao što je poznato, veličina konačnog elementa mora da bude manja od petine talasne dužine koja odgovara najvišoj frekvenciji u analizi da bi se dobili rezultati zadovoljavajuće tačnosti. Zato proračun konstrukcija primenom MKE u oblasti visokih frekvencija, koje se javljaju u akustici ili kod vibracija mašina visokih frekvencija, postaje glomazan i


INTRODUCTION
Shells are often used as elements of engineering structures, particularly in mechanical, civil, aerospace and naval engineering.The reason is that the shells can be used for a large-span structures, without intermediate supports, and because they have a favourable stiffness to weight ratio.Shells are used for the construction of reservoirs, pressure vessels, fuselage, naval vehicles, etc.During their lifetime and exploitation, these structures are subjected to various loadings under complex boundary conditions, such as violent dynamic load, that can lead to collapse of structures.Therefore, a complete knowledge of the dynamic characteristics of these structural elements is very important to ensure a safe, successful and economically feasible design.
The dynamic stiffness method (DSM), which has recently been attracting the attention of a growing number of researchers, is an alternative to the finite element method (FEM) in the vibration and buckling analysis.As known, the size of the finite elements must be less than a fifth of the wavelength corresponding to the highest frequency in the analysis, in order to obtain satisfactory results.Therefore, the analysis of structures by FEM in the high frequencies range, that occur in acoustics and vibrations of machine at high frequencies, becomes cumbersome and time consuming because of the need for condensing the finite element mesh.
Ako je debljina ljuske manja od 1/20 talasne dužine i/ili poluprečnika krivine, teorija tankih ljuski, gde su zanemarene smičuća deformacija i rotaciona inercija, generalno je prihvatljiva.U zavisnosti od pretpostavki učinjenih tokom izvođenja kinematičkih relacija, definisanja presečnih sila i uslova ravnoteže, formulisane su različite teorije tankih ljuski, čiji je pregled dat u knjizi od Leisse-а [12].U ovom radu je formulisana dinamička matrica krutosti kružne cilindrične ljuske po Flügge-ovoj teoriji.U sistemu parcijalnih diferencijalnih jednačina, kojima je definisan problem slobodnih vibracija, izvodi po vremenu su eliminisani korišćenjem spektralne dekompozicije, dok je tangencionalna koordinata eliminisana primenom rešenja u obliku Fourier-ovog reda.Na taj način je dobijen sistem od tri obične diferencijalne jednačine koje su funkcija samo The formulation of the dynamic stiffness (DS) matrix, which relates forces and displacements at the boundaries of the element, is based on the exact solutions of the free vibration problem in the frequency domain.Therefore, the stiffness matrices determined in this manner are accurate, frequency dependent, and contain among the stiffness, inertia and damping.Therefore, the elements based on these matrices are called continuous elements.By their usage, the discretization of observed domain is reduced to a minimum, i.e. one continuous element is sufficient to describe behaviour of a structural element of constant material and geometrical properties.The global stiffness matrix is obtained using the assembly procedure in the same way as the FEM, with one exception -in the DSM two dimensional elements, instead at nodes, are connected along the boundary lines.Finally, it can be concluded that the DSM uses advantages of analytical as well as numerical methods (FEM).
The dynamic stiffness matrix was first developed for one-dimensional elements.Based on the closed form solution of the free vibration of a straight beam, the dynamic stiffness matrix according to the Euler-Bernoulli and Timoshenko beam theory was formulated [1][2][3][4].Casimir [5] developed the dynamic stiffness matrix of a curved beam, whereby the differential equations were solved numerically.In the same way, Le Sourne [6] derived the dynamic stiffness matrix of a circular cylindrical shell.Continuous elements formulated in this way are called numerical continuous elements.
For a rectangular plate, the free vibration problem can be solved in a closed form only for special combinations of boundary conditions.Gorman [7] applied the method of superposition and obtained analytical solution for free vibration of a rectangular Kirchhoff plate with arbitrary boundary conditions, in the form of the Fourier series.Casimir [8] used the method of superposition, along with the method of projection for discretization of boundary conditions, and derived dynamic stiffness matrix of a rectangular Kirchhoff plate [8].In the same manner, continuous rectangular plate elements for the in-plane vibration [9], transverse vibration according to the Mindlin plate theory [10] and transverse vibration of a laminated plate using the higher-order shear deformation theory [11] have been formulated, as well.
Thin shell theory, which neglects shear deformation and rotary inertia, is generally acceptable if the shell thickness is less than 1/20 th of the wavelength and/or radius of curvature.Depending on the assumptions made in the strain-displacement relations, the force and moment resultant equations and equilibrium equations, various thin shell theories have been formulated.A comprehensive overview of the shell theories is given in the book of Leissa [12].In this paper the dynamic stiffness matrix of a circular cylindrical shell element according to the Flügge thin shell theory is formulated.In the governing partial differential equations time derivatives are eliminated by using the spectral decomposition, while the circumferential coordinate is eliminated by applying the solution in the form of Fourier series.Therefore, the governing system of three partial differential equations is transformed into the system of three ordinary differential equations, which depends only on the longitudinal coordinate.By expanding the podužne koordinate.Razvijanjem determinante ovog sistema jednačina dobijena je jedna obična diferencijalna jednačina osmog reda s konstantnim koeficijentima, čije je rešenje dobro poznato.Povezivanjem komponenti vektora sila i vektora pomeranja na krajevima ljuske formirana je dinamička matrica krutosti, koja je implementirana u za tu svrhu napisani Matlab [13] program za određivanje sopstvenih vrednosti i oblika oscilovanja kružne cilindrične ljuske.Metod dinamičke krutosti primenjen je u analizi slobodnih vibracija kružne cilindrične ljuske s različitim kombinacijama graničnih uslova, skokovitom promenom debljine, kao i s prstenastim međuosloncima.Putem numeričkih primera izvršena je verifikacija dobijenih rezultata poređenjem s dostupnim analitičkim rešenjima u literaturi, kao i s rezultatima Abaqus-a [14].
Sa u , v i w označena su komponentalna pomeranja srednje površi ljuske.determinant of this system of equations, one eighth order ordinary differential equation with constant coefficients is obtained.The solution of this equation is well known.The dynamic stiffness matrix is formulated by relating the force and displacement vectors at the ends of a circular cylindrical shell element.It is implemented in the computer program, written in Matlab [13], for computing natural frequencies and mode shapes of circular cylindrical shell assemblies.The dynamic stiffness method is applied to the free vibrations analysis of circular cylindrical shells with different combinations of boundary conditions, stepped thickness variation and intermediate ring supports.Verification of the obtained results is carried out in comparison with available analytical solutions in the literature, as well as with the results of the Abaqus [13].

FORMULATION OF THE PROBLEM
Figure (1) shows a closed circular cylindrical shell of constant thickness h, radius a and length L, where u, v and w are the displacements of the mid surface in x, ϕ and z directions, respectively.Slika 1. Geometrija i koordinatni sistem zatvorene kružne cilindrične ljuske Fig. 1.Geometry and coordinate system for a closed circular cylindrical shell

Governing differential equations for free vibration of a circular cylindrical shell according to the Flügge thin shell theory
The governing differential equations of the Flügge thin shell theory are [12]: je krutost u "ravni" ljuske, ρ je gustina mase i E je Young-ov modul elastičnosti.
Izrazi za presečne sile, u funkciji komponentalnih pomeranja, jesu [12]: In Eqs (2) ν is the Poisson's ratio, ( ) is the stiffness in the mid surface of shell, ρ is the mass density and E is the Young's modulus of elasticity.The expressions of forces and moments in terms of the displacements are [12]: Shear forces x Q and ϕ Q are determined from the equilibrium conditions: Konvencija o pozitivnim presečnim silama data je na sledećoj slici: The positive directions of the force and moment resultants are given in the next figure:

General solution of the free vibration problem
By using the method of separation of variables, the general solution of the system of Eqs. ( 1) is sought in the following form:
Izbor trigonometrijskih funkcija u 1. delu rešenja, jednačina (7), omogućava transformaciju sistema jednačina (6) u sistem tri obične diferencijalne jednačine: procedure for the 1 st part of Eq. ( 7) will be presented, since the solution for the 2 nd part can be obtained in the same way.It should be noted that the natural frequencies obtained for the 1 st and 2 nd part of Eq. ( 7) are the same, which means that for the closed circular cylindrical shell all the natural frequencies are double.

DYNAMIC STIFFNESS MATRIX K Dm
The displacement vector q and the force vector Q , contain the displacements and rotations, i.e. the forces and moments at the boundaries = 0 x and = x L , respectively: Slika 3. Komponente vektora pomeranja q i vektora sila Q

Fig. 3. Components of the displacement and force vectors
Komponente vektora pomeranja q i sila Q prikazane su na slici (3) i definisane su sledećim izrazima: The components of the vectors q and Q , shown in Fig. (3) are defined by the following expressions: 1 ˆcos cos Vektori pomeranja ˆm q i sila ˆm Q koji sadrže amplitude pomeranja i rotacija, odnosno presečnih sila za = 0 The new displacement and force vector, namely ˆm q and ˆm Q , that contain the amplitudes of displacements and rotations, i.e. forces, on the boundaries = 0 x and = x L for m th harmonic are: Veza između vektora ˆm q i vektora integracionih konstanti m C uspostavlja se preko matrice m D , odnosno veza između vektora ˆm Q i vektora integracionih konstanti se uspostavlja preko matrice m F , na sledeći način: The vector ˆm q is related to the vector of integration constants m C by the matrix m D , while the vectors ˆm Q and m C are related through the matrix m F , as follows: where the vector of integration constants 8 Ako se iz jednačine (32) izrazi m C u funkciji od ˆm q i zameni u jednačinu (33) dobija se veza između vektora ˆm If vector m C is expressed from the Eqs. (32)as a function of ˆm q and replaced in the Eq. ( 33), the relation between vectors ˆm Q and ˆm q is obtained in the following form:
When the global dynamic stiffness matrix is determined, the boundary conditions are applied by removing the rows and columns that correspond to the constrained displacements and rotations.In the numerical examples only the results for the most commonly used boundary conditions will be presented: • free F: componential displacements and rotation ϕ ψˆ are different from zero, • shear diaphragm SD: If a circular cylindrical shell with C-C boundary conditions is modelled with one continuous element, all rows and columns have to be deleted.Thus, for this combination of boundary conditions minimum two continuous elements have to be used.
The natural frequencies are obtained from the condition that ( ) = K det 0

Dm
. Since the dynamic stiffness matrix is transcendent matrix, the natural frequencies are obtained using searching techniques.In this paper instead of finding zeros of the determinant K Dm , the peaks of the expression ( ) K 1 logdet Dm [15]   are sought.Therefore, a computer program written in Matlab has been used.
When the exact natural frequencies are determined, the corresponding mode shapes are computed routinely.The mode shapes do not have true nodal lines, i.e. lines on the surface of a shell for which all displacement components are zero, but have lines along which two of displacement components are zero and the third has maximum values.In the case of axisymmetric vibration, the axial, circumferential and radial displacements are uncoupled, giving distinct nodal lines.

NUMERICAL RESULTS
The dynamic stiffness matrix based on the Flügge thin shell theory, formulated above, is verified through four numerical examples.The exact natural frequencies computed by the DSM are compared with the available analytical solution in the literature, as well as with the results from Abaqus.The adopted value of Poisson ratio

Example 5.1
In this example, the first nine natural frequencies (Hz) of a circular cylindrical shell with C-C boundary conditions and the following properties: = 210 a m and = 0.01 h m, are determined.Table 1 shows the results obtained by the DSM according to the Flügge thin shell theory, the results of Zhang [16], as well as the results obtained by Abaqus.Zhang used the Love-Timoshenko thin shell theory [17] and the wave propagation approach.The axial wave number is determined approximately as the wave number of an equivalent beam with similar boundary conditions.For the Abaqus models, two different finite elements are used: STRI3 (flat element with six nodes) and S4R (4-node, quadrilateral shell element with reduced integration and a large-strain formulation).The sign "-" in the table means that Zhang missed the fifth natural frequency.By applying a single continuous element, total agreement with the results of FEM analysis has been obtained.The total number of elements in the FE model was 50526 for S4R elements, i.e. 101052 for STRI3 elements.

Example 5.2
Table 2 shows the first dimensionless natural frequencies ) with F-F boundary conditions and given m= 1,2,3,4,5,6.The results obtained by the DSM based on the Flügge thin shell theory are compared with the exact results of Zhang [18], as well as those obtained by Abaqus.Xhang used the Goldenveizer-Novozhilov thin shell theory ( [19], [20]) and the state-space method to obtain the homogenous differential equations.In Abaqus, the circular cylindrical shell of radius a = 1 m has been modelled with 101052 STRI3 finite elements.Table 2 shows close agreement between the results obtained by Xhang, DSM and Abaqus.The relative differences between these results are negligible.In Fig. 6. the corresponding mode shapes obtained by the program written in the Matlab are presented.In this example, the application of the DSM in free vibration analysis of a stepped circular cylindrical shell has been demonstrated.The first four dimensionless natural frequencies of circular cylindrical shells with onestep thickness variation and with C-F, C-SD and C-C boundary conditions are calculated for a given m.Then results obtained by the DSM are compared with the results of Zhang [21].Zhang used the Flügge thin shell theory, the state-space method and domain decomposition method in order to satisfy continuity requirements between shell segments.In Table 3, the results are presented for shells with following boundary conditions and geometry: L are the thickness and the length of the first segment of the shell, 2 h is the thickness of the second segment and L is the total length of the shell.In Table 3 the relative differences in the percentage between the solutions are presented, which do not exceed 0.01%.  1 a E za ljusku sa SD-SD, C-C i F-F graničnim uslovima koja ima dva, odnosno tri prstenasta međuoslonca na jednakim rastojanjima.Na mestu prstenastog međuoslonca sprečeno je radijalno pomeranje, tj.

Example 5.4
In the last example, the fundamental dimensionless natural frequencies 1 a E for a shell with SD-SD, C-C and F-F boundary conditions and two, i.e. three equally spaced intermediate ring supports are computed.The radial displacement is prevented at the locations of ring supports, i.e.
The results obtained by the DSM using the Flügge theory are compared with the results of Xiang [22], who used the Goldenveizer-Novozhilov thin shell theory, the state-space technique and the domain decomposition method.Fundamental frequencies, along with the corresponding circumferential number m, are given in Table 4.The relative differences in percentage, which is also presented in

CONCLUSION
In this paper the dynamic stiffness matrix for a circular cylindrical shell element is formulated.Dynamic stiffness matrix is developed using the exact solution of the governing differential equations of free vibration according to the Flügge thin shell theory.The derived dynamic stiffness matrix is implemented in a Matlab program for calculation of natural frequencies and mode shapes of circular cylindrical shells and circular cylindrical shell assemblies.For several numerical examples natural frequencies are calculated using this program.The obtained results are validated against the available analytical results in the literature, as well as the results of the FE program Abaqus.
The analyses of the obtained results show excellent agreement with analytical solutions as well as the results obtained by Abaqus.The advantage of using the dynamic stiffness matrix is obvious: (1) for shells of constant cross section and arbitrary boundary conditions at the ends, it is sufficient to use only one DS element, except for shells clamped at both side, where it is necesje koristiti samo jedan element, izuzev kod obostrano uklještene ljuske gde je potrebno koristiti dva kontinualna elementa; (2) za ljuske promenljive debljine i za ljuske s međuosloncima potreban je samo jedan sary to apply two continuous elements; (2) for shells with stepped thickness variation and shells with intermediate ring supports it is necessary to apply one element for each segment of the shell.Contrary to the DSM in the element za svaki segment ljuske.Za razliku od MDK u MKE je potrebno mnogo elemenata (preko 100000) da bi se postigla željena tačnost.
FEM a significant number of elements (over 100000) have to be applied in numerical model in order to achieve the desired accuracy.Presented results demonstrate the main advantages of DSM in comparison with FEM -a significant reduction of the size of model and highly accurate results in a wide frequency range.

Nevenka KOLAREVIC Marija NEFOVSKA DANILOVIC Mira PETRONIJEVIC
In this paper the dynamic stiffness method is used for free vibration analysis of a circular cylindrical shell.The dynamic stiffness matrix is formulated on the base of the exact solution for free vibration of a circular cylindrical shell according to the Flügge thin shell theory.The matrix is frequency dependent and, besides the stiffness, includes inertia and damping effects.The derived dynamic stiffness matrix is implemented in the code developed in a Matlab program for computing natural frequencies and mode shapes of a circular cylindrical shell.Several numerical examples are carried out.The obtained results are validated against the results obtained by using the commercial finite element program Abaqus as well as the available analytical solutions from the literature.Key words: free vibration, circular cylindrical shell, Flügge thin shell theory, dynamic stiffness matrix,

Table 1 .
First nine natural frequencies (Hz) for a circular cylindrical shell with C-C boundary conditions:

Table 2 .
Fundamental dimensionless natural frequency

Table 3 .
First four dimensionless natural frequencies

Table 4 .
Fundamental dimensionless natural frequencies