Document Type: Research Paper
Authors
^{1} University Complex of Materials and Manufacturing Technology, Malek Ashtar University of Technology, Lavizan, Tehran, Iran
^{2} University Complex of Materials and Manufacturing Technology, Maleke Ashtar University of Technology, Lavizan, Tehran, Iran
Abstract
Keywords
Thermal Buckling and Thermal Induced free Vibration Analysis of Perforated Composite Plates: a Mathematical Model
S. Soleimanian, A. Davar, J. Eskandari Jam*, M.R. Zamani, M. Heydari Beni
University Complex of Materials and Manufacturing Technology, Malek Ashtar University of Technology, Lavizan, Tehran., Iran
KEYWORDS 

ABSTRACT 
Thermal buckling Free vibration Perforated plate Heaviside function Mathematical modeling 
This artcile is concerned with thermal buckling and thermal induced free vibration analyses of PCPs (perforated composite plates) with simply supported edges applying a mathematical model. The stiffness and density of PCP are defined locally using Heaviside distribution functions. The governing equations are derived based on CLPT. The present solution gives reasonable results in comparison with the few literatures. In order to inspect the structural behaviour of PCPs subjected to initial thermal loads, many parametric studies have been carried out. Results indicated that the presence of perforations has a significant effect on thermal buckling and thermal induced fundamental frequency. 
The increased application of light weight engineering structures has stimulated interest in the improvement of methods for the analysis of PCPs under different loadings. In order to design PCPs in thermal environments, considerations must be applied to the thermally induced vibration behavior and the possibility of thermal buckling. As the first attempt on thermal buckling of flat or initially imperfect isotropic plates, Gossard et al. [1] performed an approximate solution based on large deflection plate theory. Tauchert and Huang [2] employed the RayleighRitz method in order to solve the thermal buckling equations for symmetric angleply laminated plates. Thangaratnam applied [3] the linear theory and the finite element method (FEM) to solve the thermal buckling problem of composite laminated plates. It was concluded that fiber orientation, number of layers, aspect ratio, and edge conditions could influence the critical buckling temperature and mode shapes significantly. Sunt and Hsu [4] investigated the thermal buckling problem of symmetric crossply laminated plates applying Kirchhoff deformation theory with the inclusion of transverse shear deformation in the displacement field. The critical temperature results were acquired by the Navier solution procedure. They revealed that the discrepancy of results obtained by the formulation with or without shear deformation components is noticeable for length to thickness ratios smaller than 20. Thronton [5] demonstrated research to review the temperature distributions in thin walled structures and thermal buckling analysis methods of isotropic and composite shells and plates. Whitney [6] examined expansional strain effects in laminated plates and derived formulations for bending and thermal buckling problems. Kant and Babu [7] applied two shear deformable finite element models based on first order shear deformation theory (FSDT) and higher order shear deformation theory (HSDT) for thermal buckling of skew fiber reinforced composite and sandwich plates. They indicated that the critical temperature values increase with an increase in skew angle and it is more pronounced in thin laminates than in thick ones. Yapici [8] investigated the thermal buckling behavior of hybridcomposite angleply laminated plates with an inclined crack using FEM. Pradeep and Ganesan [9] explored thermal buckling of multilayer rectangular viscoelastic clamped sandwich plates using FEM. Duran et al. [10] acquired thermal critical buckling temperatures of composite plates with spatial varying fiber orientations using classical lamination theory (CLPT) and FEM. Li et al. [11] applied CLPT in order to investigate buckling and vibroacoustic responses of the clamped composite plates in thermal environments excited by a concentrated harmonic force. Jin et al. [12] performed digital image correlation (DIC) technique to inspect thermal buckling measurement of a circular laminated composite plate under uniform temperature distribution. Through making comparisons among experimental results and those obtained by theory and nonlinear FEM analysis, they exhibited that DIC is promising for studying thermal buckling of composite structures in diverse fields. Bouazza et al. [13] studied thermal buckling of laminated crossply plates applying a refined hyperbolic shear deformation theory. Cetkovic [14] proposed a layerwise displacement model to analyze thermal buckling problem of composite plates and showed that the critical buckling temperature increases by increasing the aspect ratio of plate. Kalita and Haldar [15] reported the free vibration analysis of rectangular isotropic plates using a nine node isoparametric plate element in conjunction with firstorder shear deformation theory contemplating the effect of rotary inertia. Perforated structures primarily have been modeled applying equivalent stiffness. A few researches [1618] have been published to analyze structures including discontinuities using exact modeling. Takabatake [16] examined the static analysis of isotropic plates with voids using the unit step function to define the structure stiffness. Takabatake used the Galerkin method to solve the differential equation of motion. Li and Cheng [17] analyzed grid stiffened composite sandwich panels with simply supported edges subjected to lateral uniform pressure. For an orthogrid stiffened plate, they considered two material regions, the cell, and the surrounding ribs. Based on this concept, they modeled the grid shape in terms of Heaviside functions, which results in the local definition for ABD matrices. The governing equations are solved by deliberating only one component of displacement w, so the solution is limited to symmetric sandwich layups. Wilson et al. [18] have operated research on elastic stability of stepped and stiffened plates. They modeled structures with variation in thickness such as single step, double stepped and latitudinal stiffened plates applying piecewise functions for thickness. The current study presents an analytical solution to thermal buckling and thermal induced free vibration analysis of PCPs taking into account uniform temperature rise and simply supported (SSSS) boundary conditions. Developing a MATLAB code, the analytical results are validated with the results obtained using ABAQUS finite element commercial software and those available in the literature. Many parametric studies have been manifested by applying the present analytical method and ABAQUS.
A rectangular perforated plate lies in (0,0,0) ≤ (x, y, z) ≤ (a, b, h) is considered as indicated in Fig. 1a. The plate consists of a repetitive pattern of rectangular voids with equal distances from each other.
Pursuant to Fig 1. it is required to define a mathematical function considering material properties in orthogonal paths. To this end, Heaviside distribution functions are introduced by Eqs. 1(a) and 1(b) [16].

