Research Article  Open Access
Framework for DynamicStiffnessBased Free Vibration Analysis of PlateLike Structures
Abstract
A framework for free vibration analysis of platelike structures is presented in the paper. Based on the previously formulated dynamic stiffness elements, FREEVIB objectoriented software in Python environment has been created. Software design and structure as well as a wide range of possible structural problems that could be analyzed using the FREEVIB are presented. Through several illustrative examples including free vibration analysis of stepped, stiffened and folded plate structures, implying isotropic or orthotropic material formulations, the efficiency and accuracy of the FREEVIB is demonstrated. The possibilities of further extensions and improvements of the software are discussed.
1. Introduction
Composite laminates have attracted great attention over past decades in a wide range of engineering fields where they are applied as structural components of aircrafts, space structures, highway and pedestrian bridges and platforms amongst others, due to its high stiffnesstoweight ratio and high fatigue and environmental resistance. Consequently, vibration and buckling analysis of such structures has become a great concern for researchers.
Analytical solutions of the free vibration problem of rectangular composite laminated plate having special types of boundary conditions have been derived based on both classical and shear deformation theories [1, 2]. Since composite laminated plates are usually assembled in a structure consisting of plates connected at different angles, with or without stiffeners and arbitrarily assigned boundary conditions, it becomes difficult or even impossible to find closedform solutions of the free vibration equations of such structures. Therefore, the finite element method (FEM) [3, 4] is most frequently used to solve the corresponding free vibration problems. However, in mid and highfrequency ranges, the computation of the free vibration response using conventional FEM requires very fine mesh of finite elements, leading to the large number of governing equations to be solved and increasing the computational cost. Consequently, the results may become inaccurate and unreliable. To overcome these issues, the dynamic stiffness method (DSM) has been applied to solve a wide range of vibration problems [5, 6].
The core of the DSM is the frequency domainbased strongform solution of the governing equations of motion derived for the corresponding elastodynamic problem. Therefore, the number of dynamic stiffness elements necessary for structural discretization is frequency independent and influenced only by the change in the geometrical and/or material properties of the structure. At the same time, using simple assembly procedure as in the FEM makes the DSM an efficient and highly accurate method capable of solving a broad variety of vibration problems.
Wittrick [7] was the first who derived twodimensional dynamic stiffness elements and corresponding dynamic stiffness matrices for the free vibration and buckling analysis of rectangular isotropic long plate assemblies with vibration and buckling modes varying sinusoidally in the longitudinal direction. His research has been extended to anisotropic plates based on which the computer program VIPASA [8] has been developed for calculating natural frequencies and critical factors of stiffened plate assemblies. Later on, based on the VIPASA, the VICON software [9] has been developed, which enabled constraints and supporting structure to be accounted for in the free vibration and buckling analysis of plate assemblies. Finally, in the VICONOPT computer code [10], the dynamic stiffness elements based on both classical plate theory (CPT) and firstorder shear deformation theory (FSDT) have been implemented and applied in the designoptimization process of stiffened plate structures. Boscolo and Banerjee [11–14] developed dynamic stiffness elements for rectangular Levytype isotropic and laminated composite plates undergoing both inplane and transverse vibration and demonstrated high accuracy and reliability of the DSM implemented in a computer code DySAP to compute free vibration characteristics of composite stiffened plate assemblies [15, 16].
Problems arising from finding the closedform solution for the most general case, regardless of the boundary conditions, have been overcome in the works of Casimir et al. [17] and Banerjee et al. [18, 19] who formulated the dynamic stiffness matrix of rectangular isotropic plate element undergoing transverse vibration based on CPT. Based on the classical laminated plate theory (CLPT), the dynamic stiffness matrix for an orthotropic plate element has been formulated recently by Ghorbel et al. [20, 21], Liu and Banerjee [22, 23], Liu et al. [24], and Papkov and Banerjee [25]. The spectral dynamic stiffness formulation for inplane modal analysis of composite plates and prismatic solids was developed by Liu [26]. In a series of contributions, dynamic stiffness matrices for a completely free rectangular plate element based on the FSDT [27] and higherorder shear deformation theory (HSDT) [28], as well as for the plate element ongoing inplane vibrations [29], have been developed by the authors. Recently, the above formulations have been broadened to the free transverse vibration analysis of sandwich [30], symmetric crossply laminated composite plates [31, 32], and composite stiffened plate assemblies and GFRP beams with or without straight through cracks [33].
The dynamic stiffness elements presented in previous authors’ work served as a basis for the development of computational framework for free vibration analysis of laminated composite platelike structures, within which the objectoriented software FREEVIB has been created. The main objective of the paper is to present the structure, capabilities, possible improvements, further extensions, and application of the FREEVIB in the vibration analysis of platelike structures.
The paper is organized as follows. Brief overview of the basic equations describing the kinematics and constitutive law of laminated composite plates are given in Sections 2.1 and 2.2. Afterwards, the key components of the dynamic stiffness method regarding the dynamic stiffness matrix development, rotation and assembly procedure, and mode shape computation are briefly explained in Section 2.3. A wide range of structural problems that could be analyzed using the FREEVIB is depicted in Section 2.4. Structure, design, and possible extensions and improvements of the FREEVIB software are explained in detail in Section 3. Through several illustrative examples, the efficiency and accuracy of the FREEVIB is presented in Section 4. Finally, concluding remarks are given in Section 5.
2. Theoretical Basis of the Dynamic Stiffness Method
The current capability of the presented computational framework implies the formulation of several dynamic stiffness elements, based on different plate theories. The overview of the dynamic stiffness method is presented in this section. For more details, the reader is referred to authors’ previous work [27–33].
2.1. Displacement Field
We consider an assembly of rectangular laminated composite plates, each having the dimensions 2a × 2b. Each plate is composed of n orthotropic layers, as shown in Figure 1(a). The overall plate thickness of the plate in the assembly is denoted as h_{i}, while the thickness of the layer of the plate is denoted as h_{ki}, where k = 1, 2, …, n denotes the layer number starting from the bottom of the laminated plate. The global coordinate system (X, Y, Z) is the Cartesian coordinate system, while the local coordinate system of each plate is located in the midplane of the laminated plate, with the z axis pointing upwards.
It is worth mentioning that the plates are connected along the boundary lines instead at nodes (like in FEM). Therefore, the plate assembly shown in Figure 1(a) is modeled in the dynamic stiffness method by using 7 dynamic stiffness elements (Figure 1(b)).
For each plate in the assembly, depending on the element type, the displacement field is defined in the local coordinate system, in an arbitrary point (x, y, z), at an arbitrary time point t, according to either CLPT, FSDT, or HSDT [1].
First, we define the displacement field of the HSDT: assuming that there is no deformation in the midplane of the laminate, the displacement field of an arbitrary point (x,y,z) at an arbitrary time point t is given as [2]:where and are the rotations about the x and y axes, respectively, are the slopes and c_{1} = 4/(3 h^{2}).
The HSDT is based on the assumption that the transverse normal is inextensible and is not perpendicular to the midplane after deformation. In the FSDT, the transverse normal remains straight after deformation (but rotated with respect to the midplane (Figure 2(b))), while in the HSDT this criterion is relaxed by expanding the inplane displacements u and as cubic functions with respect to the thickness coordinate (Figure 2(c)). This eliminates the application of the shear correction factor in the HSDT, while the shear correction factor k is required in the FSDT to satisfy the compatibility conditions of shear stresses.
(a)
(b)
(c)
The displacement fields of both FSDT and CLPT (Figure 2(a)) can be easily derived from the HSDT. Finally, the presented formulation can be easily reduced for the analysis of singlelayer isotropic and orthotropic plates, as well as for the analysis of multilayer plates with transversely isotropic layers (i.e., sandwich panels).
The linear strain fields of HSDT, FSDT, and CLPT, associated with the assumed displacement field, can be derived in the usual manner [1]. Note that HSDT and FSDT are referred as sheardeformable plate theories.
2.2. Constitutive Equations
The constitutive equations for the orthotropic layer within the laminate in the assembly, assuming linear elastic material behavior, in the local (material) coordinate system of a single layer which coincides with the fiber direction can be written as [1]where is the matrix of the reduced stiffness components of the layer for the plane stress case [1], while and denote the stress and strain vectors of the layer in the local coordinate system of a single layer. Since each laminate in the assembly is made of several orthotropic layers, having the material axes oriented arbitrarily with respect to the xyz coordinates of the laminate (Figure 1(a)), the constitutive relations must be transformed from the layer coordinate system to the laminate coordinate system using the transformation matrix T^{k}:where is the constitutive matrix for the layer in the xyz coordinate system of the considered laminate in the assembly. The stressstrain relations for the orthotropic layer in the xyz coordinates can now be computed as
For the symmetric crossply laminates considered within the proposed framework, = = = 0.
2.3. Dynamic Stiffness Formulation
Assuming there is no external loading acting on the plate, the EulerLagrange equations of motion for the HSDT, FSDT, and CLPT can be derived using Hamilton’s principle in terms of the displacement components [1]. The dynamic stiffness matrix of a laminated composite plate is obtained from the governing equations of motion through the following steps.
2.3.1. Transformation of the Governing Equations of Motion to Frequency Domain
The governing equations of motion are transformed into the frequency domain by assuming a harmonic representation of the displacement/rotation field:where are the amplitudes of the displacement/rotation in the frequency domain. Having in mind that equation (5) is valid for all angular frequencies ω in the considered frequency range, the argument ω will be omitted in further representations. After the substitution of the above transformation into the governing equations of motion, the equations of motion are transformed into the following set of partial differential equations:where L is the matrix of the differential operators defined in terms of the plate stiffness coefficients, the mass moments of inertia, and the angular frequency ω [27–29].
2.3.2. Superposition of Symmetry Contributions
Displacement/rotation amplitudes of a rectangular plate element can be presented as a superposition of four symmetry contributions: symmetricsymmetric (SS), symmetricantisymmetric (SA), antisymmetricsymmetric (AS), and antisymmetricantisymmetric (AA), according to [17–19]. By the superposition of the symmetry contributions, it is possible to analyze only one quarter of a rectangular laminated plate, which significantly reduces the size of the corresponding dynamic stiffness matrices. By using the method of separation of variables, the general solution for each symmetry contribution can be represented in the Fourier series form aswhere is the corresponding displacement/rotation function, and are the unknown displacement/rotation functions, while and are the Fourier base functions, defined depending on the symmetric or antisymmetric contribution [27–29]. In addition, it should be pointed out that Fourier base functions differ, depending upon the corresponding elastodynamic problem (inplane or transverse vibration) and applied plate theory (CLPT, FSDT, or HSDT). More details about the base functions are given in the Appendix. In practical calculations, the infinite Fourier series must be truncated. Consequently, the solution of equation (7) is approximate, and its accuracy depends on the number of terms in the general solution. The solutions for all symmetry contributions are given in [27–29].
2.3.3. Projection Method
The vector of displacement components along plate boundaries x = a and y = b is called the displacement vector . The corresponding force vector consists of force components along plate boundaries x = a and y = b, and it is denoted as . These vectors are functions of spatial variables x and y, and consequently, they cannot be related explicitly as in the case of onedimensional elements. This issue can be overcome with the aid of the projection method [17]. Therefore, instead of using the vectors and , new projection vectors and are introduced:where L denotes the length of the corresponding boundary line (2a or 2b), s is related to x or y coordinate, while is the projection matrix composed of the Fourier base functions. Now, the components of the projection vectors are the Fourier coefficients in the series expansion (see [27–33] for details). The relation between the projection vectors and and the vector of integration constants is obtained as
In order to preserve the squared form of matrices and , the number of terms in the Fourier series expansion in equation (9) should be the same as the number of terms in the general solutions given in equation (7). Since integrals given in equation (8) are computed analytically, elements of matrices and are derived in explicit form ([27–33]). Consequently, the code execution is significantly accelerated. Finally, by eliminating the vector C^{ij} from equation (9), the following relation between the projection vectors and is obtained:where is the dynamic stiffness matrix for the ij contribution (where ij = SS, SA, AS, or AA), for the transverse vibrations.
The dynamic stiffness matrix for a completely free dynamic stiffness element based on CLPT, FSDT, or HSDT, which relates the projections of the forces and displacements on the four boundary lines of the plate, is obtained by using the following expression (for details see [27–33]):where T is the transformation matrix [27–29]. The size of the dynamic stiffness matrix depends on the number of terms in the general solution M and is equal to 32M + 12 (HSDT), 24M + 8 (FSDT) and 16M + 8 (CLPT).
Finally, considering that transverse and inplane vibrations of a single plate represent two independent states (for single layer plates, symmetric sandwich panels, or symmetric crossply laminated composite plates), the dynamic stiffness matrix of the single plate can be written aswhere denotes the dynamic stiffness matrix of plate element for transverse vibration, while is the dynamic stiffness matrix of plate element undergoing inplane vibration which can be determined in the same way as the .
2.3.4. Rotation of Dynamic Stiffness Matrices to Global Coordinates
For stiffened plate assemblies where the plates (i.e., plate 1 and plate 2) are connected to each other with an arbitrary angle between them, vibrations of plate 1 cause the vibrations of the corresponding plate 2, and vice versa. Consequently, it is necessary to relate the displacement and force projection vectors and in the local coordinate system xyz of the single plate (Figure 1) and the corresponding projection vectors and in the global coordinate system XYZ of the plate assembly. This is accomplished by using the rotation matrix T_{R}:
The rotation matrix T_{R} depends on the number of terms in the general solution, the angle between the local (xyz) and global coordinate system (XYZ), and the selected dynamic stiffness element (Figure 1). The rotation matrices are given in the Appendix for the CLPT, FSDT, and HSDT. According to the established relations between the projection vectors in the local and global coordinate systems, dynamic stiffness matrix of the single plate in global coordinate system is derived as
2.3.5. Assembly Procedure and Assignment of Boundary Conditions
Dynamic stiffness matrices of single plates are assembled in the global dynamic stiffness matrix of plate assembly, using similar assembly procedure as in the conventional FEM. The procedure was demonstrated in authors’ previous works [27–33]. In the analysis, arbitrary boundary conditions can be applied by removing the rows and columns of the global dynamic stiffness matrix that correspond to the components of the constrained displacement projections. The following combinations of boundary conditions are usually used along the plate edges (Table 1).

