APPLICATION OF DYNAMIC STIFFNESS METHOD IN NUMERICAL FREE VIBRATION ANALYSIS OF STIFFENED PLATES

Primena građevinskih materijala visokih mehaničkih karakteristika, pre svega čelika, dovodi do smanjenja dimenzija (debljine) konstruktivnih pločastih elemenata. Time se postiže slična nosivost za manji utrošak materijala, kao i ušteda u ceni. Međutim, vitki elementi postaju osetljivi na izbočavanje, što uzrokuje upotrebu ploča sa ukrućenjima. Ploče sa ukrućenjima imaju široku primenu u građevinarstvu i drugim inženjerskim disciplinama, posebno pri projektovanju, na primer, mostova većih raspona i manjih poprečnih preseka, te izradi oplata brodova, konstrukcija aviona. Tokom svog radnog veka, ove konstrukcije često su izložene dinamičkom opterećenju, te je


INTRODUCTION
Application of building materials of high mechanical properties (i.e. steel), leads to the reduction of dimensions (thickness) of structural plate-like elements. This maintains the same structural capacity with reduced material cost, as well the price reduction. However, slender elements become vulnerable to buckling, leading to the use of stiffened plates.
Stiffened plates are widely used in civil engineering and other engineering fields, especially in design of longspan bridges with small cross-sections, construction of ship hulls, aircrafts, etc. During their working life, these structures are usually exposed to dynamic loading, thus the calculation of dynamic structural response is of high proračun njihovog dinamičkog odgovora veoma značajan u inženjerskoj praksi. U takvim slučajevima, neophodno je odrediti osnovne dinamičke karakteristike sistema, kao što su sopstvene frekvencije i oblici oscilovanja.
Dynamic response of thin plates can be obtained using Kirchhoff classical plate theory (CPT). In the case of thick plates, this theory fails to provide adequate results because the transverse shear deformation is neglected. Therefore, it is necessary to use Mindlin plate theory (first-order shear deformation theory -FSDT), which accounts for the constant transverse shear deformation through the plate thickness. Along with the development of different plate theories, the adequate analytical methods have been derived, too [1]. However, these methods are based on the exact solution of differential equations of motion of plates having special boundary conditions. Leissa [2] provided the comprehensive review of analytical solutions of the free vibration problem for plates of different shape, based on the Kirchhoff classical plate theory. Liew et al. [3] studied the free vibration of thick plates with arbitrary boundary conditions using the Raylegh-Ritz method. The application of the above mentioned methods is restricted to the free vibration analysis of single plates and cannot be easily extended to the analysis of more complex plate assemblies with different geometry and material properties, which are mostly used in engineering practice (stiffened plates, sandwich panels, stepped plates, etc.). In such cases, numerical methods are applied, i.e. the finite element method (FEM) [4]. It is a well-known fact that in the dynamic finite element analysis, the size of the finite element depends on the wavelength of the highest frequency of interest for the analysis. According to [5,6], the ratio between the wavelength of the highest frequency and the finite element size should be between 10 and 20, in order to obtain the results of satisfactory accuracy. Minimum number of finite elements is directly proportional to the highest considered frequency, thus for the complex structures (where higher modes of vibration must be taken into account), the number of finite elements becomes very high, increasing the overall computational time.
Lately, the dynamic stiffness method (DSM) [7][8][9][10][11][12][13][14][15] is increasingly used in the free vibration analysis of plates. DSM combines the properties of the FEM (such as physical discretization and the possibility of assembly of single elements into the unique global system) with the solution of the displacement field in the form of exact solution of differential equation of the free vibration problem. Having in mind that interpolation functions which describe the displacement field in the DSM represent the exact solution of the differential equation of motion in the frequency domain, the approximation errors are eliminated. The plate discretization is necessary only if some geometrical and/or physical discontinuity is present. This reduces both the numbers of elements in the analysis and degrees of freedom, as well as computational time and error possibility, in comparison with the FEM. In contrary to the FEM where the mass is lumped in nodes of finite elements, in the DSM the mass is continuously distributed. In addition, material damping can be easily included in the analysis using the complex elasticity modulus Ec = E(1+2iζ), where ζ is the relative damping coefficient. Using this Na taj način, omogućeno je da različiti elementi konstrukcije imaju različito prigušenje, što je još jedna od prednosti MDK u poređenju s MKE.
In the paper, the numerical model for free undamped vibration analysis of stiffened Mindlin plates with arbitrary boundary conditions, based on the DSM, has been presented. Starting from the dynamic stiffness matrices for the free vibration analysis of transverse and in-plane vibrations, rotation matrix for different plate positions (perpendicular to the reference plane), has been derived [16]. Similar procedure as in the FEM has been applied for the development of the global dynamic stiffness matrix of stiffened plate assembly. Computer program in MATLAB [17] has been developed for the free vibration analysis of plate assemblies. The validation of the presented procedure has been performed by comparison of the obtained results with those obtained by using the software package Abaqus [18].