(1a) 
(1b) 
where (x_{ci}, y_{ci}) is the location of the center of voids. The parameters c and d are the length & width of voids.
The Heaviside distribution (HD) function is presented by Eq. 2 [16].
By plotting the HD function Fig 2., it can be observed that values of one and zero are allocated to white and black regions.
The stiffness matrix of PCP layers can be given by:
where Q^{k} is the stiffness matrix of an orthotropic lamina given by Eq. 4 [6].
Fig. 1. Perforated plate: (a) isometric view, (b) normal view.
Fig. 2. The plot of the HD function.
The linear displacement field is contemplated as [6]:
Applying CLPT, equilibrium equations have been derived as [6]:
Force and moment resultants are given by [6]:
where the A, B, and D coefficients can be acquired by [6]:
And the initial uniform thermal load can be shown by [6]:
where the parameter 𝛥T refers to uniformtemperature difference.
The governing equilibrium equations of PCPs can be obtained by substitution of Eqs. (11) and (15) into Eqs. (810) as:
where the second indices are devoted to local derivatives.
Byconsidering displacement functions corresponding to SSSS edge conditions as [19]:
and applying the Galerkin method, the system of PDEs given by Eqs. (1618) lead to Eqs. (22) and (23) for thermal buckling (T(t)=1) and thermal induced free vibration (T(t)=e^{i}^{𝜔}^{t}) analyses, respectively. The parameters 𝜔 and t refer to frequency and time, respectively.
For thermal buckling and thermal induced free vibration analyses, the eigenvalue problems given by Eq. (22) and (23) [19] have been solved.
where and refer to elastic and thermal stiffness matrix coefficients, respectively. M is the mass matrix, and δ is the displacement vector. , and M matrices are expanded in Appendices (13).
In order to verify the analytical results, FEM models are performed using ABAQUS to simulate thermal buckling and thermal induced free vibration problems.
The FEM meshed model Fig 3. is produced utilizing S4R elements. For a PCP taking into account the material and geometry properties according to Tables 1 and 2, mesh convergence is reported for critical temperature difference and fundamental frequency in Table 3. As indicated in Table 3, enough value for element edge length to plate edge length ratio is achieved by 0.00225.
Fig. 3. FEM Meshed model.
Table 1. Material properties for the PCP [20]
ρ (kg/m^{3}) 
α_{2 } (1/°C) 
α_{1} _{ }(1/°C) 
υ_{12} 
G_{12 }(GPa) 
E_{2 }(GPa) 
E_{1} (GPa) 
1800 
22.1×10^{6} 
8.6×10^{6} 
0.26 
4.14 
8.27 
38.6 
Table 2. Geometry properties for the PCP
R 
h(mm) 
b(mm) 
a (mm) 
10% 
4 
200 
200 
Table 3. Mesh convergence for the PCP
element size to plate length ratio 
𝛥T_{c} (°C) 
Fundamental frequency (Hz) 
A partial view of the meshed model 
0.02 
38.228 
258.62 

