# Framework for Dynamic-Stiffness-Based Free Vibration Analysis of Plate-Like Structures

Authors: Miroslav Marjanović; Marija Nefovska-Danilović; Emilija Damnjanović
Journal: Shock and Vibration
 (2019)

---

## 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.


---

## Body

## 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 Sections2.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 Section3. 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 ith plate in the assembly is denoted as hi, while the thickness of the kth layer of the ith plate is denoted as hki, 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.Figure 1
Assembly of (a) three arbitrarily oriented multilayer plates and (b) seven dynamic stiffness elements.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 Figure1(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]:(1)ux,y,z,t=zϕyx,y,t−c1⋅z3ϕyx,y,t+∂wx,y,t∂x,vx,y,z,t=−zϕxx,y,t−c1⋅z3−ϕxx,y,t+∂wx,y,t∂y,wx,y,z,t=wx,y,t,where ϕxx,y,t and ϕyx,y,t are the rotations about the x and y axes, respectively, ∂wx,y,t/∂xand∂wx,y,t/∂y are the slopes and c1 = 4/(3 h2).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 (Figure2(b))), while in the HSDT this criterion is relaxed by expanding the in-plane displacements u and v 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.Figure 2
Deformed line element in the CLPT (a), FSDT (b), and HSDT (c).
(a)
(b)
(c)The displacement fields of both FSDT and CLPT (Figure2(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 thekth orthotropic layer within the ith 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](2)σ¯k=Qkε¯k,where Qk is the matrix of the reduced stiffness components of the kth layer for the plane stress case [1], while σ¯k and ε¯k denote the stress and strain vectors of the kth 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 Tk:(3)Q¯k=Tk−1QkTk,where Q¯k is the constitutive matrix for the kth layer in the xyz coordinate system of the considered laminate in the assembly. The stress-strain relations for the kth orthotropic layer in the xyz coordinates can now be computed as(4)σk=σxσyτxyτxzτyzk=Q¯11Q¯12Q¯1600Q¯12Q¯22Q¯2600Q¯16Q¯26Q¯6600000Q¯55Q¯45000Q¯45Q¯44kεxεyγxyγxzγyzk=Q¯kεk.For the symmetric cross-ply laminates considered within the proposed framework,Q¯16 = Q¯26 = Q¯45 = 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:(5)ux,y,t=u^x,y,ωeiωt,where u^x,y,ω are the amplitudes of the displacement/rotation ux,y,t 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:(6)Lu^x,y=0,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 elementu^x,y 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 as(7)u^ijx,y=∑m∞1Umijx1fmijy+∑m∞2Umijy2fmijxi,j=SorA,where u^ijx,y is the corresponding displacement/rotation function, 1Umijx and 2Umijy are the unknown displacement/rotation functions, while 1fmijy and 2fmijx 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 componentsu^ij along plate boundaries x = a and y = b is called the displacement vector q^ij. The corresponding force vector consists of force components along plate boundaries x = a and y = b, and it is denoted as Q^ij. 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 q^ij and Q^ij, new projection vectors q˜ij and Q˜ij are introduced:(8)q˜ij=2L∫sHijq^ijds,Q˜ij=2L∫sHij,Q^ijds,where L denotes the length of the corresponding boundary line (2a or 2b), s is related to x or y coordinate, while Hij 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 q˜ij and Q˜ij and the vector of integration constants Cij is obtained as(9)q˜ij=D˜DijCij,Q˜ij=F˜DijCij.In order to preserve the squared form of matricesD˜Dij and F˜Dij, 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 D˜Dij and F˜Dij are derived in explicit form ([27–33]). Consequently, the code execution is significantly accelerated. Finally, by eliminating the vector Cij from equation (9), the following relation between the projection vectors Q˜ij and q˜ij is obtained:(10)Q˜ij=F˜DijD˜Dij−1q˜ij=K˜Dtijq˜ij,where K˜Dtij 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]):(11)K˜Dt=12TTK˜DtSS0000K˜DtSA0000K˜DtAS0000K˜DtAAT,where T is the transformation matrix [27–29]. The size of the dynamic stiffness matrix K˜Dt 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 as(12)K˜D=K˜Dt00K˜Di,where K˜Dt denotes the dynamic stiffness matrix of plate element for transverse vibration, while K˜Di is the dynamic stiffness matrix of plate element undergoing in-plane vibration which can be determined in the same way as the K˜Dt.