2.3.6. Computation of Natural Frequencies and Mode Shapes
The dynamic stiffness matrix is square, frequencydependent matrix whose size depends on the number of terms M in the general solution. The natural frequencies can be computed from the following equation:where is the global dynamic stiffness submatrix of the plate assembly related to the unknown generalized displacement projections of the plate assembly. Since the elements of the dynamic stiffness matrix contain transcendental functions, the solutions can be obtained by applying some of the search methods. To avoid numerical difficulties when calculating the zeros of equation (15), the natural frequencies can be determined as maxima of the following expression:
The above expression is computed for all frequencies in the frequency range of interest with a frequency increment Δω. Consequently, the accuracy of the computed natural frequencies is affected only by the frequency increment. A drawback of this method is that the coincident natural frequencies arising in the case of symmetric geometry and boundary conditions can be omitted. In this case, the Wittrick–Williams algorithm [7, 8] can be applied to obtain the number of natural frequencies that are less than a trial frequency . After the natural frequencies have been computed, the mode shape that corresponds to the natural frequency ω_{i} is obtained in the usual manner.
2.4. Analyzed Problems
Using the procedure described in the previous sections, a variety of structural problems can be analyzed (Figure 3): (a) individual plates, (b) plate assembly of arbitrary shape, (c) stepped plates, (d) stiffened plates, and (e) cracked plates. In addition, different material properties of single and multilayer plates can be considered: (f) single layer isotropic or orthotropic plates, (g) symmetric sandwich panels, or (h) symmetric crossply laminated composite plates. Note that arbitrary combinations of boundary conditions can be assigned along plate edges, as shown in Section 2.3.5. Finally, no frequency limitations have been detected, allowing the same mesh to be utilized for the calculation of any frequency (in contrary to FEM).
3. Software Design
FREEVIB software is completely written in Python. It is an interpreted, interactive, objectoriented programming language, which incorporates modules, exceptions, dynamic typing, very high level dynamic data types, and classes. Python combines remarkable power with very clear syntax [34]. The language comes with a large standard library that covers areas such as string processing, software engineering, and operating system interfaces. The multiprocessing module in Python’s Standard Library enables the parallel execution of the code to notably speed up the computation with a few modifications.
3.1. Program Design
Procedureoriented programming (POP) and objectoriented programming (OOP) are two possible approaches in the design of engineering software. OOP has been widely utilized in engineering software development based on the finite element method ([35–43]). However, the dynamic stiffness method has been used so far mainly in the form of the POP. In the paper, objectoriented design is introduced to enable code encapsulation, class inheritance, and further code reusability. The class relationships of the presented software are shown in Figure 4. For the sake of simplicity, the arguments and return value types of a method and the constructor of a class are omitted. Properties and methods inherited from a super class are not listed, as well.
The input parameters are inserted as properties of InputData class. For the input, simple text file may be created in the designated format or generated using the existing graphical preprocessors such as “GiD—the personal pre and postprocessor,” developed in CIMNE, Barcelona [44] or “FreeCAD: Your Own 3D Parametric Modeler” [45]. The input parameters related to the geometry, material, and number of terms in trigonometric series (M) are transferred to the basic classes: Point, DOF, Line, Material, Lamina, and Section, as well as to the DynStiffElement class, which is the super class which implements almost all methods related to the creation of the element dynamic stiffness matrix.
Through the ObjectConstructor module, the instances of all basic classes are created based on the input parameters. Line class is created by the composition from Point and DOF class instances, while Section class is created by the composition from Lamina class instances. For all Material class instances, constitutive matrix Q^{k} is calculated by using the calcMatrixQ() method (equation (2)). Transformation matrix T^{k} and the global constitutive matrix of the lamina (equation (3)) are calculated by using the calcMatrixT() and calcMatrixQGlobal() methods of the Lamina class. Finally, all stiffness and inertia coefficients based on different plate theories and related to the instances of the Section class (equation (6)) are calculated by using the calcCoefficients() method of the Section class.
DynStiffElement class is the super class which implements almost all methods related to the creation of the dynamic stiffness matrix. In addition, CLPTElement, FSDTElement, and HSDTElement classes are inherited from this class. Their class structures are identical, but the implemented methods calcMatrixSSb(ω,M), calcMatrixSAb(ω,M), calcMatrixASb(ω,M), calcMatrixAAb(ω,M), calcMatrixTb(M), and calcMatrixTRb(M) are related to the calculation of different dynamic stiffness, transformation and rotation matrices for transverse vibration, explained in the Section 2.3.3, equations (10) and (11). The dynamic stiffness matrices of all dynamic stiffness elements are calculated according to equation (14) in global coordinates, by using the calcMatrixKD() method of the DynStiffElement class. They are stored as the properties of the DynStiffElement class instances, for all elements in the assembly.
In the FreeVibrationAnalysis class, which has only one instance which stores the data related to the performed analysis, the calculation of the dynamic stiffness matrices is performed for all elements in the assembly using the above procedure. This is performed in the forloop which loops through every frequency of interest, in the range F_{start} : F_{step} : F_{end} (input data). The global dynamic stiffness matrix associated with the unknown displacement components is created through the deriveMatrixKnn() method of the FreeVibrationAnalysis class.
As soon as the is created for every frequency of interest, the algorithm for calculation of natural frequencies is called through the deriveFrequencies() method of the FreeVibrationAnalysis class (see equations (15) and (16)). The obtained natural frequencies are listed and plotted in the output file.
The mathematical calculations are performed using numpy [46] and scipy [47] packages for scientific computing.
3.2. Possibilities for Software Improvement
3.2.1. Element Formulation
The presented code could be extended by adding new dynamic stiffness element formulations. The example is open cylindrical shell element derived by Kolarević [48]. The developed element is compatible with the already implemented element formulations by means of the degrees of freedom along element edges. In addition, the formulation of the dynamic stiffness element based on a sophisticated layerwise theory of Reddy [1] could improve the accuracy of the analysis of thick laminated composite plates. This formulation is presented by Boscolo and Banerjee [14].
3.2.2. Analysis Types
The family of implemented algorithms could be extended by adding new analysis types, such as response analysis or buckling analysis. For example, the response analysis of plate assembly using the dynamic stiffness method (see [29] for details) could be implemented by using the class relations shown in Figure 5, which implies implementation of the deriveResponse() method in the ResponseAnalysis subclass inherited from the Analysis class, followed with the additional set of input parameters (InputData class) related to the transient loading (intensity, duration, etc.).
3.2.3. Parallelization
Multicore computers and cluster computers are the common two kinds of hardware that support parallel computing [49]. With the rapid development of multicore CPU technique, developing the program on a multicore desktop computer using the multiprocessing module in Python’s Standard Library enables the parallel execution of the code. The possibility for the development of the algorithm for the parallel execution is motivated by the fact that, after the global dynamic stiffness matrix is assembled , the evaluation of equation (16) should be performed for every frequency of interest, in the range F_{start} : F_{step} : F_{end}. Having in mind that the results for every frequency are mutually independent, the parallelization of the code using the parallel forloops would increase the speed of the sequential execution by the factor ∼n_{CPU} (where n_{CPU} is the number of the CPUs in the multicore desktop computer).
4. Numerical Examples
To demonstrate the efficiency and accuracy of the FREEVIB, several examples have been presented in this section. For the FSDTbased dynamic stiffness element, shear correction factor k = 5/6 has been used in all calculations.
Example 1. The first example deals with the free vibration analysis of thick stepped concrete isotropic CSFS plate (Figure 6). Plate geometry is defined as h_{1}/a = 0.05, h_{2}/a = 0.03. The following boundary conditions have been assigned: Cclamped edge; Ssimplysupported edge; Ffree edge.
The present solution based on the CLPT, FSDT, and HSDT using M = 5 terms in the series expansion has been validated against the results of the eigenvalue analysis from the commercial software Abaqus, using 2500 S8R elements (eightnode doubly curved thick shell element with reduced integration and homogeneous shell section). The dimensionless natural frequencies are calculated for the first 6 modes. The results are elaborated in Table 2.
Obviously, FREEVIB is capable of accurately predicting natural frequencies of the thick stepped isotropic plate, when sheardeformable dynamic stiffness elements are applied. The numerical solution in Abaqus is closer to the obtained results when refined theory is applied. As expected, the average relative difference is 12.17% for the CLPT, 0.064% for the FSDT, and 0.059% for the HSDT.
Finally, first 6 mode shapes along with their corresponding frequencies are calculated and plotted for the HSDT element and compared with the Abaqus solution in Figure 7. Excellent agreement is obtained.