0.01 
37.936 
255.42 

0.005 
38.038 
255.58 

0.00225 
38.072 
255.18 

0.00125 
38.067 
255.06 
The accuracy of analytical solutions for buckling and free vibration are compared with those reported in the literature. For an isotropic plate (υ=0.3, α=7.4 ×106 1/°C) considering different values of a/h ratio, the accuracy of the present analytical and ABAQUS methods for thermal buckling are examined through comparisons with the Ref. [14] as indicated in Table 4. A maximum discrepancy of 1.4% is acquired between the present analytical and Cetkovic results. Consequently, the present analytical method can be applied reliably for isotropic plates with a/h ratios bigger than or equal to 20.
Considering a composite plate (layup: [0/90] s, E_{1}/E_{2}=20, G_{12}/E_{2}=G_{13}/E_{2}=G_{23}/E_{2}=0.5, υ_{12}=υ_{13}=υ_{23}= 0.25, α_{1}/α_{2}= 2), the second validation is conducted to check the accuracy of the present analytical solution with the analytical method developed by Bouzza [13]. As indicated in Fig 4., the ABAQUS results are closer to the Ref. [13] values in comparison with the analytical results. The maximum discrepancy between the present analytical results and Ref. [13] is about 5.9%, which corresponds to a/h=2.5.
Table 5 illustrated the results of Nondimensional frequency (𝜔*=𝜔a^{2}√ (ρh⁄ (D (2,2)))) with corresponding mode shapes for an isotropic plate with singular central perforation (υ=0.3 and c/a=0.2). Table 5 includes the results obtained by Kalita and Haldar [15] and the present study results and the maximum discrepancy is 3.08%. The present ABAQUS results are in closer agreement with the Ref. [15] since both methods are developed based on nonlinear displacement field.
In order to evaluate the thermal buckling behavior of PCPs efficiently, some case studies are demonstrated.
For examined perforated plates, the number of voids (m_{x} = n_{y}) are contemplated equal to 20, and the volume fraction of voids is defined as:
Fig. 4. Comparison of critical temperature difference results versus a/b ratio for a composite plate.
Table 4. Comparisons of critical temperature difference results (°C) for an isotropic plate (υ=0.3, α=7.4 ×10^{6} 1/°C)
Solution Method 
a/h 


20 
40 
60 
80 
100 

CLPT (2002) (Adapted from [14]) 
427.477 
106.869 
47.497 
26.717 
17.099 

Cetkovic (2016) [14] 
421.584 
106.497 
47.423 
26.694 
17.089 

Present Analytical 
427.478 
106.869 
47.498 
26.717 
17.099 

Present ABAQUS 
403.99 
104.27 
46.828 
26.435 
16.964 

Discrepancy between present analytical and Cetkovic results (2016) 
1.4% 
0.35% 
0.160% 
0.09% 
0.06% 