#### 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 vectorsq˜ and Q˜ in the local coordinate system xyz of the single plate (Figure 1) and the corresponding projection vectors q˜∗ and Q˜∗ in the global coordinate system XYZ of the plate assembly. This is accomplished by using the rotation matrix TR:(13)q˜=TR⋅q˜∗.The rotation matrixTR 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(14)K˜D∗=TRTK˜DTR.


#### 2.3.5. Assembly Procedure and Assignment of Boundary Conditions

Dynamic stiffness matrices of single platesK˜D∗ are assembled in the global dynamic stiffness matrix K˜DG 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).Table 1
Assignment of boundary conditions.
Edge type
Simply-supported edge
Clamped edge
Edge parallel to theX-axis
u = v= w = 0, ϕy = ϕz = 0, w,x = w,z = 0
u = v = w = 0, ϕx = ϕy = ϕz = 0, w,x = w,y = w,z = 0
Edge parallel to theY-axis
u = v = w = 0, ϕx = ϕz = 0, w,y = w,z = 0
Edge parallel to theZ-axis
u = v = w = 0, ϕx = ϕy = 0, w,x = w,y = 0


#### 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 termsM in the general solution. The natural frequencies can be computed from the following equation:(15)detK˜D,nnGω=0,where K˜D,nnG is the global dynamic stiffness submatrix of the plate assembly related to the unknown generalized displacement projections q˜nG of the plate assembly. Since the elements of the dynamic stiffness matrix K˜D,nnG 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:(16)gω=log1detK˜D,nnGω.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 ith 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 (Figure3): (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).Figure 3
Structural problems which can be analyzed using the proposed computational framework.


## 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 ith plate in the assembly is denoted as hi, while the thickness of the kth layer of the ith plate is denoted as hki, 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.Figure 1
Assembly of (a) three arbitrarily oriented multilayer plates and (b) seven dynamic stiffness elements.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 Figure1(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]:(1)ux,y,z,t=zϕyx,y,t−c1⋅z3ϕyx,y,t+∂wx,y,t∂x,vx,y,z,t=−zϕxx,y,t−c1⋅z3−ϕxx,y,t+∂wx,y,t∂y,wx,y,z,t=wx,y,t,where ϕxx,y,t and ϕyx,y,t are the rotations about the x and y axes, respectively, ∂wx,y,t/∂xand∂wx,y,t/∂y are the slopes and c1 = 4/(3 h2).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 (Figure2(b))), while in the HSDT this criterion is relaxed by expanding the in-plane displacements u and v 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.Figure 2
Deformed line element in the CLPT (a), FSDT (b), and HSDT (c).
(a)
(b)
(c)The displacement fields of both FSDT and CLPT (Figure2(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 thekth orthotropic layer within the ith 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](2)σ¯k=Qkε¯k,where Qk is the matrix of the reduced stiffness components of the kth layer for the plane stress case [1], while σ¯k and ε¯k denote the stress and strain vectors of the kth 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 Tk:(3)Q¯k=Tk−1QkTk,where Q¯k is the constitutive matrix for the kth layer in the xyz coordinate system of the considered laminate in the assembly. The stress-strain relations for the kth orthotropic layer in the xyz coordinates can now be computed as(4)σk=σxσyτxyτxzτyzk=Q¯11Q¯12Q¯1600Q¯12Q¯22Q¯2600Q¯16Q¯26Q¯6600000Q¯55Q¯45000Q¯45Q¯44kεxεyγxyγxzγyzk=Q¯kεk.For the symmetric cross-ply laminates considered within the proposed framework,Q¯16 = Q¯26 = Q¯45 = 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:(5)ux,y,t=u^x,y,ωeiωt,where u^x,y,ω are the amplitudes of the displacement/rotation ux,y,t 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:(6)Lu^x,y=0,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 elementu^x,y 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 as(7)u^ijx,y=∑m∞1Umijx1fmijy+∑m∞2Umijy2fmijxi,j=SorA,where u^ijx,y is the corresponding displacement/rotation function, 1Umijx and 2Umijy are the unknown displacement/rotation functions, while 1fmijy and 2fmijx 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 componentsu^ij along plate boundaries x = a and y = b is called the displacement vector q^ij. The corresponding force vector consists of force components along plate boundaries x = a and y = b, and it is denoted as Q^ij. 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 q^ij and Q^ij, new projection vectors q˜ij and Q˜ij are introduced:(8)q˜ij=2L∫sHijq^ijds,Q˜ij=2L∫sHij,Q^ijds,where L denotes the length of the corresponding boundary line (2a or 2b), s is related to x or y coordinate, while Hij 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 q˜ij and Q˜ij and the vector of integration constants Cij is obtained as(9)q˜ij=D˜DijCij,Q˜ij=F˜DijCij.In order to preserve the squared form of matricesD˜Dij and F˜Dij, 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 D˜Dij and F˜Dij are derived in explicit form ([27–33]). Consequently, the code execution is significantly accelerated. Finally, by eliminating the vector Cij from equation (9), the following relation between the projection vectors Q˜ij and q˜ij is obtained:(10)Q˜ij=F˜DijD˜Dij−1q˜ij=K˜Dtijq˜ij,where K˜Dtij 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]):(11)K˜Dt=12TTK˜DtSS0000K˜DtSA0000K˜DtAS0000K˜DtAAT,where T is the transformation matrix [27–29]. The size of the dynamic stiffness matrix K˜Dt 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 as(12)K˜D=K˜Dt00K˜Di,where K˜Dt denotes the dynamic stiffness matrix of plate element for transverse vibration, while K˜Di is the dynamic stiffness matrix of plate element undergoing in-plane vibration which can be determined in the same way as the K˜Dt.