Example 2. The second example deals with the free vibration analysis of square laminated composite plate with a single stiffener (Figure 8). The plate consists of 4 orthotropic layers in crossply symmetric stacking sequence (0/90/90/0). Boundary conditions: Cclamped edge; Ssimplysupported edge; Ffree edge. Three different h/a ratios have been considered (h/a = 0.02, h/a = 0.05, and h/a = 0.10), where h is the overall plate thickness and a is the side length of the plate.
Each layer is made of orthotropic material having the following properties: E_{2}/E_{1} = 40, G_{12}/E_{2} = G_{13}/E_{2} = 0.6, G_{23}/E_{2} = 0.5, and ν_{12} = 0.25.
The present solution based both on FSDT and HSDT with M = 5 terms in the series expansion has been applied and validated against the results of the eigenvalue analysis from the commercial software Abaqus, using 3000 S8R elements with composite shell section. The dimensionless natural frequencies are calculated for the first 6 modes and presented in Table 3. The average relative differences when compared with the FEM solution in Abaqus are presented, as well.
When applied to the free vibration analysis of laminated composite plates with a single stiffener, FREEVIB more accurately predicts natural frequencies in comparison to the Abaqus software. The Abaqus solution is in better agreement with the obtained results when thin laminate is considered, while for the thick laminates (and therefore for higher frequency range), the differences are higher varying from 4.80%–11.30% in comparison with the exact solution from FREEVIB.