PROCEDURE FOR DEVELOPMENT OF DYNAMIC STIFFNESS MATRIX OF RECTANGULAR PLATE
Development of the dynamic stiffness matrix for rectangular plate element undergoing transverse and inplane vibrations has been given in detail in [13,15]. In the paper, only the basic steps in the procedure will be presented. The procedure starts from with equations of motion of Mindlin plate element in the time domain: Pretpostavlja se da su pomeranja harmonijske funkcije frekvencije ω, odnosno: where u, v, w, φx, φy are the displacement components, h is the plate thickness, ρ is the mass density, E is the modulus of elasticity, G is the shear modulus, ν is the Poisson's coefficient, It is assumed that the displacements are the harmonic functions of frequency ω, i.e.: Na slici 1a prikazano je polje pomeranja pravougaone ploče po Mindlin-ovoj teoriji.
is the amplitude of the displacement field (u, v, w, φx or φy) in the frequency domain. Using the above assumption, equations of motion (1) are transformed from time to frequency domain. Figure 1 shows the displacement field of the rectangular plate based on Mindlin plate theory Arbitrary displacement field can be represented as a superposition of four symmetry contributions with respect to x and y coordinate axes: symmetric-symmetric (SS), symmetric-anti-symmetric (SA), anti-symmetricsymmetric (AS) and anti-symmetric-anti-symmetric (AA): Na slici 1b prikazana je simetrična deformacija ploče u odnosu na obe koordinatne ose (SS). Podelom polja pomeranja na četiri slučaja simetrije, moguće je analizirati samo jednu četvrtinu ploče, čime se znatno umanjuje red dinamičkih matrica krutosti i ubrzava proračun. Rešenje jednačina kretanja u frekventnom domenu pretpostavlja se u obliku beskonačnog Fourierovog reda u sledećem obliku: Figure 1b illustrates the double-symmetric plate deformation with respect to both coordinate axes (SS). By dividing the displacement field into four symmetry contributions, it is possible to analyze only a quarter of a plate, which significantly reduces the order of the dynamic stiffness matrices and accelerates the calculation. The solution of the equations of motion in the frequency domain is assumed in the form of an infinite Fourier's series as: xyCfxgy ω u (4) gde su: Na sličan način, vektor sila ˆi j Q na konturama ploče dobija se zamenom koordinata kontura u jednačinu (5): In a similar manner, the force vector ˆi j Q on the plate boundaries is obtained by replacing the coordinates of the boundaries in equation (5): S obzirom da su komponente vektora pomeranja i sila na konturama ploče funkcije prostornih koordinata x i y, nije moguće direktno uspostaviti vezu između tih vektorasa jedne strane i vektora integracionih konstanti C -sa druge strane. Ovaj problem može se rešiti pomoću metode projekcije koja se bazira na predstavljanju funkcija pomeranja i sila na konturi ploče u vidu Fourier-ovog reda: Since the components of the displacement and force vectors along plate boundaries are functions of spatial coordinates x and y, it is impossible to establish a direct relation between these vectors on one side, and the vector of integration constants C on the other side. This problem can be solved by using the projection method, which is based on the representation of the displacement and force functions along the plate boundary in the Fourier series form: Dinamička matrica krutosti čitave ploče % D K može se dobiti primenom transfer matrice, kao što je pokazano u radovima [13,15].