### 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 vectorsq˜ and Q˜ in the local coordinate system xyz of the single plate (Figure 1) and the corresponding projection vectors q˜∗ and Q˜∗ in the global coordinate system XYZ of the plate assembly. This is accomplished by using the rotation matrix TR:(13)q˜=TR⋅q˜∗.The rotation matrixTR 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(14)K˜D∗=TRTK˜DTR.


### 2.3.5. Assembly Procedure and Assignment of Boundary Conditions

Dynamic stiffness matrices of single platesK˜D∗ are assembled in the global dynamic stiffness matrix K˜DG 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).Table 1
Assignment of boundary conditions.
Edge type
Simply-supported edge
Clamped edge
Edge parallel to theX-axis
u = v= w = 0, ϕy = ϕz = 0, w,x = w,z = 0
u = v = w = 0, ϕx = ϕy = ϕz = 0, w,x = w,y = w,z = 0
Edge parallel to theY-axis
u = v = w = 0, ϕx = ϕz = 0, w,y = w,z = 0
Edge parallel to theZ-axis
u = v = w = 0, ϕx = ϕy = 0, w,x = w,y = 0


### 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 termsM in the general solution. The natural frequencies can be computed from the following equation:(15)detK˜D,nnGω=0,where K˜D,nnG is the global dynamic stiffness submatrix of the plate assembly related to the unknown generalized displacement projections q˜nG of the plate assembly. Since the elements of the dynamic stiffness matrix K˜D,nnG 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:(16)gω=log1detK˜D,nnGω.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 ith mode shape that corresponds to the natural frequency ωi is obtained in the usual manner.


## 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:(5)ux,y,t=u^x,y,ωeiωt,where u^x,y,ω are the amplitudes of the displacement/rotation ux,y,t 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:(6)Lu^x,y=0,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 elementu^x,y 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 as(7)u^ijx,y=∑m∞1Umijx1fmijy+∑m∞2Umijy2fmijxi,j=SorA,where u^ijx,y is the corresponding displacement/rotation function, 1Umijx and 2Umijy are the unknown displacement/rotation functions, while 1fmijy and 2fmijx 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 componentsu^ij along plate boundaries x = a and y = b is called the displacement vector q^ij. The corresponding force vector consists of force components along plate boundaries x = a and y = b, and it is denoted as Q^ij. 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 q^ij and Q^ij, new projection vectors q˜ij and Q˜ij are introduced:(8)q˜ij=2L∫sHijq^ijds,Q˜ij=2L∫sHij,Q^ijds,where L denotes the length of the corresponding boundary line (2a or 2b), s is related to x or y coordinate, while Hij 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 q˜ij and Q˜ij and the vector of integration constants Cij is obtained as(9)q˜ij=D˜DijCij,Q˜ij=F˜DijCij.In order to preserve the squared form of matricesD˜Dij and F˜Dij, 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 D˜Dij and F˜Dij are derived in explicit form ([27–33]). Consequently, the code execution is significantly accelerated. Finally, by eliminating the vector Cij from equation (9), the following relation between the projection vectors Q˜ij and q˜ij is obtained:(10)Q˜ij=F˜DijD˜Dij−1q˜ij=K˜Dtijq˜ij,where K˜Dtij 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]):(11)K˜Dt=12TTK˜DtSS0000K˜DtSA0000K˜DtAS0000K˜DtAAT,where T is the transformation matrix [27–29]. The size of the dynamic stiffness matrix K˜Dt 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 as(12)K˜D=K˜Dt00K˜Di,where K˜Dt denotes the dynamic stiffness matrix of plate element for transverse vibration, while K˜Di is the dynamic stiffness matrix of plate element undergoing in-plane vibration which can be determined in the same way as the K˜Dt.