Example 3. The final example deals with the free vibration analysis of isotropic cantilever onefold plate clamped along the edge X = 0, with h/a = 0.02 and a_{1}/a = a_{2}/a = 0.5. The plate geometry is shown in Figure 9. The study has been made for a crank angle of α = 120°C. The material properties for both plates are E = 10.92 GPa, ν = 0.3.
The convergence study of the present model is performed using M = 1, 3, or 5 terms in the series expansion and dynamic stiffness elements based both on FSDT and HSDT. The converged solution (M = 5) for both plate bending theories is validated against the results of the eigenvalue analysis from the commercial software Abaqus, using 2500 S8R elements. The dimensionless natural frequencies are calculated for the first 4 modes. An additional comparison is made with the results obtained by Liu and Huang [50] using finite elementtransfer matrix method (FETM) based on the CPT. The results are elaborated in Table 4.
Obviously, FREEVIB demonstrates high accuracy in the calculation of natural frequencies of folded plate structures. The converged solution is achieved even for a very small number of terms in series solution (M = 5).
The computational aspect of FREEVIB is very important. In the presented example, the FE mesh of 2500 S8R elements in Abaqus implies 2652 nodes with 6 DOFs per node in the global stiffness and mass matrices. On the other hand, in the FREEVIB, the assembly of 2 dynamic stiffness elements implies only 7 edges with 9 + 18M DOFs per edge. This leads to the following equation:This confirms the computational advantage of the proposed method in comparison with the conventional FEM.