STIFFENED PLATES
Transverse and in-plane vibrations of a single plate represent two independent states. Therefore, the dynamic stiffness matrix of the single plate can be written as:  (11) U slučaju ploča sa ukrućenjima, gde su ploče međusobno spojene pod pravim uglom, poprečne vibracije jedne ploče izazivaju vibracije u ravni druge ploče i obrnuto. Zbog toga je potrebno uspostaviti vezu između vektora projekcija pomeranja i sila % q i % Q u lokalnom i vektora % * q i % * Q u globalnom koordinatnom sistemu (slike 2 i 3), pomoću matrice rotacije T: For stiffened plates (where the plates are perpendicular to each other), transverse vibrations of a single plate cause the in-plane vibrations of a corresponding perpendicular plate, and vice versa. Consequently, it is necessary to establish the relation between the projection vectors of displacements % q and forces % Q in the local, and the corresponding projection vectors % * q and % * Q in the global c. s. (Figures 2 and   3). This is accomplished by using the rotation matrix T: Figure 2. Displacement components on the plate boundaries in the local coordinate system Slika 3. Komponente pomeranja na konturama ploča u globalnom koordinatnom sistemu Figure 3. Displacement components on the plate boundaries in the global coordinate system U skladu sa definisanim relacijama projekcija u lokalnom i globalnom koordinatnom sistemu, dinamička matrica krutosti ploče u globalnom koordinatnom sistemu dobija se kao:
Dinamička matrica krutosti je kvadratna, frekventno zavisna matrica čiji red zavisi od broja članova reda M usvojenog rešenja. Sopstvene frekvencije određuju se iz sledeće jednačine: Dynamic stiffness matrices of individual plates are assembled in the global dynamic stiffness matrix of the plate assembly, using similar assembly procedure as in the conventional FEM. Note that the plates are connected along the boundaries instead at nodes. In the analysis, arbitrary boundary conditions can be applied.
The dynamic stiffness matrix is the square, frequency-dependent matrix whose size depends on the number of terms M in the general solution. The natural frequencies can be determined from the following equation: gde je % * To avoid numerical difficulties when calculating the zeroes of equation (14), the natural frequencies can be determined as maxima of the following expression:

NUMERICAL EXAMPLES
Application of the dynamic stiffness method in the free vibration analysis of stiffened plates is illustrated in the following examples. In order to validate the presented method based on the DSM, a computer code in MATLAB [17] was developed for the free vibration analysis of stiffened plates with arbitrary boundary conditions. The results were compared with the results obtained using Abaqus [18].
The first example is concerned with a reinforced concrete girder of a box cross-section, (E = 31.5GPa, ν = 0.2 and ρ = 2500kg/m 3 ). The geometry and boundary conditions are presented in Figure 4a. In this case, the number of dynamic stiffness elements is equal to 9 (which is minimal number o f elements), while the number of boundary lines is 26. Boundary lines parallel to x-and z-axis are simply supported (Sx, Sz), while the boundary lines parallel to y-axis are either free (F) or clamped (C). Boundary conditions are assigned in the following way: • simply supported edge parallel to x-axis (Sx): u=v=w=φy=φz=0, • simply supported edge parallel to z-axis (Sz): u=v=w=φx=φy=0, • uklještena ivica paralelna sa y-osom (C): u=v=w=φx=φy=φz=0, • slobodna ivica paralelna sa y-osom (F): φz=0.
In order to check the convergence of the proposed method, the first ten natural frequencies f = ω/2π have been calculated using different number of terms in the general solution. The results are presented in Table 1 and compared with the numerical solutions obtained using 46250 Abaqus S4R type plate elements (element size 0.1 m). For the solutions with M = 7 terms in the series expansion, the differences ∆ from the solutions obtained in Abaqus are calculated.