## 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 vectorsq˜ and Q˜ in the local coordinate system xyz of the single plate (Figure 1) and the corresponding projection vectors q˜∗ and Q˜∗ in the global coordinate system XYZ of the plate assembly. This is accomplished by using the rotation matrix TR:(13)q˜=TR⋅q˜∗.The rotation matrixTR 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(14)K˜D∗=TRTK˜DTR.


## 2.3.5. Assembly Procedure and Assignment of Boundary Conditions

Dynamic stiffness matrices of single platesK˜D∗ are assembled in the global dynamic stiffness matrix K˜DG 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).Table 1
Assignment of boundary conditions.
Edge type
Simply-supported edge
Clamped edge
Edge parallel to theX-axis
u = v= w = 0, ϕy = ϕz = 0, w,x = w,z = 0
u = v = w = 0, ϕx = ϕy = ϕz = 0, w,x = w,y = w,z = 0
Edge parallel to theY-axis
u = v = w = 0, ϕx = ϕz = 0, w,y = w,z = 0
Edge parallel to theZ-axis
u = v = w = 0, ϕx = ϕy = 0, w,x = w,y = 0


## 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 termsM in the general solution. The natural frequencies can be computed from the following equation:(15)detK˜D,nnGω=0,where K˜D,nnG is the global dynamic stiffness submatrix of the plate assembly related to the unknown generalized displacement projections q˜nG of the plate assembly. Since the elements of the dynamic stiffness matrix K˜D,nnG 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:(16)gω=log1detK˜D,nnGω.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 ith 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 (Figure3): (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).Figure 3
Structural problems which can be analyzed using the proposed computational framework.