5. Conclusions
Based on the authors’ previous research, within a computational framework for free vibration analysis of laminated composite plate assemblies using the dynamic stiffness method, FREEVIB objectoriented software has been developed in Python. Key features of the FREEVIB including the structure, design, application, possible extensions, and improvements have been presented in the paper. The efficiency and accuracy of the FREEVIB has been demonstrated through several illustrative numerical examples.
The combination of simple assembly procedure as in the FEM with the advantages of the analytical methods, makes the FREEVIB an efficient numerical tool for free vibration analysis of a wide range of structural problems including individual plates, plate assemblies of arbitrary shape, and stepped and stiffened plates connected at arbitrary angles, with or without cracks. In addition, different material properties of single and multilayer plates can be considered: singlelayer isotropic or orthotropic plates, sandwich panels, or laminated composite plates, without any restrictions regarding the boundary conditions or the frequency limitations. In most applications, only 3–5 terms in the series expansions were sufficient to achieve high accuracy and reliability of the computed natural frequencies and mode shapes in comparison with the conventional FEM solutions. Moreover, the structural discretization is influenced only by the change in geometrical/material properties of plate assemblies.
Owing to its design and structure, FREEVIB could be easily extended to include new dynamic stiffness element formulations and analysis types allowing dynamic response and buckling analysis of plate assemblies.
Appendix
Submatrix of direction cosines between the global coordinate system (XYZ) and local coordinate system (xyz) of the single plate in the assembly can be written aswhere λ_{ij} is the cosine of the angle between the axis in the XYZ and axis in the xyz system.
Transformation matrix T_{R} for the dynamic stiffness element can be written as
The submatrix T_{R b} is the rotation matrix related to the outofplane displacement components, while T_{R i} is the rotation matrix related to the inplane displacement components. The submatrix T_{R i} is the same for all element types and can be written aswhere t_{i0} is the 4 × 36 matrix with the following nonzero terms: t_{i0}[1,1] = λ_{Xx}, t_{i0} [1,2] = λ_{Xy}, t_{i0} [1,3] = λ_{Xz}, t_{i0} [2,10] = λ_{Yx}, t_{i0} [2,11] = λ_{Yy}, t_{i0} [2,12] = λ_{Yz}. t_{i0} [3,19] = λ_{Xx}, t_{i0} [3,20] = λ_{Xy}, t_{i0} [3,21] = λ_{Xz}, t_{i0} [4,28] = λ_{Yx}, t_{i0} [4,29] = λ_{Yy}, and t_{i0} [4,30] = λ_{Yz}, while
Submatrix T_{R b} for the CLPT dynamic stiffness element can be written as follows:where t_{b0} is the following 8 × 36 matrix:while:
Submatrix T_{R b} for the FSDT dynamic stiffness element can be written as follows:where t_{b0} is the following 8 × 36 matrix:while:
Submatrix T_{R b} for the HSDT dynamic stiffness element can be written aswhere t_{b0} is the following 12 × 36 matrix:while:
General solution for inplane displacements (u, ) and rotations (Φ_{x} and Φ_{y}):
General solution for transverse displacement :
Fourier base functions:
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The financial support of the Government of the Republic of Serbia  Ministry of Education, Science and Technological Development, under the projects TR36046 and TR36048, is acknowledged.
References
 J. N. Reddy, Mechanics of laminated composite plates and shells. Theory and analysis, CRC Press, Boca Raton, FL, USA, 2nd edition, 2004.
 J. N. Reddy, “A simple higherorder theory for laminated composite plates,” Journal of Applied Mechanics, vol. 51, no. 4, pp. 745–752, 1984. View at: Publisher Site  Google Scholar
 O. C. Zienkiewicz and R. L. Taylor, The Finite Element MethodVolume 1: The Basis, ButterworthHeinemann, Oxford, UK, 5th edition, 2000.
 O. C. Zienkiewicz and R. L. Taylor, The Finite Element MethodVolume 2: Solid Mechanics, ButterworthHeinemann, Oxford, UK, 5th edition, 2000.
 J. F. Doyle, Wave Propagation in Structures: Spectral Analysis Using Fast Discrete Fourier Transforms, Springer, Berlin, Germany, 2nd edition, 1997.
 U. Lee, Spectral Element Method in Structural Dynamics, John Wiley and Sons, Hoboken, NJ, USA, 2009.
 W. H. Wittrick, “General sinusoidal stiffness matrices for buckling and vibration analyses of thin flatwalled structures,” International Journal of Mechanical Sciences, vol. 10, no. 12, pp. 949–966, 1968. View at: Publisher Site  Google Scholar
 W. H. Wittrick and F. W. Williams, “Buckling and vibration of anisotropic or isotropic plate assemblies under combined loadings,” International Journal of Mechanical Sciences, vol. 16, no. 4, pp. 209–239, 1974. View at: Publisher Site  Google Scholar
 M. S. Anderson, F. W. Williams, and C. J. Wright, “Buckling and vibration of any prismatic assembly of shear and compression loaded anisotropic plates with an arbitrary supporting structure,” International Journal of Mechanical Sciences, vol. 25, no. 8, pp. 585–596, 1983. View at: Publisher Site  Google Scholar
 D. M. McGowan and M. S. Anderson, “Development of curvedplate elements for the exact buckling analysis of composite plate assemblies including transverse shear effects,” American Institute of Aeronautics and Astronautics, Plymouth, NH, USA, 1999, Technical Report. View at: Google Scholar
 M. Boscolo and J. R. Banerjee, “Dynamic stiffness elements and their applications for plates using first order shear deformation theory,” Computers and Structures, vol. 89, pp. 395–410, 2011. View at: Publisher Site  Google Scholar
 M. Boscolo and J. R. Banerjee, “Dynamic stiffness method for exact inplane free vibration analysis of plates and plate assemblies,” Journal of Sound and Vibration, vol. 330, no. 12, pp. 2928–2936, 2011. View at: Publisher Site  Google Scholar
 F. A. Fazzolari, M. Boscolo, and J. R. Banerjee, “An exact dynamic stiffness element using a higher order shear deformation theory for free vibration analysis of composite plate assemblies,” Composite Structures, vol. 96, pp. 262–278, 2013. View at: Publisher Site  Google Scholar
 M. Boscolo and J. R. Banerjee, “Layerwise dynamic stiffness solution for free vibration analysis of laminated composite plates,” Journal of Sound and Vibration, vol. 333, no. 1, pp. 200–227, 2014. View at: Publisher Site  Google Scholar
 M. Boscolo and J. R. Banerjee, “Dynamic stiffness formulation for composite Mindlin plates for exact modal analysis of structures. Part I: Theory,” Computers and Structures, vol. 9697, pp. 61–73, 2012. View at: Publisher Site  Google Scholar
 M. Boscolo and J. R. Banerjee, “Dynamic stiffness formulation for composite Mindlin plates for exact modal analysis of structures. Part II: results and applications,” Computers and Structures, vol. 9697, pp. 74–83, 2012. View at: Publisher Site  Google Scholar
 J. B. Casimir, S. Kevorkian, and T. Vinh, “The dynamic stiffness matrix of twodimensional elements: application to Kirchhoff's plate continuous elements,” Journal of Sound and Vibration, vol. 287, no. 3, pp. 571–589, 2005. View at: Publisher Site  Google Scholar
 J. R. Banerjee, S. O. Papkov, X. Liu, and D. Kennedy, “Dynamic stiffness matrix of a rectangular plate for the general case,” Journal of Sound and Vibration, vol. 342, pp. 177–199, 2015. View at: Publisher Site  Google Scholar
 X. Liu and J. R. Banerjee, “Free vibration analysis for plates with arbitrary boundary conditions using a novel spectraldynamic stiffness method,” Computers and Structures, vol. 164, pp. 108–126, 2016. View at: Publisher Site  Google Scholar
 O. Ghorbel, J. B. Casimir, L. Hammami, I. Tawfiq, and M. Haddar, “Dynamic stiffness formulation for free orthotropic plates,” Journal of Sound and Vibration, vol. 346, pp. 361–375, 2015. View at: Publisher Site  Google Scholar
 O. Ghorbel, J. B. Casimir, L. Hammami, I. Tawfiq, and M. Haddar, “Inplane dynamic stiffness matrix for a free orthotropic plate,” Journal of Sound and Vibration, vol. 364, pp. 234–246, 2016. View at: Publisher Site  Google Scholar
 X. Liu and J. R. Banerjee, “An exact spectraldynamic stiffness method for free flexural vibration analysis of orthotropic composite plate assembliespart I: theory,” Composite Structures, vol. 132, pp. 1274–1287, 2015. View at: Publisher Site  Google Scholar
 X. Liu and J. R. Banerjee, “An exact spectraldynamic stiffness method for free flexural vibration analysis of orthotropic composite plate assembliespart II: applications,” Composite Structures, vol. 132, pp. 1288–1302, 2015. View at: Publisher Site  Google Scholar
 X. Liu, H. I. Kassem, and J. R. Banerjee, “An exact spectral dynamic stiffness theory for composite platelike structures with arbitrary nonuniform elastic supports, mass attachments and coupling constraints,” Composite Structures, vol. 142, pp. 140–154, 2016. View at: Publisher Site  Google Scholar
 S. O. Papkov and J. R. Banerjee, “A new method for free vibration and buckling analysis of rectangular orthotropic plates,” Journal of Sound and Vibration, vol. 339, pp. 342–358, 2015. View at: Publisher Site  Google Scholar
 X. Liu, “Spectral dynamic stiffness formulation for inplane modal analysis of composite plate assemblies and prismatic solids with arbitrary classical/nonclassical boundary conditions,” Composite Structures, vol. 158, pp. 262–280, 2016. View at: Publisher Site  Google Scholar
 N. Kolarevic, M. NefovskaDanilovic, and M. Petronijevic, “Dynamic stiffness elements for free vibration analysis of rectangular Mindlin plate assemblies,” Journal of Sound and Vibration, vol. 359, pp. 84–106, 2015. View at: Publisher Site  Google Scholar
 N. Kolarevic, M. Marjanović, M. NefovskaDanilovic, and M. Petronijevic, “Free vibration analysis of plate assemblies using the dynamic stiffness method based on the higher order shear deformation theory,” Journal of Sound and Vibration, vol. 364, pp. 110–132, 2016. View at: Publisher Site  Google Scholar
 M. NefovskaDanilovic and M. Petronijevic, “Inplane free vibration and response analysis of isotropic rectangular plates using the dynamic stiffness method,” Computers and Structures, vol. 152, pp. 82–95, 2015. View at: Publisher Site  Google Scholar
 M. Marjanović, N. Kolarević, M. NefovskaDanilović, and M. Petronijević, “Free vibration study of sandwich plates using a family of novel shear deformable dynamic stiffness elements: limitations and comparison with the finite element solutions,” ThinWalled Structures, vol. 107, pp. 678–694, 2016. View at: Publisher Site  Google Scholar
 M. Marjanovic, N. Kolarevic, M. NefovskaDanilovic, and M. Petronijevic, “Shear deformable dynamic stiffness elements for a free vibration analysis of composite plate assembliespart I: theory,” Composite Structures, vol. 159, pp. 728–744, 2017. View at: Publisher Site  Google Scholar
 M. Marjanovic, N. Kolarevic, M. NefovskaDanilovic, and M. Petronijevic, “Shear deformable dynamic stiffness elements for a free vibration analysis of composite plate assembliespart II: numerical examples,” Composite Structures, vol. 159, pp. 183–196, 2017. View at: Publisher Site  Google Scholar
 E. Damnjanovic, M. Marjanovic, and M. NefovskaDanilovic, “Free vibration analysis of stiffened and cracked composite plate assemblies using sheardeformable dynamic stiffness elements,” Composite Structures, vol. 180, pp. 723–740, 2017. View at: Publisher Site  Google Scholar
 Python Software Foundation, “Python language reference, version 2.7,” 2019, http://www.python.org. View at: Google Scholar
 J. R. Filho and P. R. B. Devloo, “Object oriented programming in scientific computations: the beginning of a new era,” Engineering Computations, vol. 8, no. 1, pp. 81–87, 1991. View at: Publisher Site  Google Scholar
 G. R. Miller, “An objectoriented approach to structural analysis and design,” Computers and Structures, vol. 40, no. 1, pp. 75–82, 1991. View at: Publisher Site  Google Scholar
 T. Zimmermann, Y. DuboisPèlerin, and P. Bomme, “Objectoriented finite element programming: I. Governing principles,” Computer Methods in Applied Mechanics and Engineering, vol. 98, no. 2, pp. 291–303, 1992. View at: Publisher Site  Google Scholar
 R. M. V. Pidaparti and A. V. Hudli, “Dynamic analysis of structures using objectoriented techniques,” Computers and Structures, vol. 49, no. 1, pp. 149–156, 1993. View at: Publisher Site  Google Scholar
 G. R. Miller, S. Banerjee, and K. Sribalaskandarajah, “A framework for interactive computational analysis in geomechanics,” Computers and Geotechnics, vol. 17, no. 1, pp. 17–37, 1995. View at: Publisher Site  Google Scholar
 G. C. Archer, “Objectoriented finite element analysis,” University of California at Berkeley, Berkeley, CA, USA, 1996, Ph.D. thesis. View at: Google Scholar
 F. T. McKenna, “Objectoriented finite element programming: frameworks for analysis, algorithms and parallel computing,” University of California at Berkeley, Berkeley, CA, USA, 1997, Ph.D. thesis. View at: Google Scholar
 P. Dadvand, “Framework for developing finite element code for multidisciplinary applications,” CIMNEInternational Centre for Numerical Methods in Engineering, Barcelona, Spain, 2007, Ph.D. thesis. View at: Google Scholar
 H. Qin, Z. Liu, Y. Liu, and H. Zhong, “An objectoriented MATLAB toolbox for automotive body conceptual design using distributed parallel optimization,” Advances in Engineering Software, vol. 106, pp. 19–32, 2017. View at: Publisher Site  Google Scholar
 A. Coll, R. Ribó, M. Pasenau et al., 13 User Manual, CIMNEInternational Centre for Numerical Methods in Engineering, Barcelona, Spain, 2016.
 J. Riegel, W. Mayer, and Y. van Havre, “FreeCAD (20012017),” 2019, http://www.freecadweb.org. View at: Google Scholar
 T. Oliphant, A Guide to NumPy, Trelgol Publishing, USA, 2006.
 E. Jones, T. Oliphant, P. Peterson et al., “SciPy: open source scientific tools for Python,” 2019, http://www.scipy.org. View at: Google Scholar
 N. Kolarević, “Vibrations and buckling of plates and shells using dynamic stiffness method,” Faculty of Civil Engineering, University of Belgrade, Belgrade, Serbia, 2016, Ph.D. thesis. View at: Google Scholar
 D. E. Culler, J. P. Singh, and A. Gupta, Parallel Computer Architecture: a Hardware/Software Approach, Gulf Professional Publishing, San Francisco, CA, USA, 1999.
 W. H. Liu and C. C. Huang, “Vibration analysis of folded plates,” Journal of Sound and Vibration, vol. 157, no. 1, pp. 123–137, 1992. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2019 Miroslav Marjanović et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.