. Dimensionless natural frequencies [Hz] of reinforced concrete girder of a box cross-section
The second example is concerned with the reinforced concrete ribbed slab, (E = 31.5GPa, ν = 0.2 and ρ = 2500kg/m 3 ), whose geometry and boundary conditions are presented in Figure 4b. In this case the number of dynamic stiffness elements is equal to 5, while the number of boundary lines is equal to 16. The first ten natural frequencies f = ω/2π have b een calculated using different number of terms in the general solution. The results are presented in Table 2 and compared with the FEM numerical solutions obtained using 23520 Abaqus S4R type plate element (element size 0.025 m). For the solutions with M = 7 terms in the series expansion, the differences ∆ from the solutions obtained using Abaqus are calculated. As in the previous example, the natural frequencies computed using the proposed method are in excellent agreement with the FEM (the average difference is 0.22% for S-F-S-F and 0.14% for S-C-S-C case). Na slici 5 izvršeno je upoređivanje pojedinih oblika oscilovanja i odgovarajućih sopstvenih frekvencija (tonovi 2, 4, 7 i 8) armirano-betonske korube, za kombinaciju konturnih uslova S-C-S-C i M=7 članova reda. Oblici oscilovanja sračunati su primenom originalnog programa u MATLAB-uzasnovanog na metodi dinamičke krutosti, kao i primenom komercijalnog programa Abaqus. Dobijeno je odlično poklapanje oblika oscilovanja.
To illustrate the advantage of the DSM in comparison with the FEM by means of the necessary computational time, number of degrees of freedom in the numerical model of the girder with the box cross-section, formulated using the DSM or the FEM, will be compared. It is worth mentioning that the number of degrees of freedom in the DSM depends on the number of contours and number of terms in series expansion M, while in the FEM it depends on the number of nodes and the element type (in the considered example, the finite element has 6 nodal degrees of freedom). Number of degrees of freedom in the DSM (for M=7)and the FEM is therefore:
that confirms the potential of the presented method by means of the necessary computational time.

CONCLUSION
This paper deals with the application of the dynamic stiffness method in the free vibration analysis of stiffened plates. The procedure for the development of dynamic stiffness matrix of rectangular plate, as well as rotation matrices for different plate orientations and perpendicularly oriented with respect to the reference plane, has been provided. The similar procedure as in the FEM for the development of the global dynamic stiffness matrix of the stiffened plate, has been applied. Computer sistema ploča. Sopstvene frekvencije i oblici oscilovanja određeni su za različite kombinacije konturnih uslova.

ZAHVALNICA
Autori zahvaljuju Ministarstvu prosvete, nauke i tehnološkog razvoja Republike Srbije na finansijskoj podršci u okviru projekata TR-36046 i TR-36048. frequencies and mode shapes were calculated for various combinations of boundary conditions. Verification of the results was performed by comparing the results obtained using Abaqus software package based on the FEM. The accuracy of the results obtained by the DSM depends solely on the number of terms in the series expansion of the adopted solution. Rapid convergence of the results obtained by the DSM has been detected. By using only three to five terms, the results of high accuracy have been obtained. However, for the higher modes of vibration, the number of required members in the series expansion should be increased.
Based on these results, it can be concluded that the DSM has great potential in the dynamic analysis of structures that can be extended to the analysis of plates and shells using higher-order theory, the plates connected with arbitrary angles, as well as for the structural elements made of modern composite materials.