## 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.Figure 4
Relationships of classes in FREEVIB.The input parameters are inserted as properties ofInputData 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 theObjectConstructor 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 Qk is calculated by using the calcMatrixQ() method (equation (2)). Transformation matrix Tk and the global constitutive matrix of the lamina Q¯k(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 K˜D∗ 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 theFreeVibrationAnalysis 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 Fstart : Fstep : Fend (input data). The global dynamic stiffness matrix associated with the unknown displacement components K˜D,nnGω is created through the deriveMatrixKnn() method of the FreeVibrationAnalysis class.As soon as theK˜D,nnGω 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 usingnumpy [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.).Figure 5
Example of possible software extension.


#### 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 K˜D,nnGω, the evaluation of equation (16) should be performed for every frequency of interest, in the range Fstart : Fstep : Fend. 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 ∼nCPU (where nCPU is the number of the CPUs in the multicore desktop computer).


## 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.Figure 4
Relationships of classes in FREEVIB.The input parameters are inserted as properties ofInputData 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 theObjectConstructor 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 Qk is calculated by using the calcMatrixQ() method (equation (2)). Transformation matrix Tk and the global constitutive matrix of the lamina Q¯k(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 K˜D∗ 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 theFreeVibrationAnalysis 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 Fstart : Fstep : Fend (input data). The global dynamic stiffness matrix associated with the unknown displacement components K˜D,nnGω is created through the deriveMatrixKnn() method of the FreeVibrationAnalysis class.As soon as theK˜D,nnGω 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 usingnumpy [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.).Figure 5
Example of possible software extension.


### 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 K˜D,nnGω, the evaluation of equation (16) should be performed for every frequency of interest, in the range Fstart : Fstep : Fend. 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 ∼nCPU (where nCPU is the number of the CPUs in the multicore desktop computer).


## 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.).Figure 5
Example of possible software extension.


## 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 K˜D,nnGω, the evaluation of equation (16) should be performed for every frequency of interest, in the range Fstart : Fstep : Fend. 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 ∼nCPU (where nCPU 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 factork = 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 (Figure6). Plate geometry is defined as h1/a = 0.05, h2/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 usingM = 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 Figure7. Excellent agreement is obtained.Figure 6
Geometry and boundary conditions of the isotropic stepped plate.Table 2
Dimensionless frequenciesω∗ = ωa (ρ(1 − ν2)/E)1/2 of stepped concrete isotropic CSFS plates (h1/a = 0.05, h2/a = 0.03).
Mode
CLPT
FSDT
HSDT
Abaqus
1
0.1287
0.1498
0.1498
0.1499
2
0.2319
0.2525
0.2525
0.2526
3
0.3320
0.4679
0.4680
0.4684
4
0.5518
0.5453
0.5453
0.5454
5
0.6122
0.6277
0.6277
0.6281
6
0.6373
0.7749
0.7750
0.7756Figure 7
First 6 mode shapes and corresponding natural frequencies of the isotropic stepped plate, plotted using Abaqus (rows 1 and 3) and FREEVIB (rows 2 and 4).Example 2.
The second example deals with the free vibration analysis of square laminated composite plate with a single stiffener (Figure8). 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:E2/E1 = 40, G12/E2 = G13/E2 = 0.6, G23/E2 = 0.5, and ν12 = 0.25.
The present solution based both on FSDT and HSDT withM = 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.Figure 8
Geometry and boundary conditions of laminated composite plate with a single stiffener.Table 3
Dimensionless frequenciesω∗ = ωa (ρ(1 − ν122)/E2)1/2 of laminated composite plate with a single stiffener (E2/E1 = 40, G12/E2 = G13/E2 = 0.6, G23/E2 = 0.5, and ν12 = 0.25).
Mode
h/a = 0.02
h/a = 0.05
h/a = 0.1
FSDT
HSDT
Abaqus
FSDT
HSDT
Abaqus
FSDT
HSDT
Abaqus
1
0.3197
0.3189
0.3179
0.7548
0.7448
0.7293
1.2855
1.2428
1.1705
2
0.5312
0.5297
0.5287
1.0053
0.9997
0.9870
1.4084
1.3762
1.3223
3
1.0066
0.8389
0.8346
1.9700
1.8417
1.7739
2.1554
2.1546
2.1401
4
1.0838
0.9529
0.9486
2.0863
1.9666
1.9554
2.8721
2.7650
2.5407
5
1.3453
1.3457
1.3407
2.1882
2.0333
1.9702
3.1976
2.8675
2.7107
6
1.5578
1.5583
1.5574
2.3647
2.0889
2.0832
3.3807
2.9119
2.8235
Average difference (%)
6.04
0.32
7.95
1.88
11.30
4.80Example 3.
The final example deals with the free vibration analysis of isotropic cantilever one-fold plate clamped along the edgeX = 0, with h/a = 0.02 and a1/a = a2/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 usingM = 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 + 18M DOFs per edge. This leads to the following equation:(17)NDOFFEM=2500×6=15000,NDOFDSM=7×9+18×5=693.This confirms the computational advantage of the proposed method in comparison with the conventional FEM.Figure 9
Geometry and boundary conditions of isotropic folded plate.Table 4
Dimensionless frequenciesω∗ = ωa (ρ(1-ν2)/E)1/2 of isotropic cantilever folded plate clamped along the edge X = 0 for crank angle α = 120° (h = 0.02a, a1 = a2 = 0.5a).
Mode
FSDT
HSDT
Abaqus
FETM [50]
M = 1
M = 3
M = 5
M = 1
M = 3
M = 5
1
0.0481
0.0484
0.0484
0.0479
0.0482
0.0482
0.0487
0.0491
2
0.0969
0.0967
0.0970
0.0976
0.0976
0.0977
0.0938
0.0943
3
0.1727
0.1756
0.1759
0.1723
0.1743
0.1756
0.1770
0.1787
4
0.2060
0.2068
0.2070
0.2066
0.2075
0.2075
0.2049
0.2065


## 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.



---

*Source: 1369235-2019-01-28.xml*