Table 5. Nondimensional frequency (𝜔^{*}= 𝜔a^{2} results for isotropic plate with singular central perforation (υ=0.3, c/a=1/5)
Solution Method 
Mode(1,1) 
Mode(1,2) 
Mode(2,1) 
Mode(2,2) 

19.1 
47.496 
47.496 
76.25 

Kalita (2016) [15] 

Present ABAQUS 
19.128 
47.6649 
47.6649 
76.39 

Present analytical 
18. 753 
48.961 
48.961 
77.17 

% Discrepancy, Present 
1.82 
3.08 
3.08 
1.21 

Analytical 

% Discrepancy, Present 
0.15 
0.36 
0.36 
0.18 

ABAQUS 
In order to study the effect of parameters E2/E1 and α2/α1 on thermal buckling of PCPs, the material properties are shown in Table 6. The stacking sequence is assumed to be [0]. As illustrated in Fig. 5, by increasing the values of E2/E1 and α2/α1 ratios, the critical temperature difference decreases. A maximum discrepancy of 5.82% is accomplished between ABAQUS and analytical results where E2/E1=0.2 and α2/α1=1.2.
The critical temperature difference of PCP and the same weight monolithic composite plate (SWCCP) is indicated in Table 7. The stacking sequence is considered to be [0/90] s. As the volume fraction of the voids increases, the PCP demonstrates higher resistance to thermal buckling than the SWCCP. It can be observed that a PCP can ameliorate the critical temperature difference four times higher than SWCCP at R=50%. Figs 6(a). and 6(b). indicate the critical thermal buckling modes for PCP and SWCCP, respectively.
In order to explore the effect of the E2/E1 ratio on the thermal induced fundamental frequency of PCPs, the material properties are considered as presented in Table 8.
Table 6. Material properties to study the effect of E_{2}/E_{1}and α_{2}/α_{1} ratios
_{ }α_{2 }(1/°C) 
α_{1}(1/°C) 
υ_{12} 
G_{12 }(GPa) 
E_{2} (GPa) 
E_{1} (GPa) 
variable 
10^{6} 
0.26 
15.873 
variable 
40 
Table 7. Critical temperature difference results for PCP (R=50%) and SWCCP.
R (%) 
PCP [0/90]s 
SWCCP [0/90]s 
𝛥T_{c (SWCCP)}/𝛥T_{c (PCP)} 

t (mm) 
𝛥T_{c (Analytical)} (°C) 
𝛥T_{c (ABAQUS)} (°C) 
t (mm) 
𝛥T_{c (Analytical)} (°C) 
𝛥T_{c (ABAQUS)} (°C) 
Analytical 
ABAQUS 

10 
4 
38.2243 
38.072 
3.6 
30.9689 
30.531 
1.23 
1.25 
20 
4 
38.2153 
38.216 
3.2 
24.4693 
24.168 
1.56 
1.58 
30 
4 
38.2062 
37.773 
2.8 
18.7343 
18.537 
2.04 
2.04 
40 
4 
38.1971 
36.773 
2.4 
13.764 
13.642 
2.78 
2.7 
50 
4 
38.188 
35.279 
2 
9.5583 
9.4896 
4 
3.72 
(a) (b)
Fig. 6. Thermal buckling critical mode shape obtained by modeling a quarter of the model using ABAQUS: (a) SWCCP, (b) PCP (R=50%).
The stacking sequence is assumed to be [0]. As displayed in Fig 7., by increasing the E2/E1 ratio, the thermal induced fundamental frequencies demonstrate ascending and descending behavior before and after the intersection point at the coordinate (32.28,272.67). A maximum discrepancy of 8.27% is achieved between ABAQUS and analytical results at E2/E1=0.6.
In order to study the effect of the α2/α1 ratio on the thermal induced fundamental frequency of PCPs, the material properties are considered as indicated in Table 6. The stacking sequence is assumed [0]. As portrayed in Fig. 8, by increasing the α2/α1 ratio, the thermal induced fundamental frequency is decreased. A maximum discrepancy of 9.14% is observed between ABAQUS and analytical results for T=0 and α2/α1=1.25.
Table 8. Material properties to study the effect of E_{2}/E_{1} ratio
_{ }α_{2}(1/°C) 
α_{1}(1/°C) 
υ_{12} 
G_{12}(GPa) 
E_{2} (GPa) 
E_{1} (GPa) 
10^{6} 
10^{6} 
0.26 
15.873 
variable 
40 
Fig. 7. Analytical results for the thermal induced fundamental frequency of PCPs considering different E_{2}/E_{1} ratios.
In order to study the effect of the stacking sequence on the free vibration behavior of PCPs, [0/90] s and [0/90] stacking sequences are considered. Pursuant to Fig. 9, the gap between the thermal induced fundamental frequency corresponded to the symmetric stacking sequence, and the asymmetric one is large. It can be concluded that the PCP with symmetric stacking sequence vibrates with higher frequency. Furthermore, it can be observed that the symmetric the structure buckles at higher temperature difference than the asymmetric one.
The thermal induced fundamental frequency of PCP and SWCCP is indicated in Table 7 for the void volume fraction of R=10%. The fundamental frequency of the PCP decreases with lower rate in comparison with SWCCP as the temperature difference increases
Fig. 8. Analytical results for the thermal induced fundamental frequency of PCPs considering different α_{2}/α_{1} ratios.
(Hz) 
A mathematical approach to thermal buckling and thermal induced free vibration analyses of perforated composite plates (PCPs) are described and discussed in this research. The present results are compared to the previously published researches [1315], and they were found to be in close agreement.
Pursuant to results, several concluding remarks are listed as follows:
=
=
=
=
=
=
[1] Gossard ML, Seide P, Roberts WM. Thermal Buckling of Plates, NACA TND. USA; 1995.
[2] Tauchert TR, Huang N N. Thermal Buckling of Symmetric Angleply Laminated Plates. Compos. Struct 1987; 4: 42435.
[3] Thangaratnam KR, Palaninathan, Ramachandaran J. Thermal Buckling of Composite Laminated Plates. Comput. Struct 1989; 32: 111724.
[4] Sunt LX, Hsij TR. Thermal Buckling of Laminated Composite Plates with Transverse Shear Deformation. Comput. Struct 1990; 36: 88389.
[5] Thornton EA. Thermal Buckling of Plates and Shells. Appl. Mech 1993; 46: 48506.
[6] Whitney JM. Structural Analysis of Laminated Anisotropic Plates. Lancaster, PA: Technomic Publishing Co. Inc; 1987.
[7] Kant T, Babu CS. Thermal Buckling Analysis of Skew FiberReinforced Composite and Sandwich Plates Using Shear Deformable Finite Element Models. Compos. Struct 2000; 49: 7785.
[8] Yapici A. Thermal Buckling Behavior of HybridComposite AnglePly Laminated Plates with an Inclined Crack. Mech. Compos. Mat 2005; 41: 13138.
[9] Pradeep V, Ganesan N. Thermal Buckling and Vibration Behavior of MultiLayer Rectangular Viscoelastic Sandwich Plates. J. Sound. Vib 2008; 310: 169 – 83.
[10] Duran AV, Fasanella NA, Sundararaghavan V, Waas A. Thermal Buckling of Composite Plates with Spatial Varying Fiber Orientations, Compos. Struct 2015; 124: 22835.
[11] Li X, Yu K, Han J, Song H, Zhao R. Buckling and VibroAcoustic Response of the Clamped Composite Laminated Plate in Thermal Environment. Int. J. Mech. Sci 2016; 119: 37082.
[12] Jin T, Ha NS, Le VT, Goo NS, Thermal Buckling Measurement of a Laminated Composite Plate Under a Uniform Temperature Distribution Using the Digital Image Correlation Method. Compos. Struct 2015; 123: 4202.
[13] Bouzza M, Lairedj A, Benseddiq N, Khalki S. a Refined Hyperbolic Shear Deformation Theory for Thermal Buckling Analysis of CrossPly Laminated Plates. Mech. Res. Commun 2016; 73: 11726.
[14] Cetkovic M, Thermal Buckling of Laminated Composite Plates Using Layerwise Displacement Model, Compos. Struct 2016; 142: 23853.
[15] Kalita K, Haldar S. Free Vibration Analysis of Rectangular Plates with Central Cutout. Cogent Eng 2016; 3: 1163781.
[16] Takabatake H. Static Analayses of Elastic Plates with Voids. Int. J. Solids. Struct 1991; 28: 17996.
[17] Li G, Cheng J. A Generalized Analytical Modeling of Grid Stiffened Composite Structures. Compos. Struct 2012; 94: 1117–27.
[18] Wilson A J, Rajasekaran S. Stability of All Edges Simply Supported, Stepped and Stiffened Rectangular Plate Under Biaxial Loading. Appl. Math. Modeling 2014; 38: 479–95.
[19] Jones R M. Buckling of Bars, Plates and Shells. Bull Ridge Corporation; 2006.
[20] Kaw A K. Mechanics of Composite Materials, Second Edition. New York: CRC Press, Taylor & Francis Group 2006.