#### Abstract

A framework for free vibration analysis of plate-like structures is presented in the paper. Based on the previously formulated dynamic stiffness elements, FREEVIB object-oriented 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 stiffness-to-weight 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 closed-form 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 high-frequency 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 domain-based strong-form 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 two-dimensional 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 first-order shear deformation theory (FSDT) have been implemented and applied in the design-optimization process of stiffened plate structures. Boscolo and Banerjee [11–14] developed dynamic stiffness elements for rectangular Levy-type isotropic and laminated composite plates undergoing both in-plane 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 closed-form 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 in-plane 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 higher-order shear deformation theory (HSDT) [28], as well as for the plate element ongoing in-plane 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 cross-ply 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 plate-like structures, within which the object-oriented 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 plate-like 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 2*a* × 2*b.* 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 in-plane 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 single-layer 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 shear-deformable 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 stress-strain relations for the orthotropic layer in the *xyz* coordinates can now be computed as

For the symmetric cross-ply laminates considered within the proposed framework, = = = 0.

##### 2.3. Dynamic Stiffness Formulation

Assuming there is no external loading acting on the plate, the Euler-Lagrange 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: symmetric-symmetric (SS), symmetric-antisymmetric (SA), antisymmetric-symmetric (AS), and antisymmetric-antisymmetric (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 (in-plane 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 one-dimensional 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 (2*a* or 2*b*), *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 in-plane vibrations of a single plate represent two independent states (for single layer plates, symmetric sandwich panels, or symmetric cross-ply 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 in-plane 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, frequency-dependent 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 cross-ply 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, object-oriented 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

Procedure-oriented programming (POP) and object-oriented 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, object-oriented 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 for-loop 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 for-loops 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 FSDT-based 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: C-clamped edge; S-simply-supported edge; F-free 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 (eight-node 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 shear-deformable 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 cross-ply symmetric stacking sequence (0/90/90/0). Boundary conditions: C-clamped edge; S-simply-supported edge; F-free 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 one-fold 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 element-transfer 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 + 18

*M*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 object-oriented 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: single-layer 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 out-of-plane displacement components, while **T**_{R i} is the rotation matrix related to the in-plane 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 in-plane 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 TR-36046 and TR-36048, is acknowledged.